filters.c

   1 #include "image.h"
   2 #include <stdlib.h>
   3 #include <math.h>
   4
   5
   6 /*
   7 =head1 NAME
   8
   9 filters.c - implements filters that operate on images
  10
  11 =head1 SYNOPSIS
  12
  13
  14   i_contrast(im, 0.8);
  15   i_hardinvert(im);
  16   // and more
  17
  18 =head1 DESCRIPTION
  19
  20 filters.c implements basic filters for Imager.  These filters
  21 should be accessible from the filter interface as defined in
  22 the pod for Imager.
  23
  24 =head1 FUNCTION REFERENCE
  25
  26 Some of these functions are internal.
  27
  28 =over 4
  29
  30 =cut
  31 */
  32
  33
  34
  35
  36
  37
  38
  39 /*
  40 =item i_contrast(im, intensity)
  41
  42 Scales the pixel values by the amount specified.
  43
  44   im        - image object
  45   intensity - scalefactor
  46
  47 =cut
  48 */
  49
  50 void
  51 i_contrast(i_img *im, float intensity) {
  52   int x, y;
  53   unsigned char ch;
  54   unsigned int new_color;
  55   i_color rcolor;
  56
  57   mm_log((1,"i_contrast(im %p, intensity %f)\n", im, intensity));
  58
  59   if(intensity < 0) return;
  60
  61   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
  62     i_gpix(im, x, y, &rcolor);
  63
  64     for(ch = 0; ch < im->channels; ch++) {
  65       new_color = (unsigned int) rcolor.channel[ch];
  66       new_color *= intensity;
  67
  68       if(new_color > 255) {
  69         new_color = 255;
  70       }
  71       rcolor.channel[ch] = (unsigned char) new_color;
  72     }
  73     i_ppix(im, x, y, &rcolor);
  74   }
  75 }
  76
  77
  78 /*
  79 =item i_hardinvert(im)
  80
  81 Inverts the pixel values of the input image.
  82
  83   im        - image object
  84
  85 =cut
  86 */
  87
  88 void
  89 i_hardinvert(i_img *im) {
  90   int x, y;
  91   unsigned char ch;
  92
  93   i_color rcolor;
  94
  95     mm_log((1,"i_hardinvert(im %p)\n", im));
  96
  97   for(y = 0; y < im->ysize; y++) {
  98     for(x = 0; x < im->xsize; x++) {
  99       i_gpix(im, x, y, &rcolor);
 100
 101       for(ch = 0; ch < im->channels; ch++) {
 102         rcolor.channel[ch] = 255 - rcolor.channel[ch];
 103       }
 104
 105       i_ppix(im, x, y, &rcolor);
 106     }
 107   }
 108 }
 109
 110
 111
 112 /*
 113 =item i_noise(im, amount, type)
 114
 115 Inverts the pixel values by the amount specified.
 116
 117   im     - image object
 118   amount - deviation in pixel values
 119   type   - noise individual for each channel if true
 120
 121 =cut
 122 */
 123
 124 #ifdef _MSC_VER
 125 /* random() is non-ASCII, even if it is better than rand() */
 126 #define random() rand()
 127 #endif
 128
 129 void
 130 i_noise(i_img *im, float amount, unsigned char type) {
 131   int x, y;
 132   unsigned char ch;
 133   int new_color;
 134   float damount = amount * 2;
 135   i_color rcolor;
 136   int color_inc = 0;
 137
 138   mm_log((1,"i_noise(im %p, intensity %.2f\n", im, amount));
 139
 140   if(amount < 0) return;
 141
 142   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 143     i_gpix(im, x, y, &rcolor);
 144
 145     if(type == 0) {
 146       color_inc = (amount - (damount * ((float)random() / RAND_MAX)));
 147     }
 148
 149     for(ch = 0; ch < im->channels; ch++) {
 150       new_color = (int) rcolor.channel[ch];
 151
 152       if(type != 0) {
 153         new_color += (amount - (damount * ((float)random() / RAND_MAX)));
 154       } else {
 155         new_color += color_inc;
 156       }
 157
 158       if(new_color < 0) {
 159         new_color = 0;
 160       }
 161       if(new_color > 255) {
 162         new_color = 255;
 163       }
 164
 165       rcolor.channel[ch] = (unsigned char) new_color;
 166     }
 167
 168     i_ppix(im, x, y, &rcolor);
 169   }
 170 }
 171
 172
 173 /*
 174 =item i_noise(im, amount, type)
 175
 176 Inverts the pixel values by the amount specified.
 177
 178   im     - image object
 179   amount - deviation in pixel values
 180   type   - noise individual for each channel if true
 181
 182 =cut
 183 */
 184
 185
 186 /*
 187 =item i_applyimage(im, add_im, mode)
 188
 189 Apply's an image to another image
 190
 191   im     - target image
 192   add_im - image that is applied to target
 193   mode   - what method is used in applying:
 194
 195   0   Normal
 196   1   Multiply
 197   2   Screen
 198   3   Overlay
 199   4   Soft Light
 200   5   Hard Light
 201   6   Color dodge
 202   7   Color Burn
 203   8   Darker
 204   9   Lighter
 205   10  Add
 206   11  Subtract
 207   12  Difference
 208   13  Exclusion
 209
 210 =cut
 211 */
 212
 213 void i_applyimage(i_img *im, i_img *add_im, unsigned char mode) {
 214   int x, y;
 215   int mx, my;
 216
 217   mm_log((1, "i_applyimage(im %p, add_im %p, mode %d", im, add_im, mode));
 218
 219   mx = (add_im->xsize <= im->xsize) ? add_im->xsize : add_im->xsize;
 220   my = (add_im->ysize <= im->ysize) ? add_im->ysize : add_im->ysize;
 221
 222   for(x = 0; x < mx; x++) {
 223     for(y = 0; y < my; y++) {
 224     }
 225   }
 226 }
 227
 228
 229 /*
 230 =item i_bumpmap(im, bump, channel, light_x, light_y, st)
 231
 232 Makes a bumpmap on image im using the bump image as the elevation map.
 233
 234   im      - target image
 235   bump    - image that contains the elevation info
 236   channel - to take the elevation information from
 237   light_x - x coordinate of light source
 238   light_y - y coordinate of light source
 239   st      - length of shadow
 240
 241 =cut
 242 */
 243
 244 void
 245 i_bumpmap(i_img *im, i_img *bump, int channel, int light_x, int light_y, int st) {
 246   int x, y, ch;
 247   int mx, my;
 248   i_color x1_color, y1_color, x2_color, y2_color, dst_color;
 249   double nX, nY;
 250   double tX, tY, tZ;
 251   double aX, aY, aL;
 252   double fZ;
 253   unsigned char px1, px2, py1, py2;
 254
 255   i_img new_im;
 256
 257   mm_log((1, "i_bumpmap(im %p, add_im %p, channel %d, light_x %d, light_y %d, st %d)\n",
 258           im, bump, channel, light_x, light_y, st));
 259
 260
 261   if(channel >= bump->channels) {
 262     mm_log((1, "i_bumpmap: channel = %d while bump image only has %d channels\n", channel, bump->channels));
 263     return;
 264   }
 265
 266   mx = (bump->xsize <= im->xsize) ? bump->xsize : im->xsize;
 267   my = (bump->ysize <= im->ysize) ? bump->ysize : im->ysize;
 268
 269   i_img_empty_ch(&new_im, im->xsize, im->ysize, im->channels);
 270
 271   aX = (light_x > (mx >> 1)) ? light_x : mx - light_x;
 272   aY = (light_y > (my >> 1)) ? light_y : my - light_y;
 273
 274   aL = sqrt((aX * aX) + (aY * aY));
 275
 276   for(y = 1; y < my - 1; y++) {
 277     for(x = 1; x < mx - 1; x++) {
 278       i_gpix(bump, x + st, y, &x1_color);
 279       i_gpix(bump, x, y + st, &y1_color);
 280       i_gpix(bump, x - st, y, &x2_color);
 281       i_gpix(bump, x, y - st, &y2_color);
 282
 283       i_gpix(im, x, y, &dst_color);
 284
 285       px1 = x1_color.channel[channel];
 286       py1 = y1_color.channel[channel];
 287       px2 = x2_color.channel[channel];
 288       py2 = y2_color.channel[channel];
 289
 290       nX = px1 - px2;
 291       nY = py1 - py2;
 292
 293       nX += 128;
 294       nY += 128;
 295
 296       fZ = (sqrt((nX * nX) + (nY * nY)) / aL);
 297
 298       tX = abs(x - light_x) / aL;
 299       tY = abs(y - light_y) / aL;
 300
 301       tZ = 1 - (sqrt((tX * tX) + (tY * tY)) * fZ);
 302
 303       if(tZ < 0) tZ = 0;
 304       if(tZ > 2) tZ = 2;
 305
 306       for(ch = 0; ch < im->channels; ch++)
 307         dst_color.channel[ch] = (unsigned char) (float)(dst_color.channel[ch] * tZ);
 308
 309       i_ppix(&new_im, x, y, &dst_color);
 310     }
 311   }
 312
 313   i_copyto(im, &new_im, 0, 0, (int)im->xsize, (int)im->ysize, 0, 0);
 314
 315   i_img_exorcise(&new_im);
 316 }
 317
 318
 319
 320 /*
 321 =item i_postlevels(im, levels)
 322
 323 Quantizes Images to fewer levels.
 324
 325   im      - target image
 326   levels  - number of levels
 327
 328 =cut
 329 */
 330
 331 void
 332 i_postlevels(i_img *im, int levels) {
 333   int x, y, ch;
 334   float pv;
 335   int rv;
 336   float av;
 337
 338   i_color rcolor;
 339
 340   rv = (int) ((float)(256 / levels));
 341   av = (float)levels;
 342
 343   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 344     i_gpix(im, x, y, &rcolor);
 345
 346     for(ch = 0; ch < im->channels; ch++) {
 347       pv = (((float)rcolor.channel[ch] / 255)) * av;
 348       pv = (int) ((int)pv * rv);
 349
 350       if(pv < 0) pv = 0;
 351       else if(pv > 255) pv = 255;
 352
 353       rcolor.channel[ch] = (unsigned char) pv;
 354     }
 355     i_ppix(im, x, y, &rcolor);
 356   }
 357 }
 358
 359
 360 /*
 361 =item i_mosaic(im, size)
 362
 363 Makes an image looks like a mosaic with tilesize of size
 364
 365   im      - target image
 366   size    - size of tiles
 367
 368 =cut
 369 */
 370
 371 void
 372 i_mosaic(i_img *im, int size) {
 373   int x, y, ch;
 374   int lx, ly, z;
 375   long sqrsize;
 376
 377   i_color rcolor;
 378   long col[256];
 379
 380   sqrsize = size * size;
 381
 382   for(y = 0; y < im->ysize; y += size) for(x = 0; x < im->xsize; x += size) {
 383     for(z = 0; z < 256; z++) col[z] = 0;
 384
 385     for(lx = 0; lx < size; lx++) {
 386       for(ly = 0; ly < size; ly++) {
 387         i_gpix(im, (x + lx), (y + ly), &rcolor);
 388
 389         for(ch = 0; ch < im->channels; ch++) {
 390           col[ch] += rcolor.channel[ch];
 391         }
 392       }
 393     }
 394
 395     for(ch = 0; ch < im->channels; ch++)
 396       rcolor.channel[ch] = (int) ((float)col[ch] / sqrsize);
 397
 398
 399     for(lx = 0; lx < size; lx++)
 400       for(ly = 0; ly < size; ly++)
 401       i_ppix(im, (x + lx), (y + ly), &rcolor);
 402
 403   }
 404 }
 405
 406 /*
 407 =item saturate(in)
 408
 409 Clamps the input value between 0 and 255. (internal)
 410
 411   in - input integer
 412
 413 =cut
 414 */
 415
 416 static
 417 unsigned char
 418 saturate(int in) {
 419   if (in>255) { return 255; }
 420   else if (in>0) return in;
 421   return 0;
 422 }
 423
 424
 425 /*
 426 =item i_watermark(im, wmark, tx, ty, pixdiff)
 427
 428 Applies a watermark to the target image
 429
 430   im      - target image
 431   wmark   - watermark image
 432   tx      - x coordinate of where watermark should be applied
 433   ty      - y coordinate of where watermark should be applied
 434   pixdiff - the magnitude of the watermark, controls how visible it is
 435
 436 =cut
 437 */
 438
 439 void
 440 i_watermark(i_img *im, i_img *wmark, int tx, int ty, int pixdiff) {
 441   int vx, vy, ch;
 442   i_color val, wval;
 443
 444   for(vx=0;vx<128;vx++) for(vy=0;vy<110;vy++) {
 445
 446     i_gpix(im,    tx+vx, ty+vy,&val );
 447     i_gpix(wmark, vx,    vy,   &wval);
 448
 449     for(ch=0;ch<im->channels;ch++)
 450       val.channel[ch] = saturate( val.channel[ch] + (pixdiff* (wval.channel[0]-128) )/128 );
 451
 452     i_ppix(im,tx+vx,ty+vy,&val);
 453   }
 454 }
 455
 456
 457 /*
 458 =item i_autolevels(im, lsat, usat, skew)
 459
 460 Scales and translates each color such that it fills the range completely.
 461 Skew is not implemented yet - purpose is to control the color skew that can
 462 occur when changing the contrast.
 463
 464   im   - target image
 465   lsat - fraction of pixels that will be truncated at the lower end of the spectrum
 466   usat - fraction of pixels that will be truncated at the higher end of the spectrum
 467   skew - not used yet
 468
 469 =cut
 470 */
 471
 472 void
 473 i_autolevels(i_img *im, float lsat, float usat, float skew) {
 474   i_color val;
 475   int i, x, y, rhist[256], ghist[256], bhist[256];
 476   int rsum, rmin, rmax;
 477   int gsum, gmin, gmax;
 478   int bsum, bmin, bmax;
 479   int rcl, rcu, gcl, gcu, bcl, bcu;
 480
 481   mm_log((1,"i_autolevels(im %p, lsat %f,usat %f,skew %f)\n", im, lsat,usat,skew));
 482
 483   rsum=gsum=bsum=0;
 484   for(i=0;i<256;i++) rhist[i]=ghist[i]=bhist[i] = 0;
 485   /* create histogram for each channel */
 486   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 487     i_gpix(im, x, y, &val);
 488     rhist[val.channel[0]]++;
 489     ghist[val.channel[1]]++;
 490     bhist[val.channel[2]]++;
 491   }
 492
 493   for(i=0;i<256;i++) {
 494     rsum+=rhist[i];
 495     gsum+=ghist[i];
 496     bsum+=bhist[i];
 497   }
 498
 499   rmin = gmin = bmin = 0;
 500   rmax = gmax = bmax = 255;
 501
 502   rcu = rcl = gcu = gcl = bcu = bcl = 0;
 503
 504   for(i=0; i<256; i++) {
 505     rcl += rhist[i];     if ( (rcl<rsum*lsat) ) rmin=i;
 506     rcu += rhist[255-i]; if ( (rcu<rsum*usat) ) rmax=255-i;
 507
 508     gcl += ghist[i];     if ( (gcl<gsum*lsat) ) gmin=i;
 509     gcu += ghist[255-i]; if ( (gcu<gsum*usat) ) gmax=255-i;
 510
 511     bcl += bhist[i];     if ( (bcl<bsum*lsat) ) bmin=i;
 512     bcu += bhist[255-i]; if ( (bcu<bsum*usat) ) bmax=255-i;
 513   }
 514
 515   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 516     i_gpix(im, x, y, &val);
 517     val.channel[0]=saturate((val.channel[0]-rmin)*255/(rmax-rmin));
 518     val.channel[1]=saturate((val.channel[1]-gmin)*255/(gmax-gmin));
 519     val.channel[2]=saturate((val.channel[2]-bmin)*255/(bmax-bmin));
 520     i_ppix(im, x, y, &val);
 521   }
 522 }
 523
 524 /*
 525 =item Noise(x,y)
 526
 527 Pseudo noise utility function used to generate perlin noise. (internal)
 528
 529   x - x coordinate
 530   y - y coordinate
 531
 532 =cut
 533 */
 534
 535 static
 536 float
 537 Noise(int x, int y) {
 538   int n = x + y * 57;
 539   n = (n<<13) ^ n;
 540   return ( 1.0 - ( (n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
 541 }
 542
 543 /*
 544 =item SmoothedNoise1(x,y)
 545
 546 Pseudo noise utility function used to generate perlin noise. (internal)
 547
 548   x - x coordinate
 549   y - y coordinate
 550
 551 =cut
 552 */
 553
 554 static
 555 float
 556 SmoothedNoise1(float x, float y) {
 557   float corners = ( Noise(x-1, y-1)+Noise(x+1, y-1)+Noise(x-1, y+1)+Noise(x+1, y+1) ) / 16;
 558   float sides   = ( Noise(x-1, y)  +Noise(x+1, y)  +Noise(x, y-1)  +Noise(x, y+1) ) /  8;
 559   float center  =  Noise(x, y) / 4;
 560   return corners + sides + center;
 561 }
 562
 563
 564 /*
 565 =item G_Interpolate(a, b, x)
 566
 567 Utility function used to generate perlin noise. (internal)
 568
 569 =cut
 570 */
 571
 572 static
 573 float C_Interpolate(float a, float b, float x) {
 574   /*  float ft = x * 3.1415927; */
 575   float ft = x * PI;
 576   float f = (1 - cos(ft)) * .5;
 577   return  a*(1-f) + b*f;
 578 }
 579
 580
 581 /*
 582 =item InterpolatedNoise(x, y)
 583
 584 Utility function used to generate perlin noise. (internal)
 585
 586 =cut
 587 */
 588
 589 static
 590 float
 591 InterpolatedNoise(float x, float y) {
 592
 593   int integer_X = x;
 594   float fractional_X = x - integer_X;
 595   int integer_Y = y;
 596   float fractional_Y = y - integer_Y;
 597
 598   float v1 = SmoothedNoise1(integer_X,     integer_Y);
 599   float v2 = SmoothedNoise1(integer_X + 1, integer_Y);
 600   float v3 = SmoothedNoise1(integer_X,     integer_Y + 1);
 601   float v4 = SmoothedNoise1(integer_X + 1, integer_Y + 1);
 602
 603   float i1 = C_Interpolate(v1 , v2 , fractional_X);
 604   float i2 = C_Interpolate(v3 , v4 , fractional_X);
 605
 606   return C_Interpolate(i1 , i2 , fractional_Y);
 607 }
 608
 609
 610
 611 /*
 612 =item PerlinNoise_2D(x, y)
 613
 614 Utility function used to generate perlin noise. (internal)
 615
 616 =cut
 617 */
 618
 619 static
 620 float
 621 PerlinNoise_2D(float x, float y) {
 622   int i,frequency;
 623   float amplitude;
 624   float total = 0;
 625   int Number_Of_Octaves=6;
 626   int n = Number_Of_Octaves - 1;
 627
 628   for(i=0;i<n;i++) {
 629     frequency = 2*i;
 630     amplitude = PI;
 631     total = total + InterpolatedNoise(x * frequency, y * frequency) * amplitude;
 632   }
 633
 634   return total;
 635 }
 636
 637
 638 /*
 639 =item i_radnoise(im, xo, yo, rscale, ascale)
 640
 641 Perlin-like radial noise.
 642
 643   im     - target image
 644   xo     - x coordinate of center
 645   yo     - y coordinate of center
 646   rscale - radial scale
 647   ascale - angular scale
 648
 649 =cut
 650 */
 651
 652 void
 653 i_radnoise(i_img *im, int xo, int yo, float rscale, float ascale) {
 654   int x, y, ch;
 655   i_color val;
 656   unsigned char v;
 657   float xc, yc, r;
 658   double a;
 659
 660   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 661     xc = (float)x-xo+0.5;
 662     yc = (float)y-yo+0.5;
 663     r = rscale*sqrt(xc*xc+yc*yc)+1.2;
 664     a = (PI+atan2(yc,xc))*ascale;
 665     v = saturate(128+100*(PerlinNoise_2D(a,r)));
 666     /* v=saturate(120+12*PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale));  Good soft marble */
 667     for(ch=0; ch<im->channels; ch++) val.channel[ch]=v;
 668     i_ppix(im, x, y, &val);
 669   }
 670 }
 671
 672
 673 /*
 674 =item i_turbnoise(im, xo, yo, scale)
 675
 676 Perlin-like 2d noise noise.
 677
 678   im     - target image
 679   xo     - x coordinate translation
 680   yo     - y coordinate translation
 681   scale  - scale of noise
 682
 683 =cut
 684 */
 685
 686 void
 687 i_turbnoise(i_img *im, float xo, float yo, float scale) {
 688   int x,y,ch;
 689   unsigned char v;
 690   i_color val;
 691
 692   for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
 693     /*    v=saturate(125*(1.0+PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale))); */
 694     v = saturate(120*(1.0+sin(xo+(float)x/scale+PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale))));
 695     for(ch=0; ch<im->channels; ch++) val.channel[ch] = v;
 696     i_ppix(im, x, y, &val);
 697   }
 698 }
 699
 700
 701
 702 /*
 703 =item i_gradgen(im, num, xo, yo, ival, dmeasure)
 704
 705 Gradient generating function.
 706
 707   im     - target image
 708   num    - number of points given
 709   xo     - array of x coordinates
 710   yo     - array of y coordinates
 711   ival   - array of i_color objects
 712   dmeasure - distance measure to be used.
 713     0 = Euclidean
 714     1 = Euclidean squared
 715     2 = Manhattan distance
 716
 717 =cut
 718 */
 719
 720
 721 void
 722 i_gradgen(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
 723
 724   i_color val;
 725   int p, x, y, ch;
 726   int channels = im->channels;
 727   int xsize    = im->xsize;
 728   int ysize    = im->ysize;
 729
 730   float *fdist;
 731
 732   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 733
 734   for(p = 0; p<num; p++) {
 735     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
 736     ICL_info(&ival[p]);
 737   }
 738
 739   fdist = mymalloc( sizeof(float) * num );
 740
 741   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 742     float cs = 0;
 743     float csd = 0;
 744     for(p = 0; p<num; p++) {
 745       int xd    = x-xo[p];
 746       int yd    = y-yo[p];
 747       switch (dmeasure) {
 748       case 0: /* euclidean */
 749         fdist[p]  = sqrt(xd*xd + yd*yd); /* euclidean distance */
 750         break;
 751       case 1: /* euclidean squared */
 752         fdist[p]  = xd*xd + yd*yd; /* euclidean distance */
 753         break;
 754       case 2: /* euclidean squared */
 755         fdist[p]  = max(xd*xd, yd*yd); /* manhattan distance */
 756         break;
 757       default:
 758         m_fatal(3,"i_gradgen: Unknown distance measure\n");
 759       }
 760       cs += fdist[p];
 761     }
 762
 763     csd = 1/((num-1)*cs);
 764
 765     for(p = 0; p<num; p++) fdist[p] = (cs-fdist[p])*csd;
 766
 767     for(ch = 0; ch<channels; ch++) {
 768       int tres = 0;
 769       for(p = 0; p<num; p++) tres += ival[p].channel[ch] * fdist[p];
 770       val.channel[ch] = saturate(tres);
 771     }
 772     i_ppix(im, x, y, &val);
 773   }
 774
 775 }
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792 void
 793 i_nearest_color_foo(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
 794
 795   int p, x, y;
 796   int xsize    = im->xsize;
 797   int ysize    = im->ysize;
 798
 799   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
 800
 801   for(p = 0; p<num; p++) {
 802     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
 803     ICL_info(&ival[p]);
 804   }
 805
 806   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 807     int   midx    = 0;
 808     float mindist = 0;
 809     float curdist = 0;
 810
 811     int xd        = x-xo[0];
 812     int yd        = y-yo[0];
 813
 814     switch (dmeasure) {
 815     case 0: /* euclidean */
 816       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 817       break;
 818     case 1: /* euclidean squared */
 819       mindist = xd*xd + yd*yd; /* euclidean distance */
 820       break;
 821     case 2: /* euclidean squared */
 822       mindist = max(xd*xd, yd*yd); /* manhattan distance */
 823       break;
 824     default:
 825       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 826     }
 827
 828     for(p = 1; p<num; p++) {
 829       xd = x-xo[p];
 830       yd = y-yo[p];
 831       switch (dmeasure) {
 832       case 0: /* euclidean */
 833         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 834         break;
 835       case 1: /* euclidean squared */
 836         curdist = xd*xd + yd*yd; /* euclidean distance */
 837         break;
 838       case 2: /* euclidean squared */
 839         curdist = max(xd*xd, yd*yd); /* manhattan distance */
 840         break;
 841       default:
 842         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 843       }
 844       if (curdist < mindist) {
 845         mindist = curdist;
 846         midx = p;
 847       }
 848     }
 849     i_ppix(im, x, y, &ival[midx]);
 850   }
 851 }
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863 void
 864 i_nearest_color(i_img *im, int num, int *xo, int *yo, i_color *oval, int dmeasure) {
 865   i_color *ival;
 866   float *tval;
 867   float c1, c2;
 868   i_color val;
 869   int p, x, y, ch;
 870   int xsize    = im->xsize;
 871   int ysize    = im->ysize;
 872   int *cmatch;
 873
 874   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));
 875
 876   tval   = mymalloc( sizeof(float)*num*im->channels );
 877   ival   = mymalloc( sizeof(i_color)*num );
 878   cmatch = mymalloc( sizeof(int)*num     );
 879
 880   for(p = 0; p<num; p++) {
 881     for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
 882     cmatch[p] = 0;
 883   }
 884
 885
 886   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
 887     int   midx    = 0;
 888     float mindist = 0;
 889     float curdist = 0;
 890
 891     int xd        = x-xo[0];
 892     int yd        = y-yo[0];
 893
 894     switch (dmeasure) {
 895     case 0: /* euclidean */
 896       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 897       break;
 898     case 1: /* euclidean squared */
 899       mindist = xd*xd + yd*yd; /* euclidean distance */
 900       break;
 901     case 2: /* euclidean squared */
 902       mindist = max(xd*xd, yd*yd); /* manhattan distance */
 903       break;
 904     default:
 905       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 906     }
 907
 908     for(p = 1; p<num; p++) {
 909       xd = x-xo[p];
 910       yd = y-yo[p];
 911       switch (dmeasure) {
 912       case 0: /* euclidean */
 913         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
 914         break;
 915       case 1: /* euclidean squared */
 916         curdist = xd*xd + yd*yd; /* euclidean distance */
 917         break;
 918       case 2: /* euclidean squared */
 919         curdist = max(xd*xd, yd*yd); /* manhattan distance */
 920         break;
 921       default:
 922         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
 923       }
 924       if (curdist < mindist) {
 925         mindist = curdist;
 926         midx = p;
 927       }
 928     }
 929
 930     cmatch[midx]++;
 931     i_gpix(im, x, y, &val);
 932     c2 = 1.0/(float)(cmatch[midx]);
 933     c1 = 1.0-c2;
 934
 935     //    printf("pixel [%d %d %d] c1+c2 = %f\n", val.channel[0], val.channel[1], val.channel[2], c1+c2);
 936     //    printf("cmatch = %d, c1 = %f, c2 = %f tval=[%f %f %f]\n", cmatch[midx], c1, c2, tval[midx*im->channels], tval[midx*im->channels+1], tval[midx*im->channels+2] );
 937
 938     for(ch = 0; ch<im->channels; ch++)
 939       tval[midx*im->channels + ch] = c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
 940
 941
 942   }
 943
 944   for(p = 0; p<num; p++) for(ch = 0; ch<im->channels; ch++) ival[p].channel[ch] = tval[p*im->channels + ch];
 945
 946   i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
 947 }