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