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   myfree(cmatch);
1337   myfree(ival);
1338   myfree(tval);
1339
1340   return 1;
1341 }
1342
1343 /*
1344 =item i_unsharp_mask(im, stddev, scale)
1345
1346 Perform an usharp mask, which is defined as subtracting the blurred
1347 image from double the original.
1348
1349 =cut
1350 */
1351
1352 void
1353 i_unsharp_mask(i_img *im, double stddev, double scale) {
1354   i_img *copy;
1355   i_img_dim x, y;
1356   int ch;
1357
1358   if (scale < 0)
1359     return;
1360   /* it really shouldn't ever be more than 1.0, but maybe ... */
1361   if (scale > 100)
1362     scale = 100;
1363
1364   copy = i_copy(im);
1365   i_gaussian(copy, stddev);
1366   if (im->bits == i_8_bits) {
1367     i_color *blur = mymalloc(im->xsize * sizeof(i_color)); /* checked 17feb2005 tonyc */
1368     i_color *out = mymalloc(im->xsize * sizeof(i_color)); /* checked 17feb2005 tonyc */
1369
1370     for (y = 0; y < im->ysize; ++y) {
1371       i_glin(copy, 0, copy->xsize, y, blur);
1372       i_glin(im, 0, im->xsize, y, out);
1373       for (x = 0; x < im->xsize; ++x) {
1374         for (ch = 0; ch < im->channels; ++ch) {
1375           /*int temp = out[x].channel[ch] +
1376             scale * (out[x].channel[ch] - blur[x].channel[ch]);*/
1377           int temp = out[x].channel[ch] * 2 - blur[x].channel[ch];
1378           if (temp < 0)
1379             temp = 0;
1380           else if (temp > 255)
1381             temp = 255;
1382           out[x].channel[ch] = temp;
1383         }
1384       }
1385       i_plin(im, 0, im->xsize, y, out);
1386     }
1387
1388     myfree(blur);
1389     myfree(out);
1390   }
1391   else {
1392     i_fcolor *blur = mymalloc(im->xsize * sizeof(i_fcolor)); /* checked 17feb2005 tonyc */
1393     i_fcolor *out = mymalloc(im->xsize * sizeof(i_fcolor)); /* checked 17feb2005 tonyc */
1394
1395     for (y = 0; y < im->ysize; ++y) {
1396       i_glinf(copy, 0, copy->xsize, y, blur);
1397       i_glinf(im, 0, im->xsize, y, out);
1398       for (x = 0; x < im->xsize; ++x) {
1399         for (ch = 0; ch < im->channels; ++ch) {
1400           double temp = out[x].channel[ch] +
1401             scale * (out[x].channel[ch] - blur[x].channel[ch]);
1402           if (temp < 0)
1403             temp = 0;
1404           else if (temp > 1.0)
1405             temp = 1.0;
1406           out[x].channel[ch] = temp;
1407         }
1408       }
1409       i_plinf(im, 0, im->xsize, y, out);
1410     }
1411
1412     myfree(blur);
1413     myfree(out);
1414   }
1415   i_img_destroy(copy);
1416 }
1417
1418 /*
1419 =item i_diff_image(im1, im2, mindist)
1420
1421 Creates a new image that is transparent, except where the pixel in im2
1422 is different from im1, where it is the pixel from im2.
1423
1424 The samples must differ by at least mindiff to be considered different.
1425
1426 =cut
1427 */
1428
1429 i_img *
1430 i_diff_image(i_img *im1, i_img *im2, double mindist) {
1431   i_img *out;
1432   int outchans, diffchans;
1433   i_img_dim xsize, ysize;
1434   dIMCTXim(im1);
1435
1436   i_clear_error();
1437   if (im1->channels != im2->channels) {
1438     i_push_error(0, "different number of channels");
1439     return NULL;
1440   }
1441
1442   outchans = diffchans = im1->channels;
1443   if (outchans == 1 || outchans == 3)
1444     ++outchans;
1445
1446   xsize = i_min(im1->xsize, im2->xsize);
1447   ysize = i_min(im1->ysize, im2->ysize);
1448
1449   out = i_sametype_chans(im1, xsize, ysize, outchans);
1450
1451   if (im1->bits == i_8_bits && im2->bits == i_8_bits) {
1452     i_color *line1 = mymalloc(xsize * sizeof(*line1)); /* checked 17feb2005 tonyc */
1453     i_color *line2 = mymalloc(xsize * sizeof(*line1)); /* checked 17feb2005 tonyc */
1454     i_color empty;
1455     i_img_dim x, y;
1456     int ch;
1457     int imindist = (int)mindist;
1458
1459     for (ch = 0; ch < MAXCHANNELS; ++ch)
1460       empty.channel[ch] = 0;
1461
1462     for (y = 0; y < ysize; ++y) {
1463       i_glin(im1, 0, xsize, y, line1);
1464       i_glin(im2, 0, xsize, y, line2);
1465       if (outchans != diffchans) {
1466         /* give the output an alpha channel since it doesn't have one */
1467         for (x = 0; x < xsize; ++x)
1468           line2[x].channel[diffchans] = 255;
1469       }
1470       for (x = 0; x < xsize; ++x) {
1471         int diff = 0;
1472         for (ch = 0; ch < diffchans; ++ch) {
1473           if (line1[x].channel[ch] != line2[x].channel[ch]
1474               && abs(line1[x].channel[ch] - line2[x].channel[ch]) > imindist) {
1475             diff = 1;
1476             break;
1477           }
1478         }
1479         if (!diff)
1480           line2[x] = empty;
1481       }
1482       i_plin(out, 0, xsize, y, line2);
1483     }
1484     myfree(line1);
1485     myfree(line2);
1486   }
1487   else {
1488     i_fcolor *line1 = mymalloc(xsize * sizeof(*line1)); /* checked 17feb2005 tonyc */
1489     i_fcolor *line2 = mymalloc(xsize * sizeof(*line2)); /* checked 17feb2005 tonyc */
1490     i_fcolor empty;
1491     i_img_dim x, y;
1492     int ch;
1493     double dist = mindist / 255.0;
1494
1495     for (ch = 0; ch < MAXCHANNELS; ++ch)
1496       empty.channel[ch] = 0;
1497
1498     for (y = 0; y < ysize; ++y) {
1499       i_glinf(im1, 0, xsize, y, line1);
1500       i_glinf(im2, 0, xsize, y, line2);
1501       if (outchans != diffchans) {
1502         /* give the output an alpha channel since it doesn't have one */
1503         for (x = 0; x < xsize; ++x)
1504           line2[x].channel[diffchans] = 1.0;
1505       }
1506       for (x = 0; x < xsize; ++x) {
1507         int diff = 0;
1508         for (ch = 0; ch < diffchans; ++ch) {
1509           if (line1[x].channel[ch] != line2[x].channel[ch]
1510               && fabs(line1[x].channel[ch] - line2[x].channel[ch]) > dist) {
1511             diff = 1;
1512             break;
1513           }
1514         }
1515         if (!diff)
1516           line2[x] = empty;
1517       }
1518       i_plinf(out, 0, xsize, y, line2);
1519     }
1520     myfree(line1);
1521     myfree(line2);
1522   }
1523
1524   return out;
1525 }
1526
1527 /*
1528 =item i_rgbdiff_image(im1, im2)
1529
1530 Creates a new image that is black, except where the pixel in im2 is
1531 different from im1, where it is the arithmetical difference to im2 per
1532 color.
1533
1534 =cut
1535 */
1536
1537 i_img *
1538 i_rgbdiff_image(i_img *im1, i_img *im2) {
1539   i_img *out;
1540   int outchans, diffchans;
1541   i_img_dim xsize, ysize;
1542   dIMCTXim(im1);
1543
1544   i_clear_error();
1545   if (im1->channels != im2->channels) {
1546     i_push_error(0, "different number of channels");
1547     return NULL;
1548   }
1549
1550   outchans = diffchans = im1->channels;
1551   if (outchans == 2 || outchans == 4)
1552     --outchans;
1553
1554   xsize = i_min(im1->xsize, im2->xsize);
1555   ysize = i_min(im1->ysize, im2->ysize);
1556
1557   out = i_sametype_chans(im1, xsize, ysize, outchans);
1558
1559 #code im1->bits == i_8_bits && im2->bits == i_8_bits
1560     IM_COLOR *line1 = mymalloc(xsize * sizeof(*line1));
1561     IM_COLOR *line2 = mymalloc(xsize * sizeof(*line1));
1562     i_img_dim x, y;
1563     int ch;
1564
1565     for (y = 0; y < ysize; ++y) {
1566       IM_GLIN(im1, 0, xsize, y, line1);
1567       IM_GLIN(im2, 0, xsize, y, line2);
1568       for (x = 0; x < xsize; ++x) {
1569         for (ch = 0; ch < outchans; ++ch) {
1570           line2[x].channel[ch] = IM_ABS(line1[x].channel[ch] - line2[x].channel[ch]);
1571         }
1572       }
1573       IM_PLIN(out, 0, xsize, y, line2);
1574     }
1575     myfree(line1);
1576     myfree(line2);
1577 #/code
1578
1579   return out;
1580 }
1581
1582 struct fount_state;
1583 static double linear_fount_f(double x, double y, struct fount_state *state);
1584 static double bilinear_fount_f(double x, double y, struct fount_state *state);
1585 static double radial_fount_f(double x, double y, struct fount_state *state);
1586 static double square_fount_f(double x, double y, struct fount_state *state);
1587 static double revolution_fount_f(double x, double y,
1588                                  struct fount_state *state);
1589 static double conical_fount_f(double x, double y, struct fount_state *state);
1590
1591 typedef double (*fount_func)(double, double, struct fount_state *);
1592 static fount_func fount_funcs[] =
1593 {
1594   linear_fount_f,
1595   bilinear_fount_f,
1596   radial_fount_f,
1597   square_fount_f,
1598   revolution_fount_f,
1599   conical_fount_f,
1600 };
1601
1602 static double linear_interp(double pos, i_fountain_seg *seg);
1603 static double sine_interp(double pos, i_fountain_seg *seg);
1604 static double sphereup_interp(double pos, i_fountain_seg *seg);
1605 static double spheredown_interp(double pos, i_fountain_seg *seg);
1606 typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
1607 static fount_interp fount_interps[] =
1608 {
1609   linear_interp,
1610   linear_interp,
1611   sine_interp,
1612   sphereup_interp,
1613   spheredown_interp,
1614 };
1615
1616 static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1617 static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1618 static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1619 typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
1620 static fount_cinterp fount_cinterps[] =
1621 {
1622   direct_cinterp,
1623   hue_up_cinterp,
1624   hue_down_cinterp,
1625 };
1626
1627 typedef double (*fount_repeat)(double v);
1628 static double fount_r_none(double v);
1629 static double fount_r_sawtooth(double v);
1630 static double fount_r_triangle(double v);
1631 static double fount_r_saw_both(double v);
1632 static double fount_r_tri_both(double v);
1633 static fount_repeat fount_repeats[] =
1634 {
1635   fount_r_none,
1636   fount_r_sawtooth,
1637   fount_r_triangle,
1638   fount_r_saw_both,
1639   fount_r_tri_both,
1640 };
1641
1642 static int simple_ssample(i_fcolor *out, double x, double y,
1643                            struct fount_state *state);
1644 static int random_ssample(i_fcolor *out, double x, double y,
1645                            struct fount_state *state);
1646 static int circle_ssample(i_fcolor *out, double x, double y,
1647                            struct fount_state *state);
1648 typedef int (*fount_ssample)(i_fcolor *out, double x, double y,
1649                               struct fount_state *state);
1650 static fount_ssample fount_ssamples[] =
1651 {
1652   NULL,
1653   simple_ssample,
1654   random_ssample,
1655   circle_ssample,
1656 };
1657
1658 static int
1659 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
1660
1661 /*
1662   Keep state information used by each type of fountain fill
1663 */
1664 struct fount_state {
1665   /* precalculated for the equation of the line perpendicular to the line AB */
1666   double lA, lB, lC;
1667   double AB;
1668   double sqrtA2B2;
1669   double mult;
1670   double cos;
1671   double sin;
1672   double theta;
1673   i_img_dim xa, ya;
1674   void *ssample_data;
1675   fount_func ffunc;
1676   fount_repeat rpfunc;
1677   fount_ssample ssfunc;
1678   double parm;
1679   i_fountain_seg *segs;
1680   int count;
1681 };
1682
1683 static void
1684 fount_init_state(struct fount_state *state, double xa, double ya,
1685                  double xb, double yb, i_fountain_type type,
1686                  i_fountain_repeat repeat, int combine, int super_sample,
1687                  double ssample_param, int count, i_fountain_seg *segs);
1688
1689 static void
1690 fount_finish_state(struct fount_state *state);
1691
1692 #define EPSILON (1e-6)
1693
1694 /*
1695 =item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1696
1697 Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1698
1699 I<type> controls how the reference points are used:
1700
1701 =over
1702
1703 =item i_ft_linear
1704
1705 linear, where A is 0 and B is 1.
1706
1707 =item i_ft_bilinear
1708
1709 linear in both directions from A.
1710
1711 =item i_ft_radial
1712
1713 circular, where A is the centre of the fill, and B is a point
1714 on the radius.
1715
1716 =item i_ft_radial_square
1717
1718 where A is the centre of the fill and B is the centre of
1719 one side of the square.
1720
1721 =item i_ft_revolution
1722
1723 where A is the centre of the fill and B defines the 0/1.0
1724 angle of the fill.
1725
1726 =item i_ft_conical
1727
1728 similar to i_ft_revolution, except that the revolution goes in both
1729 directions
1730
1731 =back
1732
1733 I<repeat> can be one of:
1734
1735 =over
1736
1737 =item i_fr_none
1738
1739 values < 0 are treated as zero, values > 1 are treated as 1.
1740
1741 =item i_fr_sawtooth
1742
1743 negative values are treated as 0, positive values are modulo 1.0
1744
1745 =item i_fr_triangle
1746
1747 negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1748 mod 1.0), otherwise the same as for sawtooth.
1749
1750 =item i_fr_saw_both
1751
1752 like i_fr_sawtooth, except that the sawtooth pattern repeats into
1753 negative values.
1754
1755 =item i_fr_tri_both
1756
1757 Like i_fr_triangle, except that negative values are handled as their
1758 absolute values.
1759
1760 =back
1761
1762 If combine is non-zero then non-opaque values are combined with the
1763 underlying color.
1764
1765 I<super_sample> controls super sampling, if any.  At some point I'll
1766 probably add a adaptive super-sampler.  Current possible values are:
1767
1768 =over
1769
1770 =item i_fts_none
1771
1772 No super-sampling is done.
1773
1774 =item i_fts_grid
1775
1776 A square grid of points withing the pixel are sampled.
1777
1778 =item i_fts_random
1779
1780 Random points within the pixel are sampled.
1781
1782 =item i_fts_circle
1783
1784 Points on the radius of a circle are sampled.  This produces fairly
1785 good results, but is fairly slow since sin() and cos() are evaluated
1786 for each point.
1787
1788 =back
1789
1790 I<ssample_param> is intended to be roughly the number of points
1791 sampled within the pixel.
1792
1793 I<count> and I<segs> define the segments of the fill.
1794
1795 =cut
1796
1797 */
1798
1799 int
1800 i_fountain(i_img *im, double xa, double ya, double xb, double yb,
1801            i_fountain_type type, i_fountain_repeat repeat,
1802            int combine, int super_sample, double ssample_param,
1803            int count, i_fountain_seg *segs) {
1804   struct fount_state state;
1805   i_img_dim x, y;
1806   i_fcolor *line = NULL;
1807   i_fcolor *work = NULL;
1808   size_t line_bytes;
1809   i_fill_combine_f combine_func = NULL;
1810   i_fill_combinef_f combinef_func = NULL;
1811   dIMCTXim(im);
1812
1813   i_clear_error();
1814
1815   /* i_fountain() allocates floating colors even for 8-bit images,
1816      so we need to do this check */
1817   line_bytes = sizeof(i_fcolor) * im->xsize;
1818   if (line_bytes / sizeof(i_fcolor) != im->xsize) {
1819     i_push_error(0, "integer overflow calculating memory allocation");
1820     return 0;
1821   }
1822
1823   line = mymalloc(line_bytes); /* checked 17feb2005 tonyc */
1824
1825   i_get_combine(combine, &combine_func, &combinef_func);
1826   if (combinef_func) {
1827     work = mymalloc(line_bytes); /* checked 17feb2005 tonyc */
1828   }
1829
1830   fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
1831                    super_sample, ssample_param, count, segs);
1832
1833   for (y = 0; y < im->ysize; ++y) {
1834     i_glinf(im, 0, im->xsize, y, line);
1835     for (x = 0; x < im->xsize; ++x) {
1836       i_fcolor c;
1837       int got_one;
1838       if (super_sample == i_fts_none)
1839         got_one = fount_getat(&c, x, y, &state);
1840       else
1841         got_one = state.ssfunc(&c, x, y, &state);
1842       if (got_one) {
1843         if (combinef_func)
1844           work[x] = c;
1845         else
1846           line[x] = c;
1847       }
1848     }
1849     if (combinef_func)
1850       combinef_func(line, work, im->channels, im->xsize);
1851     i_plinf(im, 0, im->xsize, y, line);
1852   }
1853   fount_finish_state(&state);
1854   myfree(work);
1855   myfree(line);
1856
1857   return 1;
1858 }
1859
1860 typedef struct {
1861   i_fill_t base;
1862   struct fount_state state;
1863 } i_fill_fountain_t;
1864
1865 static void
1866 fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width, int channels,
1867             i_fcolor *data);
1868 static void
1869 fount_fill_destroy(i_fill_t *fill);
1870
1871 static i_fill_fountain_t
1872 fount_fill_proto =
1873   {
1874     {
1875       NULL,
1876       fill_fountf,
1877       fount_fill_destroy
1878     }
1879   };
1880
1881
1882 /*
1883 =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>)
1884
1885 =category Fills
1886 =synopsis fill = i_new_fill_fount(0, 0, 100, 100, i_ft_linear, i_ft_linear,
1887 =synopsis                         i_fr_triangle, 0, i_fts_grid, 9, 1, segs);
1888
1889
1890 Creates a new general fill which fills with a fountain fill.
1891
1892 =cut
1893 */
1894
1895 i_fill_t *
1896 i_new_fill_fount(double xa, double ya, double xb, double yb,
1897                  i_fountain_type type, i_fountain_repeat repeat,
1898                  int combine, int super_sample, double ssample_param,
1899                  int count, i_fountain_seg *segs) {
1900   i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1901
1902   *fill = fount_fill_proto;
1903   if (combine)
1904     i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1905   else {
1906     fill->base.combine = NULL;
1907     fill->base.combinef = NULL;
1908   }
1909   fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1910                    super_sample, ssample_param, count, segs);
1911
1912   return &fill->base;
1913 }
1914
1915 /*
1916 =back
1917
1918 =head1 INTERNAL FUNCTIONS
1919
1920 =over
1921
1922 =item fount_init_state(...)
1923
1924 Used by both the fountain fill filter and the fountain fill.
1925
1926 =cut
1927 */
1928
1929 static void
1930 fount_init_state(struct fount_state *state, double xa, double ya,
1931                  double xb, double yb, i_fountain_type type,
1932                  i_fountain_repeat repeat, int combine, int super_sample,
1933                  double ssample_param, int count, i_fountain_seg *segs) {
1934   int i, j;
1935   size_t bytes;
1936   i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count); /* checked 2jul06 - duplicating original */
1937   /*int have_alpha = im->channels == 2 || im->channels == 4;*/
1938
1939   memset(state, 0, sizeof(*state));
1940   /* we keep a local copy that we can adjust for speed */
1941   for (i = 0; i < count; ++i) {
1942     i_fountain_seg *seg = my_segs + i;
1943
1944     *seg = segs[i];
1945     if (seg->type < 0 || seg->type >= i_fst_end)
1946       seg->type = i_fst_linear;
1947     if (seg->color < 0 || seg->color >= i_fc_end)
1948       seg->color = i_fc_direct;
1949     if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1950       /* so we don't have to translate to HSV on each request, do it here */
1951       for (j = 0; j < 2; ++j) {
1952         i_rgb_to_hsvf(seg->c+j);
1953       }
1954       if (seg->color == i_fc_hue_up) {
1955         if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1956           seg->c[1].channel[0] += 1.0;
1957       }
1958       else {
1959         if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1960           seg->c[0].channel[0] += 1.0;
1961       }
1962     }
1963     /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1964            seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1965            seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1966            seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1967            seg->c[1].channel[3], seg->type, seg->color);*/
1968
1969   }
1970
1971   /* initialize each engine */
1972   /* these are so common ... */
1973   state->lA = xb - xa;
1974   state->lB = yb - ya;
1975   state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1976   state->xa = xa;
1977   state->ya = ya;
1978   switch (type) {
1979   default:
1980     type = i_ft_linear; /* make the invalid value valid */
1981   case i_ft_linear:
1982   case i_ft_bilinear:
1983     state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1984     state->mult = 1;
1985     state->mult = 1/linear_fount_f(xb, yb, state);
1986     break;
1987
1988   case i_ft_radial:
1989     state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1990                              + (double)(yb-ya)*(yb-ya));
1991     break;
1992
1993   case i_ft_radial_square:
1994     state->cos = state->lA / state->AB;
1995     state->sin = state->lB / state->AB;
1996     state->mult = 1.0 / state->AB;
1997     break;
1998
1999   case i_ft_revolution:
2000     state->theta = atan2(yb-ya, xb-xa);
2001     state->mult = 1.0 / (PI * 2);
2002     break;
2003
2004   case i_ft_conical:
2005     state->theta = atan2(yb-ya, xb-xa);
2006     state->mult = 1.0 / PI;
2007     break;
2008   }
2009   state->ffunc = fount_funcs[type];
2010   if (super_sample < 0
2011       || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
2012     super_sample = 0;
2013   }
2014   state->ssample_data = NULL;
2015   switch (super_sample) {
2016   case i_fts_grid:
2017     ssample_param = floor(0.5 + sqrt(ssample_param));
2018     bytes = ssample_param * ssample_param * sizeof(i_fcolor);
2019     if (bytes / sizeof(i_fcolor) == ssample_param * ssample_param) {
2020       state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param); /* checked 1jul06 tonyc */
2021     }
2022     else {
2023       super_sample = i_fts_none;
2024     }
2025     break;
2026
2027   case i_fts_random:
2028   case i_fts_circle:
2029     ssample_param = floor(0.5+ssample_param);
2030     if (im_size_t_max / sizeof(i_fcolor) > ssample_param) {
2031       bytes = sizeof(i_fcolor) * ssample_param;
2032       state->ssample_data = mymalloc(bytes);
2033     }
2034     else {
2035       dIMCTX;
2036       im_log((aIMCTX, 1,"size_t overflow calculating super-sample array size for random or circl"));
2037       super_sample = i_fts_none;
2038     }
2039     break;
2040   }
2041   state->parm = ssample_param;
2042   state->ssfunc = fount_ssamples[super_sample];
2043   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
2044     repeat = 0;
2045   state->rpfunc = fount_repeats[repeat];
2046   state->segs = my_segs;
2047   state->count = count;
2048 }
2049
2050 static void
2051 fount_finish_state(struct fount_state *state) {
2052   if (state->ssample_data)
2053     myfree(state->ssample_data);
2054   myfree(state->segs);
2055 }
2056
2057
2058 /*
2059 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
2060
2061 Evaluates the fountain fill at the given point.
2062
2063 This is called by both the non-super-sampling and super-sampling code.
2064
2065 You might think that it would make sense to sample the fill parameter
2066 instead, and combine those, but this breaks badly.
2067
2068 =cut
2069 */
2070
2071 static int
2072 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
2073   double v = (state->rpfunc)((state->ffunc)(x, y, state));
2074   int i;
2075
2076   i = 0;
2077   while (i < state->count
2078          && (v < state->segs[i].start || v > state->segs[i].end)) {
2079     ++i;
2080   }
2081   if (i < state->count) {
2082     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
2083     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
2084     return 1;
2085   }
2086   else
2087     return 0;
2088 }
2089
2090 /*
2091 =item linear_fount_f(x, y, state)
2092
2093 Calculate the fill parameter for a linear fountain fill.
2094
2095 Uses the point to line distance function, with some precalculation
2096 done in i_fountain().
2097
2098 =cut
2099 */
2100 static double
2101 linear_fount_f(double x, double y, struct fount_state *state) {
2102   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
2103 }
2104
2105 /*
2106 =item bilinear_fount_f(x, y, state)
2107
2108 Calculate the fill parameter for a bi-linear fountain fill.
2109
2110 =cut
2111 */
2112 static double
2113 bilinear_fount_f(double x, double y, struct fount_state *state) {
2114   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
2115 }
2116
2117 /*
2118 =item radial_fount_f(x, y, state)
2119
2120 Calculate the fill parameter for a radial fountain fill.
2121
2122 Simply uses the distance function.
2123
2124 =cut
2125  */
2126 static double
2127 radial_fount_f(double x, double y, struct fount_state *state) {
2128   return sqrt((double)(state->xa-x)*(state->xa-x)
2129               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
2130 }
2131
2132 /*
2133 =item square_fount_f(x, y, state)
2134
2135 Calculate the fill parameter for a square fountain fill.
2136
2137 Works by rotating the reference co-ordinate around the centre of the
2138 square.
2139
2140 =cut
2141 */
2142 static double
2143 square_fount_f(double x, double y, struct fount_state *state) {
2144   i_img_dim xc, yc; /* centred on A */
2145   double xt, yt; /* rotated by theta */
2146   xc = x - state->xa;
2147   yc = y - state->ya;
2148   xt = fabs(xc * state->cos + yc * state->sin);
2149   yt = fabs(-xc * state->sin + yc * state->cos);
2150   return (xt > yt ? xt : yt) * state->mult;
2151 }
2152
2153 /*
2154 =item revolution_fount_f(x, y, state)
2155
2156 Calculates the fill parameter for the revolution fountain fill.
2157
2158 =cut
2159 */
2160 static double
2161 revolution_fount_f(double x, double y, struct fount_state *state) {
2162   double angle = atan2(y - state->ya, x - state->xa);
2163
2164   angle -= state->theta;
2165   if (angle < 0) {
2166     angle = fmod(angle+ PI * 4, PI*2);
2167   }
2168
2169   return angle * state->mult;
2170 }
2171
2172 /*
2173 =item conical_fount_f(x, y, state)
2174
2175 Calculates the fill parameter for the conical fountain fill.
2176
2177 =cut
2178 */
2179 static double
2180 conical_fount_f(double x, double y, struct fount_state *state) {
2181   double angle = atan2(y - state->ya, x - state->xa);
2182
2183   angle -= state->theta;
2184   if (angle < -PI)
2185     angle += PI * 2;
2186   else if (angle > PI)
2187     angle -= PI * 2;
2188
2189   return fabs(angle) * state->mult;
2190 }
2191
2192 /*
2193 =item linear_interp(pos, seg)
2194
2195 Calculates linear interpolation on the fill parameter.  Breaks the
2196 segment into 2 regions based in the I<middle> value.
2197
2198 =cut
2199 */
2200 static double
2201 linear_interp(double pos, i_fountain_seg *seg) {
2202   if (pos < seg->middle) {
2203     double len = seg->middle - seg->start;
2204     if (len < EPSILON)
2205       return 0.0;
2206     else
2207       return (pos - seg->start) / len / 2;
2208   }
2209   else {
2210     double len = seg->end - seg->middle;
2211     if (len < EPSILON)
2212       return 1.0;
2213     else
2214       return 0.5 + (pos - seg->middle) / len / 2;
2215   }
2216 }
2217
2218 /*
2219 =item sine_interp(pos, seg)
2220
2221 Calculates sine function interpolation on the fill parameter.
2222
2223 =cut
2224 */
2225 static double
2226 sine_interp(double pos, i_fountain_seg *seg) {
2227   /* I wonder if there's a simple way to smooth the transition for this */
2228   double work = linear_interp(pos, seg);
2229
2230   return (1-cos(work * PI))/2;
2231 }
2232
2233 /*
2234 =item sphereup_interp(pos, seg)
2235
2236 Calculates spherical interpolation on the fill parameter, with the cusp
2237 at the low-end.
2238
2239 =cut
2240 */
2241 static double
2242 sphereup_interp(double pos, i_fountain_seg *seg) {
2243   double work = linear_interp(pos, seg);
2244
2245   return sqrt(1.0 - (1-work) * (1-work));
2246 }
2247
2248 /*
2249 =item spheredown_interp(pos, seg)
2250
2251 Calculates spherical interpolation on the fill parameter, with the cusp
2252 at the high-end.
2253
2254 =cut
2255 */
2256 static double
2257 spheredown_interp(double pos, i_fountain_seg *seg) {
2258   double work = linear_interp(pos, seg);
2259
2260   return 1-sqrt(1.0 - work * work);
2261 }
2262
2263 /*
2264 =item direct_cinterp(out, pos, seg)
2265
2266 Calculates the fountain color based on direct scaling of the channels
2267 of the color channels.
2268
2269 =cut
2270 */
2271 static void
2272 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2273   int ch;
2274   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2275     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2276       + seg->c[1].channel[ch] * pos;
2277   }
2278 }
2279
2280 /*
2281 =item hue_up_cinterp(put, pos, seg)
2282
2283 Calculates the fountain color based on scaling a HSV value.  The hue
2284 increases as the fill parameter increases.
2285
2286 =cut
2287 */
2288 static void
2289 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2290   int ch;
2291   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2292     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2293       + seg->c[1].channel[ch] * pos;
2294   }
2295   i_hsv_to_rgbf(out);
2296 }
2297
2298 /*
2299 =item hue_down_cinterp(put, pos, seg)
2300
2301 Calculates the fountain color based on scaling a HSV value.  The hue
2302 decreases as the fill parameter increases.
2303
2304 =cut
2305 */
2306 static void
2307 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2308   int ch;
2309   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2310     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2311       + seg->c[1].channel[ch] * pos;
2312   }
2313   i_hsv_to_rgbf(out);
2314 }
2315
2316 /*
2317 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2318
2319 Simple grid-based super-sampling.
2320
2321 =cut
2322 */
2323 static int
2324 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
2325   i_fcolor *work = state->ssample_data;
2326   i_img_dim dx, dy;
2327   int grid = state->parm;
2328   double base = -0.5 + 0.5 / grid;
2329   double step = 1.0 / grid;
2330   int ch, i;
2331   int samp_count = 0;
2332
2333   for (dx = 0; dx < grid; ++dx) {
2334     for (dy = 0; dy < grid; ++dy) {
2335       if (fount_getat(work+samp_count, x + base + step * dx,
2336                       y + base + step * dy, state)) {
2337         ++samp_count;
2338       }
2339     }
2340   }
2341   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2342     out->channel[ch] = 0;
2343     for (i = 0; i < samp_count; ++i) {
2344       out->channel[ch] += work[i].channel[ch];
2345     }
2346     /* we divide by 4 rather than samp_count since if there's only one valid
2347        sample it should be mostly transparent */
2348     out->channel[ch] /= grid * grid;
2349   }
2350   return samp_count;
2351 }
2352
2353 /*
2354 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2355
2356 Random super-sampling.
2357
2358 =cut
2359 */
2360 static int
2361 random_ssample(i_fcolor *out, double x, double y,
2362                struct fount_state *state) {
2363   i_fcolor *work = state->ssample_data;
2364   int i, ch;
2365   int maxsamples = state->parm;
2366   double rand_scale = 1.0 / RAND_MAX;
2367   int samp_count = 0;
2368   for (i = 0; i < maxsamples; ++i) {
2369     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
2370                     y - 0.5 + rand() * rand_scale, state)) {
2371       ++samp_count;
2372     }
2373   }
2374   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2375     out->channel[ch] = 0;
2376     for (i = 0; i < samp_count; ++i) {
2377       out->channel[ch] += work[i].channel[ch];
2378     }
2379     /* we divide by maxsamples rather than samp_count since if there's
2380        only one valid sample it should be mostly transparent */
2381     out->channel[ch] /= maxsamples;
2382   }
2383   return samp_count;
2384 }
2385
2386 /*
2387 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2388
2389 Super-sampling around the circumference of a circle.
2390
2391 I considered saving the sin()/cos() values and transforming step-size
2392 around the circle, but that's inaccurate, though it may not matter
2393 much.
2394
2395 =cut
2396  */
2397 static int
2398 circle_ssample(i_fcolor *out, double x, double y,
2399                struct fount_state *state) {
2400   i_fcolor *work = state->ssample_data;
2401   int i, ch;
2402   int maxsamples = state->parm;
2403   double angle = 2 * PI / maxsamples;
2404   double radius = 0.3; /* semi-random */
2405   int samp_count = 0;
2406   for (i = 0; i < maxsamples; ++i) {
2407     if (fount_getat(work+samp_count, x + radius * cos(angle * i),
2408                     y + radius * sin(angle * i), state)) {
2409       ++samp_count;
2410     }
2411   }
2412   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2413     out->channel[ch] = 0;
2414     for (i = 0; i < samp_count; ++i) {
2415       out->channel[ch] += work[i].channel[ch];
2416     }
2417     /* we divide by maxsamples rather than samp_count since if there's
2418        only one valid sample it should be mostly transparent */
2419     out->channel[ch] /= maxsamples;
2420   }
2421   return samp_count;
2422 }
2423
2424 /*
2425 =item fount_r_none(v)
2426
2427 Implements no repeats.  Simply clamps the fill value.
2428
2429 =cut
2430 */
2431 static double
2432 fount_r_none(double v) {
2433   return v < 0 ? 0 : v > 1 ? 1 : v;
2434 }
2435
2436 /*
2437 =item fount_r_sawtooth(v)
2438
2439 Implements sawtooth repeats.  Clamps negative values and uses fmod()
2440 on others.
2441
2442 =cut
2443 */
2444 static double
2445 fount_r_sawtooth(double v) {
2446   return v < 0 ? 0 : fmod(v, 1.0);
2447 }
2448
2449 /*
2450 =item fount_r_triangle(v)
2451
2452 Implements triangle repeats.  Clamps negative values, uses fmod to get
2453 a range 0 through 2 and then adjusts values > 1.
2454
2455 =cut
2456 */
2457 static double
2458 fount_r_triangle(double v) {
2459   if (v < 0)
2460     return 0;
2461   else {
2462     v = fmod(v, 2.0);
2463     return v > 1.0 ? 2.0 - v : v;
2464   }
2465 }
2466
2467 /*
2468 =item fount_r_saw_both(v)
2469
2470 Implements sawtooth repeats in the both postive and negative directions.
2471
2472 Adjusts the value to be postive and then just uses fmod().
2473
2474 =cut
2475 */
2476 static double
2477 fount_r_saw_both(double v) {
2478   if (v < 0)
2479     v += 1+(int)(-v);
2480   return fmod(v, 1.0);
2481 }
2482
2483 /*
2484 =item fount_r_tri_both(v)
2485
2486 Implements triangle repeats in the both postive and negative directions.
2487
2488 Uses fmod on the absolute value, and then adjusts values > 1.
2489
2490 =cut
2491 */
2492 static double
2493 fount_r_tri_both(double v) {
2494   v = fmod(fabs(v), 2.0);
2495   return v > 1.0 ? 2.0 - v : v;
2496 }
2497
2498 /*
2499 =item fill_fountf(fill, x, y, width, channels, data)
2500
2501 The fill function for fountain fills.
2502
2503 =cut
2504 */
2505 static void
2506 fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width,
2507             int channels, i_fcolor *data) {
2508   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2509
2510   while (width--) {
2511     i_fcolor c;
2512     int got_one;
2513
2514     if (f->state.ssfunc)
2515       got_one = f->state.ssfunc(&c, x, y, &f->state);
2516     else
2517       got_one = fount_getat(&c, x, y, &f->state);
2518
2519     if (got_one)
2520       *data++ = c;
2521
2522     ++x;
2523   }
2524 }
2525
2526 /*
2527 =item fount_fill_destroy(fill)
2528
2529 =cut
2530 */
2531 static void
2532 fount_fill_destroy(i_fill_t *fill) {
2533   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2534   fount_finish_state(&f->state);
2535 }
2536
2537 /*
2538 =back
2539
2540 =head1 AUTHOR
2541
2542 Arnar M. Hrafnkelsson <addi@umich.edu>
2543
2544 Tony Cook <tony@develop-help.com> (i_fountain())
2545
2546 =head1 SEE ALSO
2547
2548 Imager(3)
2549
2550 =cut
2551 */