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