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   i_fcolor *work = NULL;
1147   int ch;
1148   i_fountain_seg *my_segs;
1149   i_fill_combine_f combine_func = NULL;
1150   i_fill_combinef_f combinef_func = NULL;
1151
1152   i_get_combine(combine, &combine_func, &combinef_func);
1153   if (combinef_func)
1154     work = mymalloc(sizeof(i_fcolor) * im->xsize);
1155
1156   fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
1157                    super_sample, ssample_param, count, segs);
1158   my_segs = state.segs;
1159
1160   for (y = 0; y < im->ysize; ++y) {
1161     i_glinf(im, 0, im->xsize, y, line);
1162     for (x = 0; x < im->xsize; ++x) {
1163       i_fcolor c;
1164       int got_one;
1165       double v;
1166       if (super_sample == i_fts_none)
1167         got_one = fount_getat(&c, x, y, &state);
1168       else
1169         got_one = state.ssfunc(&c, x, y, &state);
1170       if (got_one) {
1171         if (combine)
1172           work[x] = c;
1173         else
1174           line[x] = c;
1175       }
1176     }
1177     if (combine)
1178       combinef_func(line, work, im->channels, im->xsize);
1179     i_plinf(im, 0, im->xsize, y, line);
1180   }
1181   fount_finish_state(&state);
1182   myfree(line);
1183 }
1184
1185 typedef struct {
1186   i_fill_t base;
1187   struct fount_state state;
1188 } i_fill_fountain_t;
1189
1190 static void
1191 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
1192             i_fcolor *data, i_fcolor *work);
1193 static void
1194 fount_fill_destroy(i_fill_t *fill);
1195
1196 /*
1197 =item i_new_fount(xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1198
1199 Creates a new general fill which fills with a fountain fill.
1200
1201 =cut
1202 */
1203
1204 i_fill_t *
1205 i_new_fill_fount(double xa, double ya, double xb, double yb,
1206                  i_fountain_type type, i_fountain_repeat repeat,
1207                  int combine, int super_sample, double ssample_param,
1208                  int count, i_fountain_seg *segs) {
1209   i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1210
1211   fill->base.fill_with_color = NULL;
1212   fill->base.fill_with_fcolor = fill_fountf;
1213   fill->base.destroy = fount_fill_destroy;
1214   if (combine)
1215     i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1216   else {
1217     fill->base.combine = NULL;
1218     fill->base.combinef = NULL;
1219   }
1220   fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1221                    super_sample, ssample_param, count, segs);
1222
1223   return &fill->base;
1224 }
1225
1226 /*
1227 =back
1228
1229 =head1 INTERNAL FUNCTIONS
1230
1231 =over
1232
1233 =item fount_init_state(...)
1234
1235 Used by both the fountain fill filter and the fountain fill.
1236
1237 =cut
1238 */
1239
1240 static void
1241 fount_init_state(struct fount_state *state, double xa, double ya,
1242                  double xb, double yb, i_fountain_type type,
1243                  i_fountain_repeat repeat, int combine, int super_sample,
1244                  double ssample_param, int count, i_fountain_seg *segs) {
1245   int i, j;
1246   i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count);
1247   /*int have_alpha = im->channels == 2 || im->channels == 4;*/
1248   int ch;
1249
1250   memset(state, 0, sizeof(*state));
1251   /* we keep a local copy that we can adjust for speed */
1252   for (i = 0; i < count; ++i) {
1253     i_fountain_seg *seg = my_segs + i;
1254
1255     *seg = segs[i];
1256     if (seg->type < 0 || type >= i_ft_end)
1257       seg->type = i_ft_linear;
1258     if (seg->color < 0 || seg->color >= i_fc_end)
1259       seg->color = i_fc_direct;
1260     if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1261       /* so we don't have to translate to HSV on each request, do it here */
1262       for (j = 0; j < 2; ++j) {
1263         i_rgb_to_hsvf(seg->c+j);
1264       }
1265       if (seg->color == i_fc_hue_up) {
1266         if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1267           seg->c[1].channel[0] += 1.0;
1268       }
1269       else {
1270         if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1271           seg->c[0].channel[0] += 1.0;
1272       }
1273     }
1274     /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1275            seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1276            seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1277            seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1278            seg->c[1].channel[3], seg->type, seg->color);*/
1279
1280   }
1281
1282   /* initialize each engine */
1283   /* these are so common ... */
1284   state->lA = xb - xa;
1285   state->lB = yb - ya;
1286   state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1287   state->xa = xa;
1288   state->ya = ya;
1289   switch (type) {
1290   default:
1291     type = i_ft_linear; /* make the invalid value valid */
1292   case i_ft_linear:
1293   case i_ft_bilinear:
1294     state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1295     state->mult = 1;
1296     state->mult = 1/linear_fount_f(xb, yb, state);
1297     break;
1298
1299   case i_ft_radial:
1300     state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1301                              + (double)(yb-ya)*(yb-ya));
1302     break;
1303
1304   case i_ft_radial_square:
1305     state->cos = state->lA / state->AB;
1306     state->sin = state->lB / state->AB;
1307     state->mult = 1.0 / state->AB;
1308     break;
1309
1310   case i_ft_revolution:
1311     state->theta = atan2(yb-ya, xb-xa);
1312     state->mult = 1.0 / (PI * 2);
1313     break;
1314
1315   case i_ft_conical:
1316     state->theta = atan2(yb-ya, xb-xa);
1317     state->mult = 1.0 / PI;
1318     break;
1319   }
1320   state->ffunc = fount_funcs[type];
1321   if (super_sample < 0
1322       || super_sample >= (sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
1323     super_sample = 0;
1324   }
1325   state->ssample_data = NULL;
1326   switch (super_sample) {
1327   case i_fts_grid:
1328     ssample_param = floor(0.5 + sqrt(ssample_param));
1329     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param);
1330     break;
1331
1332   case i_fts_random:
1333   case i_fts_circle:
1334     ssample_param = floor(0.5+ssample_param);
1335     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
1336     break;
1337   }
1338   state->parm = ssample_param;
1339   state->ssfunc = fount_ssamples[super_sample];
1340   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1341     repeat = 0;
1342   state->rpfunc = fount_repeats[repeat];
1343   state->segs = my_segs;
1344   state->count = count;
1345 }
1346
1347 static void
1348 fount_finish_state(struct fount_state *state) {
1349   if (state->ssample_data)
1350     myfree(state->ssample_data);
1351   myfree(state->segs);
1352 }
1353
1354
1355 /*
1356 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
1357
1358 Evaluates the fountain fill at the given point.
1359
1360 This is called by both the non-super-sampling and super-sampling code.
1361
1362 You might think that it would make sense to sample the fill parameter
1363 instead, and combine those, but this breaks badly.
1364
1365 =cut
1366 */
1367
1368 static int
1369 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
1370   double v = (state->rpfunc)((state->ffunc)(x, y, state));
1371   int i;
1372
1373   i = 0;
1374   while (i < state->count
1375          && (v < state->segs[i].start || v > state->segs[i].end)) {
1376     ++i;
1377   }
1378   if (i < state->count) {
1379     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
1380     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
1381     return 1;
1382   }
1383   else
1384     return 0;
1385 }
1386
1387 /*
1388 =item linear_fount_f(x, y, state)
1389
1390 Calculate the fill parameter for a linear fountain fill.
1391
1392 Uses the point to line distance function, with some precalculation
1393 done in i_fountain().
1394
1395 =cut
1396 */
1397 static double
1398 linear_fount_f(double x, double y, struct fount_state *state) {
1399   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
1400 }
1401
1402 /*
1403 =item bilinear_fount_f(x, y, state)
1404
1405 Calculate the fill parameter for a bi-linear fountain fill.
1406
1407 =cut
1408 */
1409 static double
1410 bilinear_fount_f(double x, double y, struct fount_state *state) {
1411   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
1412 }
1413
1414 /*
1415 =item radial_fount_f(x, y, state)
1416
1417 Calculate the fill parameter for a radial fountain fill.
1418
1419 Simply uses the distance function.
1420
1421 =cut
1422  */
1423 static double
1424 radial_fount_f(double x, double y, struct fount_state *state) {
1425   return sqrt((double)(state->xa-x)*(state->xa-x)
1426               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
1427 }
1428
1429 /*
1430 =item square_fount_f(x, y, state)
1431
1432 Calculate the fill parameter for a square fountain fill.
1433
1434 Works by rotating the reference co-ordinate around the centre of the
1435 square.
1436
1437 =cut
1438 */
1439 static double
1440 square_fount_f(double x, double y, struct fount_state *state) {
1441   int xc, yc; /* centred on A */
1442   double xt, yt; /* rotated by theta */
1443   xc = x - state->xa;
1444   yc = y - state->ya;
1445   xt = fabs(xc * state->cos + yc * state->sin);
1446   yt = fabs(-xc * state->sin + yc * state->cos);
1447   return (xt > yt ? xt : yt) * state->mult;
1448 }
1449
1450 /*
1451 =item revolution_fount_f(x, y, state)
1452
1453 Calculates the fill parameter for the revolution fountain fill.
1454
1455 =cut
1456 */
1457 static double
1458 revolution_fount_f(double x, double y, struct fount_state *state) {
1459   double angle = atan2(y - state->ya, x - state->xa);
1460
1461   angle -= state->theta;
1462   if (angle < 0) {
1463     angle = fmod(angle+ PI * 4, PI*2);
1464   }
1465
1466   return angle * state->mult;
1467 }
1468
1469 /*
1470 =item conical_fount_f(x, y, state)
1471
1472 Calculates the fill parameter for the conical fountain fill.
1473
1474 =cut
1475 */
1476 static double
1477 conical_fount_f(double x, double y, struct fount_state *state) {
1478   double angle = atan2(y - state->ya, x - state->xa);
1479
1480   angle -= state->theta;
1481   if (angle < -PI)
1482     angle += PI * 2;
1483   else if (angle > PI)
1484     angle -= PI * 2;
1485
1486   return fabs(angle) * state->mult;
1487 }
1488
1489 /*
1490 =item linear_interp(pos, seg)
1491
1492 Calculates linear interpolation on the fill parameter.  Breaks the
1493 segment into 2 regions based in the I<middle> value.
1494
1495 =cut
1496 */
1497 static double
1498 linear_interp(double pos, i_fountain_seg *seg) {
1499   if (pos < seg->middle) {
1500     double len = seg->middle - seg->start;
1501     if (len < EPSILON)
1502       return 0.0;
1503     else
1504       return (pos - seg->start) / len / 2;
1505   }
1506   else {
1507     double len = seg->end - seg->middle;
1508     if (len < EPSILON)
1509       return 1.0;
1510     else
1511       return 0.5 + (pos - seg->middle) / len / 2;
1512   }
1513 }
1514
1515 /*
1516 =item sine_interp(pos, seg)
1517
1518 Calculates sine function interpolation on the fill parameter.
1519
1520 =cut
1521 */
1522 static double
1523 sine_interp(double pos, i_fountain_seg *seg) {
1524   /* I wonder if there's a simple way to smooth the transition for this */
1525   double work = linear_interp(pos, seg);
1526
1527   return (1-cos(work * PI))/2;
1528 }
1529
1530 /*
1531 =item sphereup_interp(pos, seg)
1532
1533 Calculates spherical interpolation on the fill parameter, with the cusp
1534 at the low-end.
1535
1536 =cut
1537 */
1538 static double
1539 sphereup_interp(double pos, i_fountain_seg *seg) {
1540   double work = linear_interp(pos, seg);
1541
1542   return sqrt(1.0 - (1-work) * (1-work));
1543 }
1544
1545 /*
1546 =item spheredown_interp(pos, seg)
1547
1548 Calculates spherical interpolation on the fill parameter, with the cusp
1549 at the high-end.
1550
1551 =cut
1552 */
1553 static double
1554 spheredown_interp(double pos, i_fountain_seg *seg) {
1555   double work = linear_interp(pos, seg);
1556
1557   return 1-sqrt(1.0 - work * work);
1558 }
1559
1560 /*
1561 =item direct_cinterp(out, pos, seg)
1562
1563 Calculates the fountain color based on direct scaling of the channels
1564 of the color channels.
1565
1566 =cut
1567 */
1568 static void
1569 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1570   int ch;
1571   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1572     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1573       + seg->c[1].channel[ch] * pos;
1574   }
1575 }
1576
1577 /*
1578 =item hue_up_cinterp(put, pos, seg)
1579
1580 Calculates the fountain color based on scaling a HSV value.  The hue
1581 increases as the fill parameter increases.
1582
1583 =cut
1584 */
1585 static void
1586 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1587   int ch;
1588   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1589     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1590       + seg->c[1].channel[ch] * pos;
1591   }
1592   i_hsv_to_rgbf(out);
1593 }
1594
1595 /*
1596 =item hue_down_cinterp(put, pos, seg)
1597
1598 Calculates the fountain color based on scaling a HSV value.  The hue
1599 decreases as the fill parameter increases.
1600
1601 =cut
1602 */
1603 static void
1604 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1605   int ch;
1606   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1607     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1608       + seg->c[1].channel[ch] * pos;
1609   }
1610   i_hsv_to_rgbf(out);
1611 }
1612
1613 /*
1614 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1615
1616 Simple grid-based super-sampling.
1617
1618 =cut
1619 */
1620 static int
1621 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
1622   i_fcolor *work = state->ssample_data;
1623   int dx, dy;
1624   int grid = state->parm;
1625   double base = -0.5 + 0.5 / grid;
1626   double step = 1.0 / grid;
1627   int ch, i;
1628   int samp_count = 0;
1629
1630   for (dx = 0; dx < grid; ++dx) {
1631     for (dy = 0; dy < grid; ++dy) {
1632       if (fount_getat(work+samp_count, x + base + step * dx,
1633                       y + base + step * dy, state)) {
1634         ++samp_count;
1635       }
1636     }
1637   }
1638   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1639     out->channel[ch] = 0;
1640     for (i = 0; i < samp_count; ++i) {
1641       out->channel[ch] += work[i].channel[ch];
1642     }
1643     /* we divide by 4 rather than samp_count since if there's only one valid
1644        sample it should be mostly transparent */
1645     out->channel[ch] /= grid * grid;
1646   }
1647   return samp_count;
1648 }
1649
1650 /*
1651 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1652
1653 Random super-sampling.
1654
1655 =cut
1656 */
1657 static int
1658 random_ssample(i_fcolor *out, double x, double y,
1659                struct fount_state *state) {
1660   i_fcolor *work = state->ssample_data;
1661   int i, ch;
1662   int maxsamples = state->parm;
1663   double rand_scale = 1.0 / RAND_MAX;
1664   int samp_count = 0;
1665   for (i = 0; i < maxsamples; ++i) {
1666     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
1667                     y - 0.5 + rand() * rand_scale, state)) {
1668       ++samp_count;
1669     }
1670   }
1671   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1672     out->channel[ch] = 0;
1673     for (i = 0; i < samp_count; ++i) {
1674       out->channel[ch] += work[i].channel[ch];
1675     }
1676     /* we divide by maxsamples rather than samp_count since if there's
1677        only one valid sample it should be mostly transparent */
1678     out->channel[ch] /= maxsamples;
1679   }
1680   return samp_count;
1681 }
1682
1683 /*
1684 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1685
1686 Super-sampling around the circumference of a circle.
1687
1688 I considered saving the sin()/cos() values and transforming step-size
1689 around the circle, but that's inaccurate, though it may not matter
1690 much.
1691
1692 =cut
1693  */
1694 static int
1695 circle_ssample(i_fcolor *out, double x, double y,
1696                struct fount_state *state) {
1697   i_fcolor *work = state->ssample_data;
1698   int i, ch;
1699   int maxsamples = state->parm;
1700   double angle = 2 * PI / maxsamples;
1701   double radius = 0.3; /* semi-random */
1702   int samp_count = 0;
1703   for (i = 0; i < maxsamples; ++i) {
1704     if (fount_getat(work+samp_count, x + radius * cos(angle * i),
1705                     y + radius * sin(angle * i), state)) {
1706       ++samp_count;
1707     }
1708   }
1709   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1710     out->channel[ch] = 0;
1711     for (i = 0; i < samp_count; ++i) {
1712       out->channel[ch] += work[i].channel[ch];
1713     }
1714     /* we divide by maxsamples rather than samp_count since if there's
1715        only one valid sample it should be mostly transparent */
1716     out->channel[ch] /= maxsamples;
1717   }
1718   return samp_count;
1719 }
1720
1721 /*
1722 =item fount_r_none(v)
1723
1724 Implements no repeats.  Simply clamps the fill value.
1725
1726 =cut
1727 */
1728 static double
1729 fount_r_none(double v) {
1730   return v < 0 ? 0 : v > 1 ? 1 : v;
1731 }
1732
1733 /*
1734 =item fount_r_sawtooth(v)
1735
1736 Implements sawtooth repeats.  Clamps negative values and uses fmod()
1737 on others.
1738
1739 =cut
1740 */
1741 static double
1742 fount_r_sawtooth(double v) {
1743   return v < 0 ? 0 : fmod(v, 1.0);
1744 }
1745
1746 /*
1747 =item fount_r_triangle(v)
1748
1749 Implements triangle repeats.  Clamps negative values, uses fmod to get
1750 a range 0 through 2 and then adjusts values > 1.
1751
1752 =cut
1753 */
1754 static double
1755 fount_r_triangle(double v) {
1756   if (v < 0)
1757     return 0;
1758   else {
1759     v = fmod(v, 2.0);
1760     return v > 1.0 ? 2.0 - v : v;
1761   }
1762 }
1763
1764 /*
1765 =item fount_r_saw_both(v)
1766
1767 Implements sawtooth repeats in the both postive and negative directions.
1768
1769 Adjusts the value to be postive and then just uses fmod().
1770
1771 =cut
1772 */
1773 static double
1774 fount_r_saw_both(double v) {
1775   if (v < 0)
1776     v += 1+(int)(-v);
1777   return fmod(v, 1.0);
1778 }
1779
1780 /*
1781 =item fount_r_tri_both(v)
1782
1783 Implements triangle repeats in the both postive and negative directions.
1784
1785 Uses fmod on the absolute value, and then adjusts values > 1.
1786
1787 =cut
1788 */
1789 static double
1790 fount_r_tri_both(double v) {
1791   v = fmod(fabs(v), 2.0);
1792   return v > 1.0 ? 2.0 - v : v;
1793 }
1794
1795 /*
1796 =item fill_fountf(fill, x, y, width, channels, data)
1797
1798 The fill function for fountain fills.
1799
1800 =cut
1801 */
1802 static void
1803 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
1804             i_fcolor *data, i_fcolor *work) {
1805   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
1806
1807   if (fill->combinef) {
1808     i_fcolor *wstart = work;
1809     int count = width;
1810
1811     while (width--) {
1812       i_fcolor c;
1813       int got_one;
1814       double v;
1815       if (f->state.ssfunc)
1816         got_one = f->state.ssfunc(&c, x, y, &f->state);
1817       else
1818         got_one = fount_getat(&c, x, y, &f->state);
1819
1820       *work++ = c;
1821
1822       ++x;
1823     }
1824     (fill->combinef)(data, wstart, channels, count);
1825   }
1826   else {
1827     while (width--) {
1828       i_fcolor c;
1829       int got_one;
1830       double v;
1831       if (f->state.ssfunc)
1832         got_one = f->state.ssfunc(&c, x, y, &f->state);
1833       else
1834         got_one = fount_getat(&c, x, y, &f->state);
1835
1836       *data++ = c;
1837
1838       ++x;
1839     }
1840   }
1841 }
1842
1843 /*
1844 =item fount_fill_destroy(fill)
1845
1846 =cut
1847 */
1848 static void
1849 fount_fill_destroy(i_fill_t *fill) {
1850   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
1851   fount_finish_state(&f->state);
1852 }
1853
1854 /*
1855 =back
1856
1857 =head1 AUTHOR
1858
1859 Arnar M. Hrafnkelsson <addi@umich.edu>
1860
1861 Tony Cook <tony@develop-help.com> (i_fountain())
1862
1863 =head1 SEE ALSO
1864
1865 Imager(3)
1866
1867 =cut
1868 */