filters.c

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