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