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