added OO interfaces for some filters
[imager.git] / 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 void
778 i_nearest_color_foo(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
779
780   int p, x, y;
781   int xsize    = im->xsize;
782   int ysize    = im->ysize;
783
784   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
785   
786   for(p = 0; p<num; p++) {
787     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
788     ICL_info(&ival[p]);
789   }
790
791   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
792     int   midx    = 0;
793     float mindist = 0;
794     float curdist = 0;
795
796     int xd        = x-xo[0];
797     int yd        = y-yo[0];
798
799     switch (dmeasure) {
800     case 0: /* euclidean */
801       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
802       break;
803     case 1: /* euclidean squared */
804       mindist = xd*xd + yd*yd; /* euclidean distance */
805       break;
806     case 2: /* euclidean squared */
807       mindist = max(xd*xd, yd*yd); /* manhattan distance */
808       break;
809     default:
810       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
811     }
812
813     for(p = 1; p<num; p++) {
814       xd = x-xo[p];
815       yd = y-yo[p];
816       switch (dmeasure) {
817       case 0: /* euclidean */
818         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
819         break;
820       case 1: /* euclidean squared */
821         curdist = xd*xd + yd*yd; /* euclidean distance */
822         break;
823       case 2: /* euclidean squared */
824         curdist = max(xd*xd, yd*yd); /* manhattan distance */
825         break;
826       default:
827         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
828       }
829       if (curdist < mindist) {
830         mindist = curdist;
831         midx = p;
832       }
833     }
834     i_ppix(im, x, y, &ival[midx]); 
835   }
836 }
837
838
839
840
841
842
843
844
845
846
847
848 void
849 i_nearest_color(i_img *im, int num, int *xo, int *yo, i_color *oval, int dmeasure) {
850   i_color *ival;
851   float *tval;
852   float c1, c2;
853   i_color val;
854   int p, x, y, ch;
855   int xsize    = im->xsize;
856   int ysize    = im->ysize;
857   int *cmatch;
858
859   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));
860
861   tval   = mymalloc( sizeof(float)*num*im->channels );
862   ival   = mymalloc( sizeof(i_color)*num );
863   cmatch = mymalloc( sizeof(int)*num     );
864
865   for(p = 0; p<num; p++) {
866     for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
867     cmatch[p] = 0;
868   }
869
870   
871   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
872     int   midx    = 0;
873     float mindist = 0;
874     float curdist = 0;
875     
876     int xd        = x-xo[0];
877     int yd        = y-yo[0];
878
879     switch (dmeasure) {
880     case 0: /* euclidean */
881       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
882       break;
883     case 1: /* euclidean squared */
884       mindist = xd*xd + yd*yd; /* euclidean distance */
885       break;
886     case 2: /* euclidean squared */
887       mindist = max(xd*xd, yd*yd); /* manhattan distance */
888       break;
889     default:
890       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
891     }
892     
893     for(p = 1; p<num; p++) {
894       xd = x-xo[p];
895       yd = y-yo[p];
896       switch (dmeasure) {
897       case 0: /* euclidean */
898         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
899         break;
900       case 1: /* euclidean squared */
901         curdist = xd*xd + yd*yd; /* euclidean distance */
902         break;
903       case 2: /* euclidean squared */
904         curdist = max(xd*xd, yd*yd); /* manhattan distance */
905         break;
906       default:
907         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
908       }
909       if (curdist < mindist) {
910         mindist = curdist;
911         midx = p;
912       }
913     }
914
915     cmatch[midx]++;
916     i_gpix(im, x, y, &val);
917     c2 = 1.0/(float)(cmatch[midx]);
918     c1 = 1.0-c2;
919     
920     for(ch = 0; ch<im->channels; ch++) 
921       tval[midx*im->channels + ch] = c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
922   
923     
924   }
925
926   for(p = 0; p<num; p++) for(ch = 0; ch<im->channels; ch++) ival[p].channel[ch] = tval[p*im->channels + ch];
927
928   i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
929 }