filters.im

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