filters.c

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