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