correct an old bug link
[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     if (im_size_t_max / sizeof(i_fcolor) > ssample_param) {
1976       bytes = sizeof(i_fcolor) * ssample_param;
1977       state->ssample_data = mymalloc(bytes);
1978     }
1979     else {
1980       dIMCTX;
1981       im_log((aIMCTX, 1,"size_t overflow calculating super-sample array size for random or circl"));
1982       super_sample = i_fts_none;
1983     }
1984     break;
1985   }
1986   state->parm = ssample_param;
1987   state->ssfunc = fount_ssamples[super_sample];
1988   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1989     repeat = 0;
1990   state->rpfunc = fount_repeats[repeat];
1991   state->segs = my_segs;
1992   state->count = count;
1993 }
1994
1995 static void
1996 fount_finish_state(struct fount_state *state) {
1997   if (state->ssample_data)
1998     myfree(state->ssample_data);
1999   myfree(state->segs);
2000 }
2001
2002
2003 /*
2004 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
2005
2006 Evaluates the fountain fill at the given point.
2007
2008 This is called by both the non-super-sampling and super-sampling code.
2009
2010 You might think that it would make sense to sample the fill parameter
2011 instead, and combine those, but this breaks badly.
2012
2013 =cut
2014 */
2015
2016 static int
2017 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
2018   double v = (state->rpfunc)((state->ffunc)(x, y, state));
2019   int i;
2020
2021   i = 0;
2022   while (i < state->count 
2023          && (v < state->segs[i].start || v > state->segs[i].end)) {
2024     ++i;
2025   }
2026   if (i < state->count) {
2027     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
2028     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
2029     return 1;
2030   }
2031   else
2032     return 0;
2033 }
2034
2035 /*
2036 =item linear_fount_f(x, y, state)
2037
2038 Calculate the fill parameter for a linear fountain fill.
2039
2040 Uses the point to line distance function, with some precalculation
2041 done in i_fountain().
2042
2043 =cut
2044 */
2045 static double
2046 linear_fount_f(double x, double y, struct fount_state *state) {
2047   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
2048 }
2049
2050 /*
2051 =item bilinear_fount_f(x, y, state)
2052
2053 Calculate the fill parameter for a bi-linear fountain fill.
2054
2055 =cut
2056 */
2057 static double
2058 bilinear_fount_f(double x, double y, struct fount_state *state) {
2059   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
2060 }
2061
2062 /*
2063 =item radial_fount_f(x, y, state)
2064
2065 Calculate the fill parameter for a radial fountain fill.
2066
2067 Simply uses the distance function.
2068
2069 =cut
2070  */
2071 static double
2072 radial_fount_f(double x, double y, struct fount_state *state) {
2073   return sqrt((double)(state->xa-x)*(state->xa-x) 
2074               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
2075 }
2076
2077 /*
2078 =item square_fount_f(x, y, state)
2079
2080 Calculate the fill parameter for a square fountain fill.
2081
2082 Works by rotating the reference co-ordinate around the centre of the
2083 square.
2084
2085 =cut
2086 */
2087 static double
2088 square_fount_f(double x, double y, struct fount_state *state) {
2089   i_img_dim xc, yc; /* centred on A */
2090   double xt, yt; /* rotated by theta */
2091   xc = x - state->xa;
2092   yc = y - state->ya;
2093   xt = fabs(xc * state->cos + yc * state->sin);
2094   yt = fabs(-xc * state->sin + yc * state->cos);
2095   return (xt > yt ? xt : yt) * state->mult;
2096 }
2097
2098 /*
2099 =item revolution_fount_f(x, y, state)
2100
2101 Calculates the fill parameter for the revolution fountain fill.
2102
2103 =cut
2104 */
2105 static double
2106 revolution_fount_f(double x, double y, struct fount_state *state) {
2107   double angle = atan2(y - state->ya, x - state->xa);
2108   
2109   angle -= state->theta;
2110   if (angle < 0) {
2111     angle = fmod(angle+ PI * 4, PI*2);
2112   }
2113
2114   return angle * state->mult;
2115 }
2116
2117 /*
2118 =item conical_fount_f(x, y, state)
2119
2120 Calculates the fill parameter for the conical fountain fill.
2121
2122 =cut
2123 */
2124 static double
2125 conical_fount_f(double x, double y, struct fount_state *state) {
2126   double angle = atan2(y - state->ya, x - state->xa);
2127   
2128   angle -= state->theta;
2129   if (angle < -PI)
2130     angle += PI * 2;
2131   else if (angle > PI) 
2132     angle -= PI * 2;
2133
2134   return fabs(angle) * state->mult;
2135 }
2136
2137 /*
2138 =item linear_interp(pos, seg)
2139
2140 Calculates linear interpolation on the fill parameter.  Breaks the
2141 segment into 2 regions based in the I<middle> value.
2142
2143 =cut
2144 */
2145 static double
2146 linear_interp(double pos, i_fountain_seg *seg) {
2147   if (pos < seg->middle) {
2148     double len = seg->middle - seg->start;
2149     if (len < EPSILON)
2150       return 0.0;
2151     else
2152       return (pos - seg->start) / len / 2;
2153   }
2154   else {
2155     double len = seg->end - seg->middle;
2156     if (len < EPSILON)
2157       return 1.0;
2158     else
2159       return 0.5 + (pos - seg->middle) / len / 2;
2160   }
2161 }
2162
2163 /*
2164 =item sine_interp(pos, seg)
2165
2166 Calculates sine function interpolation on the fill parameter.
2167
2168 =cut
2169 */
2170 static double
2171 sine_interp(double pos, i_fountain_seg *seg) {
2172   /* I wonder if there's a simple way to smooth the transition for this */
2173   double work = linear_interp(pos, seg);
2174
2175   return (1-cos(work * PI))/2;
2176 }
2177
2178 /*
2179 =item sphereup_interp(pos, seg)
2180
2181 Calculates spherical interpolation on the fill parameter, with the cusp 
2182 at the low-end.
2183
2184 =cut
2185 */
2186 static double
2187 sphereup_interp(double pos, i_fountain_seg *seg) {
2188   double work = linear_interp(pos, seg);
2189
2190   return sqrt(1.0 - (1-work) * (1-work));
2191 }
2192
2193 /*
2194 =item spheredown_interp(pos, seg)
2195
2196 Calculates spherical interpolation on the fill parameter, with the cusp 
2197 at the high-end.
2198
2199 =cut
2200 */
2201 static double
2202 spheredown_interp(double pos, i_fountain_seg *seg) {
2203   double work = linear_interp(pos, seg);
2204
2205   return 1-sqrt(1.0 - work * work);
2206 }
2207
2208 /*
2209 =item direct_cinterp(out, pos, seg)
2210
2211 Calculates the fountain color based on direct scaling of the channels
2212 of the color channels.
2213
2214 =cut
2215 */
2216 static void
2217 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2218   int ch;
2219   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2220     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
2221       + seg->c[1].channel[ch] * pos;
2222   }
2223 }
2224
2225 /*
2226 =item hue_up_cinterp(put, pos, seg)
2227
2228 Calculates the fountain color based on scaling a HSV value.  The hue
2229 increases as the fill parameter increases.
2230
2231 =cut
2232 */
2233 static void
2234 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2235   int ch;
2236   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2237     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
2238       + seg->c[1].channel[ch] * pos;
2239   }
2240   i_hsv_to_rgbf(out);
2241 }
2242
2243 /*
2244 =item hue_down_cinterp(put, pos, seg)
2245
2246 Calculates the fountain color based on scaling a HSV value.  The hue
2247 decreases as the fill parameter increases.
2248
2249 =cut
2250 */
2251 static void
2252 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2253   int ch;
2254   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2255     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
2256       + seg->c[1].channel[ch] * pos;
2257   }
2258   i_hsv_to_rgbf(out);
2259 }
2260
2261 /*
2262 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2263
2264 Simple grid-based super-sampling.
2265
2266 =cut
2267 */
2268 static int
2269 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
2270   i_fcolor *work = state->ssample_data;
2271   i_img_dim dx, dy;
2272   int grid = state->parm;
2273   double base = -0.5 + 0.5 / grid;
2274   double step = 1.0 / grid;
2275   int ch, i;
2276   int samp_count = 0;
2277
2278   for (dx = 0; dx < grid; ++dx) {
2279     for (dy = 0; dy < grid; ++dy) {
2280       if (fount_getat(work+samp_count, x + base + step * dx, 
2281                       y + base + step * dy, state)) {
2282         ++samp_count;
2283       }
2284     }
2285   }
2286   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2287     out->channel[ch] = 0;
2288     for (i = 0; i < samp_count; ++i) {
2289       out->channel[ch] += work[i].channel[ch];
2290     }
2291     /* we divide by 4 rather than samp_count since if there's only one valid
2292        sample it should be mostly transparent */
2293     out->channel[ch] /= grid * grid;
2294   }
2295   return samp_count;
2296 }
2297
2298 /*
2299 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2300
2301 Random super-sampling.
2302
2303 =cut
2304 */
2305 static int
2306 random_ssample(i_fcolor *out, double x, double y, 
2307                struct fount_state *state) {
2308   i_fcolor *work = state->ssample_data;
2309   int i, ch;
2310   int maxsamples = state->parm;
2311   double rand_scale = 1.0 / RAND_MAX;
2312   int samp_count = 0;
2313   for (i = 0; i < maxsamples; ++i) {
2314     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale, 
2315                     y - 0.5 + rand() * rand_scale, state)) {
2316       ++samp_count;
2317     }
2318   }
2319   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2320     out->channel[ch] = 0;
2321     for (i = 0; i < samp_count; ++i) {
2322       out->channel[ch] += work[i].channel[ch];
2323     }
2324     /* we divide by maxsamples rather than samp_count since if there's
2325        only one valid sample it should be mostly transparent */
2326     out->channel[ch] /= maxsamples;
2327   }
2328   return samp_count;
2329 }
2330
2331 /*
2332 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2333
2334 Super-sampling around the circumference of a circle.
2335
2336 I considered saving the sin()/cos() values and transforming step-size
2337 around the circle, but that's inaccurate, though it may not matter
2338 much.
2339
2340 =cut
2341  */
2342 static int
2343 circle_ssample(i_fcolor *out, double x, double y, 
2344                struct fount_state *state) {
2345   i_fcolor *work = state->ssample_data;
2346   int i, ch;
2347   int maxsamples = state->parm;
2348   double angle = 2 * PI / maxsamples;
2349   double radius = 0.3; /* semi-random */
2350   int samp_count = 0;
2351   for (i = 0; i < maxsamples; ++i) {
2352     if (fount_getat(work+samp_count, x + radius * cos(angle * i), 
2353                     y + radius * sin(angle * i), state)) {
2354       ++samp_count;
2355     }
2356   }
2357   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2358     out->channel[ch] = 0;
2359     for (i = 0; i < samp_count; ++i) {
2360       out->channel[ch] += work[i].channel[ch];
2361     }
2362     /* we divide by maxsamples rather than samp_count since if there's
2363        only one valid sample it should be mostly transparent */
2364     out->channel[ch] /= maxsamples;
2365   }
2366   return samp_count;
2367 }
2368
2369 /*
2370 =item fount_r_none(v)
2371
2372 Implements no repeats.  Simply clamps the fill value.
2373
2374 =cut
2375 */
2376 static double
2377 fount_r_none(double v) {
2378   return v < 0 ? 0 : v > 1 ? 1 : v;
2379 }
2380
2381 /*
2382 =item fount_r_sawtooth(v)
2383
2384 Implements sawtooth repeats.  Clamps negative values and uses fmod()
2385 on others.
2386
2387 =cut
2388 */
2389 static double
2390 fount_r_sawtooth(double v) {
2391   return v < 0 ? 0 : fmod(v, 1.0);
2392 }
2393
2394 /*
2395 =item fount_r_triangle(v)
2396
2397 Implements triangle repeats.  Clamps negative values, uses fmod to get
2398 a range 0 through 2 and then adjusts values > 1.
2399
2400 =cut
2401 */
2402 static double
2403 fount_r_triangle(double v) {
2404   if (v < 0)
2405     return 0;
2406   else {
2407     v = fmod(v, 2.0);
2408     return v > 1.0 ? 2.0 - v : v;
2409   }
2410 }
2411
2412 /*
2413 =item fount_r_saw_both(v)
2414
2415 Implements sawtooth repeats in the both postive and negative directions.
2416
2417 Adjusts the value to be postive and then just uses fmod().
2418
2419 =cut
2420 */
2421 static double
2422 fount_r_saw_both(double v) {
2423   if (v < 0)
2424     v += 1+(int)(-v);
2425   return fmod(v, 1.0);
2426 }
2427
2428 /*
2429 =item fount_r_tri_both(v)
2430
2431 Implements triangle repeats in the both postive and negative directions.
2432
2433 Uses fmod on the absolute value, and then adjusts values > 1.
2434
2435 =cut
2436 */
2437 static double
2438 fount_r_tri_both(double v) {
2439   v = fmod(fabs(v), 2.0);
2440   return v > 1.0 ? 2.0 - v : v;
2441 }
2442
2443 /*
2444 =item fill_fountf(fill, x, y, width, channels, data)
2445
2446 The fill function for fountain fills.
2447
2448 =cut
2449 */
2450 static void
2451 fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width,
2452             int channels, i_fcolor *data) {
2453   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2454   
2455   while (width--) {
2456     i_fcolor c;
2457     int got_one;
2458     
2459     if (f->state.ssfunc)
2460       got_one = f->state.ssfunc(&c, x, y, &f->state);
2461     else
2462       got_one = fount_getat(&c, x, y, &f->state);
2463
2464     if (got_one)
2465       *data++ = c;
2466     
2467     ++x;
2468   }
2469 }
2470
2471 /*
2472 =item fount_fill_destroy(fill)
2473
2474 =cut
2475 */
2476 static void
2477 fount_fill_destroy(i_fill_t *fill) {
2478   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2479   fount_finish_state(&f->state);
2480 }
2481
2482 /*
2483 =back
2484
2485 =head1 AUTHOR
2486
2487 Arnar M. Hrafnkelsson <addi@umich.edu>
2488
2489 Tony Cook <tony@develop-help.com> (i_fountain())
2490
2491 =head1 SEE ALSO
2492
2493 Imager(3)
2494
2495 =cut
2496 */