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