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