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