filters.c

   1 #include "image.h"
   2 #include <stdlib.h>
   3 #include <math.h>
   4
   5
   6 /*
   7 =head1 NAME
   8
   9 filters.c - implements filters that operate on images
  10
  11 =head1 SYNOPSIS
  12
  13
  14   i_contrast(im, 0.8);
  15   i_hardinvert(im);
  16   // and more
  17
  18 =head1 DESCRIPTION
  19
  20 filters.c implements basic filters for Imager.  These filters
  21 should be accessible from the filter interface as defined in
  22 the pod for Imager.
  23
  24 =head1 FUNCTION REFERENCE
  25
  26 Some of these functions are internal.
  27
  28 =over 4
  29
  30 =cut
  31 */
  32
  33
  34
  35
  36
  37
  38
  39 /*
  40 =item i_contrast(im, intensity)
  41
  42 Scales the pixel values by the amount specified.
  43
  44   im        - image object
  45   intensity - scalefactor
  46
  47 =cut
  48 */
  49
  50 void
  51 i_contrast(i_img *im, float intensity) {
  52   int x, y;
  53   unsigned char ch;
  54   unsigned int new_color;
  55   i_color rcolor;
  56
  57   mm_log((1,"i_contrast(im %p, intensity %f)\n", im, intensity));
  58
  59   if(intensity < 0) return;
  60
  61   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
  62     i_gpix(im, x, y, &rcolor);
  63
  64     for(ch = 0; ch < im->channels; ch++) {
  65       new_color = (unsigned int) rcolor.channel[ch];
  66       new_color *= intensity;
  67
  68       if(new_color > 255) {
  69         new_color = 255;
  70       }
  71       rcolor.channel[ch] = (unsigned char) new_color;
  72     }
  73     i_ppix(im, x, y, &rcolor);
  74   }
  75 }
  76
  77
  78 /*
  79 =item i_hardinvert(im)
  80
  81 Inverts the pixel values of the input image.
  82
  83   im        - image object
  84
  85 =cut
  86 */
  87
  88 void
  89 i_hardinvert(i_img *im) {
  90   int x, y;
  91   unsigned char ch;
  92
  93   i_color rcolor;
  94
  95     mm_log((1,"i_hardinvert(im %p)\n", im));
  96
  97   for(y = 0; y < im->ysize; y++) {
  98     for(x = 0; x < im->xsize; x++) {
  99       i_gpix(im, x, y, &rcolor);
 100
 101       for(ch = 0; ch < im->channels; ch++) {
 102         rcolor.channel[ch] = 255 - rcolor.channel[ch];
 103       }
 104
 105       i_ppix(im, x, y, &rcolor);
 106     }
 107   }
 108 }
 109
 110
 111
 112 /*
 113 =item i_noise(im, amount, type)
 114
 115 Inverts the pixel values by the amount specified.
 116
 117   im     - image object
 118   amount - deviation in pixel values
 119   type   - noise individual for each channel if true
 120
 121 =cut
 122 */
 123
 124 #ifdef _MSC_VER
 125 /* random() is non-ASCII, even if it is better than rand() */
 126 #define random() rand()
 127 #endif
 128
 129 void
 130 i_noise(i_img *im, float amount, unsigned char type) {
 131   int x, y;
 132   unsigned char ch;
 133   int new_color;
 134   float damount = amount * 2;
 135   i_color rcolor;
 136   int color_inc = 0;
 137
 138   mm_log((1,"i_noise(im %p, intensity %.2f\n", im, amount));
 139
 140   if(amount < 0) return;
 141
 142   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 143     i_gpix(im, x, y, &rcolor);
 144
 145     if(type == 0) {
 146       color_inc = (amount - (damount * ((float)random() / RAND_MAX)));
 147     }
 148
 149     for(ch = 0; ch < im->channels; ch++) {
 150       new_color = (int) rcolor.channel[ch];
 151
 152       if(type != 0) {
 153         new_color += (amount - (damount * ((float)random() / RAND_MAX)));
 154       } else {
 155         new_color += color_inc;
 156       }
 157
 158       if(new_color < 0) {
 159         new_color = 0;
 160       }
 161       if(new_color > 255) {
 162         new_color = 255;
 163       }
 164
 165       rcolor.channel[ch] = (unsigned char) new_color;
 166     }
 167
 168     i_ppix(im, x, y, &rcolor);
 169   }
 170 }
 171
 172
 173 /*
 174 =item i_noise(im, amount, type)
 175
 176 Inverts the pixel values by the amount specified.
 177
 178   im     - image object
 179   amount - deviation in pixel values
 180   type   - noise individual for each channel if true
 181
 182 =cut
 183 */
 184
 185
 186 /*
 187 =item i_applyimage(im, add_im, mode)
 188
 189 Apply's an image to another image
 190
 191   im     - target image
 192   add_im - image that is applied to target
 193   mode   - what method is used in applying:
 194
 195   0   Normal
 196   1   Multiply
 197   2   Screen
 198   3   Overlay
 199   4   Soft Light
 200   5   Hard Light
 201   6   Color dodge
 202   7   Color Burn
 203   8   Darker
 204   9   Lighter
 205   10  Add
 206   11  Subtract
 207   12  Difference
 208   13  Exclusion
 209
 210 =cut
 211 */
 212
 213 void i_applyimage(i_img *im, i_img *add_im, unsigned char mode) {
 214   int x, y;
 215   int mx, my;
 216
 217   mm_log((1, "i_applyimage(im %p, add_im %p, mode %d", im, add_im, mode));
 218
 219   mx = (add_im->xsize <= im->xsize) ? add_im->xsize : add_im->xsize;
 220   my = (add_im->ysize <= im->ysize) ? add_im->ysize : add_im->ysize;
 221
 222   for(x = 0; x < mx; x++) {
 223     for(y = 0; y < my; y++) {
 224     }
 225   }
 226 }
 227
 228
 229 /*
 230 =item i_bumpmap(im, bump, channel, light_x, light_y, st)
 231
 232 Makes a bumpmap on image im using the bump image as the elevation map.
 233
 234   im      - target image
 235   bump    - image that contains the elevation info
 236   channel - to take the elevation information from
 237   light_x - x coordinate of light source
 238   light_y - y coordinate of light source
 239   st      - length of shadow
 240
 241 =cut
 242 */
 243
 244 void
 245 i_bumpmap(i_img *im, i_img *bump, int channel, int light_x, int light_y, int st) {
 246   int x, y, ch;
 247   int mx, my;
 248   i_color x1_color, y1_color, x2_color, y2_color, dst_color;
 249   double nX, nY;
 250   double tX, tY, tZ;
 251   double aX, aY, aL;
 252   double fZ;
 253   unsigned char px1, px2, py1, py2;
 254
 255   i_img new_im;
 256
 257   mm_log((1, "i_bumpmap(im %p, add_im %p, channel %d, light_x %d, light_y %d, st %d)\n",
 258           im, bump, channel, light_x, light_y, st));
 259
 260
 261   if(channel >= bump->channels) {
 262     mm_log((1, "i_bumpmap: channel = %d while bump image only has %d channels\n", channel, bump->channels));
 263     return;
 264   }
 265
 266   mx = (bump->xsize <= im->xsize) ? bump->xsize : im->xsize;
 267   my = (bump->ysize <= im->ysize) ? bump->ysize : im->ysize;
 268
 269   i_img_empty_ch(&new_im, im->xsize, im->ysize, im->channels);
 270
 271   aX = (light_x > (mx >> 1)) ? light_x : mx - light_x;
 272   aY = (light_y > (my >> 1)) ? light_y : my - light_y;
 273
 274   aL = sqrt((aX * aX) + (aY * aY));
 275
 276   for(y = 1; y < my - 1; y++) {
 277     for(x = 1; x < mx - 1; x++) {
 278       i_gpix(bump, x + st, y, &x1_color);
 279       i_gpix(bump, x, y + st, &y1_color);
 280       i_gpix(bump, x - st, y, &x2_color);
 281       i_gpix(bump, x, y - st, &y2_color);
 282
 283       i_gpix(im, x, y, &dst_color);
 284
 285       px1 = x1_color.channel[channel];
 286       py1 = y1_color.channel[channel];
 287       px2 = x2_color.channel[channel];
 288       py2 = y2_color.channel[channel];
 289
 290       nX = px1 - px2;
 291       nY = py1 - py2;
 292
 293       nX += 128;
 294       nY += 128;
 295
 296       fZ = (sqrt((nX * nX) + (nY * nY)) / aL);
 297
 298       tX = abs(x - light_x) / aL;
 299       tY = abs(y - light_y) / aL;
 300
 301       tZ = 1 - (sqrt((tX * tX) + (tY * tY)) * fZ);
 302
 303       if(tZ < 0) tZ = 0;
 304       if(tZ > 2) tZ = 2;
 305
 306       for(ch = 0; ch < im->channels; ch++)
 307         dst_color.channel[ch] = (unsigned char) (float)(dst_color.channel[ch] * tZ);
 308
 309       i_ppix(&new_im, x, y, &dst_color);
 310     }
 311   }
 312
 313   i_copyto(im, &new_im, 0, 0, (int)im->xsize, (int)im->ysize, 0, 0);
 314
 315   i_img_exorcise(&new_im);
 316 }
 317
 318
 319
 320 /*
 321 =item i_postlevels(im, levels)
 322
 323 Quantizes Images to fewer levels.
 324
 325   im      - target image
 326   levels  - number of levels
 327
 328 =cut
 329 */
 330
 331 void
 332 i_postlevels(i_img *im, int levels) {
 333   int x, y, ch;
 334   float pv;
 335   int rv;
 336   float av;
 337
 338   i_color rcolor;
 339
 340   rv = (int) ((float)(256 / levels));
 341   av = (float)levels;
 342
 343   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 344     i_gpix(im, x, y, &rcolor);
 345
 346     for(ch = 0; ch < im->channels; ch++) {
 347       pv = (((float)rcolor.channel[ch] / 255)) * av;
 348       pv = (int) ((int)pv * rv);
 349
 350       if(pv < 0) pv = 0;
 351       else if(pv > 255) pv = 255;
 352
 353       rcolor.channel[ch] = (unsigned char) pv;
 354     }
 355     i_ppix(im, x, y, &rcolor);
 356   }
 357 }
 358
 359
 360 /*
 361 =item i_mosaic(im, size)
 362
 363 Makes an image looks like a mosaic with tilesize of size
 364
 365   im      - target image
 366   size    - size of tiles
 367
 368 =cut
 369 */
 370
 371 void
 372 i_mosaic(i_img *im, int size) {
 373   int x, y, ch;
 374   int lx, ly, z;
 375   long sqrsize;
 376
 377   i_color rcolor;
 378   long col[256];
 379
 380   sqrsize = size * size;
 381
 382   for(y = 0; y < im->ysize; y += size) for(x = 0; x < im->xsize; x += size) {
 383     for(z = 0; z < 256; z++) col[z] = 0;
 384
 385     for(lx = 0; lx < size; lx++) {
 386       for(ly = 0; ly < size; ly++) {
 387         i_gpix(im, (x + lx), (y + ly), &rcolor);
 388
 389         for(ch = 0; ch < im->channels; ch++) {
 390           col[ch] += rcolor.channel[ch];
 391         }
 392       }
 393     }
 394
 395     for(ch = 0; ch < im->channels; ch++)
 396       rcolor.channel[ch] = (int) ((float)col[ch] / sqrsize);
 397
 398
 399     for(lx = 0; lx < size; lx++)
 400       for(ly = 0; ly < size; ly++)
 401       i_ppix(im, (x + lx), (y + ly), &rcolor);
 402
 403   }
 404 }
 405
 406 /*
 407 =item saturate(in)
 408
 409 Clamps the input value between 0 and 255. (internal)
 410
 411   in - input integer
 412
 413 =cut
 414 */
 415
 416 static
 417 unsigned char
 418 saturate(int in) {
 419   if (in>255) { return 255; }
 420   else if (in>0) return in;
 421   return 0;
 422 }
 423
 424
 425 /*
 426 =item i_watermark(im, wmark, tx, ty, pixdiff)
 427
 428 Applies a watermark to the target image
 429
 430   im      - target image
 431   wmark   - watermark image
 432   tx      - x coordinate of where watermark should be applied
 433   ty      - y coordinate of where watermark should be applied
 434   pixdiff - the magnitude of the watermark, controls how visible it is
 435
 436 =cut
 437 */
 438
 439 void
 440 i_watermark(i_img *im, i_img *wmark, int tx, int ty, int pixdiff) {
 441   int vx, vy, ch;
 442   i_color val, wval;
 443
 444   for(vx=0;vx<128;vx++) for(vy=0;vy<110;vy++) {
 445
 446     i_gpix(im,    tx+vx, ty+vy,&val );
 447     i_gpix(wmark, vx,    vy,   &wval);
 448
 449     for(ch=0;ch<im->channels;ch++)
 450       val.channel[ch] = saturate( val.channel[ch] + (pixdiff* (wval.channel[0]-128) )/128 );
 451
 452     i_ppix(im,tx+vx,ty+vy,&val);
 453   }
 454 }
 455
 456
 457 /*
 458 =item i_autolevels(im, lsat, usat, skew)
 459
 460 Scales and translates each color such that it fills the range completely.
 461 Skew is not implemented yet - purpose is to control the color skew that can
 462 occur when changing the contrast.
 463
 464   im   - target image
 465   lsat - fraction of pixels that will be truncated at the lower end of the spectrum
 466   usat - fraction of pixels that will be truncated at the higher end of the spectrum
 467   skew - not used yet
 468
 469 =cut
 470 */
 471
 472 void
 473 i_autolevels(i_img *im, float lsat, float usat, float skew) {
 474   i_color val;
 475   int i, x, y, rhist[256], ghist[256], bhist[256];
 476   int rsum, rmin, rmax;
 477   int gsum, gmin, gmax;
 478   int bsum, bmin, bmax;
 479   int rcl, rcu, gcl, gcu, bcl, bcu;
 480
 481   mm_log((1,"i_autolevels(im %p, lsat %f,usat %f,skew %f)\n", im, lsat,usat,skew));
 482
 483   rsum=gsum=bsum=0;
 484   for(i=0;i<256;i++) rhist[i]=ghist[i]=bhist[i] = 0;
 485   /* create histogram for each channel */
 486   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 487     i_gpix(im, x, y, &val);
 488     rhist[val.channel[0]]++;
 489     ghist[val.channel[1]]++;
 490     bhist[val.channel[2]]++;
 491   }
 492
 493   for(i=0;i<256;i++) {
 494     rsum+=rhist[i];
 495     gsum+=ghist[i];
 496     bsum+=bhist[i];
 497   }
 498
 499   rmin = gmin = bmin = 0;
 500   rmax = gmax = bmax = 255;
 501
 502   rcu = rcl = gcu = gcl = bcu = bcl = 0;
 503
 504   for(i=0; i<256; i++) {
 505     rcl += rhist[i];     if ( (rcl<rsum*lsat) ) rmin=i;
 506     rcu += rhist[255-i]; if ( (rcu<rsum*usat) ) rmax=255-i;
 507
 508     gcl += ghist[i];     if ( (gcl<gsum*lsat) ) gmin=i;
 509     gcu += ghist[255-i]; if ( (gcu<gsum*usat) ) gmax=255-i;
 510
 511     bcl += bhist[i];     if ( (bcl<bsum*lsat) ) bmin=i;
 512     bcu += bhist[255-i]; if ( (bcu<bsum*usat) ) bmax=255-i;
 513   }
 514
 515   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 516     i_gpix(im, x, y, &val);
 517     val.channel[0]=saturate((val.channel[0]-rmin)*255/(rmax-rmin));
 518     val.channel[1]=saturate((val.channel[1]-gmin)*255/(gmax-gmin));
 519     val.channel[2]=saturate((val.channel[2]-bmin)*255/(bmax-bmin));
 520     i_ppix(im, x, y, &val);
 521   }
 522 }
 523
 524 /*
 525 =item Noise(x,y)
 526
 527 Pseudo noise utility function used to generate perlin noise. (internal)
 528
 529   x - x coordinate
 530   y - y coordinate
 531
 532 =cut
 533 */
 534
 535 static
 536 float
 537 Noise(int x, int y) {
 538   int n = x + y * 57;
 539   n = (n<<13) ^ n;
 540   return ( 1.0 - ( (n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
 541 }
 542
 543 /*
 544 =item SmoothedNoise1(x,y)
 545
 546 Pseudo noise utility function used to generate perlin noise. (internal)
 547
 548   x - x coordinate
 549   y - y coordinate
 550
 551 =cut
 552 */
 553
 554 static
 555 float
 556 SmoothedNoise1(float x, float y) {
 557   float corners = ( Noise(x-1, y-1)+Noise(x+1, y-1)+Noise(x-1, y+1)+Noise(x+1, y+1) ) / 16;
 558   float sides   = ( Noise(x-1, y)  +Noise(x+1, y)  +Noise(x, y-1)  +Noise(x, y+1) ) /  8;
 559   float center  =  Noise(x, y) / 4;
 560   return corners + sides + center;
 561 }
 562
 563
 564 /*
 565 =item G_Interpolate(a, b, x)
 566
 567 Utility function used to generate perlin noise. (internal)
 568
 569 =cut
 570 */
 571
 572 static
 573 float C_Interpolate(float a, float b, float x) {
 574   /*  float ft = x * 3.1415927; */
 575   float ft = x * PI;
 576   float f = (1 - cos(ft)) * .5;
 577   return  a*(1-f) + b*f;
 578 }
 579
 580
 581 /*
 582 =item InterpolatedNoise(x, y)
 583
 584 Utility function used to generate perlin noise. (internal)
 585
 586 =cut
 587 */
 588
 589 static
 590 float
 591 InterpolatedNoise(float x, float y) {
 592
 593   int integer_X = x;
 594   float fractional_X = x - integer_X;
 595   int integer_Y = y;
 596   float fractional_Y = y - integer_Y;
 597
 598   float v1 = SmoothedNoise1(integer_X,     integer_Y);
 599   float v2 = SmoothedNoise1(integer_X + 1, integer_Y);
 600   float v3 = SmoothedNoise1(integer_X,     integer_Y + 1);
 601   float v4 = SmoothedNoise1(integer_X + 1, integer_Y + 1);
 602
 603   float i1 = C_Interpolate(v1 , v2 , fractional_X);
 604   float i2 = C_Interpolate(v3 , v4 , fractional_X);
 605
 606   return C_Interpolate(i1 , i2 , fractional_Y);
 607 }
 608
 609
 610
 611 /*
 612 =item PerlinNoise_2D(x, y)
 613
 614 Utility function used to generate perlin noise. (internal)
 615
 616 =cut
 617 */
 618
 619 static
 620 float
 621 PerlinNoise_2D(float x, float y) {
 622   int i,frequency;
 623   float amplitude;
 624   float total = 0;
 625   int Number_Of_Octaves=6;
 626   int n = Number_Of_Octaves - 1;
 627
 628   for(i=0;i<n;i++) {
 629     frequency = 2*i;
 630     amplitude = PI;
 631     total = total + InterpolatedNoise(x * frequency, y * frequency) * amplitude;
 632   }
 633
 634   return total;
 635 }
 636
 637
 638 /*
 639 =item i_radnoise(im, xo, yo, rscale, ascale)
 640
 641 Perlin-like radial noise.
 642
 643   im     - target image
 644   xo     - x coordinate of center
 645   yo     - y coordinate of center
 646   rscale - radial scale
 647   ascale - angular scale
 648
 649 =cut
 650 */
 651
 652 void
 653 i_radnoise(i_img *im, int xo, int yo, float rscale, float ascale) {
 654   int x, y, ch;
 655   i_color val;
 656   unsigned char v;
 657   float xc, yc, r;
 658   double a;
 659
 660   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 661     xc = (float)x-xo+0.5;
 662     yc = (float)y-yo+0.5;
 663     r = rscale*sqrt(xc*xc+yc*yc)+1.2;
 664     a = (PI+atan2(yc,xc))*ascale;
 665     v = saturate(128+100*(PerlinNoise_2D(a,r)));
 666     /* v=saturate(120+12*PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale));  Good soft marble */
 667     for(ch=0; ch<im->channels; ch++) val.channel[ch]=v;
 668     i_ppix(im, x, y, &val);
 669   }
 670 }
 671
 672
 673 /*
 674 =item i_turbnoise(im, xo, yo, scale)
 675
 676 Perlin-like 2d noise noise.
 677
 678   im     - target image
 679   xo     - x coordinate translation
 680   yo     - y coordinate translation
 681   scale  - scale of noise
 682
 683 =cut
 684 */
 685
 686 void
 687 i_turbnoise(i_img *im, float xo, float yo, float scale) {
 688   int x,y,ch;
 689   unsigned char v;
 690   i_color val;
 691
 692   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 693     /*    v=saturate(125*(1.0+PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale))); */
 694     v = saturate(120*(1.0+sin(xo+(float)x/scale+PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale))));
 695     for(ch=0; ch<im->channels; ch++) val.channel[ch] = v;
 696     i_ppix(im, x, y, &val);
 697   }
 698 }
 699
 700
 701
 702 /*
 703 =item i_gradgen(im, num, xo, yo, ival, dmeasure)
 704
 705 Gradient generating function.
 706
 707   im     - target image
 708   num    - number of points given
 709   xo     - array of x coordinates
 710   yo     - array of y coordinates
 711   ival   - array of i_color objects
 712   dmeasure - distance measure to be used.
 713     0 = Euclidean
 714     1 = Euclidean squared
 715     2 = Manhattan distance
 716
 717 =cut
 718 */
 719
 720
 721 void
 722 i_gradgen(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
 723
 724   i_color val;
 725   int p, x, y, ch;
 726   int channels = im->channels;
 727   int xsize    = im->xsize;
 728   int ysize    = im->ysize;
 729
 730   float *fdist;
 731
 732   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 733
 734   for(p = 0; p<num; p++) {
 735     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
 736     ICL_info(&ival[p]);
 737   }
 738
 739   fdist = mymalloc( sizeof(float) * num );
 740
 741   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 742     float cs = 0;
 743     float csd = 0;
 744     for(p = 0; p<num; p++) {
 745       int xd    = x-xo[p];
 746       int yd    = y-yo[p];
 747       switch (dmeasure) {
 748       case 0: /* euclidean */
 749         fdist[p]  = sqrt(xd*xd + yd*yd); /* euclidean distance */
 750         break;
 751       case 1: /* euclidean squared */
 752         fdist[p]  = xd*xd + yd*yd; /* euclidean distance */
 753         break;
 754       case 2: /* euclidean squared */
 755         fdist[p]  = max(xd*xd, yd*yd); /* manhattan distance */
 756         break;
 757       default:
 758         m_fatal(3,"i_gradgen: Unknown distance measure\n");
 759       }
 760       cs += fdist[p];
 761     }
 762
 763     csd = 1/((num-1)*cs);
 764
 765     for(p = 0; p<num; p++) fdist[p] = (cs-fdist[p])*csd;
 766
 767     for(ch = 0; ch<channels; ch++) {
 768       int tres = 0;
 769       for(p = 0; p<num; p++) tres += ival[p].channel[ch] * fdist[p];
 770       val.channel[ch] = saturate(tres);
 771     }
 772     i_ppix(im, x, y, &val);
 773   }
 774
 775 }
 776
 777 void
 778 i_nearest_color_foo(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
 779
 780   int p, x, y;
 781   int xsize    = im->xsize;
 782   int ysize    = im->ysize;
 783
 784   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 785
 786   for(p = 0; p<num; p++) {
 787     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
 788     ICL_info(&ival[p]);
 789   }
 790
 791   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 792     int   midx    = 0;
 793     float mindist = 0;
 794     float curdist = 0;
 795
 796     int xd        = x-xo[0];
 797     int yd        = y-yo[0];
 798
 799     switch (dmeasure) {
 800     case 0: /* euclidean */
 801       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 802       break;
 803     case 1: /* euclidean squared */
 804       mindist = xd*xd + yd*yd; /* euclidean distance */
 805       break;
 806     case 2: /* euclidean squared */
 807       mindist = max(xd*xd, yd*yd); /* manhattan distance */
 808       break;
 809     default:
 810       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 811     }
 812
 813     for(p = 1; p<num; p++) {
 814       xd = x-xo[p];
 815       yd = y-yo[p];
 816       switch (dmeasure) {
 817       case 0: /* euclidean */
 818         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 819         break;
 820       case 1: /* euclidean squared */
 821         curdist = xd*xd + yd*yd; /* euclidean distance */
 822         break;
 823       case 2: /* euclidean squared */
 824         curdist = max(xd*xd, yd*yd); /* manhattan distance */
 825         break;
 826       default:
 827         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 828       }
 829       if (curdist < mindist) {
 830         mindist = curdist;
 831         midx = p;
 832       }
 833     }
 834     i_ppix(im, x, y, &ival[midx]);
 835   }
 836 }
 837
 838 void
 839 i_nearest_color(i_img *im, int num, int *xo, int *yo, i_color *oval, int dmeasure) {
 840   i_color *ival;
 841   float *tval;
 842   float c1, c2;
 843   i_color val;
 844   int p, x, y, ch;
 845   int xsize    = im->xsize;
 846   int ysize    = im->ysize;
 847   int *cmatch;
 848
 849   mm_log((1,"i_nearest_color(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, oval, dmeasure));
 850
 851   tval   = mymalloc( sizeof(float)*num*im->channels );
 852   ival   = mymalloc( sizeof(i_color)*num );
 853   cmatch = mymalloc( sizeof(int)*num     );
 854
 855   for(p = 0; p<num; p++) {
 856     for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
 857     cmatch[p] = 0;
 858   }
 859
 860
 861   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 862     int   midx    = 0;
 863     float mindist = 0;
 864     float curdist = 0;
 865
 866     int xd        = x-xo[0];
 867     int yd        = y-yo[0];
 868
 869     switch (dmeasure) {
 870     case 0: /* euclidean */
 871       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 872       break;
 873     case 1: /* euclidean squared */
 874       mindist = xd*xd + yd*yd; /* euclidean distance */
 875       break;
 876     case 2: /* euclidean squared */
 877       mindist = max(xd*xd, yd*yd); /* manhattan distance */
 878       break;
 879     default:
 880       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 881     }
 882
 883     for(p = 1; p<num; p++) {
 884       xd = x-xo[p];
 885       yd = y-yo[p];
 886       switch (dmeasure) {
 887       case 0: /* euclidean */
 888         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 889         break;
 890       case 1: /* euclidean squared */
 891         curdist = xd*xd + yd*yd; /* euclidean distance */
 892         break;
 893       case 2: /* euclidean squared */
 894         curdist = max(xd*xd, yd*yd); /* manhattan distance */
 895         break;
 896       default:
 897         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 898       }
 899       if (curdist < mindist) {
 900         mindist = curdist;
 901         midx = p;
 902       }
 903     }
 904
 905     cmatch[midx]++;
 906     i_gpix(im, x, y, &val);
 907     c2 = 1.0/(float)(cmatch[midx]);
 908     c1 = 1.0-c2;
 909
 910     for(ch = 0; ch<im->channels; ch++)
 911       tval[midx*im->channels + ch] = c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
 912
 913
 914   }
 915
 916   for(p = 0; p<num; p++) for(ch = 0; ch<im->channels; ch++) ival[p].channel[ch] = tval[p*im->channels + ch];
 917
 918   i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
 919 }
 920
 921 struct fount_state;
 922 static double linear_fount_f(double x, double y, struct fount_state *state);
 923 static double bilinear_fount_f(double x, double y, struct fount_state *state);
 924 static double radial_fount_f(double x, double y, struct fount_state *state);
 925 static double square_fount_f(double x, double y, struct fount_state *state);
 926 static double revolution_fount_f(double x, double y,
 927                                  struct fount_state *state);
 928 static double conical_fount_f(double x, double y, struct fount_state *state);
 929
 930 typedef double (*fount_func)(double, double, struct fount_state *);
 931 static fount_func fount_funcs[] =
 932 {
 933   linear_fount_f,
 934   bilinear_fount_f,
 935   radial_fount_f,
 936   square_fount_f,
 937   revolution_fount_f,
 938   conical_fount_f,
 939 };
 940
 941 static double linear_interp(double pos, i_fountain_seg *seg);
 942 static double sine_interp(double pos, i_fountain_seg *seg);
 943 static double sphereup_interp(double pos, i_fountain_seg *seg);
 944 static double spheredown_interp(double pos, i_fountain_seg *seg);
 945 typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
 946 static fount_interp fount_interps[] =
 947 {
 948   linear_interp,
 949   linear_interp,
 950   sine_interp,
 951   sphereup_interp,
 952   spheredown_interp,
 953 };
 954
 955 static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
 956 static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
 957 static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
 958 typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
 959 static fount_cinterp fount_cinterps[] =
 960 {
 961   direct_cinterp,
 962   hue_up_cinterp,
 963   hue_down_cinterp,
 964 };
 965
 966 typedef double (*fount_repeat)(double v);
 967 static double fount_r_none(double v);
 968 static double fount_r_sawtooth(double v);
 969 static double fount_r_triangle(double v);
 970 static double fount_r_saw_both(double v);
 971 static double fount_r_tri_both(double v);
 972 static fount_repeat fount_repeats[] =
 973 {
 974   fount_r_none,
 975   fount_r_sawtooth,
 976   fount_r_triangle,
 977   fount_r_saw_both,
 978   fount_r_tri_both,
 979 };
 980
 981 static int simple_ssample(i_fcolor *out, double x, double y,
 982                            struct fount_state *state);
 983 static int random_ssample(i_fcolor *out, double x, double y,
 984                            struct fount_state *state);
 985 static int circle_ssample(i_fcolor *out, double x, double y,
 986                            struct fount_state *state);
 987 typedef int (*fount_ssample)(i_fcolor *out, double x, double y,
 988                               struct fount_state *state);
 989 static fount_ssample fount_ssamples[] =
 990 {
 991   NULL,
 992   simple_ssample,
 993   random_ssample,
 994   circle_ssample,
 995 };
 996
 997 static int
 998 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
 999
1000 /*
1001   Keep state information used by each type of fountain fill
1002 */
1003 struct fount_state {
1004   /* precalculated for the equation of the line perpendicular to the line AB */
1005   double lA, lB, lC;
1006   double AB;
1007   double sqrtA2B2;
1008   double mult;
1009   double cos;
1010   double sin;
1011   double theta;
1012   int xa, ya;
1013   void *ssample_data;
1014   fount_func ffunc;
1015   fount_repeat rpfunc;
1016   fount_ssample ssfunc;
1017   double parm;
1018   i_fountain_seg *segs;
1019   int count;
1020 };
1021
1022 static void
1023 fount_init_state(struct fount_state *state, double xa, double ya,
1024                  double xb, double yb, i_fountain_type type,
1025                  i_fountain_repeat repeat, int combine, int super_sample,
1026                  double ssample_param, int count, i_fountain_seg *segs);
1027
1028 static void
1029 fount_finish_state(struct fount_state *state);
1030
1031 #define EPSILON (1e-6)
1032
1033 /*
1034 =item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1035
1036 Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1037
1038 I<type> controls how the reference points are used:
1039
1040 =over
1041
1042 =item i_ft_linear
1043
1044 linear, where A is 0 and B is 1.
1045
1046 =item i_ft_bilinear
1047
1048 linear in both directions from A.
1049
1050 =item i_ft_radial
1051
1052 circular, where A is the centre of the fill, and B is a point
1053 on the radius.
1054
1055 =item i_ft_radial_square
1056
1057 where A is the centre of the fill and B is the centre of
1058 one side of the square.
1059
1060 =item i_ft_revolution
1061
1062 where A is the centre of the fill and B defines the 0/1.0
1063 angle of the fill.
1064
1065 =item i_ft_conical
1066
1067 similar to i_ft_revolution, except that the revolution goes in both
1068 directions
1069
1070 =back
1071
1072 I<repeat> can be one of:
1073
1074 =over
1075
1076 =item i_fr_none
1077
1078 values < 0 are treated as zero, values > 1 are treated as 1.
1079
1080 =item i_fr_sawtooth
1081
1082 negative values are treated as 0, positive values are modulo 1.0
1083
1084 =item i_fr_triangle
1085
1086 negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1087 mod 1.0), otherwise the same as for sawtooth.
1088
1089 =item i_fr_saw_both
1090
1091 like i_fr_sawtooth, except that the sawtooth pattern repeats into
1092 negative values.
1093
1094 =item i_fr_tri_both
1095
1096 Like i_fr_triangle, except that negative values are handled as their
1097 absolute values.
1098
1099 =back
1100
1101 If combine is non-zero then non-opaque values are combined with the
1102 underlying color.
1103
1104 I<super_sample> controls super sampling, if any.  At some point I'll
1105 probably add a adaptive super-sampler.  Current possible values are:
1106
1107 =over
1108
1109 =item i_fts_none
1110
1111 No super-sampling is done.
1112
1113 =item i_fts_grid
1114
1115 A square grid of points withing the pixel are sampled.
1116
1117 =item i_fts_random
1118
1119 Random points within the pixel are sampled.
1120
1121 =item i_fts_circle
1122
1123 Points on the radius of a circle are sampled.  This produces fairly
1124 good results, but is fairly slow since sin() and cos() are evaluated
1125 for each point.
1126
1127 =back
1128
1129 I<ssample_param> is intended to be roughly the number of points
1130 sampled within the pixel.
1131
1132 I<count> and I<segs> define the segments of the fill.
1133
1134 =cut
1135
1136 */
1137
1138 void
1139 i_fountain(i_img *im, double xa, double ya, double xb, double yb,
1140            i_fountain_type type, i_fountain_repeat repeat,
1141            int combine, int super_sample, double ssample_param,
1142            int count, i_fountain_seg *segs) {
1143   struct fount_state state;
1144   int x, y;
1145   i_fcolor *line = mymalloc(sizeof(i_fcolor) * im->xsize);
1146   int ch;
1147   i_fountain_seg *my_segs;
1148
1149   fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
1150                    super_sample, ssample_param, count, segs);
1151   my_segs = state.segs;
1152
1153   for (y = 0; y < im->ysize; ++y) {
1154     i_glinf(im, 0, im->xsize, y, line);
1155     for (x = 0; x < im->xsize; ++x) {
1156       i_fcolor c;
1157       int got_one;
1158       double v;
1159       if (super_sample == i_fts_none)
1160         got_one = fount_getat(&c, x, y, &state);
1161       else
1162         got_one = state.ssfunc(&c, x, y, &state);
1163       if (got_one) {
1164         if (combine) {
1165           for (ch = 0; ch < im->channels; ++ch) {
1166             line[x].channel[ch] = line[x].channel[ch] * (1.0 - c.channel[3])
1167               + c.channel[ch] * c.channel[3];
1168           }
1169         }
1170         else
1171           line[x] = c;
1172       }
1173     }
1174     i_plinf(im, 0, im->xsize, y, line);
1175   }
1176   fount_finish_state(&state);
1177   myfree(line);
1178 }
1179
1180 typedef struct {
1181   i_fill_t base;
1182   struct fount_state state;
1183 } i_fill_fountain_t;
1184
1185 static void
1186 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
1187             i_fcolor *data);
1188 static void
1189 fount_fill_destroy(i_fill_t *fill);
1190
1191 /*
1192 =item i_new_fount(xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1193
1194 =cut
1195 */
1196
1197 i_fill_t *
1198 i_new_fill_fount(double xa, double ya, double xb, double yb,
1199                  i_fountain_type type, i_fountain_repeat repeat,
1200                  int combine, int super_sample, double ssample_param,
1201                  int count, i_fountain_seg *segs) {
1202   i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1203
1204   fill->base.fill_with_color = NULL;
1205   fill->base.fill_with_fcolor = fill_fountf;
1206   fill->base.destroy = fount_fill_destroy;
1207   fill->base.combines = combine;
1208   fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1209                    super_sample, ssample_param, count, segs);
1210
1211   return &fill->base;
1212 }
1213
1214 /*
1215 =back
1216
1217 =head1 INTERNAL FUNCTIONS
1218
1219 =over
1220
1221 =item fount_init_state(...)
1222
1223 Used by both the fountain fill filter and the fountain fill.
1224
1225 =cut
1226 */
1227
1228 static void
1229 fount_init_state(struct fount_state *state, double xa, double ya,
1230                  double xb, double yb, i_fountain_type type,
1231                  i_fountain_repeat repeat, int combine, int super_sample,
1232                  double ssample_param, int count, i_fountain_seg *segs) {
1233   int i, j;
1234   i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count);
1235   /*int have_alpha = im->channels == 2 || im->channels == 4;*/
1236   int ch;
1237
1238   memset(state, 0, sizeof(*state));
1239   /* we keep a local copy that we can adjust for speed */
1240   for (i = 0; i < count; ++i) {
1241     i_fountain_seg *seg = my_segs + i;
1242
1243     *seg = segs[i];
1244     if (seg->type < 0 || type >= i_ft_end)
1245       seg->type = i_ft_linear;
1246     if (seg->color < 0 || seg->color >= i_fc_end)
1247       seg->color = i_fc_direct;
1248     if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1249       /* so we don't have to translate to HSV on each request, do it here */
1250       for (j = 0; j < 2; ++j) {
1251         i_rgb_to_hsvf(seg->c+j);
1252       }
1253       if (seg->color == i_fc_hue_up) {
1254         if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1255           seg->c[1].channel[0] += 1.0;
1256       }
1257       else {
1258         if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1259           seg->c[0].channel[0] += 1.0;
1260       }
1261     }
1262     /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1263            seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1264            seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1265            seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1266            seg->c[1].channel[3], seg->type, seg->color);*/
1267
1268   }
1269
1270   /* initialize each engine */
1271   /* these are so common ... */
1272   state->lA = xb - xa;
1273   state->lB = yb - ya;
1274   state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1275   state->xa = xa;
1276   state->ya = ya;
1277   switch (type) {
1278   default:
1279     type = i_ft_linear; /* make the invalid value valid */
1280   case i_ft_linear:
1281   case i_ft_bilinear:
1282     state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1283     state->mult = 1;
1284     state->mult = 1/linear_fount_f(xb, yb, state);
1285     break;
1286
1287   case i_ft_radial:
1288     state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1289                              + (double)(yb-ya)*(yb-ya));
1290     break;
1291
1292   case i_ft_radial_square:
1293     state->cos = state->lA / state->AB;
1294     state->sin = state->lB / state->AB;
1295     state->mult = 1.0 / state->AB;
1296     break;
1297
1298   case i_ft_revolution:
1299     state->theta = atan2(yb-ya, xb-xa);
1300     state->mult = 1.0 / (PI * 2);
1301     break;
1302
1303   case i_ft_conical:
1304     state->theta = atan2(yb-ya, xb-xa);
1305     state->mult = 1.0 / PI;
1306     break;
1307   }
1308   state->ffunc = fount_funcs[type];
1309   if (super_sample < 0
1310       || super_sample >= (sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
1311     super_sample = 0;
1312   }
1313   state->ssample_data = NULL;
1314   switch (super_sample) {
1315   case i_fts_grid:
1316     ssample_param = floor(0.5 + sqrt(ssample_param));
1317     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param);
1318     break;
1319
1320   case i_fts_random:
1321   case i_fts_circle:
1322     ssample_param = floor(0.5+ssample_param);
1323     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
1324     break;
1325   }
1326   state->parm = ssample_param;
1327   state->ssfunc = fount_ssamples[super_sample];
1328   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1329     repeat = 0;
1330   state->rpfunc = fount_repeats[repeat];
1331   state->segs = my_segs;
1332   state->count = count;
1333 }
1334
1335 static void
1336 fount_finish_state(struct fount_state *state) {
1337   if (state->ssample_data)
1338     myfree(state->ssample_data);
1339   myfree(state->segs);
1340 }
1341
1342
1343 /*
1344 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
1345
1346 Evaluates the fountain fill at the given point.
1347
1348 This is called by both the non-super-sampling and super-sampling code.
1349
1350 You might think that it would make sense to sample the fill parameter
1351 instead, and combine those, but this breaks badly.
1352
1353 =cut
1354 */
1355
1356 static int
1357 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
1358   double v = (state->rpfunc)((state->ffunc)(x, y, state));
1359   int i;
1360
1361   i = 0;
1362   while (i < state->count
1363          && (v < state->segs[i].start || v > state->segs[i].end)) {
1364     ++i;
1365   }
1366   if (i < state->count) {
1367     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
1368     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
1369     return 1;
1370   }
1371   else
1372     return 0;
1373 }
1374
1375 /*
1376 =item linear_fount_f(x, y, state)
1377
1378 Calculate the fill parameter for a linear fountain fill.
1379
1380 Uses the point to line distance function, with some precalculation
1381 done in i_fountain().
1382
1383 =cut
1384 */
1385 static double
1386 linear_fount_f(double x, double y, struct fount_state *state) {
1387   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
1388 }
1389
1390 /*
1391 =item bilinear_fount_f(x, y, state)
1392
1393 Calculate the fill parameter for a bi-linear fountain fill.
1394
1395 =cut
1396 */
1397 static double
1398 bilinear_fount_f(double x, double y, struct fount_state *state) {
1399   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
1400 }
1401
1402 /*
1403 =item radial_fount_f(x, y, state)
1404
1405 Calculate the fill parameter for a radial fountain fill.
1406
1407 Simply uses the distance function.
1408
1409 =cut
1410  */
1411 static double
1412 radial_fount_f(double x, double y, struct fount_state *state) {
1413   return sqrt((double)(state->xa-x)*(state->xa-x)
1414               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
1415 }
1416
1417 /*
1418 =item square_fount_f(x, y, state)
1419
1420 Calculate the fill parameter for a square fountain fill.
1421
1422 Works by rotating the reference co-ordinate around the centre of the
1423 square.
1424
1425 =cut
1426 */
1427 static double
1428 square_fount_f(double x, double y, struct fount_state *state) {
1429   int xc, yc; /* centred on A */
1430   double xt, yt; /* rotated by theta */
1431   xc = x - state->xa;
1432   yc = y - state->ya;
1433   xt = fabs(xc * state->cos + yc * state->sin);
1434   yt = fabs(-xc * state->sin + yc * state->cos);
1435   return (xt > yt ? xt : yt) * state->mult;
1436 }
1437
1438 /*
1439 =item revolution_fount_f(x, y, state)
1440
1441 Calculates the fill parameter for the revolution fountain fill.
1442
1443 =cut
1444 */
1445 static double
1446 revolution_fount_f(double x, double y, struct fount_state *state) {
1447   double angle = atan2(y - state->ya, x - state->xa);
1448
1449   angle -= state->theta;
1450   if (angle < 0) {
1451     angle = fmod(angle+ PI * 4, PI*2);
1452   }
1453
1454   return angle * state->mult;
1455 }
1456
1457 /*
1458 =item conical_fount_f(x, y, state)
1459
1460 Calculates the fill parameter for the conical fountain fill.
1461
1462 =cut
1463 */
1464 static double
1465 conical_fount_f(double x, double y, struct fount_state *state) {
1466   double angle = atan2(y - state->ya, x - state->xa);
1467
1468   angle -= state->theta;
1469   if (angle < -PI)
1470     angle += PI * 2;
1471   else if (angle > PI)
1472     angle -= PI * 2;
1473
1474   return fabs(angle) * state->mult;
1475 }
1476
1477 /*
1478 =item linear_interp(pos, seg)
1479
1480 Calculates linear interpolation on the fill parameter.  Breaks the
1481 segment into 2 regions based in the I<middle> value.
1482
1483 =cut
1484 */
1485 static double
1486 linear_interp(double pos, i_fountain_seg *seg) {
1487   if (pos < seg->middle) {
1488     double len = seg->middle - seg->start;
1489     if (len < EPSILON)
1490       return 0.0;
1491     else
1492       return (pos - seg->start) / len / 2;
1493   }
1494   else {
1495     double len = seg->end - seg->middle;
1496     if (len < EPSILON)
1497       return 1.0;
1498     else
1499       return 0.5 + (pos - seg->middle) / len / 2;
1500   }
1501 }
1502
1503 /*
1504 =item sine_interp(pos, seg)
1505
1506 Calculates sine function interpolation on the fill parameter.
1507
1508 =cut
1509 */
1510 static double
1511 sine_interp(double pos, i_fountain_seg *seg) {
1512   /* I wonder if there's a simple way to smooth the transition for this */
1513   double work = linear_interp(pos, seg);
1514
1515   return (1-cos(work * PI))/2;
1516 }
1517
1518 /*
1519 =item sphereup_interp(pos, seg)
1520
1521 Calculates spherical interpolation on the fill parameter, with the cusp
1522 at the low-end.
1523
1524 =cut
1525 */
1526 static double
1527 sphereup_interp(double pos, i_fountain_seg *seg) {
1528   double work = linear_interp(pos, seg);
1529
1530   return sqrt(1.0 - (1-work) * (1-work));
1531 }
1532
1533 /*
1534 =item spheredown_interp(pos, seg)
1535
1536 Calculates spherical interpolation on the fill parameter, with the cusp
1537 at the high-end.
1538
1539 =cut
1540 */
1541 static double
1542 spheredown_interp(double pos, i_fountain_seg *seg) {
1543   double work = linear_interp(pos, seg);
1544
1545   return 1-sqrt(1.0 - work * work);
1546 }
1547
1548 /*
1549 =item direct_cinterp(out, pos, seg)
1550
1551 Calculates the fountain color based on direct scaling of the channels
1552 of the color channels.
1553
1554 =cut
1555 */
1556 static void
1557 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1558   int ch;
1559   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1560     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1561       + seg->c[1].channel[ch] * pos;
1562   }
1563 }
1564
1565 /*
1566 =item hue_up_cinterp(put, pos, seg)
1567
1568 Calculates the fountain color based on scaling a HSV value.  The hue
1569 increases as the fill parameter increases.
1570
1571 =cut
1572 */
1573 static void
1574 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1575   int ch;
1576   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1577     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1578       + seg->c[1].channel[ch] * pos;
1579   }
1580   i_hsv_to_rgbf(out);
1581 }
1582
1583 /*
1584 =item hue_down_cinterp(put, pos, seg)
1585
1586 Calculates the fountain color based on scaling a HSV value.  The hue
1587 decreases as the fill parameter increases.
1588
1589 =cut
1590 */
1591 static void
1592 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1593   int ch;
1594   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1595     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1596       + seg->c[1].channel[ch] * pos;
1597   }
1598   i_hsv_to_rgbf(out);
1599 }
1600
1601 /*
1602 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1603
1604 Simple grid-based super-sampling.
1605
1606 =cut
1607 */
1608 static int
1609 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
1610   i_fcolor *work = state->ssample_data;
1611   int dx, dy;
1612   int grid = state->parm;
1613   double base = -0.5 + 0.5 / grid;
1614   double step = 1.0 / grid;
1615   int ch, i;
1616   int samp_count = 0;
1617
1618   for (dx = 0; dx < grid; ++dx) {
1619     for (dy = 0; dy < grid; ++dy) {
1620       if (fount_getat(work+samp_count, x + base + step * dx,
1621                       y + base + step * dy, state)) {
1622         ++samp_count;
1623       }
1624     }
1625   }
1626   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1627     out->channel[ch] = 0;
1628     for (i = 0; i < samp_count; ++i) {
1629       out->channel[ch] += work[i].channel[ch];
1630     }
1631     /* we divide by 4 rather than samp_count since if there's only one valid
1632        sample it should be mostly transparent */
1633     out->channel[ch] /= grid * grid;
1634   }
1635   return samp_count;
1636 }
1637
1638 /*
1639 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1640
1641 Random super-sampling.
1642
1643 =cut
1644 */
1645 static int
1646 random_ssample(i_fcolor *out, double x, double y,
1647                struct fount_state *state) {
1648   i_fcolor *work = state->ssample_data;
1649   int i, ch;
1650   int maxsamples = state->parm;
1651   double rand_scale = 1.0 / RAND_MAX;
1652   int samp_count = 0;
1653   for (i = 0; i < maxsamples; ++i) {
1654     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
1655                     y - 0.5 + rand() * rand_scale, state)) {
1656       ++samp_count;
1657     }
1658   }
1659   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1660     out->channel[ch] = 0;
1661     for (i = 0; i < samp_count; ++i) {
1662       out->channel[ch] += work[i].channel[ch];
1663     }
1664     /* we divide by maxsamples rather than samp_count since if there's
1665        only one valid sample it should be mostly transparent */
1666     out->channel[ch] /= maxsamples;
1667   }
1668   return samp_count;
1669 }
1670
1671 /*
1672 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1673
1674 Super-sampling around the circumference of a circle.
1675
1676 I considered saving the sin()/cos() values and transforming step-size
1677 around the circle, but that's inaccurate, though it may not matter
1678 much.
1679
1680 =cut
1681  */
1682 static int
1683 circle_ssample(i_fcolor *out, double x, double y,
1684                struct fount_state *state) {
1685   i_fcolor *work = state->ssample_data;
1686   int i, ch;
1687   int maxsamples = state->parm;
1688   double angle = 2 * PI / maxsamples;
1689   double radius = 0.3; /* semi-random */
1690   int samp_count = 0;
1691   for (i = 0; i < maxsamples; ++i) {
1692     if (fount_getat(work+samp_count, x + radius * cos(angle * i),
1693                     y + radius * sin(angle * i), state)) {
1694       ++samp_count;
1695     }
1696   }
1697   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1698     out->channel[ch] = 0;
1699     for (i = 0; i < samp_count; ++i) {
1700       out->channel[ch] += work[i].channel[ch];
1701     }
1702     /* we divide by maxsamples rather than samp_count since if there's
1703        only one valid sample it should be mostly transparent */
1704     out->channel[ch] /= maxsamples;
1705   }
1706   return samp_count;
1707 }
1708
1709 /*
1710 =item fount_r_none(v)
1711
1712 Implements no repeats.  Simply clamps the fill value.
1713
1714 =cut
1715 */
1716 static double
1717 fount_r_none(double v) {
1718   return v < 0 ? 0 : v > 1 ? 1 : v;
1719 }
1720
1721 /*
1722 =item fount_r_sawtooth(v)
1723
1724 Implements sawtooth repeats.  Clamps negative values and uses fmod()
1725 on others.
1726
1727 =cut
1728 */
1729 static double
1730 fount_r_sawtooth(double v) {
1731   return v < 0 ? 0 : fmod(v, 1.0);
1732 }
1733
1734 /*
1735 =item fount_r_triangle(v)
1736
1737 Implements triangle repeats.  Clamps negative values, uses fmod to get
1738 a range 0 through 2 and then adjusts values > 1.
1739
1740 =cut
1741 */
1742 static double
1743 fount_r_triangle(double v) {
1744   if (v < 0)
1745     return 0;
1746   else {
1747     v = fmod(v, 2.0);
1748     return v > 1.0 ? 2.0 - v : v;
1749   }
1750 }
1751
1752 /*
1753 =item fount_r_saw_both(v)
1754
1755 Implements sawtooth repeats in the both postive and negative directions.
1756
1757 Adjusts the value to be postive and then just uses fmod().
1758
1759 =cut
1760 */
1761 static double
1762 fount_r_saw_both(double v) {
1763   if (v < 0)
1764     v += 1+(int)(-v);
1765   return fmod(v, 1.0);
1766 }
1767
1768 /*
1769 =item fount_r_tri_both(v)
1770
1771 Implements triangle repeats in the both postive and negative directions.
1772
1773 Uses fmod on the absolute value, and then adjusts values > 1.
1774
1775 =cut
1776 */
1777 static double
1778 fount_r_tri_both(double v) {
1779   v = fmod(fabs(v), 2.0);
1780   return v > 1.0 ? 2.0 - v : v;
1781 }
1782
1783 /*
1784 =item fill_fountf(fill, x, y, width, channels, data)
1785
1786 The fill function for fountain fills.
1787
1788 =cut
1789 */
1790 static void
1791 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
1792             i_fcolor *data) {
1793   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
1794   int ch;
1795
1796   while (width--) {
1797     i_fcolor c;
1798     int got_one;
1799     double v;
1800     if (f->state.ssfunc)
1801       got_one = f->state.ssfunc(&c, x, y, &f->state);
1802     else
1803       got_one = fount_getat(&c, x, y, &f->state);
1804
1805     if (got_one) {
1806       if (f->base.combines) {
1807         for (ch = 0; ch < channels; ++ch) {
1808           data->channel[ch] = data->channel[ch] * (1.0 - c.channel[3])
1809             + c.channel[ch] * c.channel[3];
1810         }
1811       }
1812       else
1813         *data = c;
1814     }
1815
1816     ++x;
1817     ++data;
1818   }
1819 }
1820
1821 /*
1822 =item fount_fill_destroy(fill)
1823
1824 =cut
1825 */
1826 static void
1827 fount_fill_destroy(i_fill_t *fill) {
1828   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
1829   fount_finish_state(&f->state);
1830 }
1831
1832 /*
1833 =back
1834
1835 =head1 AUTHOR
1836
1837 Arnar M. Hrafnkelsson <addi@umich.edu>
1838
1839 Tony Cook <tony@develop-help.com> (i_fountain())
1840
1841 =head1 SEE ALSO
1842
1843 Imager(3)
1844
1845 =cut
1846 */