prevent an unsigned overflow in FT1 has_chars() implementation
[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
1775   fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine, 
1776                    super_sample, ssample_param, count, segs);
1777
1778   for (y = 0; y < im->ysize; ++y) {
1779     i_glinf(im, 0, im->xsize, y, line);
1780     for (x = 0; x < im->xsize; ++x) {
1781       i_fcolor c;
1782       int got_one;
1783       if (super_sample == i_fts_none)
1784         got_one = fount_getat(&c, x, y, &state);
1785       else
1786         got_one = state.ssfunc(&c, x, y, &state);
1787       if (got_one) {
1788         if (combinef_func)
1789           work[x] = c;
1790         else 
1791           line[x] = c;
1792       }
1793     }
1794     if (combinef_func)
1795       combinef_func(line, work, im->channels, im->xsize);
1796     i_plinf(im, 0, im->xsize, y, line);
1797   }
1798   fount_finish_state(&state);
1799   myfree(work);
1800   myfree(line);
1801
1802   return 1;
1803 }
1804
1805 typedef struct {
1806   i_fill_t base;
1807   struct fount_state state;
1808 } i_fill_fountain_t;
1809
1810 static void
1811 fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width, int channels, 
1812             i_fcolor *data);
1813 static void
1814 fount_fill_destroy(i_fill_t *fill);
1815
1816 static i_fill_fountain_t
1817 fount_fill_proto =
1818   {
1819     {
1820       NULL,
1821       fill_fountf,
1822       fount_fill_destroy
1823     }
1824   };
1825
1826
1827 /*
1828 =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>)
1829
1830 =category Fills
1831 =synopsis fill = i_new_fill_fount(0, 0, 100, 100, i_ft_linear, i_ft_linear, 
1832 =synopsis                         i_fr_triangle, 0, i_fts_grid, 9, 1, segs);
1833
1834
1835 Creates a new general fill which fills with a fountain fill.
1836
1837 =cut
1838 */
1839
1840 i_fill_t *
1841 i_new_fill_fount(double xa, double ya, double xb, double yb, 
1842                  i_fountain_type type, i_fountain_repeat repeat, 
1843                  int combine, int super_sample, double ssample_param, 
1844                  int count, i_fountain_seg *segs) {
1845   i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1846   
1847   *fill = fount_fill_proto;
1848   if (combine)
1849     i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1850   else {
1851     fill->base.combine = NULL;
1852     fill->base.combinef = NULL;
1853   }
1854   fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine, 
1855                    super_sample, ssample_param, count, segs);
1856
1857   return &fill->base;
1858 }
1859
1860 /*
1861 =back
1862
1863 =head1 INTERNAL FUNCTIONS
1864
1865 =over
1866
1867 =item fount_init_state(...)
1868
1869 Used by both the fountain fill filter and the fountain fill.
1870
1871 =cut
1872 */
1873
1874 static void
1875 fount_init_state(struct fount_state *state, double xa, double ya, 
1876                  double xb, double yb, i_fountain_type type, 
1877                  i_fountain_repeat repeat, int combine, int super_sample, 
1878                  double ssample_param, int count, i_fountain_seg *segs) {
1879   int i, j;
1880   size_t bytes;
1881   i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count); /* checked 2jul06 - duplicating original */
1882   /*int have_alpha = im->channels == 2 || im->channels == 4;*/
1883   
1884   memset(state, 0, sizeof(*state));
1885   /* we keep a local copy that we can adjust for speed */
1886   for (i = 0; i < count; ++i) {
1887     i_fountain_seg *seg = my_segs + i;
1888
1889     *seg = segs[i];
1890     if (seg->type < 0 || seg->type >= i_fst_end)
1891       seg->type = i_fst_linear;
1892     if (seg->color < 0 || seg->color >= i_fc_end)
1893       seg->color = i_fc_direct;
1894     if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1895       /* so we don't have to translate to HSV on each request, do it here */
1896       for (j = 0; j < 2; ++j) {
1897         i_rgb_to_hsvf(seg->c+j);
1898       }
1899       if (seg->color == i_fc_hue_up) {
1900         if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1901           seg->c[1].channel[0] += 1.0;
1902       }
1903       else {
1904         if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1905           seg->c[0].channel[0] += 1.0;
1906       }
1907     }
1908     /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n", 
1909            seg->start, seg->middle, seg->end, seg->c[0].channel[0], 
1910            seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1911            seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2], 
1912            seg->c[1].channel[3], seg->type, seg->color);*/
1913            
1914   }
1915
1916   /* initialize each engine */
1917   /* these are so common ... */
1918   state->lA = xb - xa;
1919   state->lB = yb - ya;
1920   state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1921   state->xa = xa;
1922   state->ya = ya;
1923   switch (type) {
1924   default:
1925     type = i_ft_linear; /* make the invalid value valid */
1926   case i_ft_linear:
1927   case i_ft_bilinear:
1928     state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1929     state->mult = 1;
1930     state->mult = 1/linear_fount_f(xb, yb, state);
1931     break;
1932
1933   case i_ft_radial:
1934     state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa) 
1935                              + (double)(yb-ya)*(yb-ya));
1936     break;
1937
1938   case i_ft_radial_square:
1939     state->cos = state->lA / state->AB;
1940     state->sin = state->lB / state->AB;
1941     state->mult = 1.0 / state->AB;
1942     break;
1943
1944   case i_ft_revolution:
1945     state->theta = atan2(yb-ya, xb-xa);
1946     state->mult = 1.0 / (PI * 2);
1947     break;
1948
1949   case i_ft_conical:
1950     state->theta = atan2(yb-ya, xb-xa);
1951     state->mult = 1.0 / PI;
1952     break;
1953   }
1954   state->ffunc = fount_funcs[type];
1955   if (super_sample < 0 
1956       || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
1957     super_sample = 0;
1958   }
1959   state->ssample_data = NULL;
1960   switch (super_sample) {
1961   case i_fts_grid:
1962     ssample_param = floor(0.5 + sqrt(ssample_param));
1963     bytes = ssample_param * ssample_param * sizeof(i_fcolor);
1964     if (bytes / sizeof(i_fcolor) == ssample_param * ssample_param) {
1965       state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param); /* checked 1jul06 tonyc */
1966     }
1967     else {
1968       super_sample = i_fts_none;
1969     }
1970     break;
1971
1972   case i_fts_random:
1973   case i_fts_circle:
1974     ssample_param = floor(0.5+ssample_param);
1975     bytes = sizeof(i_fcolor) * ssample_param;
1976     if (bytes / sizeof(i_fcolor) == ssample_param) {
1977       state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
1978     }
1979     else {
1980       super_sample = i_fts_none;
1981     }
1982     break;
1983   }
1984   state->parm = ssample_param;
1985   state->ssfunc = fount_ssamples[super_sample];
1986   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1987     repeat = 0;
1988   state->rpfunc = fount_repeats[repeat];
1989   state->segs = my_segs;
1990   state->count = count;
1991 }
1992
1993 static void
1994 fount_finish_state(struct fount_state *state) {
1995   if (state->ssample_data)
1996     myfree(state->ssample_data);
1997   myfree(state->segs);
1998 }
1999
2000
2001 /*
2002 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
2003
2004 Evaluates the fountain fill at the given point.
2005
2006 This is called by both the non-super-sampling and super-sampling code.
2007
2008 You might think that it would make sense to sample the fill parameter
2009 instead, and combine those, but this breaks badly.
2010
2011 =cut
2012 */
2013
2014 static int
2015 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
2016   double v = (state->rpfunc)((state->ffunc)(x, y, state));
2017   int i;
2018
2019   i = 0;
2020   while (i < state->count 
2021          && (v < state->segs[i].start || v > state->segs[i].end)) {
2022     ++i;
2023   }
2024   if (i < state->count) {
2025     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
2026     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
2027     return 1;
2028   }
2029   else
2030     return 0;
2031 }
2032
2033 /*
2034 =item linear_fount_f(x, y, state)
2035
2036 Calculate the fill parameter for a linear fountain fill.
2037
2038 Uses the point to line distance function, with some precalculation
2039 done in i_fountain().
2040
2041 =cut
2042 */
2043 static double
2044 linear_fount_f(double x, double y, struct fount_state *state) {
2045   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
2046 }
2047
2048 /*
2049 =item bilinear_fount_f(x, y, state)
2050
2051 Calculate the fill parameter for a bi-linear fountain fill.
2052
2053 =cut
2054 */
2055 static double
2056 bilinear_fount_f(double x, double y, struct fount_state *state) {
2057   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
2058 }
2059
2060 /*
2061 =item radial_fount_f(x, y, state)
2062
2063 Calculate the fill parameter for a radial fountain fill.
2064
2065 Simply uses the distance function.
2066
2067 =cut
2068  */
2069 static double
2070 radial_fount_f(double x, double y, struct fount_state *state) {
2071   return sqrt((double)(state->xa-x)*(state->xa-x) 
2072               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
2073 }
2074
2075 /*
2076 =item square_fount_f(x, y, state)
2077
2078 Calculate the fill parameter for a square fountain fill.
2079
2080 Works by rotating the reference co-ordinate around the centre of the
2081 square.
2082
2083 =cut
2084 */
2085 static double
2086 square_fount_f(double x, double y, struct fount_state *state) {
2087   i_img_dim xc, yc; /* centred on A */
2088   double xt, yt; /* rotated by theta */
2089   xc = x - state->xa;
2090   yc = y - state->ya;
2091   xt = fabs(xc * state->cos + yc * state->sin);
2092   yt = fabs(-xc * state->sin + yc * state->cos);
2093   return (xt > yt ? xt : yt) * state->mult;
2094 }
2095
2096 /*
2097 =item revolution_fount_f(x, y, state)
2098
2099 Calculates the fill parameter for the revolution fountain fill.
2100
2101 =cut
2102 */
2103 static double
2104 revolution_fount_f(double x, double y, struct fount_state *state) {
2105   double angle = atan2(y - state->ya, x - state->xa);
2106   
2107   angle -= state->theta;
2108   if (angle < 0) {
2109     angle = fmod(angle+ PI * 4, PI*2);
2110   }
2111
2112   return angle * state->mult;
2113 }
2114
2115 /*
2116 =item conical_fount_f(x, y, state)
2117
2118 Calculates the fill parameter for the conical fountain fill.
2119
2120 =cut
2121 */
2122 static double
2123 conical_fount_f(double x, double y, struct fount_state *state) {
2124   double angle = atan2(y - state->ya, x - state->xa);
2125   
2126   angle -= state->theta;
2127   if (angle < -PI)
2128     angle += PI * 2;
2129   else if (angle > PI) 
2130     angle -= PI * 2;
2131
2132   return fabs(angle) * state->mult;
2133 }
2134
2135 /*
2136 =item linear_interp(pos, seg)
2137
2138 Calculates linear interpolation on the fill parameter.  Breaks the
2139 segment into 2 regions based in the I<middle> value.
2140
2141 =cut
2142 */
2143 static double
2144 linear_interp(double pos, i_fountain_seg *seg) {
2145   if (pos < seg->middle) {
2146     double len = seg->middle - seg->start;
2147     if (len < EPSILON)
2148       return 0.0;
2149     else
2150       return (pos - seg->start) / len / 2;
2151   }
2152   else {
2153     double len = seg->end - seg->middle;
2154     if (len < EPSILON)
2155       return 1.0;
2156     else
2157       return 0.5 + (pos - seg->middle) / len / 2;
2158   }
2159 }
2160
2161 /*
2162 =item sine_interp(pos, seg)
2163
2164 Calculates sine function interpolation on the fill parameter.
2165
2166 =cut
2167 */
2168 static double
2169 sine_interp(double pos, i_fountain_seg *seg) {
2170   /* I wonder if there's a simple way to smooth the transition for this */
2171   double work = linear_interp(pos, seg);
2172
2173   return (1-cos(work * PI))/2;
2174 }
2175
2176 /*
2177 =item sphereup_interp(pos, seg)
2178
2179 Calculates spherical interpolation on the fill parameter, with the cusp 
2180 at the low-end.
2181
2182 =cut
2183 */
2184 static double
2185 sphereup_interp(double pos, i_fountain_seg *seg) {
2186   double work = linear_interp(pos, seg);
2187
2188   return sqrt(1.0 - (1-work) * (1-work));
2189 }
2190
2191 /*
2192 =item spheredown_interp(pos, seg)
2193
2194 Calculates spherical interpolation on the fill parameter, with the cusp 
2195 at the high-end.
2196
2197 =cut
2198 */
2199 static double
2200 spheredown_interp(double pos, i_fountain_seg *seg) {
2201   double work = linear_interp(pos, seg);
2202
2203   return 1-sqrt(1.0 - work * work);
2204 }
2205
2206 /*
2207 =item direct_cinterp(out, pos, seg)
2208
2209 Calculates the fountain color based on direct scaling of the channels
2210 of the color channels.
2211
2212 =cut
2213 */
2214 static void
2215 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2216   int ch;
2217   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2218     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
2219       + seg->c[1].channel[ch] * pos;
2220   }
2221 }
2222
2223 /*
2224 =item hue_up_cinterp(put, pos, seg)
2225
2226 Calculates the fountain color based on scaling a HSV value.  The hue
2227 increases as the fill parameter increases.
2228
2229 =cut
2230 */
2231 static void
2232 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2233   int ch;
2234   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2235     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
2236       + seg->c[1].channel[ch] * pos;
2237   }
2238   i_hsv_to_rgbf(out);
2239 }
2240
2241 /*
2242 =item hue_down_cinterp(put, pos, seg)
2243
2244 Calculates the fountain color based on scaling a HSV value.  The hue
2245 decreases as the fill parameter increases.
2246
2247 =cut
2248 */
2249 static void
2250 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2251   int ch;
2252   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2253     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
2254       + seg->c[1].channel[ch] * pos;
2255   }
2256   i_hsv_to_rgbf(out);
2257 }
2258
2259 /*
2260 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2261
2262 Simple grid-based super-sampling.
2263
2264 =cut
2265 */
2266 static int
2267 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
2268   i_fcolor *work = state->ssample_data;
2269   i_img_dim dx, dy;
2270   int grid = state->parm;
2271   double base = -0.5 + 0.5 / grid;
2272   double step = 1.0 / grid;
2273   int ch, i;
2274   int samp_count = 0;
2275
2276   for (dx = 0; dx < grid; ++dx) {
2277     for (dy = 0; dy < grid; ++dy) {
2278       if (fount_getat(work+samp_count, x + base + step * dx, 
2279                       y + base + step * dy, state)) {
2280         ++samp_count;
2281       }
2282     }
2283   }
2284   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2285     out->channel[ch] = 0;
2286     for (i = 0; i < samp_count; ++i) {
2287       out->channel[ch] += work[i].channel[ch];
2288     }
2289     /* we divide by 4 rather than samp_count since if there's only one valid
2290        sample it should be mostly transparent */
2291     out->channel[ch] /= grid * grid;
2292   }
2293   return samp_count;
2294 }
2295
2296 /*
2297 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2298
2299 Random super-sampling.
2300
2301 =cut
2302 */
2303 static int
2304 random_ssample(i_fcolor *out, double x, double y, 
2305                struct fount_state *state) {
2306   i_fcolor *work = state->ssample_data;
2307   int i, ch;
2308   int maxsamples = state->parm;
2309   double rand_scale = 1.0 / RAND_MAX;
2310   int samp_count = 0;
2311   for (i = 0; i < maxsamples; ++i) {
2312     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale, 
2313                     y - 0.5 + rand() * rand_scale, state)) {
2314       ++samp_count;
2315     }
2316   }
2317   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2318     out->channel[ch] = 0;
2319     for (i = 0; i < samp_count; ++i) {
2320       out->channel[ch] += work[i].channel[ch];
2321     }
2322     /* we divide by maxsamples rather than samp_count since if there's
2323        only one valid sample it should be mostly transparent */
2324     out->channel[ch] /= maxsamples;
2325   }
2326   return samp_count;
2327 }
2328
2329 /*
2330 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2331
2332 Super-sampling around the circumference of a circle.
2333
2334 I considered saving the sin()/cos() values and transforming step-size
2335 around the circle, but that's inaccurate, though it may not matter
2336 much.
2337
2338 =cut
2339  */
2340 static int
2341 circle_ssample(i_fcolor *out, double x, double y, 
2342                struct fount_state *state) {
2343   i_fcolor *work = state->ssample_data;
2344   int i, ch;
2345   int maxsamples = state->parm;
2346   double angle = 2 * PI / maxsamples;
2347   double radius = 0.3; /* semi-random */
2348   int samp_count = 0;
2349   for (i = 0; i < maxsamples; ++i) {
2350     if (fount_getat(work+samp_count, x + radius * cos(angle * i), 
2351                     y + radius * sin(angle * i), state)) {
2352       ++samp_count;
2353     }
2354   }
2355   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2356     out->channel[ch] = 0;
2357     for (i = 0; i < samp_count; ++i) {
2358       out->channel[ch] += work[i].channel[ch];
2359     }
2360     /* we divide by maxsamples rather than samp_count since if there's
2361        only one valid sample it should be mostly transparent */
2362     out->channel[ch] /= maxsamples;
2363   }
2364   return samp_count;
2365 }
2366
2367 /*
2368 =item fount_r_none(v)
2369
2370 Implements no repeats.  Simply clamps the fill value.
2371
2372 =cut
2373 */
2374 static double
2375 fount_r_none(double v) {
2376   return v < 0 ? 0 : v > 1 ? 1 : v;
2377 }
2378
2379 /*
2380 =item fount_r_sawtooth(v)
2381
2382 Implements sawtooth repeats.  Clamps negative values and uses fmod()
2383 on others.
2384
2385 =cut
2386 */
2387 static double
2388 fount_r_sawtooth(double v) {
2389   return v < 0 ? 0 : fmod(v, 1.0);
2390 }
2391
2392 /*
2393 =item fount_r_triangle(v)
2394
2395 Implements triangle repeats.  Clamps negative values, uses fmod to get
2396 a range 0 through 2 and then adjusts values > 1.
2397
2398 =cut
2399 */
2400 static double
2401 fount_r_triangle(double v) {
2402   if (v < 0)
2403     return 0;
2404   else {
2405     v = fmod(v, 2.0);
2406     return v > 1.0 ? 2.0 - v : v;
2407   }
2408 }
2409
2410 /*
2411 =item fount_r_saw_both(v)
2412
2413 Implements sawtooth repeats in the both postive and negative directions.
2414
2415 Adjusts the value to be postive and then just uses fmod().
2416
2417 =cut
2418 */
2419 static double
2420 fount_r_saw_both(double v) {
2421   if (v < 0)
2422     v += 1+(int)(-v);
2423   return fmod(v, 1.0);
2424 }
2425
2426 /*
2427 =item fount_r_tri_both(v)
2428
2429 Implements triangle repeats in the both postive and negative directions.
2430
2431 Uses fmod on the absolute value, and then adjusts values > 1.
2432
2433 =cut
2434 */
2435 static double
2436 fount_r_tri_both(double v) {
2437   v = fmod(fabs(v), 2.0);
2438   return v > 1.0 ? 2.0 - v : v;
2439 }
2440
2441 /*
2442 =item fill_fountf(fill, x, y, width, channels, data)
2443
2444 The fill function for fountain fills.
2445
2446 =cut
2447 */
2448 static void
2449 fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width,
2450             int channels, i_fcolor *data) {
2451   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2452   
2453   while (width--) {
2454     i_fcolor c;
2455     int got_one;
2456     
2457     if (f->state.ssfunc)
2458       got_one = f->state.ssfunc(&c, x, y, &f->state);
2459     else
2460       got_one = fount_getat(&c, x, y, &f->state);
2461
2462     if (got_one)
2463       *data++ = c;
2464     
2465     ++x;
2466   }
2467 }
2468
2469 /*
2470 =item fount_fill_destroy(fill)
2471
2472 =cut
2473 */
2474 static void
2475 fount_fill_destroy(i_fill_t *fill) {
2476   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2477   fount_finish_state(&f->state);
2478 }
2479
2480 /*
2481 =back
2482
2483 =head1 AUTHOR
2484
2485 Arnar M. Hrafnkelsson <addi@umich.edu>
2486
2487 Tony Cook <tony@develop-help.com> (i_fountain())
2488
2489 =head1 SEE ALSO
2490
2491 Imager(3)
2492
2493 =cut
2494 */