filters.c

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