filters.im

   1 #include "imager.h"
   2 #include "imageri.h"
   3 #include <stdlib.h>
   4 #include <math.h>
   5
   6
   7 /*
   8 =head1 NAME
   9
  10 filters.im - implements filters that operate on images
  11
  12 =head1 SYNOPSIS
  13
  14
  15   i_contrast(im, 0.8);
  16   i_hardinvert(im);
  17   i_hardinvertall(im);
  18   i_unsharp_mask(im, 2.0, 1.0);
  19   ... and more
  20
  21 =head1 DESCRIPTION
  22
  23 filters.c implements basic filters for Imager.  These filters
  24 should be accessible from the filter interface as defined in
  25 the pod for Imager.
  26
  27 =head1 FUNCTION REFERENCE
  28
  29 Some of these functions are internal.
  30
  31 =over
  32
  33 =cut
  34 */
  35
  36
  37
  38
  39 /*
  40 =item saturate(in)
  41
  42 Clamps the input value between 0 and 255. (internal)
  43
  44   in - input integer
  45
  46 =cut
  47 */
  48
  49 static
  50 unsigned char
  51 saturate(int in) {
  52   if (in>255) { return 255; }
  53   else if (in>0) return in;
  54   return 0;
  55 }
  56
  57
  58
  59 /*
  60 =item i_contrast(im, intensity)
  61
  62 Scales the pixel values by the amount specified.
  63
  64   im        - image object
  65   intensity - scalefactor
  66
  67 =cut
  68 */
  69
  70 void
  71 i_contrast(i_img *im, float intensity) {
  72   i_img_dim x, y;
  73   unsigned char ch;
  74   unsigned int new_color;
  75   i_color rcolor;
  76
  77   mm_log((1,"i_contrast(im %p, intensity %f)\n", im, intensity));
  78
  79   if(intensity < 0) return;
  80
  81   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
  82     i_gpix(im, x, y, &rcolor);
  83
  84     for(ch = 0; ch < im->channels; ch++) {
  85       new_color = (unsigned int) rcolor.channel[ch];
  86       new_color *= intensity;
  87
  88       if(new_color > 255) {
  89         new_color = 255;
  90       }
  91       rcolor.channel[ch] = (unsigned char) new_color;
  92     }
  93     i_ppix(im, x, y, &rcolor);
  94   }
  95 }
  96
  97
  98 static int
  99 s_hardinvert_low(i_img *im, int all) {
 100   i_img_dim x, y;
 101   int ch;
 102   int invert_channels = all ? im->channels : i_img_color_channels(im);
 103
 104   mm_log((1,"i_hardinvert(im %p)\n", im));
 105
 106 #code im->bits <= 8
 107   IM_COLOR *row, *entry;
 108
 109   /* always rooms to allocate a single line of i_color */
 110   row = mymalloc(sizeof(IM_COLOR) * im->xsize); /* checked 17feb2005 tonyc */
 111
 112   for(y = 0; y < im->ysize; y++) {
 113     IM_GLIN(im, 0, im->xsize, y, row);
 114     entry = row;
 115     for(x = 0; x < im->xsize; x++) {
 116       for(ch = 0; ch < invert_channels; ch++) {
 117         entry->channel[ch] = IM_SAMPLE_MAX - entry->channel[ch];
 118       }
 119       ++entry;
 120     }
 121     IM_PLIN(im, 0, im->xsize, y, row);
 122   }
 123   myfree(row);
 124 #/code
 125
 126   return 1;
 127 }
 128
 129 /*
 130 =item i_hardinvert(im)
 131
 132 Inverts the color channels of the input image.
 133
 134   im        - image object
 135
 136 =cut
 137 */
 138
 139 void
 140 i_hardinvert(i_img *im) {
 141   s_hardinvert_low(im, 0);
 142 }
 143
 144 /*
 145 =item i_hardinvertall(im)
 146
 147 Inverts all channels of the input image.
 148
 149   im        - image object
 150
 151 =cut
 152 */
 153
 154 void
 155 i_hardinvertall(i_img *im) {
 156   s_hardinvert_low(im, 1);
 157 }
 158
 159 /*
 160 =item i_noise(im, amount, type)
 161
 162 Adjusts the sample values randomly by the amount specified.
 163
 164 If type is 0, adjust all channels in a pixel by the same (random)
 165 amount amount, if non-zero adjust each sample independently.
 166
 167   im     - image object
 168   amount - deviation in pixel values
 169   type   - noise individual for each channel if true
 170
 171 =cut
 172 */
 173
 174 #ifdef WIN32
 175 /* random() is non-ASCII, even if it is better than rand() */
 176 #define random() rand()
 177 #endif
 178
 179 void
 180 i_noise(i_img *im, float amount, unsigned char type) {
 181   i_img_dim x, y;
 182   unsigned char ch;
 183   int new_color;
 184   float damount = amount * 2;
 185   i_color rcolor;
 186   int color_inc = 0;
 187
 188   mm_log((1,"i_noise(im %p, intensity %.2f\n", im, amount));
 189
 190   if(amount < 0) return;
 191
 192   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 193     i_gpix(im, x, y, &rcolor);
 194
 195     if(type == 0) {
 196       color_inc = (amount - (damount * ((float)random() / RAND_MAX)));
 197     }
 198
 199     for(ch = 0; ch < im->channels; ch++) {
 200       new_color = (int) rcolor.channel[ch];
 201
 202       if(type != 0) {
 203         new_color += (amount - (damount * ((float)random() / RAND_MAX)));
 204       } else {
 205         new_color += color_inc;
 206       }
 207
 208       if(new_color < 0) {
 209         new_color = 0;
 210       }
 211       if(new_color > 255) {
 212         new_color = 255;
 213       }
 214
 215       rcolor.channel[ch] = (unsigned char) new_color;
 216     }
 217
 218     i_ppix(im, x, y, &rcolor);
 219   }
 220 }
 221
 222 /*
 223 =item i_bumpmap(im, bump, channel, light_x, light_y, st)
 224
 225 Makes a bumpmap on image im using the bump image as the elevation map.
 226
 227   im      - target image
 228   bump    - image that contains the elevation info
 229   channel - to take the elevation information from
 230   light_x - x coordinate of light source
 231   light_y - y coordinate of light source
 232   st      - length of shadow
 233
 234 =cut
 235 */
 236
 237 void
 238 i_bumpmap(i_img *im, i_img *bump, int channel, i_img_dim light_x, i_img_dim light_y, i_img_dim st) {
 239   i_img_dim x, y;
 240   int ch;
 241   i_img_dim mx, my;
 242   i_color x1_color, y1_color, x2_color, y2_color, dst_color;
 243   double nX, nY;
 244   double tX, tY, tZ;
 245   double aX, aY, aL;
 246   double fZ;
 247   unsigned char px1, px2, py1, py2;
 248
 249   i_img new_im;
 250
 251   mm_log((1, "i_bumpmap(im %p, add_im %p, channel %d, light(" i_DFp "), st %" i_DF ")\n",
 252           im, bump, channel, i_DFcp(light_x, light_y), i_DFc(st)));
 253
 254
 255   if(channel >= bump->channels) {
 256     mm_log((1, "i_bumpmap: channel = %d while bump image only has %d channels\n", channel, bump->channels));
 257     return;
 258   }
 259
 260   mx = (bump->xsize <= im->xsize) ? bump->xsize : im->xsize;
 261   my = (bump->ysize <= im->ysize) ? bump->ysize : im->ysize;
 262
 263   i_img_empty_ch(&new_im, im->xsize, im->ysize, im->channels);
 264
 265   aX = (light_x > (mx >> 1)) ? light_x : mx - light_x;
 266   aY = (light_y > (my >> 1)) ? light_y : my - light_y;
 267
 268   aL = sqrt((aX * aX) + (aY * aY));
 269
 270   for(y = 1; y < my - 1; y++) {
 271     for(x = 1; x < mx - 1; x++) {
 272       i_gpix(bump, x + st, y, &x1_color);
 273       i_gpix(bump, x, y + st, &y1_color);
 274       i_gpix(bump, x - st, y, &x2_color);
 275       i_gpix(bump, x, y - st, &y2_color);
 276
 277       i_gpix(im, x, y, &dst_color);
 278
 279       px1 = x1_color.channel[channel];
 280       py1 = y1_color.channel[channel];
 281       px2 = x2_color.channel[channel];
 282       py2 = y2_color.channel[channel];
 283
 284       nX = px1 - px2;
 285       nY = py1 - py2;
 286
 287       nX += 128;
 288       nY += 128;
 289
 290       fZ = (sqrt((nX * nX) + (nY * nY)) / aL);
 291
 292       tX = i_abs(x - light_x) / aL;
 293       tY = i_abs(y - light_y) / aL;
 294
 295       tZ = 1 - (sqrt((tX * tX) + (tY * tY)) * fZ);
 296
 297       if(tZ < 0) tZ = 0;
 298       if(tZ > 2) tZ = 2;
 299
 300       for(ch = 0; ch < im->channels; ch++)
 301         dst_color.channel[ch] = (unsigned char) (float)(dst_color.channel[ch] * tZ);
 302
 303       i_ppix(&new_im, x, y, &dst_color);
 304     }
 305   }
 306
 307   i_copyto(im, &new_im, 0, 0, im->xsize, im->ysize, 0, 0);
 308
 309   i_img_exorcise(&new_im);
 310 }
 311
 312
 313
 314
 315 typedef struct {
 316   double x,y,z;
 317 } fvec;
 318
 319
 320 static
 321 float
 322 dotp(fvec *a, fvec *b) {
 323   return a->x*b->x+a->y*b->y+a->z*b->z;
 324 }
 325
 326 static
 327 void
 328 normalize(fvec *a) {
 329   double d = sqrt(dotp(a,a));
 330   a->x /= d;
 331   a->y /= d;
 332   a->z /= d;
 333 }
 334
 335
 336 /*
 337   positive directions:
 338
 339   x - right,
 340   y - down
 341   z - out of the plane
 342
 343   I = Ia + Ip*( cd*Scol(N.L) + cs*(R.V)^n )
 344
 345   Here, the variables are:
 346
 347   * Ia   - ambient colour
 348   * Ip   - intensity of the point light source
 349   * cd   - diffuse coefficient
 350   * Scol - surface colour
 351   * cs   - specular coefficient
 352   * n    - objects shinyness
 353   * N    - normal vector
 354   * L    - lighting vector
 355   * R    - reflection vector
 356   * V    - vision vector
 357
 358   static void fvec_dump(fvec *x) {
 359     printf("(%.2f %.2f %.2f)", x->x, x->y, x->z);
 360   }
 361 */
 362
 363 /* XXX: Should these return a code for success? */
 364
 365
 366
 367
 368 /*
 369 =item i_bumpmap_complex(im, bump, channel, tx, ty, Lx, Ly, Lz, Ip, cd, cs, n, Ia, Il, Is)
 370
 371 Makes a bumpmap on image im using the bump image as the elevation map.
 372
 373   im      - target image
 374   bump    - image that contains the elevation info
 375   channel - to take the elevation information from
 376   tx      - shift in x direction of where to start applying bumpmap
 377   ty      - shift in y direction of where to start applying bumpmap
 378   Lx      - x position/direction of light
 379   Ly      - y position/direction of light
 380   Lz      - z position/direction of light
 381   Ip      - light intensity
 382   cd      - diffuse coefficient
 383   cs      - specular coefficient
 384   n       - surface shinyness
 385   Ia      - ambient colour
 386   Il      - light colour
 387   Is      - specular colour
 388
 389 if z<0 then the L is taken to be the direction the light is shining in.  Otherwise
 390 the L is taken to be the position of the Light, Relative to the image.
 391
 392 =cut
 393 */
 394
 395
 396 void
 397 i_bumpmap_complex(i_img *im,
 398                   i_img *bump,
 399                   int channel,
 400                   i_img_dim tx,
 401                   i_img_dim ty,
 402                   double Lx,
 403                   double Ly,
 404                   double Lz,
 405                   float cd,
 406                   float cs,
 407                   float n,
 408                   i_color *Ia,
 409                   i_color *Il,
 410                   i_color *Is) {
 411   i_img new_im;
 412
 413   int inflight;
 414   i_img_dim x, y;
 415   int ch;
 416   i_img_dim mx, Mx, my, My;
 417
 418   float cdc[MAXCHANNELS];
 419   float csc[MAXCHANNELS];
 420
 421   i_color x1_color, y1_color, x2_color, y2_color;
 422
 423   i_color Scol;   /* Surface colour       */
 424
 425   fvec L;         /* Light vector */
 426   fvec N;         /* surface normal       */
 427   fvec R;         /* Reflection vector    */
 428   fvec V;         /* Vision vector        */
 429
 430   mm_log((1, "i_bumpmap_complex(im %p, bump %p, channel %d, t(" i_DFp
 431           "), Lx %.2f, Ly %.2f, Lz %.2f, cd %.2f, cs %.2f, n %.2f, Ia %p, Il %p, Is %p)\n",
 432           im, bump, channel, i_DFcp(tx, ty), Lx, Ly, Lz, cd, cs, n, Ia, Il, Is));
 433
 434   if (channel >= bump->channels) {
 435     mm_log((1, "i_bumpmap_complex: channel = %d while bump image only has %d channels\n", channel, bump->channels));
 436     return;
 437   }
 438
 439   for(ch=0; ch<im->channels; ch++) {
 440     cdc[ch] = (float)Il->channel[ch]*cd/255.f;
 441     csc[ch] = (float)Is->channel[ch]*cs/255.f;
 442   }
 443
 444   mx = 1;
 445   my = 1;
 446   Mx = bump->xsize-1;
 447   My = bump->ysize-1;
 448
 449   V.x = 0;
 450   V.y = 0;
 451   V.z = 1;
 452
 453   if (Lz < 0) { /* Light specifies a direction vector, reverse it to get the vector from surface to light */
 454     L.x = -Lx;
 455     L.y = -Ly;
 456     L.z = -Lz;
 457     normalize(&L);
 458   } else {      /* Light is the position of the light source */
 459     inflight = 0;
 460     L.x = -0.2;
 461     L.y = -0.4;
 462     L.z =  1;
 463     normalize(&L);
 464   }
 465
 466   i_img_empty_ch(&new_im, im->xsize, im->ysize, im->channels);
 467
 468   for(y = 0; y < im->ysize; y++) {
 469     for(x = 0; x < im->xsize; x++) {
 470       double dp1, dp2;
 471       double dx = 0, dy = 0;
 472
 473       /* Calculate surface normal */
 474       if (mx<x && x<Mx && my<y && y<My) {
 475         i_gpix(bump, x + 1, y,     &x1_color);
 476         i_gpix(bump, x - 1, y,     &x2_color);
 477         i_gpix(bump, x,     y + 1, &y1_color);
 478         i_gpix(bump, x,     y - 1, &y2_color);
 479         dx = x2_color.channel[channel] - x1_color.channel[channel];
 480         dy = y2_color.channel[channel] - y1_color.channel[channel];
 481       } else {
 482         dx = 0;
 483         dy = 0;
 484       }
 485       N.x = -dx * 0.015;
 486       N.y = -dy * 0.015;
 487       N.z = 1;
 488       normalize(&N);
 489
 490       /* Calculate Light vector if needed */
 491       if (Lz>=0) {
 492         L.x = Lx - x;
 493         L.y = Ly - y;
 494         L.z = Lz;
 495         normalize(&L);
 496       }
 497
 498       dp1 = dotp(&L,&N);
 499       R.x = -L.x + 2*dp1*N.x;
 500       R.y = -L.y + 2*dp1*N.y;
 501       R.z = -L.z + 2*dp1*N.z;
 502
 503       dp2 = dotp(&R,&V);
 504
 505       dp1 = dp1<0 ?0 : dp1;
 506       dp2 = pow(dp2<0 ?0 : dp2,n);
 507
 508       i_gpix(im, x, y, &Scol);
 509
 510       for(ch = 0; ch < im->channels; ch++)
 511         Scol.channel[ch] =
 512           saturate( Ia->channel[ch] + cdc[ch]*Scol.channel[ch]*dp1 + csc[ch]*dp2 );
 513
 514       i_ppix(&new_im, x, y, &Scol);
 515     }
 516   }
 517
 518   i_copyto(im, &new_im, 0, 0, im->xsize, im->ysize, 0, 0);
 519   i_img_exorcise(&new_im);
 520 }
 521
 522
 523 /*
 524 =item i_postlevels(im, levels)
 525
 526 Quantizes Images to fewer levels.
 527
 528   im      - target image
 529   levels  - number of levels
 530
 531 =cut
 532 */
 533
 534 void
 535 i_postlevels(i_img *im, int levels) {
 536   i_img_dim x, y;
 537   int ch;
 538   float pv;
 539   int rv;
 540   float av;
 541
 542   i_color rcolor;
 543
 544   rv = (int) ((float)(256 / levels));
 545   av = (float)levels;
 546
 547   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 548     i_gpix(im, x, y, &rcolor);
 549
 550     for(ch = 0; ch < im->channels; ch++) {
 551       pv = (((float)rcolor.channel[ch] / 255)) * av;
 552       pv = (int) ((int)pv * rv);
 553
 554       if(pv < 0) pv = 0;
 555       else if(pv > 255) pv = 255;
 556
 557       rcolor.channel[ch] = (unsigned char) pv;
 558     }
 559     i_ppix(im, x, y, &rcolor);
 560   }
 561 }
 562
 563
 564 /*
 565 =item i_mosaic(im, size)
 566
 567 Makes an image looks like a mosaic with tilesize of size
 568
 569   im      - target image
 570   size    - size of tiles
 571
 572 =cut
 573 */
 574
 575 void
 576 i_mosaic(i_img *im, i_img_dim size) {
 577   i_img_dim x, y;
 578   int ch, z;
 579   i_img_dim lx, ly;
 580   long sqrsize;
 581
 582   i_color rcolor;
 583   long col[256];
 584
 585   sqrsize = size * size;
 586
 587   for(y = 0; y < im->ysize; y += size) for(x = 0; x < im->xsize; x += size) {
 588     for(z = 0; z < 256; z++) col[z] = 0;
 589
 590     for(lx = 0; lx < size; lx++) {
 591       for(ly = 0; ly < size; ly++) {
 592         i_gpix(im, (x + lx), (y + ly), &rcolor);
 593
 594         for(ch = 0; ch < im->channels; ch++) {
 595           col[ch] += rcolor.channel[ch];
 596         }
 597       }
 598     }
 599
 600     for(ch = 0; ch < im->channels; ch++)
 601       rcolor.channel[ch] = (int) ((float)col[ch] / sqrsize);
 602
 603
 604     for(lx = 0; lx < size; lx++)
 605       for(ly = 0; ly < size; ly++)
 606       i_ppix(im, (x + lx), (y + ly), &rcolor);
 607
 608   }
 609 }
 610
 611
 612 /*
 613 =item i_watermark(im, wmark, tx, ty, pixdiff)
 614
 615 Applies a watermark to the target image
 616
 617   im      - target image
 618   wmark   - watermark image
 619   tx      - x coordinate of where watermark should be applied
 620   ty      - y coordinate of where watermark should be applied
 621   pixdiff - the magnitude of the watermark, controls how visible it is
 622
 623 =cut
 624 */
 625
 626 void
 627 i_watermark(i_img *im, i_img *wmark, i_img_dim tx, i_img_dim ty, int pixdiff) {
 628   i_img_dim vx, vy;
 629   int ch;
 630   i_color val, wval;
 631
 632         i_img_dim mx = wmark->xsize;
 633         i_img_dim my = wmark->ysize;
 634
 635   for(vx=0;vx<mx;vx++) for(vy=0;vy<my;vy++) {
 636
 637     i_gpix(im,    tx+vx, ty+vy,&val );
 638     i_gpix(wmark, vx,    vy,   &wval);
 639
 640     for(ch=0;ch<im->channels;ch++)
 641       val.channel[ch] = saturate( val.channel[ch] + (pixdiff* (wval.channel[0]-128) )/128 );
 642
 643     i_ppix(im,tx+vx,ty+vy,&val);
 644   }
 645 }
 646
 647
 648 /*
 649 =item i_autolevels(im, lsat, usat, skew)
 650
 651 Scales and translates each color such that it fills the range completely.
 652 Skew is not implemented yet - purpose is to control the color skew that can
 653 occur when changing the contrast.
 654
 655   im   - target image
 656   lsat - fraction of pixels that will be truncated at the lower end of the spectrum
 657   usat - fraction of pixels that will be truncated at the higher end of the spectrum
 658   skew - not used yet
 659
 660 =cut
 661 */
 662
 663 void
 664 i_autolevels(i_img *im, float lsat, float usat, float skew) {
 665   i_color val;
 666   i_img_dim i, x, y, rhist[256], ghist[256], bhist[256];
 667   i_img_dim rsum, rmin, rmax;
 668   i_img_dim gsum, gmin, gmax;
 669   i_img_dim bsum, bmin, bmax;
 670   i_img_dim rcl, rcu, gcl, gcu, bcl, bcu;
 671
 672   mm_log((1,"i_autolevels(im %p, lsat %f,usat %f,skew %f)\n", im, lsat,usat,skew));
 673
 674   rsum=gsum=bsum=0;
 675   for(i=0;i<256;i++) rhist[i]=ghist[i]=bhist[i] = 0;
 676   /* create histogram for each channel */
 677   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 678     i_gpix(im, x, y, &val);
 679     rhist[val.channel[0]]++;
 680     ghist[val.channel[1]]++;
 681     bhist[val.channel[2]]++;
 682   }
 683
 684   for(i=0;i<256;i++) {
 685     rsum+=rhist[i];
 686     gsum+=ghist[i];
 687     bsum+=bhist[i];
 688   }
 689
 690   rmin = gmin = bmin = 0;
 691   rmax = gmax = bmax = 255;
 692
 693   rcu = rcl = gcu = gcl = bcu = bcl = 0;
 694
 695   for(i=0; i<256; i++) {
 696     rcl += rhist[i];     if ( (rcl<rsum*lsat) ) rmin=i;
 697     rcu += rhist[255-i]; if ( (rcu<rsum*usat) ) rmax=255-i;
 698
 699     gcl += ghist[i];     if ( (gcl<gsum*lsat) ) gmin=i;
 700     gcu += ghist[255-i]; if ( (gcu<gsum*usat) ) gmax=255-i;
 701
 702     bcl += bhist[i];     if ( (bcl<bsum*lsat) ) bmin=i;
 703     bcu += bhist[255-i]; if ( (bcu<bsum*usat) ) bmax=255-i;
 704   }
 705
 706   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 707     i_gpix(im, x, y, &val);
 708     val.channel[0]=saturate((val.channel[0]-rmin)*255/(rmax-rmin));
 709     val.channel[1]=saturate((val.channel[1]-gmin)*255/(gmax-gmin));
 710     val.channel[2]=saturate((val.channel[2]-bmin)*255/(bmax-bmin));
 711     i_ppix(im, x, y, &val);
 712   }
 713 }
 714
 715 /*
 716 =item Noise(x,y)
 717
 718 Pseudo noise utility function used to generate perlin noise. (internal)
 719
 720   x - x coordinate
 721   y - y coordinate
 722
 723 =cut
 724 */
 725
 726 static
 727 double
 728 Noise(i_img_dim x, i_img_dim y) {
 729   i_img_dim n = x + y * 57;
 730   n = (n<<13) ^ n;
 731   return ( 1.0 - ( (n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
 732 }
 733
 734 /*
 735 =item SmoothedNoise1(x,y)
 736
 737 Pseudo noise utility function used to generate perlin noise. (internal)
 738
 739   x - x coordinate
 740   y - y coordinate
 741
 742 =cut
 743 */
 744
 745 static
 746 double
 747 SmoothedNoise1(double x, double y) {
 748   double corners = ( Noise(x-1, y-1)+Noise(x+1, y-1)+Noise(x-1, y+1)+Noise(x+1, y+1) ) / 16;
 749   double sides   = ( Noise(x-1, y)  +Noise(x+1, y)  +Noise(x, y-1)  +Noise(x, y+1) ) /  8;
 750   double center  =  Noise(x, y) / 4;
 751   return corners + sides + center;
 752 }
 753
 754
 755 /*
 756 =item G_Interpolate(a, b, x)
 757
 758 Utility function used to generate perlin noise. (internal)
 759
 760 =cut
 761 */
 762
 763 static
 764 double
 765 C_Interpolate(double a, double b, double x) {
 766   /*  float ft = x * 3.1415927; */
 767   double ft = x * PI;
 768   double f = (1 - cos(ft)) * .5;
 769   return  a*(1-f) + b*f;
 770 }
 771
 772
 773 /*
 774 =item InterpolatedNoise(x, y)
 775
 776 Utility function used to generate perlin noise. (internal)
 777
 778 =cut
 779 */
 780
 781 static
 782 double
 783 InterpolatedNoise(double x, double y) {
 784
 785   i_img_dim integer_X = x;
 786   double fractional_X = x - integer_X;
 787   i_img_dim integer_Y = y;
 788   double fractional_Y = y - integer_Y;
 789
 790   double v1 = SmoothedNoise1(integer_X,     integer_Y);
 791   double v2 = SmoothedNoise1(integer_X + 1, integer_Y);
 792   double v3 = SmoothedNoise1(integer_X,     integer_Y + 1);
 793   double v4 = SmoothedNoise1(integer_X + 1, integer_Y + 1);
 794
 795   double i1 = C_Interpolate(v1 , v2 , fractional_X);
 796   double i2 = C_Interpolate(v3 , v4 , fractional_X);
 797
 798   return C_Interpolate(i1 , i2 , fractional_Y);
 799 }
 800
 801
 802
 803 /*
 804 =item PerlinNoise_2D(x, y)
 805
 806 Utility function used to generate perlin noise. (internal)
 807
 808 =cut
 809 */
 810
 811 static
 812 float
 813 PerlinNoise_2D(float x, float y) {
 814   int i,frequency;
 815   double amplitude;
 816   double total = 0;
 817   int Number_Of_Octaves=6;
 818   int n = Number_Of_Octaves - 1;
 819
 820   for(i=0;i<n;i++) {
 821     frequency = 2*i;
 822     amplitude = PI;
 823     total = total + InterpolatedNoise(x * frequency, y * frequency) * amplitude;
 824   }
 825
 826   return total;
 827 }
 828
 829
 830 /*
 831 =item i_radnoise(im, xo, yo, rscale, ascale)
 832
 833 Perlin-like radial noise.
 834
 835   im     - target image
 836   xo     - x coordinate of center
 837   yo     - y coordinate of center
 838   rscale - radial scale
 839   ascale - angular scale
 840
 841 =cut
 842 */
 843
 844 void
 845 i_radnoise(i_img *im, i_img_dim xo, i_img_dim yo, double rscale, double ascale) {
 846   i_img_dim x, y;
 847   int ch;
 848   i_color val;
 849   unsigned char v;
 850   double xc, yc, r;
 851   double a;
 852
 853   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 854     xc = (double)x-xo+0.5;
 855     yc = (double)y-yo+0.5;
 856     r = rscale*sqrt(xc*xc+yc*yc)+1.2;
 857     a = (PI+atan2(yc,xc))*ascale;
 858     v = saturate(128+100*(PerlinNoise_2D(a,r)));
 859     /* v=saturate(120+12*PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale));  Good soft marble */
 860     for(ch=0; ch<im->channels; ch++) val.channel[ch]=v;
 861     i_ppix(im, x, y, &val);
 862   }
 863 }
 864
 865
 866 /*
 867 =item i_turbnoise(im, xo, yo, scale)
 868
 869 Perlin-like 2d noise noise.
 870
 871   im     - target image
 872   xo     - x coordinate translation
 873   yo     - y coordinate translation
 874   scale  - scale of noise
 875
 876 =cut
 877 */
 878
 879 void
 880 i_turbnoise(i_img *im, double xo, double yo, double scale) {
 881   i_img_dim x,y;
 882   int ch;
 883   unsigned char v;
 884   i_color val;
 885
 886   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 887     /*    v=saturate(125*(1.0+PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale))); */
 888     v = saturate(120*(1.0+sin(xo+(double)x/scale+PerlinNoise_2D(xo+(double)x/scale,yo+(float)y/scale))));
 889     for(ch=0; ch<im->channels; ch++) val.channel[ch] = v;
 890     i_ppix(im, x, y, &val);
 891   }
 892 }
 893
 894
 895
 896 /*
 897 =item i_gradgen(im, num, xo, yo, ival, dmeasure)
 898
 899 Gradient generating function.
 900
 901   im     - target image
 902   num    - number of points given
 903   xo     - array of x coordinates
 904   yo     - array of y coordinates
 905   ival   - array of i_color objects
 906   dmeasure - distance measure to be used.
 907     0 = Euclidean
 908     1 = Euclidean squared
 909     2 = Manhattan distance
 910
 911 =cut
 912 */
 913
 914
 915 void
 916 i_gradgen(i_img *im, int num, i_img_dim *xo, i_img_dim *yo, i_color *ival, int dmeasure) {
 917
 918   i_color val;
 919   int p, ch;
 920   i_img_dim x, y;
 921   int channels = im->channels;
 922   i_img_dim xsize    = im->xsize;
 923   i_img_dim ysize    = im->ysize;
 924   size_t bytes;
 925
 926   double *fdist;
 927
 928   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 929
 930   for(p = 0; p<num; p++) {
 931     mm_log((1,"i_gradgen: p%d(" i_DFp ")\n", p, i_DFcp(xo[p], yo[p])));
 932     ICL_info(&ival[p]);
 933   }
 934
 935   /* on the systems I have sizeof(float) == sizeof(int) and thus
 936      this would be same size as the arrays xo and yo point at, but this
 937      may not be true for other systems
 938
 939      since the arrays here are caller controlled, I assume that on
 940      overflow is a programming error rather than an end-user error, so
 941      calling exit() is justified.
 942   */
 943   bytes = sizeof(double) * num;
 944   if (bytes / num != sizeof(double)) {
 945     fprintf(stderr, "integer overflow calculating memory allocation");
 946     exit(1);
 947   }
 948   fdist = mymalloc( bytes ); /* checked 14jul05 tonyc */
 949
 950   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 951     double cs = 0;
 952     double csd = 0;
 953     for(p = 0; p<num; p++) {
 954       i_img_dim xd    = x-xo[p];
 955       i_img_dim yd    = y-yo[p];
 956       switch (dmeasure) {
 957       case 0: /* euclidean */
 958         fdist[p]  = sqrt(xd*xd + yd*yd); /* euclidean distance */
 959         break;
 960       case 1: /* euclidean squared */
 961         fdist[p]  = xd*xd + yd*yd; /* euclidean distance */
 962         break;
 963       case 2: /* euclidean squared */
 964         fdist[p]  = i_max(xd*xd, yd*yd); /* manhattan distance */
 965         break;
 966       default:
 967         i_fatal(3,"i_gradgen: Unknown distance measure\n");
 968       }
 969       cs += fdist[p];
 970     }
 971
 972     csd = 1/((num-1)*cs);
 973
 974     for(p = 0; p<num; p++) fdist[p] = (cs-fdist[p])*csd;
 975
 976     for(ch = 0; ch<channels; ch++) {
 977       int tres = 0;
 978       for(p = 0; p<num; p++) tres += ival[p].channel[ch] * fdist[p];
 979       val.channel[ch] = saturate(tres);
 980     }
 981     i_ppix(im, x, y, &val);
 982   }
 983   myfree(fdist);
 984
 985 }
 986
 987 void
 988 i_nearest_color_foo(i_img *im, int num, i_img_dim *xo, i_img_dim *yo, i_color *ival, int dmeasure) {
 989
 990   int p;
 991   i_img_dim x, y;
 992   i_img_dim xsize    = im->xsize;
 993   i_img_dim ysize    = im->ysize;
 994
 995   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 996
 997   for(p = 0; p<num; p++) {
 998     mm_log((1,"i_gradgen: p%d(" i_DFp ")\n", p, i_DFcp(xo[p], yo[p])));
 999     ICL_info(&ival[p]);
1000   }
1001
1002   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1003     int    midx    = 0;
1004     double mindist = 0;
1005     double curdist = 0;
1006
1007     i_img_dim xd        = x-xo[0];
1008     i_img_dim yd        = y-yo[0];
1009
1010     switch (dmeasure) {
1011     case 0: /* euclidean */
1012       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1013       break;
1014     case 1: /* euclidean squared */
1015       mindist = xd*xd + yd*yd; /* euclidean distance */
1016       break;
1017     case 2: /* euclidean squared */
1018       mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
1019       break;
1020     default:
1021       i_fatal(3,"i_nearest_color: Unknown distance measure\n");
1022     }
1023
1024     for(p = 1; p<num; p++) {
1025       xd = x-xo[p];
1026       yd = y-yo[p];
1027       switch (dmeasure) {
1028       case 0: /* euclidean */
1029         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1030         break;
1031       case 1: /* euclidean squared */
1032         curdist = xd*xd + yd*yd; /* euclidean distance */
1033         break;
1034       case 2: /* euclidean squared */
1035         curdist = i_max(xd*xd, yd*yd); /* manhattan distance */
1036         break;
1037       default:
1038         i_fatal(3,"i_nearest_color: Unknown distance measure\n");
1039       }
1040       if (curdist < mindist) {
1041         mindist = curdist;
1042         midx = p;
1043       }
1044     }
1045     i_ppix(im, x, y, &ival[midx]);
1046   }
1047 }
1048
1049 /*
1050 =item i_nearest_color(im, num, xo, yo, oval, dmeasure)
1051
1052 This wasn't document - quoth Addi:
1053
1054   An arty type of filter
1055
1056 FIXME: check IRC logs for actual text.
1057
1058 Inputs:
1059
1060 =over
1061
1062 =item *
1063
1064 i_img *im - image to render on.
1065
1066 =item *
1067
1068 int num - number of points/colors in xo, yo, oval
1069
1070 =item *
1071
1072 i_img_dim *xo - array of I<num> x positions
1073
1074 =item *
1075
1076 i_img_dim *yo - array of I<num> y positions
1077
1078 =item *
1079
1080 i_color *oval - array of I<num> colors
1081
1082 xo, yo, oval correspond to each other, the point xo[i], yo[i] has a
1083 color something like oval[i], at least closer to that color than other
1084 points.
1085
1086 =item *
1087
1088 int dmeasure - how we measure the distance from some point P(x,y) to
1089 any (xo[i], yo[i]).
1090
1091 Valid values are:
1092
1093 =over
1094
1095 =item 0
1096
1097 euclidean distance: sqrt((x2-x1)**2 + (y2-y1)**2)
1098
1099 =item 1
1100
1101 square of euclidean distance: ((x2-x1)**2 + (y2-y1)**2)
1102
1103 =item 2
1104
1105 manhattan distance: max((y2-y1)**2, (x2-x1)**2)
1106
1107 =back
1108
1109 An invalid value causes an error exit (the program is aborted).
1110
1111 =back
1112
1113 =cut
1114  */
1115
1116 int
1117 i_nearest_color(i_img *im, int num, i_img_dim *xo, i_img_dim *yo, i_color *oval, int dmeasure) {
1118   i_color *ival;
1119   float *tval;
1120   double c1, c2;
1121   i_color val;
1122   int p, ch;
1123   i_img_dim x, y;
1124   i_img_dim xsize    = im->xsize;
1125   i_img_dim ysize    = im->ysize;
1126   int *cmatch;
1127   size_t ival_bytes, tval_bytes;
1128
1129   mm_log((1,"i_nearest_color(im %p, num %d, xo %p, yo %p, oval %p, dmeasure %d)\n", im, num, xo, yo, oval, dmeasure));
1130
1131   i_clear_error();
1132
1133   if (num <= 0) {
1134     i_push_error(0, "no points supplied to nearest_color filter");
1135     return 0;
1136   }
1137
1138   if (dmeasure < 0 || dmeasure > i_dmeasure_limit) {
1139     i_push_error(0, "distance measure invalid");
1140     return 0;
1141   }
1142
1143   tval_bytes = sizeof(float)*num*im->channels;
1144   if (tval_bytes / num != sizeof(float) * im->channels) {
1145     i_push_error(0, "integer overflow calculating memory allocation");
1146     return 0;
1147   }
1148   ival_bytes  = sizeof(i_color) * num;
1149   if (ival_bytes / sizeof(i_color) != num) {
1150     i_push_error(0, "integer overflow calculating memory allocation");
1151     return 0;
1152   }
1153   tval   = mymalloc( tval_bytes ); /* checked 17feb2005 tonyc */
1154   ival   = mymalloc( ival_bytes ); /* checked 17feb2005 tonyc */
1155   cmatch = mymalloc( sizeof(int)*num     ); /* checked 17feb2005 tonyc */
1156
1157   for(p = 0; p<num; p++) {
1158     for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
1159     cmatch[p] = 0;
1160   }
1161
1162
1163   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1164     int   midx    = 0;
1165     double mindist = 0;
1166     double curdist = 0;
1167
1168     i_img_dim xd        = x-xo[0];
1169     i_img_dim yd        = y-yo[0];
1170
1171     switch (dmeasure) {
1172     case 0: /* euclidean */
1173       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1174       break;
1175     case 1: /* euclidean squared */
1176       mindist = xd*xd + yd*yd; /* euclidean distance */
1177       break;
1178     case 2: /* manhatten distance */
1179       mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
1180       break;
1181     default:
1182       i_fatal(3,"i_nearest_color: Unknown distance measure\n");
1183     }
1184
1185     for(p = 1; p<num; p++) {
1186       xd = x-xo[p];
1187       yd = y-yo[p];
1188       switch (dmeasure) {
1189       case 0: /* euclidean */
1190         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1191         break;
1192       case 1: /* euclidean squared */
1193         curdist = xd*xd + yd*yd; /* euclidean distance */
1194         break;
1195       case 2: /* euclidean squared */
1196         curdist = i_max(xd*xd, yd*yd); /* manhattan distance */
1197         break;
1198       default:
1199         i_fatal(3,"i_nearest_color: Unknown distance measure\n");
1200       }
1201       if (curdist < mindist) {
1202         mindist = curdist;
1203         midx = p;
1204       }
1205     }
1206
1207     cmatch[midx]++;
1208     i_gpix(im, x, y, &val);
1209     c2 = 1.0/(float)(cmatch[midx]);
1210     c1 = 1.0-c2;
1211
1212     for(ch = 0; ch<im->channels; ch++)
1213       tval[midx*im->channels + ch] =
1214         c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
1215
1216   }
1217
1218   for(p = 0; p<num; p++) {
1219     for(ch = 0; ch<im->channels; ch++)
1220       ival[p].channel[ch] = tval[p*im->channels + ch];
1221
1222     /* avoid uninitialized value messages from valgrind */
1223     while (ch < MAXCHANNELS)
1224       ival[p].channel[ch++] = 0;
1225   }
1226
1227   i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
1228
1229   return 1;
1230 }
1231
1232 /*
1233 =item i_unsharp_mask(im, stddev, scale)
1234
1235 Perform an usharp mask, which is defined as subtracting the blurred
1236 image from double the original.
1237
1238 =cut
1239 */
1240
1241 void
1242 i_unsharp_mask(i_img *im, double stddev, double scale) {
1243   i_img *copy;
1244   i_img_dim x, y;
1245   int ch;
1246
1247   if (scale < 0)
1248     return;
1249   /* it really shouldn't ever be more than 1.0, but maybe ... */
1250   if (scale > 100)
1251     scale = 100;
1252
1253   copy = i_copy(im);
1254   i_gaussian(copy, stddev);
1255   if (im->bits == i_8_bits) {
1256     i_color *blur = mymalloc(im->xsize * sizeof(i_color)); /* checked 17feb2005 tonyc */
1257     i_color *out = mymalloc(im->xsize * sizeof(i_color)); /* checked 17feb2005 tonyc */
1258
1259     for (y = 0; y < im->ysize; ++y) {
1260       i_glin(copy, 0, copy->xsize, y, blur);
1261       i_glin(im, 0, im->xsize, y, out);
1262       for (x = 0; x < im->xsize; ++x) {
1263         for (ch = 0; ch < im->channels; ++ch) {
1264           /*int temp = out[x].channel[ch] +
1265             scale * (out[x].channel[ch] - blur[x].channel[ch]);*/
1266           int temp = out[x].channel[ch] * 2 - blur[x].channel[ch];
1267           if (temp < 0)
1268             temp = 0;
1269           else if (temp > 255)
1270             temp = 255;
1271           out[x].channel[ch] = temp;
1272         }
1273       }
1274       i_plin(im, 0, im->xsize, y, out);
1275     }
1276
1277     myfree(blur);
1278     myfree(out);
1279   }
1280   else {
1281     i_fcolor *blur = mymalloc(im->xsize * sizeof(i_fcolor)); /* checked 17feb2005 tonyc */
1282     i_fcolor *out = mymalloc(im->xsize * sizeof(i_fcolor)); /* checked 17feb2005 tonyc */
1283
1284     for (y = 0; y < im->ysize; ++y) {
1285       i_glinf(copy, 0, copy->xsize, y, blur);
1286       i_glinf(im, 0, im->xsize, y, out);
1287       for (x = 0; x < im->xsize; ++x) {
1288         for (ch = 0; ch < im->channels; ++ch) {
1289           double temp = out[x].channel[ch] +
1290             scale * (out[x].channel[ch] - blur[x].channel[ch]);
1291           if (temp < 0)
1292             temp = 0;
1293           else if (temp > 1.0)
1294             temp = 1.0;
1295           out[x].channel[ch] = temp;
1296         }
1297       }
1298       i_plinf(im, 0, im->xsize, y, out);
1299     }
1300
1301     myfree(blur);
1302     myfree(out);
1303   }
1304   i_img_destroy(copy);
1305 }
1306
1307 /*
1308 =item i_diff_image(im1, im2, mindist)
1309
1310 Creates a new image that is transparent, except where the pixel in im2
1311 is different from im1, where it is the pixel from im2.
1312
1313 The samples must differ by at least mindiff to be considered different.
1314
1315 =cut
1316 */
1317
1318 i_img *
1319 i_diff_image(i_img *im1, i_img *im2, double mindist) {
1320   i_img *out;
1321   int outchans, diffchans;
1322   i_img_dim xsize, ysize;
1323
1324   i_clear_error();
1325   if (im1->channels != im2->channels) {
1326     i_push_error(0, "different number of channels");
1327     return NULL;
1328   }
1329
1330   outchans = diffchans = im1->channels;
1331   if (outchans == 1 || outchans == 3)
1332     ++outchans;
1333
1334   xsize = i_min(im1->xsize, im2->xsize);
1335   ysize = i_min(im1->ysize, im2->ysize);
1336
1337   out = i_sametype_chans(im1, xsize, ysize, outchans);
1338
1339   if (im1->bits == i_8_bits && im2->bits == i_8_bits) {
1340     i_color *line1 = mymalloc(xsize * sizeof(*line1)); /* checked 17feb2005 tonyc */
1341     i_color *line2 = mymalloc(xsize * sizeof(*line1)); /* checked 17feb2005 tonyc */
1342     i_color empty;
1343     i_img_dim x, y;
1344     int ch;
1345     int imindist = (int)mindist;
1346
1347     for (ch = 0; ch < MAXCHANNELS; ++ch)
1348       empty.channel[ch] = 0;
1349
1350     for (y = 0; y < ysize; ++y) {
1351       i_glin(im1, 0, xsize, y, line1);
1352       i_glin(im2, 0, xsize, y, line2);
1353       if (outchans != diffchans) {
1354         /* give the output an alpha channel since it doesn't have one */
1355         for (x = 0; x < xsize; ++x)
1356           line2[x].channel[diffchans] = 255;
1357       }
1358       for (x = 0; x < xsize; ++x) {
1359         int diff = 0;
1360         for (ch = 0; ch < diffchans; ++ch) {
1361           if (line1[x].channel[ch] != line2[x].channel[ch]
1362               && abs(line1[x].channel[ch] - line2[x].channel[ch]) > imindist) {
1363             diff = 1;
1364             break;
1365           }
1366         }
1367         if (!diff)
1368           line2[x] = empty;
1369       }
1370       i_plin(out, 0, xsize, y, line2);
1371     }
1372     myfree(line1);
1373     myfree(line2);
1374   }
1375   else {
1376     i_fcolor *line1 = mymalloc(xsize * sizeof(*line1)); /* checked 17feb2005 tonyc */
1377     i_fcolor *line2 = mymalloc(xsize * sizeof(*line2)); /* checked 17feb2005 tonyc */
1378     i_fcolor empty;
1379     i_img_dim x, y;
1380     int ch;
1381     double dist = mindist / 255.0;
1382
1383     for (ch = 0; ch < MAXCHANNELS; ++ch)
1384       empty.channel[ch] = 0;
1385
1386     for (y = 0; y < ysize; ++y) {
1387       i_glinf(im1, 0, xsize, y, line1);
1388       i_glinf(im2, 0, xsize, y, line2);
1389       if (outchans != diffchans) {
1390         /* give the output an alpha channel since it doesn't have one */
1391         for (x = 0; x < xsize; ++x)
1392           line2[x].channel[diffchans] = 1.0;
1393       }
1394       for (x = 0; x < xsize; ++x) {
1395         int diff = 0;
1396         for (ch = 0; ch < diffchans; ++ch) {
1397           if (line1[x].channel[ch] != line2[x].channel[ch]
1398               && fabs(line1[x].channel[ch] - line2[x].channel[ch]) > dist) {
1399             diff = 1;
1400             break;
1401           }
1402         }
1403         if (!diff)
1404           line2[x] = empty;
1405       }
1406       i_plinf(out, 0, xsize, y, line2);
1407     }
1408     myfree(line1);
1409     myfree(line2);
1410   }
1411
1412   return out;
1413 }
1414
1415 struct fount_state;
1416 static double linear_fount_f(double x, double y, struct fount_state *state);
1417 static double bilinear_fount_f(double x, double y, struct fount_state *state);
1418 static double radial_fount_f(double x, double y, struct fount_state *state);
1419 static double square_fount_f(double x, double y, struct fount_state *state);
1420 static double revolution_fount_f(double x, double y,
1421                                  struct fount_state *state);
1422 static double conical_fount_f(double x, double y, struct fount_state *state);
1423
1424 typedef double (*fount_func)(double, double, struct fount_state *);
1425 static fount_func fount_funcs[] =
1426 {
1427   linear_fount_f,
1428   bilinear_fount_f,
1429   radial_fount_f,
1430   square_fount_f,
1431   revolution_fount_f,
1432   conical_fount_f,
1433 };
1434
1435 static double linear_interp(double pos, i_fountain_seg *seg);
1436 static double sine_interp(double pos, i_fountain_seg *seg);
1437 static double sphereup_interp(double pos, i_fountain_seg *seg);
1438 static double spheredown_interp(double pos, i_fountain_seg *seg);
1439 typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
1440 static fount_interp fount_interps[] =
1441 {
1442   linear_interp,
1443   linear_interp,
1444   sine_interp,
1445   sphereup_interp,
1446   spheredown_interp,
1447 };
1448
1449 static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1450 static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1451 static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1452 typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
1453 static fount_cinterp fount_cinterps[] =
1454 {
1455   direct_cinterp,
1456   hue_up_cinterp,
1457   hue_down_cinterp,
1458 };
1459
1460 typedef double (*fount_repeat)(double v);
1461 static double fount_r_none(double v);
1462 static double fount_r_sawtooth(double v);
1463 static double fount_r_triangle(double v);
1464 static double fount_r_saw_both(double v);
1465 static double fount_r_tri_both(double v);
1466 static fount_repeat fount_repeats[] =
1467 {
1468   fount_r_none,
1469   fount_r_sawtooth,
1470   fount_r_triangle,
1471   fount_r_saw_both,
1472   fount_r_tri_both,
1473 };
1474
1475 static int simple_ssample(i_fcolor *out, double x, double y,
1476                            struct fount_state *state);
1477 static int random_ssample(i_fcolor *out, double x, double y,
1478                            struct fount_state *state);
1479 static int circle_ssample(i_fcolor *out, double x, double y,
1480                            struct fount_state *state);
1481 typedef int (*fount_ssample)(i_fcolor *out, double x, double y,
1482                               struct fount_state *state);
1483 static fount_ssample fount_ssamples[] =
1484 {
1485   NULL,
1486   simple_ssample,
1487   random_ssample,
1488   circle_ssample,
1489 };
1490
1491 static int
1492 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
1493
1494 /*
1495   Keep state information used by each type of fountain fill
1496 */
1497 struct fount_state {
1498   /* precalculated for the equation of the line perpendicular to the line AB */
1499   double lA, lB, lC;
1500   double AB;
1501   double sqrtA2B2;
1502   double mult;
1503   double cos;
1504   double sin;
1505   double theta;
1506   i_img_dim xa, ya;
1507   void *ssample_data;
1508   fount_func ffunc;
1509   fount_repeat rpfunc;
1510   fount_ssample ssfunc;
1511   double parm;
1512   i_fountain_seg *segs;
1513   int count;
1514 };
1515
1516 static void
1517 fount_init_state(struct fount_state *state, double xa, double ya,
1518                  double xb, double yb, i_fountain_type type,
1519                  i_fountain_repeat repeat, int combine, int super_sample,
1520                  double ssample_param, int count, i_fountain_seg *segs);
1521
1522 static void
1523 fount_finish_state(struct fount_state *state);
1524
1525 #define EPSILON (1e-6)
1526
1527 /*
1528 =item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1529
1530 Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1531
1532 I<type> controls how the reference points are used:
1533
1534 =over
1535
1536 =item i_ft_linear
1537
1538 linear, where A is 0 and B is 1.
1539
1540 =item i_ft_bilinear
1541
1542 linear in both directions from A.
1543
1544 =item i_ft_radial
1545
1546 circular, where A is the centre of the fill, and B is a point
1547 on the radius.
1548
1549 =item i_ft_radial_square
1550
1551 where A is the centre of the fill and B is the centre of
1552 one side of the square.
1553
1554 =item i_ft_revolution
1555
1556 where A is the centre of the fill and B defines the 0/1.0
1557 angle of the fill.
1558
1559 =item i_ft_conical
1560
1561 similar to i_ft_revolution, except that the revolution goes in both
1562 directions
1563
1564 =back
1565
1566 I<repeat> can be one of:
1567
1568 =over
1569
1570 =item i_fr_none
1571
1572 values < 0 are treated as zero, values > 1 are treated as 1.
1573
1574 =item i_fr_sawtooth
1575
1576 negative values are treated as 0, positive values are modulo 1.0
1577
1578 =item i_fr_triangle
1579
1580 negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1581 mod 1.0), otherwise the same as for sawtooth.
1582
1583 =item i_fr_saw_both
1584
1585 like i_fr_sawtooth, except that the sawtooth pattern repeats into
1586 negative values.
1587
1588 =item i_fr_tri_both
1589
1590 Like i_fr_triangle, except that negative values are handled as their
1591 absolute values.
1592
1593 =back
1594
1595 If combine is non-zero then non-opaque values are combined with the
1596 underlying color.
1597
1598 I<super_sample> controls super sampling, if any.  At some point I'll
1599 probably add a adaptive super-sampler.  Current possible values are:
1600
1601 =over
1602
1603 =item i_fts_none
1604
1605 No super-sampling is done.
1606
1607 =item i_fts_grid
1608
1609 A square grid of points withing the pixel are sampled.
1610
1611 =item i_fts_random
1612
1613 Random points within the pixel are sampled.
1614
1615 =item i_fts_circle
1616
1617 Points on the radius of a circle are sampled.  This produces fairly
1618 good results, but is fairly slow since sin() and cos() are evaluated
1619 for each point.
1620
1621 =back
1622
1623 I<ssample_param> is intended to be roughly the number of points
1624 sampled within the pixel.
1625
1626 I<count> and I<segs> define the segments of the fill.
1627
1628 =cut
1629
1630 */
1631
1632 int
1633 i_fountain(i_img *im, double xa, double ya, double xb, double yb,
1634            i_fountain_type type, i_fountain_repeat repeat,
1635            int combine, int super_sample, double ssample_param,
1636            int count, i_fountain_seg *segs) {
1637   struct fount_state state;
1638   i_img_dim x, y;
1639   i_fcolor *line = NULL;
1640   i_fcolor *work = NULL;
1641   size_t line_bytes;
1642   i_fountain_seg *my_segs;
1643   i_fill_combine_f combine_func = NULL;
1644   i_fill_combinef_f combinef_func = NULL;
1645
1646   i_clear_error();
1647
1648   /* i_fountain() allocates floating colors even for 8-bit images,
1649      so we need to do this check */
1650   line_bytes = sizeof(i_fcolor) * im->xsize;
1651   if (line_bytes / sizeof(i_fcolor) != im->xsize) {
1652     i_push_error(0, "integer overflow calculating memory allocation");
1653     return 0;
1654   }
1655
1656   line = mymalloc(line_bytes); /* checked 17feb2005 tonyc */
1657
1658   i_get_combine(combine, &combine_func, &combinef_func);
1659   if (combinef_func)
1660     work = mymalloc(line_bytes); /* checked 17feb2005 tonyc */
1661
1662   fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
1663                    super_sample, ssample_param, count, segs);
1664   my_segs = state.segs;
1665
1666   for (y = 0; y < im->ysize; ++y) {
1667     i_glinf(im, 0, im->xsize, y, line);
1668     for (x = 0; x < im->xsize; ++x) {
1669       i_fcolor c;
1670       int got_one;
1671       if (super_sample == i_fts_none)
1672         got_one = fount_getat(&c, x, y, &state);
1673       else
1674         got_one = state.ssfunc(&c, x, y, &state);
1675       if (got_one) {
1676         if (combine)
1677           work[x] = c;
1678         else
1679           line[x] = c;
1680       }
1681     }
1682     if (combine)
1683       combinef_func(line, work, im->channels, im->xsize);
1684     i_plinf(im, 0, im->xsize, y, line);
1685   }
1686   fount_finish_state(&state);
1687   if (work) myfree(work);
1688   myfree(line);
1689
1690   return 1;
1691 }
1692
1693 typedef struct {
1694   i_fill_t base;
1695   struct fount_state state;
1696 } i_fill_fountain_t;
1697
1698 static void
1699 fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width, int channels,
1700             i_fcolor *data);
1701 static void
1702 fount_fill_destroy(i_fill_t *fill);
1703
1704 static i_fill_fountain_t
1705 fount_fill_proto =
1706   {
1707     {
1708       NULL,
1709       fill_fountf,
1710       fount_fill_destroy
1711     }
1712   };
1713
1714
1715 /*
1716 =item i_new_fill_fount(C<xa>, C<ya>, C<xb>, C<yb>, C<type>, C<repeat>, C<combine>, C<super_sample>, C<ssample_param>, C<count>, C<segs>)
1717
1718 =category Fills
1719 =synopsis fill = i_new_fill_fount(0, 0, 100, 100, i_ft_linear, i_ft_linear,
1720 =synopsis                         i_fr_triangle, 0, i_fts_grid, 9, 1, segs);
1721
1722
1723 Creates a new general fill which fills with a fountain fill.
1724
1725 =cut
1726 */
1727
1728 i_fill_t *
1729 i_new_fill_fount(double xa, double ya, double xb, double yb,
1730                  i_fountain_type type, i_fountain_repeat repeat,
1731                  int combine, int super_sample, double ssample_param,
1732                  int count, i_fountain_seg *segs) {
1733   i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1734
1735   *fill = fount_fill_proto;
1736   if (combine)
1737     i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1738   else {
1739     fill->base.combine = NULL;
1740     fill->base.combinef = NULL;
1741   }
1742   fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1743                    super_sample, ssample_param, count, segs);
1744
1745   return &fill->base;
1746 }
1747
1748 /*
1749 =back
1750
1751 =head1 INTERNAL FUNCTIONS
1752
1753 =over
1754
1755 =item fount_init_state(...)
1756
1757 Used by both the fountain fill filter and the fountain fill.
1758
1759 =cut
1760 */
1761
1762 static void
1763 fount_init_state(struct fount_state *state, double xa, double ya,
1764                  double xb, double yb, i_fountain_type type,
1765                  i_fountain_repeat repeat, int combine, int super_sample,
1766                  double ssample_param, int count, i_fountain_seg *segs) {
1767   int i, j;
1768   size_t bytes;
1769   i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count); /* checked 2jul06 - duplicating original */
1770   /*int have_alpha = im->channels == 2 || im->channels == 4;*/
1771
1772   memset(state, 0, sizeof(*state));
1773   /* we keep a local copy that we can adjust for speed */
1774   for (i = 0; i < count; ++i) {
1775     i_fountain_seg *seg = my_segs + i;
1776
1777     *seg = segs[i];
1778     if (seg->type < 0 || seg->type >= i_fst_end)
1779       seg->type = i_fst_linear;
1780     if (seg->color < 0 || seg->color >= i_fc_end)
1781       seg->color = i_fc_direct;
1782     if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1783       /* so we don't have to translate to HSV on each request, do it here */
1784       for (j = 0; j < 2; ++j) {
1785         i_rgb_to_hsvf(seg->c+j);
1786       }
1787       if (seg->color == i_fc_hue_up) {
1788         if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1789           seg->c[1].channel[0] += 1.0;
1790       }
1791       else {
1792         if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1793           seg->c[0].channel[0] += 1.0;
1794       }
1795     }
1796     /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1797            seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1798            seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1799            seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1800            seg->c[1].channel[3], seg->type, seg->color);*/
1801
1802   }
1803
1804   /* initialize each engine */
1805   /* these are so common ... */
1806   state->lA = xb - xa;
1807   state->lB = yb - ya;
1808   state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1809   state->xa = xa;
1810   state->ya = ya;
1811   switch (type) {
1812   default:
1813     type = i_ft_linear; /* make the invalid value valid */
1814   case i_ft_linear:
1815   case i_ft_bilinear:
1816     state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1817     state->mult = 1;
1818     state->mult = 1/linear_fount_f(xb, yb, state);
1819     break;
1820
1821   case i_ft_radial:
1822     state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1823                              + (double)(yb-ya)*(yb-ya));
1824     break;
1825
1826   case i_ft_radial_square:
1827     state->cos = state->lA / state->AB;
1828     state->sin = state->lB / state->AB;
1829     state->mult = 1.0 / state->AB;
1830     break;
1831
1832   case i_ft_revolution:
1833     state->theta = atan2(yb-ya, xb-xa);
1834     state->mult = 1.0 / (PI * 2);
1835     break;
1836
1837   case i_ft_conical:
1838     state->theta = atan2(yb-ya, xb-xa);
1839     state->mult = 1.0 / PI;
1840     break;
1841   }
1842   state->ffunc = fount_funcs[type];
1843   if (super_sample < 0
1844       || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
1845     super_sample = 0;
1846   }
1847   state->ssample_data = NULL;
1848   switch (super_sample) {
1849   case i_fts_grid:
1850     ssample_param = floor(0.5 + sqrt(ssample_param));
1851     bytes = ssample_param * ssample_param * sizeof(i_fcolor);
1852     if (bytes / sizeof(i_fcolor) == ssample_param * ssample_param) {
1853       state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param); /* checked 1jul06 tonyc */
1854     }
1855     else {
1856       super_sample = i_fts_none;
1857     }
1858     break;
1859
1860   case i_fts_random:
1861   case i_fts_circle:
1862     ssample_param = floor(0.5+ssample_param);
1863     bytes = sizeof(i_fcolor) * ssample_param;
1864     if (bytes / sizeof(i_fcolor) == ssample_param) {
1865       state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
1866     }
1867     else {
1868       super_sample = i_fts_none;
1869     }
1870     break;
1871   }
1872   state->parm = ssample_param;
1873   state->ssfunc = fount_ssamples[super_sample];
1874   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1875     repeat = 0;
1876   state->rpfunc = fount_repeats[repeat];
1877   state->segs = my_segs;
1878   state->count = count;
1879 }
1880
1881 static void
1882 fount_finish_state(struct fount_state *state) {
1883   if (state->ssample_data)
1884     myfree(state->ssample_data);
1885   myfree(state->segs);
1886 }
1887
1888
1889 /*
1890 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
1891
1892 Evaluates the fountain fill at the given point.
1893
1894 This is called by both the non-super-sampling and super-sampling code.
1895
1896 You might think that it would make sense to sample the fill parameter
1897 instead, and combine those, but this breaks badly.
1898
1899 =cut
1900 */
1901
1902 static int
1903 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
1904   double v = (state->rpfunc)((state->ffunc)(x, y, state));
1905   int i;
1906
1907   i = 0;
1908   while (i < state->count
1909          && (v < state->segs[i].start || v > state->segs[i].end)) {
1910     ++i;
1911   }
1912   if (i < state->count) {
1913     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
1914     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
1915     return 1;
1916   }
1917   else
1918     return 0;
1919 }
1920
1921 /*
1922 =item linear_fount_f(x, y, state)
1923
1924 Calculate the fill parameter for a linear fountain fill.
1925
1926 Uses the point to line distance function, with some precalculation
1927 done in i_fountain().
1928
1929 =cut
1930 */
1931 static double
1932 linear_fount_f(double x, double y, struct fount_state *state) {
1933   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
1934 }
1935
1936 /*
1937 =item bilinear_fount_f(x, y, state)
1938
1939 Calculate the fill parameter for a bi-linear fountain fill.
1940
1941 =cut
1942 */
1943 static double
1944 bilinear_fount_f(double x, double y, struct fount_state *state) {
1945   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
1946 }
1947
1948 /*
1949 =item radial_fount_f(x, y, state)
1950
1951 Calculate the fill parameter for a radial fountain fill.
1952
1953 Simply uses the distance function.
1954
1955 =cut
1956  */
1957 static double
1958 radial_fount_f(double x, double y, struct fount_state *state) {
1959   return sqrt((double)(state->xa-x)*(state->xa-x)
1960               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
1961 }
1962
1963 /*
1964 =item square_fount_f(x, y, state)
1965
1966 Calculate the fill parameter for a square fountain fill.
1967
1968 Works by rotating the reference co-ordinate around the centre of the
1969 square.
1970
1971 =cut
1972 */
1973 static double
1974 square_fount_f(double x, double y, struct fount_state *state) {
1975   i_img_dim xc, yc; /* centred on A */
1976   double xt, yt; /* rotated by theta */
1977   xc = x - state->xa;
1978   yc = y - state->ya;
1979   xt = fabs(xc * state->cos + yc * state->sin);
1980   yt = fabs(-xc * state->sin + yc * state->cos);
1981   return (xt > yt ? xt : yt) * state->mult;
1982 }
1983
1984 /*
1985 =item revolution_fount_f(x, y, state)
1986
1987 Calculates the fill parameter for the revolution fountain fill.
1988
1989 =cut
1990 */
1991 static double
1992 revolution_fount_f(double x, double y, struct fount_state *state) {
1993   double angle = atan2(y - state->ya, x - state->xa);
1994
1995   angle -= state->theta;
1996   if (angle < 0) {
1997     angle = fmod(angle+ PI * 4, PI*2);
1998   }
1999
2000   return angle * state->mult;
2001 }
2002
2003 /*
2004 =item conical_fount_f(x, y, state)
2005
2006 Calculates the fill parameter for the conical fountain fill.
2007
2008 =cut
2009 */
2010 static double
2011 conical_fount_f(double x, double y, struct fount_state *state) {
2012   double angle = atan2(y - state->ya, x - state->xa);
2013
2014   angle -= state->theta;
2015   if (angle < -PI)
2016     angle += PI * 2;
2017   else if (angle > PI)
2018     angle -= PI * 2;
2019
2020   return fabs(angle) * state->mult;
2021 }
2022
2023 /*
2024 =item linear_interp(pos, seg)
2025
2026 Calculates linear interpolation on the fill parameter.  Breaks the
2027 segment into 2 regions based in the I<middle> value.
2028
2029 =cut
2030 */
2031 static double
2032 linear_interp(double pos, i_fountain_seg *seg) {
2033   if (pos < seg->middle) {
2034     double len = seg->middle - seg->start;
2035     if (len < EPSILON)
2036       return 0.0;
2037     else
2038       return (pos - seg->start) / len / 2;
2039   }
2040   else {
2041     double len = seg->end - seg->middle;
2042     if (len < EPSILON)
2043       return 1.0;
2044     else
2045       return 0.5 + (pos - seg->middle) / len / 2;
2046   }
2047 }
2048
2049 /*
2050 =item sine_interp(pos, seg)
2051
2052 Calculates sine function interpolation on the fill parameter.
2053
2054 =cut
2055 */
2056 static double
2057 sine_interp(double pos, i_fountain_seg *seg) {
2058   /* I wonder if there's a simple way to smooth the transition for this */
2059   double work = linear_interp(pos, seg);
2060
2061   return (1-cos(work * PI))/2;
2062 }
2063
2064 /*
2065 =item sphereup_interp(pos, seg)
2066
2067 Calculates spherical interpolation on the fill parameter, with the cusp
2068 at the low-end.
2069
2070 =cut
2071 */
2072 static double
2073 sphereup_interp(double pos, i_fountain_seg *seg) {
2074   double work = linear_interp(pos, seg);
2075
2076   return sqrt(1.0 - (1-work) * (1-work));
2077 }
2078
2079 /*
2080 =item spheredown_interp(pos, seg)
2081
2082 Calculates spherical interpolation on the fill parameter, with the cusp
2083 at the high-end.
2084
2085 =cut
2086 */
2087 static double
2088 spheredown_interp(double pos, i_fountain_seg *seg) {
2089   double work = linear_interp(pos, seg);
2090
2091   return 1-sqrt(1.0 - work * work);
2092 }
2093
2094 /*
2095 =item direct_cinterp(out, pos, seg)
2096
2097 Calculates the fountain color based on direct scaling of the channels
2098 of the color channels.
2099
2100 =cut
2101 */
2102 static void
2103 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2104   int ch;
2105   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2106     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2107       + seg->c[1].channel[ch] * pos;
2108   }
2109 }
2110
2111 /*
2112 =item hue_up_cinterp(put, pos, seg)
2113
2114 Calculates the fountain color based on scaling a HSV value.  The hue
2115 increases as the fill parameter increases.
2116
2117 =cut
2118 */
2119 static void
2120 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2121   int ch;
2122   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2123     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2124       + seg->c[1].channel[ch] * pos;
2125   }
2126   i_hsv_to_rgbf(out);
2127 }
2128
2129 /*
2130 =item hue_down_cinterp(put, pos, seg)
2131
2132 Calculates the fountain color based on scaling a HSV value.  The hue
2133 decreases as the fill parameter increases.
2134
2135 =cut
2136 */
2137 static void
2138 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2139   int ch;
2140   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2141     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2142       + seg->c[1].channel[ch] * pos;
2143   }
2144   i_hsv_to_rgbf(out);
2145 }
2146
2147 /*
2148 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2149
2150 Simple grid-based super-sampling.
2151
2152 =cut
2153 */
2154 static int
2155 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
2156   i_fcolor *work = state->ssample_data;
2157   i_img_dim dx, dy;
2158   int grid = state->parm;
2159   double base = -0.5 + 0.5 / grid;
2160   double step = 1.0 / grid;
2161   int ch, i;
2162   int samp_count = 0;
2163
2164   for (dx = 0; dx < grid; ++dx) {
2165     for (dy = 0; dy < grid; ++dy) {
2166       if (fount_getat(work+samp_count, x + base + step * dx,
2167                       y + base + step * dy, state)) {
2168         ++samp_count;
2169       }
2170     }
2171   }
2172   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2173     out->channel[ch] = 0;
2174     for (i = 0; i < samp_count; ++i) {
2175       out->channel[ch] += work[i].channel[ch];
2176     }
2177     /* we divide by 4 rather than samp_count since if there's only one valid
2178        sample it should be mostly transparent */
2179     out->channel[ch] /= grid * grid;
2180   }
2181   return samp_count;
2182 }
2183
2184 /*
2185 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2186
2187 Random super-sampling.
2188
2189 =cut
2190 */
2191 static int
2192 random_ssample(i_fcolor *out, double x, double y,
2193                struct fount_state *state) {
2194   i_fcolor *work = state->ssample_data;
2195   int i, ch;
2196   int maxsamples = state->parm;
2197   double rand_scale = 1.0 / RAND_MAX;
2198   int samp_count = 0;
2199   for (i = 0; i < maxsamples; ++i) {
2200     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
2201                     y - 0.5 + rand() * rand_scale, state)) {
2202       ++samp_count;
2203     }
2204   }
2205   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2206     out->channel[ch] = 0;
2207     for (i = 0; i < samp_count; ++i) {
2208       out->channel[ch] += work[i].channel[ch];
2209     }
2210     /* we divide by maxsamples rather than samp_count since if there's
2211        only one valid sample it should be mostly transparent */
2212     out->channel[ch] /= maxsamples;
2213   }
2214   return samp_count;
2215 }
2216
2217 /*
2218 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2219
2220 Super-sampling around the circumference of a circle.
2221
2222 I considered saving the sin()/cos() values and transforming step-size
2223 around the circle, but that's inaccurate, though it may not matter
2224 much.
2225
2226 =cut
2227  */
2228 static int
2229 circle_ssample(i_fcolor *out, double x, double y,
2230                struct fount_state *state) {
2231   i_fcolor *work = state->ssample_data;
2232   int i, ch;
2233   int maxsamples = state->parm;
2234   double angle = 2 * PI / maxsamples;
2235   double radius = 0.3; /* semi-random */
2236   int samp_count = 0;
2237   for (i = 0; i < maxsamples; ++i) {
2238     if (fount_getat(work+samp_count, x + radius * cos(angle * i),
2239                     y + radius * sin(angle * i), state)) {
2240       ++samp_count;
2241     }
2242   }
2243   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2244     out->channel[ch] = 0;
2245     for (i = 0; i < samp_count; ++i) {
2246       out->channel[ch] += work[i].channel[ch];
2247     }
2248     /* we divide by maxsamples rather than samp_count since if there's
2249        only one valid sample it should be mostly transparent */
2250     out->channel[ch] /= maxsamples;
2251   }
2252   return samp_count;
2253 }
2254
2255 /*
2256 =item fount_r_none(v)
2257
2258 Implements no repeats.  Simply clamps the fill value.
2259
2260 =cut
2261 */
2262 static double
2263 fount_r_none(double v) {
2264   return v < 0 ? 0 : v > 1 ? 1 : v;
2265 }
2266
2267 /*
2268 =item fount_r_sawtooth(v)
2269
2270 Implements sawtooth repeats.  Clamps negative values and uses fmod()
2271 on others.
2272
2273 =cut
2274 */
2275 static double
2276 fount_r_sawtooth(double v) {
2277   return v < 0 ? 0 : fmod(v, 1.0);
2278 }
2279
2280 /*
2281 =item fount_r_triangle(v)
2282
2283 Implements triangle repeats.  Clamps negative values, uses fmod to get
2284 a range 0 through 2 and then adjusts values > 1.
2285
2286 =cut
2287 */
2288 static double
2289 fount_r_triangle(double v) {
2290   if (v < 0)
2291     return 0;
2292   else {
2293     v = fmod(v, 2.0);
2294     return v > 1.0 ? 2.0 - v : v;
2295   }
2296 }
2297
2298 /*
2299 =item fount_r_saw_both(v)
2300
2301 Implements sawtooth repeats in the both postive and negative directions.
2302
2303 Adjusts the value to be postive and then just uses fmod().
2304
2305 =cut
2306 */
2307 static double
2308 fount_r_saw_both(double v) {
2309   if (v < 0)
2310     v += 1+(int)(-v);
2311   return fmod(v, 1.0);
2312 }
2313
2314 /*
2315 =item fount_r_tri_both(v)
2316
2317 Implements triangle repeats in the both postive and negative directions.
2318
2319 Uses fmod on the absolute value, and then adjusts values > 1.
2320
2321 =cut
2322 */
2323 static double
2324 fount_r_tri_both(double v) {
2325   v = fmod(fabs(v), 2.0);
2326   return v > 1.0 ? 2.0 - v : v;
2327 }
2328
2329 /*
2330 =item fill_fountf(fill, x, y, width, channels, data)
2331
2332 The fill function for fountain fills.
2333
2334 =cut
2335 */
2336 static void
2337 fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width,
2338             int channels, i_fcolor *data) {
2339   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2340
2341   while (width--) {
2342     i_fcolor c;
2343     int got_one;
2344
2345     if (f->state.ssfunc)
2346       got_one = f->state.ssfunc(&c, x, y, &f->state);
2347     else
2348       got_one = fount_getat(&c, x, y, &f->state);
2349
2350     *data++ = c;
2351
2352     ++x;
2353   }
2354 }
2355
2356 /*
2357 =item fount_fill_destroy(fill)
2358
2359 =cut
2360 */
2361 static void
2362 fount_fill_destroy(i_fill_t *fill) {
2363   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2364   fount_finish_state(&f->state);
2365 }
2366
2367 /*
2368 =back
2369
2370 =head1 AUTHOR
2371
2372 Arnar M. Hrafnkelsson <addi@umich.edu>
2373
2374 Tony Cook <tony@develop-help.com> (i_fountain())
2375
2376 =head1 SEE ALSO
2377
2378 Imager(3)
2379
2380 =cut
2381 */