filters.c

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