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 _MSC_VER
 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
 941   float *fdist;
 942
 943   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 944
 945   for(p = 0; p<num; p++) {
 946     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
 947     ICL_info(&ival[p]);
 948   }
 949
 950   fdist = mymalloc( sizeof(float) * num );
 951
 952   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 953     float cs = 0;
 954     float csd = 0;
 955     for(p = 0; p<num; p++) {
 956       int xd    = x-xo[p];
 957       int yd    = y-yo[p];
 958       switch (dmeasure) {
 959       case 0: /* euclidean */
 960         fdist[p]  = sqrt(xd*xd + yd*yd); /* euclidean distance */
 961         break;
 962       case 1: /* euclidean squared */
 963         fdist[p]  = xd*xd + yd*yd; /* euclidean distance */
 964         break;
 965       case 2: /* euclidean squared */
 966         fdist[p]  = i_max(xd*xd, yd*yd); /* manhattan distance */
 967         break;
 968       default:
 969         m_fatal(3,"i_gradgen: Unknown distance measure\n");
 970       }
 971       cs += fdist[p];
 972     }
 973
 974     csd = 1/((num-1)*cs);
 975
 976     for(p = 0; p<num; p++) fdist[p] = (cs-fdist[p])*csd;
 977
 978     for(ch = 0; ch<channels; ch++) {
 979       int tres = 0;
 980       for(p = 0; p<num; p++) tres += ival[p].channel[ch] * fdist[p];
 981       val.channel[ch] = saturate(tres);
 982     }
 983     i_ppix(im, x, y, &val);
 984   }
 985   myfree(fdist);
 986
 987 }
 988
 989 void
 990 i_nearest_color_foo(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
 991
 992   int p, x, y;
 993   int xsize    = im->xsize;
 994   int ysize    = im->ysize;
 995
 996   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 997
 998   for(p = 0; p<num; p++) {
 999     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
1000     ICL_info(&ival[p]);
1001   }
1002
1003   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1004     int   midx    = 0;
1005     float mindist = 0;
1006     float curdist = 0;
1007
1008     int xd        = x-xo[0];
1009     int yd        = y-yo[0];
1010
1011     switch (dmeasure) {
1012     case 0: /* euclidean */
1013       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1014       break;
1015     case 1: /* euclidean squared */
1016       mindist = xd*xd + yd*yd; /* euclidean distance */
1017       break;
1018     case 2: /* euclidean squared */
1019       mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
1020       break;
1021     default:
1022       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1023     }
1024
1025     for(p = 1; p<num; p++) {
1026       xd = x-xo[p];
1027       yd = y-yo[p];
1028       switch (dmeasure) {
1029       case 0: /* euclidean */
1030         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1031         break;
1032       case 1: /* euclidean squared */
1033         curdist = xd*xd + yd*yd; /* euclidean distance */
1034         break;
1035       case 2: /* euclidean squared */
1036         curdist = 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       if (curdist < mindist) {
1042         mindist = curdist;
1043         midx = p;
1044       }
1045     }
1046     i_ppix(im, x, y, &ival[midx]);
1047   }
1048 }
1049
1050 void
1051 i_nearest_color(i_img *im, int num, int *xo, int *yo, i_color *oval, int dmeasure) {
1052   i_color *ival;
1053   float *tval;
1054   float c1, c2;
1055   i_color val;
1056   int p, x, y, ch;
1057   int xsize    = im->xsize;
1058   int ysize    = im->ysize;
1059   int *cmatch;
1060
1061   mm_log((1,"i_nearest_color(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, oval, dmeasure));
1062
1063   tval   = mymalloc( sizeof(float)*num*im->channels );
1064   ival   = mymalloc( sizeof(i_color)*num );
1065   cmatch = mymalloc( sizeof(int)*num     );
1066
1067   for(p = 0; p<num; p++) {
1068     for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
1069     cmatch[p] = 0;
1070   }
1071
1072
1073   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1074     int   midx    = 0;
1075     float mindist = 0;
1076     float curdist = 0;
1077
1078     int xd        = x-xo[0];
1079     int yd        = y-yo[0];
1080
1081     switch (dmeasure) {
1082     case 0: /* euclidean */
1083       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1084       break;
1085     case 1: /* euclidean squared */
1086       mindist = xd*xd + yd*yd; /* euclidean distance */
1087       break;
1088     case 2: /* euclidean squared */
1089       mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
1090       break;
1091     default:
1092       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1093     }
1094
1095     for(p = 1; p<num; p++) {
1096       xd = x-xo[p];
1097       yd = y-yo[p];
1098       switch (dmeasure) {
1099       case 0: /* euclidean */
1100         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1101         break;
1102       case 1: /* euclidean squared */
1103         curdist = xd*xd + yd*yd; /* euclidean distance */
1104         break;
1105       case 2: /* euclidean squared */
1106         curdist = i_max(xd*xd, yd*yd); /* manhattan distance */
1107         break;
1108       default:
1109         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1110       }
1111       if (curdist < mindist) {
1112         mindist = curdist;
1113         midx = p;
1114       }
1115     }
1116
1117     cmatch[midx]++;
1118     i_gpix(im, x, y, &val);
1119     c2 = 1.0/(float)(cmatch[midx]);
1120     c1 = 1.0-c2;
1121
1122     for(ch = 0; ch<im->channels; ch++)
1123       tval[midx*im->channels + ch] = c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
1124
1125
1126   }
1127
1128   for(p = 0; p<num; p++) for(ch = 0; ch<im->channels; ch++) ival[p].channel[ch] = tval[p*im->channels + ch];
1129
1130   i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
1131 }
1132
1133 /*
1134 =item i_unsharp_mask(im, stddev, scale)
1135
1136 Perform an usharp mask, which is defined as subtracting the blurred
1137 image from double the original.
1138
1139 =cut
1140 */
1141 void i_unsharp_mask(i_img *im, double stddev, double scale) {
1142   i_img copy;
1143   int x, y, ch;
1144
1145   if (scale < 0)
1146     return;
1147   /* it really shouldn't ever be more than 1.0, but maybe ... */
1148   if (scale > 100)
1149     scale = 100;
1150
1151   i_copy(&copy, im);
1152   i_gaussian(&copy, stddev);
1153   if (im->bits == i_8_bits) {
1154     i_color *blur = mymalloc(im->xsize * sizeof(i_color) * 2);
1155     i_color *out = blur + im->xsize;
1156
1157     for (y = 0; y < im->ysize; ++y) {
1158       i_glin(&copy, 0, copy.xsize, y, blur);
1159       i_glin(im, 0, im->xsize, y, out);
1160       for (x = 0; x < im->xsize; ++x) {
1161         for (ch = 0; ch < im->channels; ++ch) {
1162           /*int temp = out[x].channel[ch] +
1163             scale * (out[x].channel[ch] - blur[x].channel[ch]);*/
1164           int temp = out[x].channel[ch] * 2 - blur[x].channel[ch];
1165           if (temp < 0)
1166             temp = 0;
1167           else if (temp > 255)
1168             temp = 255;
1169           out[x].channel[ch] = temp;
1170         }
1171       }
1172       i_plin(im, 0, im->xsize, y, out);
1173     }
1174
1175     myfree(blur);
1176   }
1177   else {
1178     i_fcolor *blur = mymalloc(im->xsize * sizeof(i_fcolor) * 2);
1179     i_fcolor *out = blur + im->xsize;
1180
1181     for (y = 0; y < im->ysize; ++y) {
1182       i_glinf(&copy, 0, copy.xsize, y, blur);
1183       i_glinf(im, 0, im->xsize, y, out);
1184       for (x = 0; x < im->xsize; ++x) {
1185         for (ch = 0; ch < im->channels; ++ch) {
1186           double temp = out[x].channel[ch] +
1187             scale * (out[x].channel[ch] - blur[x].channel[ch]);
1188           if (temp < 0)
1189             temp = 0;
1190           else if (temp > 1.0)
1191             temp = 1.0;
1192           out[x].channel[ch] = temp;
1193         }
1194       }
1195       i_plinf(im, 0, im->xsize, y, out);
1196     }
1197
1198     myfree(blur);
1199   }
1200   i_img_exorcise(&copy);
1201 }
1202
1203 /*
1204 =item i_diff_image(im1, im2, mindiff)
1205
1206 Creates a new image that is transparent, except where the pixel in im2
1207 is different from im1, where it is the pixel from im2.
1208
1209 The samples must differ by at least mindiff to be considered different.
1210
1211 =cut
1212 */
1213
1214 i_img *
1215 i_diff_image(i_img *im1, i_img *im2, int mindiff) {
1216   i_img *out;
1217   int outchans, diffchans;
1218   int xsize, ysize;
1219   i_img temp;
1220
1221   i_clear_error();
1222   if (im1->channels != im2->channels) {
1223     i_push_error(0, "different number of channels");
1224     return NULL;
1225   }
1226
1227   outchans = diffchans = im1->channels;
1228   if (outchans == 1 || outchans == 3)
1229     ++outchans;
1230
1231   xsize = i_min(im1->xsize, im2->xsize);
1232   ysize = i_min(im1->ysize, im2->ysize);
1233
1234   out = i_sametype_chans(im1, xsize, ysize, outchans);
1235
1236   if (im1->bits == i_8_bits && im2->bits == i_8_bits) {
1237     i_color *line1 = mymalloc(2 * xsize * sizeof(*line1));
1238     i_color *line2 = line1 + xsize;
1239     i_color empty;
1240     int x, y, ch;
1241
1242     for (ch = 0; ch < MAXCHANNELS; ++ch)
1243       empty.channel[ch] = 0;
1244
1245     for (y = 0; y < ysize; ++y) {
1246       i_glin(im1, 0, xsize, y, line1);
1247       i_glin(im2, 0, xsize, y, line2);
1248       if (outchans != diffchans) {
1249         /* give the output an alpha channel since it doesn't have one */
1250         for (x = 0; x < xsize; ++x)
1251           line2[x].channel[diffchans] = 255;
1252       }
1253       for (x = 0; x < xsize; ++x) {
1254         int diff = 0;
1255         for (ch = 0; ch < diffchans; ++ch) {
1256           if (line1[x].channel[ch] != line2[x].channel[ch]
1257               && abs(line1[x].channel[ch] - line2[x].channel[ch]) > mindiff) {
1258             diff = 1;
1259             break;
1260           }
1261         }
1262         if (!diff)
1263           line2[x] = empty;
1264       }
1265       i_plin(out, 0, xsize, y, line2);
1266     }
1267     myfree(line1);
1268   }
1269   else {
1270     i_fcolor *line1 = mymalloc(2 * xsize * sizeof(*line1));
1271     i_fcolor *line2 = line1 + xsize;
1272     i_fcolor empty;
1273     int x, y, ch;
1274     double dist = mindiff / 255;
1275
1276     for (ch = 0; ch < MAXCHANNELS; ++ch)
1277       empty.channel[ch] = 0;
1278
1279     for (y = 0; y < ysize; ++y) {
1280       i_glinf(im1, 0, xsize, y, line1);
1281       i_glinf(im2, 0, xsize, y, line2);
1282       if (outchans != diffchans) {
1283         /* give the output an alpha channel since it doesn't have one */
1284         for (x = 0; x < xsize; ++x)
1285           line2[x].channel[diffchans] = 1.0;
1286       }
1287       for (x = 0; x < xsize; ++x) {
1288         int diff = 0;
1289         for (ch = 0; ch < diffchans; ++ch) {
1290           if (line1[x].channel[ch] != line2[x].channel[ch]
1291               && abs(line1[x].channel[ch] - line2[x].channel[ch]) > dist) {
1292             diff = 1;
1293             break;
1294           }
1295         }
1296         if (!diff)
1297           line2[x] = empty;
1298       }
1299       i_plinf(out, 0, xsize, y, line2);
1300     }
1301     myfree(line1);
1302   }
1303
1304   return out;
1305 }
1306
1307 struct fount_state;
1308 static double linear_fount_f(double x, double y, struct fount_state *state);
1309 static double bilinear_fount_f(double x, double y, struct fount_state *state);
1310 static double radial_fount_f(double x, double y, struct fount_state *state);
1311 static double square_fount_f(double x, double y, struct fount_state *state);
1312 static double revolution_fount_f(double x, double y,
1313                                  struct fount_state *state);
1314 static double conical_fount_f(double x, double y, struct fount_state *state);
1315
1316 typedef double (*fount_func)(double, double, struct fount_state *);
1317 static fount_func fount_funcs[] =
1318 {
1319   linear_fount_f,
1320   bilinear_fount_f,
1321   radial_fount_f,
1322   square_fount_f,
1323   revolution_fount_f,
1324   conical_fount_f,
1325 };
1326
1327 static double linear_interp(double pos, i_fountain_seg *seg);
1328 static double sine_interp(double pos, i_fountain_seg *seg);
1329 static double sphereup_interp(double pos, i_fountain_seg *seg);
1330 static double spheredown_interp(double pos, i_fountain_seg *seg);
1331 typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
1332 static fount_interp fount_interps[] =
1333 {
1334   linear_interp,
1335   linear_interp,
1336   sine_interp,
1337   sphereup_interp,
1338   spheredown_interp,
1339 };
1340
1341 static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1342 static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1343 static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1344 typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
1345 static fount_cinterp fount_cinterps[] =
1346 {
1347   direct_cinterp,
1348   hue_up_cinterp,
1349   hue_down_cinterp,
1350 };
1351
1352 typedef double (*fount_repeat)(double v);
1353 static double fount_r_none(double v);
1354 static double fount_r_sawtooth(double v);
1355 static double fount_r_triangle(double v);
1356 static double fount_r_saw_both(double v);
1357 static double fount_r_tri_both(double v);
1358 static fount_repeat fount_repeats[] =
1359 {
1360   fount_r_none,
1361   fount_r_sawtooth,
1362   fount_r_triangle,
1363   fount_r_saw_both,
1364   fount_r_tri_both,
1365 };
1366
1367 static int simple_ssample(i_fcolor *out, double x, double y,
1368                            struct fount_state *state);
1369 static int random_ssample(i_fcolor *out, double x, double y,
1370                            struct fount_state *state);
1371 static int circle_ssample(i_fcolor *out, double x, double y,
1372                            struct fount_state *state);
1373 typedef int (*fount_ssample)(i_fcolor *out, double x, double y,
1374                               struct fount_state *state);
1375 static fount_ssample fount_ssamples[] =
1376 {
1377   NULL,
1378   simple_ssample,
1379   random_ssample,
1380   circle_ssample,
1381 };
1382
1383 static int
1384 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
1385
1386 /*
1387   Keep state information used by each type of fountain fill
1388 */
1389 struct fount_state {
1390   /* precalculated for the equation of the line perpendicular to the line AB */
1391   double lA, lB, lC;
1392   double AB;
1393   double sqrtA2B2;
1394   double mult;
1395   double cos;
1396   double sin;
1397   double theta;
1398   int xa, ya;
1399   void *ssample_data;
1400   fount_func ffunc;
1401   fount_repeat rpfunc;
1402   fount_ssample ssfunc;
1403   double parm;
1404   i_fountain_seg *segs;
1405   int count;
1406 };
1407
1408 static void
1409 fount_init_state(struct fount_state *state, double xa, double ya,
1410                  double xb, double yb, i_fountain_type type,
1411                  i_fountain_repeat repeat, int combine, int super_sample,
1412                  double ssample_param, int count, i_fountain_seg *segs);
1413
1414 static void
1415 fount_finish_state(struct fount_state *state);
1416
1417 #define EPSILON (1e-6)
1418
1419 /*
1420 =item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1421
1422 Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1423
1424 I<type> controls how the reference points are used:
1425
1426 =over
1427
1428 =item i_ft_linear
1429
1430 linear, where A is 0 and B is 1.
1431
1432 =item i_ft_bilinear
1433
1434 linear in both directions from A.
1435
1436 =item i_ft_radial
1437
1438 circular, where A is the centre of the fill, and B is a point
1439 on the radius.
1440
1441 =item i_ft_radial_square
1442
1443 where A is the centre of the fill and B is the centre of
1444 one side of the square.
1445
1446 =item i_ft_revolution
1447
1448 where A is the centre of the fill and B defines the 0/1.0
1449 angle of the fill.
1450
1451 =item i_ft_conical
1452
1453 similar to i_ft_revolution, except that the revolution goes in both
1454 directions
1455
1456 =back
1457
1458 I<repeat> can be one of:
1459
1460 =over
1461
1462 =item i_fr_none
1463
1464 values < 0 are treated as zero, values > 1 are treated as 1.
1465
1466 =item i_fr_sawtooth
1467
1468 negative values are treated as 0, positive values are modulo 1.0
1469
1470 =item i_fr_triangle
1471
1472 negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1473 mod 1.0), otherwise the same as for sawtooth.
1474
1475 =item i_fr_saw_both
1476
1477 like i_fr_sawtooth, except that the sawtooth pattern repeats into
1478 negative values.
1479
1480 =item i_fr_tri_both
1481
1482 Like i_fr_triangle, except that negative values are handled as their
1483 absolute values.
1484
1485 =back
1486
1487 If combine is non-zero then non-opaque values are combined with the
1488 underlying color.
1489
1490 I<super_sample> controls super sampling, if any.  At some point I'll
1491 probably add a adaptive super-sampler.  Current possible values are:
1492
1493 =over
1494
1495 =item i_fts_none
1496
1497 No super-sampling is done.
1498
1499 =item i_fts_grid
1500
1501 A square grid of points withing the pixel are sampled.
1502
1503 =item i_fts_random
1504
1505 Random points within the pixel are sampled.
1506
1507 =item i_fts_circle
1508
1509 Points on the radius of a circle are sampled.  This produces fairly
1510 good results, but is fairly slow since sin() and cos() are evaluated
1511 for each point.
1512
1513 =back
1514
1515 I<ssample_param> is intended to be roughly the number of points
1516 sampled within the pixel.
1517
1518 I<count> and I<segs> define the segments of the fill.
1519
1520 =cut
1521
1522 */
1523
1524 void
1525 i_fountain(i_img *im, double xa, double ya, double xb, double yb,
1526            i_fountain_type type, i_fountain_repeat repeat,
1527            int combine, int super_sample, double ssample_param,
1528            int count, i_fountain_seg *segs) {
1529   struct fount_state state;
1530   int x, y;
1531   i_fcolor *line = mymalloc(sizeof(i_fcolor) * im->xsize);
1532   i_fcolor *work = NULL;
1533
1534   i_fountain_seg *my_segs;
1535   i_fill_combine_f combine_func = NULL;
1536   i_fill_combinef_f combinef_func = NULL;
1537
1538   i_get_combine(combine, &combine_func, &combinef_func);
1539   if (combinef_func)
1540     work = mymalloc(sizeof(i_fcolor) * im->xsize);
1541
1542   fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
1543                    super_sample, ssample_param, count, segs);
1544   my_segs = state.segs;
1545
1546   for (y = 0; y < im->ysize; ++y) {
1547     i_glinf(im, 0, im->xsize, y, line);
1548     for (x = 0; x < im->xsize; ++x) {
1549       i_fcolor c;
1550       int got_one;
1551       if (super_sample == i_fts_none)
1552         got_one = fount_getat(&c, x, y, &state);
1553       else
1554         got_one = state.ssfunc(&c, x, y, &state);
1555       if (got_one) {
1556         if (combine)
1557           work[x] = c;
1558         else
1559           line[x] = c;
1560       }
1561     }
1562     if (combine)
1563       combinef_func(line, work, im->channels, im->xsize);
1564     i_plinf(im, 0, im->xsize, y, line);
1565   }
1566   fount_finish_state(&state);
1567   if (work) myfree(work);
1568   myfree(line);
1569 }
1570
1571 typedef struct {
1572   i_fill_t base;
1573   struct fount_state state;
1574 } i_fill_fountain_t;
1575
1576 static void
1577 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
1578             i_fcolor *data);
1579 static void
1580 fount_fill_destroy(i_fill_t *fill);
1581
1582 /*
1583 =item i_new_fount(xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1584
1585 Creates a new general fill which fills with a fountain fill.
1586
1587 =cut
1588 */
1589
1590 i_fill_t *
1591 i_new_fill_fount(double xa, double ya, double xb, double yb,
1592                  i_fountain_type type, i_fountain_repeat repeat,
1593                  int combine, int super_sample, double ssample_param,
1594                  int count, i_fountain_seg *segs) {
1595   i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1596
1597   fill->base.fill_with_color = NULL;
1598   fill->base.fill_with_fcolor = fill_fountf;
1599   fill->base.destroy = fount_fill_destroy;
1600   if (combine)
1601     i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1602   else {
1603     fill->base.combine = NULL;
1604     fill->base.combinef = NULL;
1605   }
1606   fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1607                    super_sample, ssample_param, count, segs);
1608
1609   return &fill->base;
1610 }
1611
1612 /*
1613 =back
1614
1615 =head1 INTERNAL FUNCTIONS
1616
1617 =over
1618
1619 =item fount_init_state(...)
1620
1621 Used by both the fountain fill filter and the fountain fill.
1622
1623 =cut
1624 */
1625
1626 static void
1627 fount_init_state(struct fount_state *state, double xa, double ya,
1628                  double xb, double yb, i_fountain_type type,
1629                  i_fountain_repeat repeat, int combine, int super_sample,
1630                  double ssample_param, int count, i_fountain_seg *segs) {
1631   int i, j;
1632   i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count);
1633   /*int have_alpha = im->channels == 2 || im->channels == 4;*/
1634
1635   memset(state, 0, sizeof(*state));
1636   /* we keep a local copy that we can adjust for speed */
1637   for (i = 0; i < count; ++i) {
1638     i_fountain_seg *seg = my_segs + i;
1639
1640     *seg = segs[i];
1641     if (seg->type < 0 || seg->type >= i_fst_end)
1642       seg->type = i_fst_linear;
1643     if (seg->color < 0 || seg->color >= i_fc_end)
1644       seg->color = i_fc_direct;
1645     if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1646       /* so we don't have to translate to HSV on each request, do it here */
1647       for (j = 0; j < 2; ++j) {
1648         i_rgb_to_hsvf(seg->c+j);
1649       }
1650       if (seg->color == i_fc_hue_up) {
1651         if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1652           seg->c[1].channel[0] += 1.0;
1653       }
1654       else {
1655         if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1656           seg->c[0].channel[0] += 1.0;
1657       }
1658     }
1659     /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1660            seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1661            seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1662            seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1663            seg->c[1].channel[3], seg->type, seg->color);*/
1664
1665   }
1666
1667   /* initialize each engine */
1668   /* these are so common ... */
1669   state->lA = xb - xa;
1670   state->lB = yb - ya;
1671   state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1672   state->xa = xa;
1673   state->ya = ya;
1674   switch (type) {
1675   default:
1676     type = i_ft_linear; /* make the invalid value valid */
1677   case i_ft_linear:
1678   case i_ft_bilinear:
1679     state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1680     state->mult = 1;
1681     state->mult = 1/linear_fount_f(xb, yb, state);
1682     break;
1683
1684   case i_ft_radial:
1685     state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1686                              + (double)(yb-ya)*(yb-ya));
1687     break;
1688
1689   case i_ft_radial_square:
1690     state->cos = state->lA / state->AB;
1691     state->sin = state->lB / state->AB;
1692     state->mult = 1.0 / state->AB;
1693     break;
1694
1695   case i_ft_revolution:
1696     state->theta = atan2(yb-ya, xb-xa);
1697     state->mult = 1.0 / (PI * 2);
1698     break;
1699
1700   case i_ft_conical:
1701     state->theta = atan2(yb-ya, xb-xa);
1702     state->mult = 1.0 / PI;
1703     break;
1704   }
1705   state->ffunc = fount_funcs[type];
1706   if (super_sample < 0
1707       || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
1708     super_sample = 0;
1709   }
1710   state->ssample_data = NULL;
1711   switch (super_sample) {
1712   case i_fts_grid:
1713     ssample_param = floor(0.5 + sqrt(ssample_param));
1714     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param);
1715     break;
1716
1717   case i_fts_random:
1718   case i_fts_circle:
1719     ssample_param = floor(0.5+ssample_param);
1720     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
1721     break;
1722   }
1723   state->parm = ssample_param;
1724   state->ssfunc = fount_ssamples[super_sample];
1725   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1726     repeat = 0;
1727   state->rpfunc = fount_repeats[repeat];
1728   state->segs = my_segs;
1729   state->count = count;
1730 }
1731
1732 static void
1733 fount_finish_state(struct fount_state *state) {
1734   if (state->ssample_data)
1735     myfree(state->ssample_data);
1736   myfree(state->segs);
1737 }
1738
1739
1740 /*
1741 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
1742
1743 Evaluates the fountain fill at the given point.
1744
1745 This is called by both the non-super-sampling and super-sampling code.
1746
1747 You might think that it would make sense to sample the fill parameter
1748 instead, and combine those, but this breaks badly.
1749
1750 =cut
1751 */
1752
1753 static int
1754 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
1755   double v = (state->rpfunc)((state->ffunc)(x, y, state));
1756   int i;
1757
1758   i = 0;
1759   while (i < state->count
1760          && (v < state->segs[i].start || v > state->segs[i].end)) {
1761     ++i;
1762   }
1763   if (i < state->count) {
1764     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
1765     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
1766     return 1;
1767   }
1768   else
1769     return 0;
1770 }
1771
1772 /*
1773 =item linear_fount_f(x, y, state)
1774
1775 Calculate the fill parameter for a linear fountain fill.
1776
1777 Uses the point to line distance function, with some precalculation
1778 done in i_fountain().
1779
1780 =cut
1781 */
1782 static double
1783 linear_fount_f(double x, double y, struct fount_state *state) {
1784   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
1785 }
1786
1787 /*
1788 =item bilinear_fount_f(x, y, state)
1789
1790 Calculate the fill parameter for a bi-linear fountain fill.
1791
1792 =cut
1793 */
1794 static double
1795 bilinear_fount_f(double x, double y, struct fount_state *state) {
1796   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
1797 }
1798
1799 /*
1800 =item radial_fount_f(x, y, state)
1801
1802 Calculate the fill parameter for a radial fountain fill.
1803
1804 Simply uses the distance function.
1805
1806 =cut
1807  */
1808 static double
1809 radial_fount_f(double x, double y, struct fount_state *state) {
1810   return sqrt((double)(state->xa-x)*(state->xa-x)
1811               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
1812 }
1813
1814 /*
1815 =item square_fount_f(x, y, state)
1816
1817 Calculate the fill parameter for a square fountain fill.
1818
1819 Works by rotating the reference co-ordinate around the centre of the
1820 square.
1821
1822 =cut
1823 */
1824 static double
1825 square_fount_f(double x, double y, struct fount_state *state) {
1826   int xc, yc; /* centred on A */
1827   double xt, yt; /* rotated by theta */
1828   xc = x - state->xa;
1829   yc = y - state->ya;
1830   xt = fabs(xc * state->cos + yc * state->sin);
1831   yt = fabs(-xc * state->sin + yc * state->cos);
1832   return (xt > yt ? xt : yt) * state->mult;
1833 }
1834
1835 /*
1836 =item revolution_fount_f(x, y, state)
1837
1838 Calculates the fill parameter for the revolution fountain fill.
1839
1840 =cut
1841 */
1842 static double
1843 revolution_fount_f(double x, double y, struct fount_state *state) {
1844   double angle = atan2(y - state->ya, x - state->xa);
1845
1846   angle -= state->theta;
1847   if (angle < 0) {
1848     angle = fmod(angle+ PI * 4, PI*2);
1849   }
1850
1851   return angle * state->mult;
1852 }
1853
1854 /*
1855 =item conical_fount_f(x, y, state)
1856
1857 Calculates the fill parameter for the conical fountain fill.
1858
1859 =cut
1860 */
1861 static double
1862 conical_fount_f(double x, double y, struct fount_state *state) {
1863   double angle = atan2(y - state->ya, x - state->xa);
1864
1865   angle -= state->theta;
1866   if (angle < -PI)
1867     angle += PI * 2;
1868   else if (angle > PI)
1869     angle -= PI * 2;
1870
1871   return fabs(angle) * state->mult;
1872 }
1873
1874 /*
1875 =item linear_interp(pos, seg)
1876
1877 Calculates linear interpolation on the fill parameter.  Breaks the
1878 segment into 2 regions based in the I<middle> value.
1879
1880 =cut
1881 */
1882 static double
1883 linear_interp(double pos, i_fountain_seg *seg) {
1884   if (pos < seg->middle) {
1885     double len = seg->middle - seg->start;
1886     if (len < EPSILON)
1887       return 0.0;
1888     else
1889       return (pos - seg->start) / len / 2;
1890   }
1891   else {
1892     double len = seg->end - seg->middle;
1893     if (len < EPSILON)
1894       return 1.0;
1895     else
1896       return 0.5 + (pos - seg->middle) / len / 2;
1897   }
1898 }
1899
1900 /*
1901 =item sine_interp(pos, seg)
1902
1903 Calculates sine function interpolation on the fill parameter.
1904
1905 =cut
1906 */
1907 static double
1908 sine_interp(double pos, i_fountain_seg *seg) {
1909   /* I wonder if there's a simple way to smooth the transition for this */
1910   double work = linear_interp(pos, seg);
1911
1912   return (1-cos(work * PI))/2;
1913 }
1914
1915 /*
1916 =item sphereup_interp(pos, seg)
1917
1918 Calculates spherical interpolation on the fill parameter, with the cusp
1919 at the low-end.
1920
1921 =cut
1922 */
1923 static double
1924 sphereup_interp(double pos, i_fountain_seg *seg) {
1925   double work = linear_interp(pos, seg);
1926
1927   return sqrt(1.0 - (1-work) * (1-work));
1928 }
1929
1930 /*
1931 =item spheredown_interp(pos, seg)
1932
1933 Calculates spherical interpolation on the fill parameter, with the cusp
1934 at the high-end.
1935
1936 =cut
1937 */
1938 static double
1939 spheredown_interp(double pos, i_fountain_seg *seg) {
1940   double work = linear_interp(pos, seg);
1941
1942   return 1-sqrt(1.0 - work * work);
1943 }
1944
1945 /*
1946 =item direct_cinterp(out, pos, seg)
1947
1948 Calculates the fountain color based on direct scaling of the channels
1949 of the color channels.
1950
1951 =cut
1952 */
1953 static void
1954 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1955   int ch;
1956   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1957     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1958       + seg->c[1].channel[ch] * pos;
1959   }
1960 }
1961
1962 /*
1963 =item hue_up_cinterp(put, pos, seg)
1964
1965 Calculates the fountain color based on scaling a HSV value.  The hue
1966 increases as the fill parameter increases.
1967
1968 =cut
1969 */
1970 static void
1971 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1972   int ch;
1973   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1974     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1975       + seg->c[1].channel[ch] * pos;
1976   }
1977   i_hsv_to_rgbf(out);
1978 }
1979
1980 /*
1981 =item hue_down_cinterp(put, pos, seg)
1982
1983 Calculates the fountain color based on scaling a HSV value.  The hue
1984 decreases as the fill parameter increases.
1985
1986 =cut
1987 */
1988 static void
1989 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1990   int ch;
1991   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1992     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1993       + seg->c[1].channel[ch] * pos;
1994   }
1995   i_hsv_to_rgbf(out);
1996 }
1997
1998 /*
1999 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2000
2001 Simple grid-based super-sampling.
2002
2003 =cut
2004 */
2005 static int
2006 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
2007   i_fcolor *work = state->ssample_data;
2008   int dx, dy;
2009   int grid = state->parm;
2010   double base = -0.5 + 0.5 / grid;
2011   double step = 1.0 / grid;
2012   int ch, i;
2013   int samp_count = 0;
2014
2015   for (dx = 0; dx < grid; ++dx) {
2016     for (dy = 0; dy < grid; ++dy) {
2017       if (fount_getat(work+samp_count, x + base + step * dx,
2018                       y + base + step * dy, state)) {
2019         ++samp_count;
2020       }
2021     }
2022   }
2023   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2024     out->channel[ch] = 0;
2025     for (i = 0; i < samp_count; ++i) {
2026       out->channel[ch] += work[i].channel[ch];
2027     }
2028     /* we divide by 4 rather than samp_count since if there's only one valid
2029        sample it should be mostly transparent */
2030     out->channel[ch] /= grid * grid;
2031   }
2032   return samp_count;
2033 }
2034
2035 /*
2036 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2037
2038 Random super-sampling.
2039
2040 =cut
2041 */
2042 static int
2043 random_ssample(i_fcolor *out, double x, double y,
2044                struct fount_state *state) {
2045   i_fcolor *work = state->ssample_data;
2046   int i, ch;
2047   int maxsamples = state->parm;
2048   double rand_scale = 1.0 / RAND_MAX;
2049   int samp_count = 0;
2050   for (i = 0; i < maxsamples; ++i) {
2051     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
2052                     y - 0.5 + rand() * rand_scale, state)) {
2053       ++samp_count;
2054     }
2055   }
2056   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2057     out->channel[ch] = 0;
2058     for (i = 0; i < samp_count; ++i) {
2059       out->channel[ch] += work[i].channel[ch];
2060     }
2061     /* we divide by maxsamples rather than samp_count since if there's
2062        only one valid sample it should be mostly transparent */
2063     out->channel[ch] /= maxsamples;
2064   }
2065   return samp_count;
2066 }
2067
2068 /*
2069 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2070
2071 Super-sampling around the circumference of a circle.
2072
2073 I considered saving the sin()/cos() values and transforming step-size
2074 around the circle, but that's inaccurate, though it may not matter
2075 much.
2076
2077 =cut
2078  */
2079 static int
2080 circle_ssample(i_fcolor *out, double x, double y,
2081                struct fount_state *state) {
2082   i_fcolor *work = state->ssample_data;
2083   int i, ch;
2084   int maxsamples = state->parm;
2085   double angle = 2 * PI / maxsamples;
2086   double radius = 0.3; /* semi-random */
2087   int samp_count = 0;
2088   for (i = 0; i < maxsamples; ++i) {
2089     if (fount_getat(work+samp_count, x + radius * cos(angle * i),
2090                     y + radius * sin(angle * i), state)) {
2091       ++samp_count;
2092     }
2093   }
2094   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2095     out->channel[ch] = 0;
2096     for (i = 0; i < samp_count; ++i) {
2097       out->channel[ch] += work[i].channel[ch];
2098     }
2099     /* we divide by maxsamples rather than samp_count since if there's
2100        only one valid sample it should be mostly transparent */
2101     out->channel[ch] /= maxsamples;
2102   }
2103   return samp_count;
2104 }
2105
2106 /*
2107 =item fount_r_none(v)
2108
2109 Implements no repeats.  Simply clamps the fill value.
2110
2111 =cut
2112 */
2113 static double
2114 fount_r_none(double v) {
2115   return v < 0 ? 0 : v > 1 ? 1 : v;
2116 }
2117
2118 /*
2119 =item fount_r_sawtooth(v)
2120
2121 Implements sawtooth repeats.  Clamps negative values and uses fmod()
2122 on others.
2123
2124 =cut
2125 */
2126 static double
2127 fount_r_sawtooth(double v) {
2128   return v < 0 ? 0 : fmod(v, 1.0);
2129 }
2130
2131 /*
2132 =item fount_r_triangle(v)
2133
2134 Implements triangle repeats.  Clamps negative values, uses fmod to get
2135 a range 0 through 2 and then adjusts values > 1.
2136
2137 =cut
2138 */
2139 static double
2140 fount_r_triangle(double v) {
2141   if (v < 0)
2142     return 0;
2143   else {
2144     v = fmod(v, 2.0);
2145     return v > 1.0 ? 2.0 - v : v;
2146   }
2147 }
2148
2149 /*
2150 =item fount_r_saw_both(v)
2151
2152 Implements sawtooth repeats in the both postive and negative directions.
2153
2154 Adjusts the value to be postive and then just uses fmod().
2155
2156 =cut
2157 */
2158 static double
2159 fount_r_saw_both(double v) {
2160   if (v < 0)
2161     v += 1+(int)(-v);
2162   return fmod(v, 1.0);
2163 }
2164
2165 /*
2166 =item fount_r_tri_both(v)
2167
2168 Implements triangle repeats in the both postive and negative directions.
2169
2170 Uses fmod on the absolute value, and then adjusts values > 1.
2171
2172 =cut
2173 */
2174 static double
2175 fount_r_tri_both(double v) {
2176   v = fmod(fabs(v), 2.0);
2177   return v > 1.0 ? 2.0 - v : v;
2178 }
2179
2180 /*
2181 =item fill_fountf(fill, x, y, width, channels, data)
2182
2183 The fill function for fountain fills.
2184
2185 =cut
2186 */
2187 static void
2188 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
2189             i_fcolor *data) {
2190   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2191
2192   while (width--) {
2193     i_fcolor c;
2194     int got_one;
2195
2196     if (f->state.ssfunc)
2197       got_one = f->state.ssfunc(&c, x, y, &f->state);
2198     else
2199       got_one = fount_getat(&c, x, y, &f->state);
2200
2201     *data++ = c;
2202
2203     ++x;
2204   }
2205 }
2206
2207 /*
2208 =item fount_fill_destroy(fill)
2209
2210 =cut
2211 */
2212 static void
2213 fount_fill_destroy(i_fill_t *fill) {
2214   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2215   fount_finish_state(&f->state);
2216 }
2217
2218 /*
2219 =back
2220
2221 =head1 AUTHOR
2222
2223 Arnar M. Hrafnkelsson <addi@umich.edu>
2224
2225 Tony Cook <tony@develop-help.com> (i_fountain())
2226
2227 =head1 SEE ALSO
2228
2229 Imager(3)
2230
2231 =cut
2232 */