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