hide or rename any symbols that are likely to conflict with other
[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]  = i_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   myfree(fdist);
983   
984 }
985
986 void
987 i_nearest_color_foo(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
988
989   int p, x, y;
990   int xsize    = im->xsize;
991   int ysize    = im->ysize;
992
993   mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
994   
995   for(p = 0; p<num; p++) {
996     mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
997     ICL_info(&ival[p]);
998   }
999
1000   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1001     int   midx    = 0;
1002     float mindist = 0;
1003     float curdist = 0;
1004
1005     int xd        = x-xo[0];
1006     int yd        = y-yo[0];
1007
1008     switch (dmeasure) {
1009     case 0: /* euclidean */
1010       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1011       break;
1012     case 1: /* euclidean squared */
1013       mindist = xd*xd + yd*yd; /* euclidean distance */
1014       break;
1015     case 2: /* euclidean squared */
1016       mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
1017       break;
1018     default:
1019       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1020     }
1021
1022     for(p = 1; p<num; p++) {
1023       xd = x-xo[p];
1024       yd = y-yo[p];
1025       switch (dmeasure) {
1026       case 0: /* euclidean */
1027         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1028         break;
1029       case 1: /* euclidean squared */
1030         curdist = xd*xd + yd*yd; /* euclidean distance */
1031         break;
1032       case 2: /* euclidean squared */
1033         curdist = i_max(xd*xd, yd*yd); /* manhattan distance */
1034         break;
1035       default:
1036         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1037       }
1038       if (curdist < mindist) {
1039         mindist = curdist;
1040         midx = p;
1041       }
1042     }
1043     i_ppix(im, x, y, &ival[midx]); 
1044   }
1045 }
1046
1047 void
1048 i_nearest_color(i_img *im, int num, int *xo, int *yo, i_color *oval, int dmeasure) {
1049   i_color *ival;
1050   float *tval;
1051   float c1, c2;
1052   i_color val;
1053   int p, x, y, ch;
1054   int xsize    = im->xsize;
1055   int ysize    = im->ysize;
1056   int *cmatch;
1057
1058   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));
1059
1060   tval   = mymalloc( sizeof(float)*num*im->channels );
1061   ival   = mymalloc( sizeof(i_color)*num );
1062   cmatch = mymalloc( sizeof(int)*num     );
1063
1064   for(p = 0; p<num; p++) {
1065     for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
1066     cmatch[p] = 0;
1067   }
1068
1069   
1070   for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1071     int   midx    = 0;
1072     float mindist = 0;
1073     float curdist = 0;
1074     
1075     int xd        = x-xo[0];
1076     int yd        = y-yo[0];
1077
1078     switch (dmeasure) {
1079     case 0: /* euclidean */
1080       mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1081       break;
1082     case 1: /* euclidean squared */
1083       mindist = xd*xd + yd*yd; /* euclidean distance */
1084       break;
1085     case 2: /* euclidean squared */
1086       mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
1087       break;
1088     default:
1089       m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1090     }
1091     
1092     for(p = 1; p<num; p++) {
1093       xd = x-xo[p];
1094       yd = y-yo[p];
1095       switch (dmeasure) {
1096       case 0: /* euclidean */
1097         curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1098         break;
1099       case 1: /* euclidean squared */
1100         curdist = xd*xd + yd*yd; /* euclidean distance */
1101         break;
1102       case 2: /* euclidean squared */
1103         curdist = i_max(xd*xd, yd*yd); /* manhattan distance */
1104         break;
1105       default:
1106         m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1107       }
1108       if (curdist < mindist) {
1109         mindist = curdist;
1110         midx = p;
1111       }
1112     }
1113
1114     cmatch[midx]++;
1115     i_gpix(im, x, y, &val);
1116     c2 = 1.0/(float)(cmatch[midx]);
1117     c1 = 1.0-c2;
1118     
1119     for(ch = 0; ch<im->channels; ch++) 
1120       tval[midx*im->channels + ch] = c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
1121   
1122     
1123   }
1124
1125   for(p = 0; p<num; p++) for(ch = 0; ch<im->channels; ch++) ival[p].channel[ch] = tval[p*im->channels + ch];
1126
1127   i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
1128 }
1129
1130 /*
1131 =item i_unsharp_mask(im, stddev, scale)
1132
1133 Perform an usharp mask, which is defined as subtracting the blurred
1134 image from double the original.
1135
1136 =cut
1137 */
1138 void i_unsharp_mask(i_img *im, double stddev, double scale) {
1139   i_img copy;
1140   int x, y, ch;
1141
1142   if (scale < 0)
1143     return;
1144   /* it really shouldn't ever be more than 1.0, but maybe ... */
1145   if (scale > 100)
1146     scale = 100;
1147
1148   i_copy(&copy, im);
1149   i_gaussian(&copy, stddev);
1150   if (im->bits == i_8_bits) {
1151     i_color *blur = mymalloc(im->xsize * sizeof(i_color) * 2);
1152     i_color *out = blur + im->xsize;
1153
1154     for (y = 0; y < im->ysize; ++y) {
1155       i_glin(&copy, 0, copy.xsize, y, blur);
1156       i_glin(im, 0, im->xsize, y, out);
1157       for (x = 0; x < im->xsize; ++x) {
1158         for (ch = 0; ch < im->channels; ++ch) {
1159           /*int temp = out[x].channel[ch] + 
1160             scale * (out[x].channel[ch] - blur[x].channel[ch]);*/
1161           int temp = out[x].channel[ch] * 2 - blur[x].channel[ch];
1162           if (temp < 0)
1163             temp = 0;
1164           else if (temp > 255)
1165             temp = 255;
1166           out[x].channel[ch] = temp;
1167         }
1168       }
1169       i_plin(im, 0, im->xsize, y, out);
1170     }
1171
1172     myfree(blur);
1173   }
1174   else {
1175     i_fcolor *blur = mymalloc(im->xsize * sizeof(i_fcolor) * 2);
1176     i_fcolor *out = blur + im->xsize;
1177
1178     for (y = 0; y < im->ysize; ++y) {
1179       i_glinf(&copy, 0, copy.xsize, y, blur);
1180       i_glinf(im, 0, im->xsize, y, out);
1181       for (x = 0; x < im->xsize; ++x) {
1182         for (ch = 0; ch < im->channels; ++ch) {
1183           double temp = out[x].channel[ch] +
1184             scale * (out[x].channel[ch] - blur[x].channel[ch]);
1185           if (temp < 0)
1186             temp = 0;
1187           else if (temp > 1.0)
1188             temp = 1.0;
1189           out[x].channel[ch] = temp;
1190         }
1191       }
1192       i_plinf(im, 0, im->xsize, y, out);
1193     }
1194
1195     myfree(blur);
1196   }
1197   i_img_exorcise(&copy);
1198 }
1199
1200 /*
1201 =item i_diff_image(im1, im2, mindiff)
1202
1203 Creates a new image that is transparent, except where the pixel in im2
1204 is different from im1, where it is the pixel from im2.
1205
1206 The samples must differ by at least mindiff to be considered different.
1207
1208 =cut
1209 */
1210
1211 i_img *
1212 i_diff_image(i_img *im1, i_img *im2, int mindiff) {
1213   i_img *out;
1214   int outchans, diffchans;
1215   int xsize, ysize;
1216   i_img temp;
1217
1218   i_clear_error();
1219   if (im1->channels != im2->channels) {
1220     i_push_error(0, "different number of channels");
1221     return NULL;
1222   }
1223
1224   outchans = diffchans = im1->channels;
1225   if (outchans == 1 || outchans == 3)
1226     ++outchans;
1227
1228   xsize = i_min(im1->xsize, im2->xsize);
1229   ysize = i_min(im1->ysize, im2->ysize);
1230
1231   out = i_sametype_chans(im1, xsize, ysize, outchans);
1232   
1233   if (im1->bits == i_8_bits && im2->bits == i_8_bits) {
1234     i_color *line1 = mymalloc(2 * xsize * sizeof(*line1));
1235     i_color *line2 = line1 + xsize;
1236     i_color empty;
1237     int x, y, ch;
1238
1239     for (ch = 0; ch < MAXCHANNELS; ++ch)
1240       empty.channel[ch] = 0;
1241
1242     for (y = 0; y < ysize; ++y) {
1243       i_glin(im1, 0, xsize, y, line1);
1244       i_glin(im2, 0, xsize, y, line2);
1245       for (x = 0; x < xsize; ++x) {
1246         int diff = 0;
1247         for (ch = 0; ch < diffchans; ++ch) {
1248           if (line1[x].channel[ch] != line2[x].channel[ch]
1249               && abs(line1[x].channel[ch] - line2[x].channel[ch]) > mindiff) {
1250             diff = 1;
1251             break;
1252           }
1253         }
1254         if (!diff)
1255           line2[x] = empty;
1256       }
1257       i_plin(out, 0, xsize, y, line2);
1258     }
1259     myfree(line1);
1260   }
1261   else {
1262     i_fcolor *line1 = mymalloc(2 * xsize * sizeof(*line1));
1263     i_fcolor *line2 = line1 + xsize;
1264     i_fcolor empty;
1265     int x, y, ch;
1266     double dist = mindiff / 255;
1267
1268     for (ch = 0; ch < MAXCHANNELS; ++ch)
1269       empty.channel[ch] = 0;
1270
1271     for (y = 0; y < ysize; ++y) {
1272       i_glinf(im1, 0, xsize, y, line1);
1273       i_glinf(im2, 0, xsize, y, line2);
1274       for (x = 0; x < xsize; ++x) {
1275         int diff = 0;
1276         for (ch = 0; ch < diffchans; ++ch) {
1277           if (line1[x].channel[ch] != line2[x].channel[ch]
1278               && abs(line1[x].channel[ch] - line2[x].channel[ch]) > dist) {
1279             diff = 1;
1280             break;
1281           }
1282         }
1283         if (!diff)
1284           line2[x] = empty;
1285       }
1286       i_plinf(out, 0, xsize, y, line2);
1287     }
1288     myfree(line1);
1289   }
1290
1291   return out;
1292 }
1293
1294 struct fount_state;
1295 static double linear_fount_f(double x, double y, struct fount_state *state);
1296 static double bilinear_fount_f(double x, double y, struct fount_state *state);
1297 static double radial_fount_f(double x, double y, struct fount_state *state);
1298 static double square_fount_f(double x, double y, struct fount_state *state);
1299 static double revolution_fount_f(double x, double y, 
1300                                  struct fount_state *state);
1301 static double conical_fount_f(double x, double y, struct fount_state *state);
1302
1303 typedef double (*fount_func)(double, double, struct fount_state *);
1304 static fount_func fount_funcs[] =
1305 {
1306   linear_fount_f,
1307   bilinear_fount_f,
1308   radial_fount_f,
1309   square_fount_f,
1310   revolution_fount_f,
1311   conical_fount_f,
1312 };
1313
1314 static double linear_interp(double pos, i_fountain_seg *seg);
1315 static double sine_interp(double pos, i_fountain_seg *seg);
1316 static double sphereup_interp(double pos, i_fountain_seg *seg);
1317 static double spheredown_interp(double pos, i_fountain_seg *seg);
1318 typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
1319 static fount_interp fount_interps[] =
1320 {
1321   linear_interp,
1322   linear_interp,
1323   sine_interp,
1324   sphereup_interp,
1325   spheredown_interp,
1326 };
1327
1328 static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1329 static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1330 static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1331 typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
1332 static fount_cinterp fount_cinterps[] =
1333 {
1334   direct_cinterp,
1335   hue_up_cinterp,
1336   hue_down_cinterp,
1337 };
1338
1339 typedef double (*fount_repeat)(double v);
1340 static double fount_r_none(double v);
1341 static double fount_r_sawtooth(double v);
1342 static double fount_r_triangle(double v);
1343 static double fount_r_saw_both(double v);
1344 static double fount_r_tri_both(double v);
1345 static fount_repeat fount_repeats[] =
1346 {
1347   fount_r_none,
1348   fount_r_sawtooth,
1349   fount_r_triangle,
1350   fount_r_saw_both,
1351   fount_r_tri_both,
1352 };
1353
1354 static int simple_ssample(i_fcolor *out, double x, double y, 
1355                            struct fount_state *state);
1356 static int random_ssample(i_fcolor *out, double x, double y, 
1357                            struct fount_state *state);
1358 static int circle_ssample(i_fcolor *out, double x, double y, 
1359                            struct fount_state *state);
1360 typedef int (*fount_ssample)(i_fcolor *out, double x, double y, 
1361                               struct fount_state *state);
1362 static fount_ssample fount_ssamples[] =
1363 {
1364   NULL,
1365   simple_ssample,
1366   random_ssample,
1367   circle_ssample,
1368 };
1369
1370 static int
1371 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
1372
1373 /*
1374   Keep state information used by each type of fountain fill
1375 */
1376 struct fount_state {
1377   /* precalculated for the equation of the line perpendicular to the line AB */
1378   double lA, lB, lC;
1379   double AB;
1380   double sqrtA2B2;
1381   double mult;
1382   double cos;
1383   double sin;
1384   double theta;
1385   int xa, ya;
1386   void *ssample_data;
1387   fount_func ffunc;
1388   fount_repeat rpfunc;
1389   fount_ssample ssfunc;
1390   double parm;
1391   i_fountain_seg *segs;
1392   int count;
1393 };
1394
1395 static void
1396 fount_init_state(struct fount_state *state, double xa, double ya, 
1397                  double xb, double yb, i_fountain_type type, 
1398                  i_fountain_repeat repeat, int combine, int super_sample, 
1399                  double ssample_param, int count, i_fountain_seg *segs);
1400
1401 static void
1402 fount_finish_state(struct fount_state *state);
1403
1404 #define EPSILON (1e-6)
1405
1406 /*
1407 =item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1408
1409 Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1410
1411 I<type> controls how the reference points are used:
1412
1413 =over
1414
1415 =item i_ft_linear
1416
1417 linear, where A is 0 and B is 1.
1418
1419 =item i_ft_bilinear
1420
1421 linear in both directions from A.
1422
1423 =item i_ft_radial
1424
1425 circular, where A is the centre of the fill, and B is a point
1426 on the radius.
1427
1428 =item i_ft_radial_square
1429
1430 where A is the centre of the fill and B is the centre of
1431 one side of the square.
1432
1433 =item i_ft_revolution
1434
1435 where A is the centre of the fill and B defines the 0/1.0
1436 angle of the fill.
1437
1438 =item i_ft_conical
1439
1440 similar to i_ft_revolution, except that the revolution goes in both
1441 directions
1442
1443 =back
1444
1445 I<repeat> can be one of:
1446
1447 =over
1448
1449 =item i_fr_none
1450
1451 values < 0 are treated as zero, values > 1 are treated as 1.
1452
1453 =item i_fr_sawtooth
1454
1455 negative values are treated as 0, positive values are modulo 1.0
1456
1457 =item i_fr_triangle
1458
1459 negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1460 mod 1.0), otherwise the same as for sawtooth.
1461
1462 =item i_fr_saw_both
1463
1464 like i_fr_sawtooth, except that the sawtooth pattern repeats into
1465 negative values.
1466
1467 =item i_fr_tri_both
1468
1469 Like i_fr_triangle, except that negative values are handled as their
1470 absolute values.
1471
1472 =back
1473
1474 If combine is non-zero then non-opaque values are combined with the
1475 underlying color.
1476
1477 I<super_sample> controls super sampling, if any.  At some point I'll
1478 probably add a adaptive super-sampler.  Current possible values are:
1479
1480 =over
1481
1482 =item i_fts_none
1483
1484 No super-sampling is done.
1485
1486 =item i_fts_grid
1487
1488 A square grid of points withing the pixel are sampled.
1489
1490 =item i_fts_random
1491
1492 Random points within the pixel are sampled.
1493
1494 =item i_fts_circle
1495
1496 Points on the radius of a circle are sampled.  This produces fairly
1497 good results, but is fairly slow since sin() and cos() are evaluated
1498 for each point.
1499
1500 =back
1501
1502 I<ssample_param> is intended to be roughly the number of points
1503 sampled within the pixel.
1504
1505 I<count> and I<segs> define the segments of the fill.
1506
1507 =cut
1508
1509 */
1510
1511 void
1512 i_fountain(i_img *im, double xa, double ya, double xb, double yb, 
1513            i_fountain_type type, i_fountain_repeat repeat, 
1514            int combine, int super_sample, double ssample_param, 
1515            int count, i_fountain_seg *segs) {
1516   struct fount_state state;
1517   int x, y;
1518   i_fcolor *line = mymalloc(sizeof(i_fcolor) * im->xsize);
1519   i_fcolor *work = NULL;
1520
1521   i_fountain_seg *my_segs;
1522   i_fill_combine_f combine_func = NULL;
1523   i_fill_combinef_f combinef_func = NULL;
1524
1525   i_get_combine(combine, &combine_func, &combinef_func);
1526   if (combinef_func)
1527     work = mymalloc(sizeof(i_fcolor) * im->xsize);
1528
1529   fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine, 
1530                    super_sample, ssample_param, count, segs);
1531   my_segs = state.segs;
1532
1533   for (y = 0; y < im->ysize; ++y) {
1534     i_glinf(im, 0, im->xsize, y, line);
1535     for (x = 0; x < im->xsize; ++x) {
1536       i_fcolor c;
1537       int got_one;
1538       if (super_sample == i_fts_none)
1539         got_one = fount_getat(&c, x, y, &state);
1540       else
1541         got_one = state.ssfunc(&c, x, y, &state);
1542       if (got_one) {
1543         if (combine)
1544           work[x] = c;
1545         else 
1546           line[x] = c;
1547       }
1548     }
1549     if (combine)
1550       combinef_func(line, work, im->channels, im->xsize);
1551     i_plinf(im, 0, im->xsize, y, line);
1552   }
1553   fount_finish_state(&state);
1554   if (work) myfree(work);
1555   myfree(line);
1556 }
1557
1558 typedef struct {
1559   i_fill_t base;
1560   struct fount_state state;
1561 } i_fill_fountain_t;
1562
1563 static void
1564 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels, 
1565             i_fcolor *data);
1566 static void
1567 fount_fill_destroy(i_fill_t *fill);
1568
1569 /*
1570 =item i_new_fount(xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1571
1572 Creates a new general fill which fills with a fountain fill.
1573
1574 =cut
1575 */
1576
1577 i_fill_t *
1578 i_new_fill_fount(double xa, double ya, double xb, double yb, 
1579                  i_fountain_type type, i_fountain_repeat repeat, 
1580                  int combine, int super_sample, double ssample_param, 
1581                  int count, i_fountain_seg *segs) {
1582   i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1583   
1584   fill->base.fill_with_color = NULL;
1585   fill->base.fill_with_fcolor = fill_fountf;
1586   fill->base.destroy = fount_fill_destroy;
1587   if (combine)
1588     i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1589   else {
1590     fill->base.combine = NULL;
1591     fill->base.combinef = NULL;
1592   }
1593   fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine, 
1594                    super_sample, ssample_param, count, segs);
1595
1596   return &fill->base;
1597 }
1598
1599 /*
1600 =back
1601
1602 =head1 INTERNAL FUNCTIONS
1603
1604 =over
1605
1606 =item fount_init_state(...)
1607
1608 Used by both the fountain fill filter and the fountain fill.
1609
1610 =cut
1611 */
1612
1613 static void
1614 fount_init_state(struct fount_state *state, double xa, double ya, 
1615                  double xb, double yb, i_fountain_type type, 
1616                  i_fountain_repeat repeat, int combine, int super_sample, 
1617                  double ssample_param, int count, i_fountain_seg *segs) {
1618   int i, j;
1619   i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count);
1620   /*int have_alpha = im->channels == 2 || im->channels == 4;*/
1621   
1622   memset(state, 0, sizeof(*state));
1623   /* we keep a local copy that we can adjust for speed */
1624   for (i = 0; i < count; ++i) {
1625     i_fountain_seg *seg = my_segs + i;
1626
1627     *seg = segs[i];
1628     if (seg->type < 0 || seg->type >= i_fst_end)
1629       seg->type = i_fst_linear;
1630     if (seg->color < 0 || seg->color >= i_fc_end)
1631       seg->color = i_fc_direct;
1632     if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1633       /* so we don't have to translate to HSV on each request, do it here */
1634       for (j = 0; j < 2; ++j) {
1635         i_rgb_to_hsvf(seg->c+j);
1636       }
1637       if (seg->color == i_fc_hue_up) {
1638         if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1639           seg->c[1].channel[0] += 1.0;
1640       }
1641       else {
1642         if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1643           seg->c[0].channel[0] += 1.0;
1644       }
1645     }
1646     /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n", 
1647            seg->start, seg->middle, seg->end, seg->c[0].channel[0], 
1648            seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1649            seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2], 
1650            seg->c[1].channel[3], seg->type, seg->color);*/
1651            
1652   }
1653
1654   /* initialize each engine */
1655   /* these are so common ... */
1656   state->lA = xb - xa;
1657   state->lB = yb - ya;
1658   state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1659   state->xa = xa;
1660   state->ya = ya;
1661   switch (type) {
1662   default:
1663     type = i_ft_linear; /* make the invalid value valid */
1664   case i_ft_linear:
1665   case i_ft_bilinear:
1666     state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1667     state->mult = 1;
1668     state->mult = 1/linear_fount_f(xb, yb, state);
1669     break;
1670
1671   case i_ft_radial:
1672     state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa) 
1673                              + (double)(yb-ya)*(yb-ya));
1674     break;
1675
1676   case i_ft_radial_square:
1677     state->cos = state->lA / state->AB;
1678     state->sin = state->lB / state->AB;
1679     state->mult = 1.0 / state->AB;
1680     break;
1681
1682   case i_ft_revolution:
1683     state->theta = atan2(yb-ya, xb-xa);
1684     state->mult = 1.0 / (PI * 2);
1685     break;
1686
1687   case i_ft_conical:
1688     state->theta = atan2(yb-ya, xb-xa);
1689     state->mult = 1.0 / PI;
1690     break;
1691   }
1692   state->ffunc = fount_funcs[type];
1693   if (super_sample < 0 
1694       || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
1695     super_sample = 0;
1696   }
1697   state->ssample_data = NULL;
1698   switch (super_sample) {
1699   case i_fts_grid:
1700     ssample_param = floor(0.5 + sqrt(ssample_param));
1701     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param);
1702     break;
1703
1704   case i_fts_random:
1705   case i_fts_circle:
1706     ssample_param = floor(0.5+ssample_param);
1707     state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
1708     break;
1709   }
1710   state->parm = ssample_param;
1711   state->ssfunc = fount_ssamples[super_sample];
1712   if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1713     repeat = 0;
1714   state->rpfunc = fount_repeats[repeat];
1715   state->segs = my_segs;
1716   state->count = count;
1717 }
1718
1719 static void
1720 fount_finish_state(struct fount_state *state) {
1721   if (state->ssample_data)
1722     myfree(state->ssample_data);
1723   myfree(state->segs);
1724 }
1725
1726
1727 /*
1728 =item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
1729
1730 Evaluates the fountain fill at the given point.
1731
1732 This is called by both the non-super-sampling and super-sampling code.
1733
1734 You might think that it would make sense to sample the fill parameter
1735 instead, and combine those, but this breaks badly.
1736
1737 =cut
1738 */
1739
1740 static int
1741 fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
1742   double v = (state->rpfunc)((state->ffunc)(x, y, state));
1743   int i;
1744
1745   i = 0;
1746   while (i < state->count 
1747          && (v < state->segs[i].start || v > state->segs[i].end)) {
1748     ++i;
1749   }
1750   if (i < state->count) {
1751     v = (fount_interps[state->segs[i].type])(v, state->segs+i);
1752     (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
1753     return 1;
1754   }
1755   else
1756     return 0;
1757 }
1758
1759 /*
1760 =item linear_fount_f(x, y, state)
1761
1762 Calculate the fill parameter for a linear fountain fill.
1763
1764 Uses the point to line distance function, with some precalculation
1765 done in i_fountain().
1766
1767 =cut
1768 */
1769 static double
1770 linear_fount_f(double x, double y, struct fount_state *state) {
1771   return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
1772 }
1773
1774 /*
1775 =item bilinear_fount_f(x, y, state)
1776
1777 Calculate the fill parameter for a bi-linear fountain fill.
1778
1779 =cut
1780 */
1781 static double
1782 bilinear_fount_f(double x, double y, struct fount_state *state) {
1783   return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
1784 }
1785
1786 /*
1787 =item radial_fount_f(x, y, state)
1788
1789 Calculate the fill parameter for a radial fountain fill.
1790
1791 Simply uses the distance function.
1792
1793 =cut
1794  */
1795 static double
1796 radial_fount_f(double x, double y, struct fount_state *state) {
1797   return sqrt((double)(state->xa-x)*(state->xa-x) 
1798               + (double)(state->ya-y)*(state->ya-y)) * state->mult;
1799 }
1800
1801 /*
1802 =item square_fount_f(x, y, state)
1803
1804 Calculate the fill parameter for a square fountain fill.
1805
1806 Works by rotating the reference co-ordinate around the centre of the
1807 square.
1808
1809 =cut
1810 */
1811 static double
1812 square_fount_f(double x, double y, struct fount_state *state) {
1813   int xc, yc; /* centred on A */
1814   double xt, yt; /* rotated by theta */
1815   xc = x - state->xa;
1816   yc = y - state->ya;
1817   xt = fabs(xc * state->cos + yc * state->sin);
1818   yt = fabs(-xc * state->sin + yc * state->cos);
1819   return (xt > yt ? xt : yt) * state->mult;
1820 }
1821
1822 /*
1823 =item revolution_fount_f(x, y, state)
1824
1825 Calculates the fill parameter for the revolution fountain fill.
1826
1827 =cut
1828 */
1829 static double
1830 revolution_fount_f(double x, double y, struct fount_state *state) {
1831   double angle = atan2(y - state->ya, x - state->xa);
1832   
1833   angle -= state->theta;
1834   if (angle < 0) {
1835     angle = fmod(angle+ PI * 4, PI*2);
1836   }
1837
1838   return angle * state->mult;
1839 }
1840
1841 /*
1842 =item conical_fount_f(x, y, state)
1843
1844 Calculates the fill parameter for the conical fountain fill.
1845
1846 =cut
1847 */
1848 static double
1849 conical_fount_f(double x, double y, struct fount_state *state) {
1850   double angle = atan2(y - state->ya, x - state->xa);
1851   
1852   angle -= state->theta;
1853   if (angle < -PI)
1854     angle += PI * 2;
1855   else if (angle > PI) 
1856     angle -= PI * 2;
1857
1858   return fabs(angle) * state->mult;
1859 }
1860
1861 /*
1862 =item linear_interp(pos, seg)
1863
1864 Calculates linear interpolation on the fill parameter.  Breaks the
1865 segment into 2 regions based in the I<middle> value.
1866
1867 =cut
1868 */
1869 static double
1870 linear_interp(double pos, i_fountain_seg *seg) {
1871   if (pos < seg->middle) {
1872     double len = seg->middle - seg->start;
1873     if (len < EPSILON)
1874       return 0.0;
1875     else
1876       return (pos - seg->start) / len / 2;
1877   }
1878   else {
1879     double len = seg->end - seg->middle;
1880     if (len < EPSILON)
1881       return 1.0;
1882     else
1883       return 0.5 + (pos - seg->middle) / len / 2;
1884   }
1885 }
1886
1887 /*
1888 =item sine_interp(pos, seg)
1889
1890 Calculates sine function interpolation on the fill parameter.
1891
1892 =cut
1893 */
1894 static double
1895 sine_interp(double pos, i_fountain_seg *seg) {
1896   /* I wonder if there's a simple way to smooth the transition for this */
1897   double work = linear_interp(pos, seg);
1898
1899   return (1-cos(work * PI))/2;
1900 }
1901
1902 /*
1903 =item sphereup_interp(pos, seg)
1904
1905 Calculates spherical interpolation on the fill parameter, with the cusp 
1906 at the low-end.
1907
1908 =cut
1909 */
1910 static double
1911 sphereup_interp(double pos, i_fountain_seg *seg) {
1912   double work = linear_interp(pos, seg);
1913
1914   return sqrt(1.0 - (1-work) * (1-work));
1915 }
1916
1917 /*
1918 =item spheredown_interp(pos, seg)
1919
1920 Calculates spherical interpolation on the fill parameter, with the cusp 
1921 at the high-end.
1922
1923 =cut
1924 */
1925 static double
1926 spheredown_interp(double pos, i_fountain_seg *seg) {
1927   double work = linear_interp(pos, seg);
1928
1929   return 1-sqrt(1.0 - work * work);
1930 }
1931
1932 /*
1933 =item direct_cinterp(out, pos, seg)
1934
1935 Calculates the fountain color based on direct scaling of the channels
1936 of the color channels.
1937
1938 =cut
1939 */
1940 static void
1941 direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1942   int ch;
1943   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1944     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
1945       + seg->c[1].channel[ch] * pos;
1946   }
1947 }
1948
1949 /*
1950 =item hue_up_cinterp(put, pos, seg)
1951
1952 Calculates the fountain color based on scaling a HSV value.  The hue
1953 increases as the fill parameter increases.
1954
1955 =cut
1956 */
1957 static void
1958 hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1959   int ch;
1960   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1961     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
1962       + seg->c[1].channel[ch] * pos;
1963   }
1964   i_hsv_to_rgbf(out);
1965 }
1966
1967 /*
1968 =item hue_down_cinterp(put, pos, seg)
1969
1970 Calculates the fountain color based on scaling a HSV value.  The hue
1971 decreases as the fill parameter increases.
1972
1973 =cut
1974 */
1975 static void
1976 hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1977   int ch;
1978   for (ch = 0; ch < MAXCHANNELS; ++ch) {
1979     out->channel[ch] = seg->c[0].channel[ch] * (1 - pos) 
1980       + seg->c[1].channel[ch] * pos;
1981   }
1982   i_hsv_to_rgbf(out);
1983 }
1984
1985 /*
1986 =item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1987
1988 Simple grid-based super-sampling.
1989
1990 =cut
1991 */
1992 static int
1993 simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
1994   i_fcolor *work = state->ssample_data;
1995   int dx, dy;
1996   int grid = state->parm;
1997   double base = -0.5 + 0.5 / grid;
1998   double step = 1.0 / grid;
1999   int ch, i;
2000   int samp_count = 0;
2001
2002   for (dx = 0; dx < grid; ++dx) {
2003     for (dy = 0; dy < grid; ++dy) {
2004       if (fount_getat(work+samp_count, x + base + step * dx, 
2005                       y + base + step * dy, state)) {
2006         ++samp_count;
2007       }
2008     }
2009   }
2010   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2011     out->channel[ch] = 0;
2012     for (i = 0; i < samp_count; ++i) {
2013       out->channel[ch] += work[i].channel[ch];
2014     }
2015     /* we divide by 4 rather than samp_count since if there's only one valid
2016        sample it should be mostly transparent */
2017     out->channel[ch] /= grid * grid;
2018   }
2019   return samp_count;
2020 }
2021
2022 /*
2023 =item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2024
2025 Random super-sampling.
2026
2027 =cut
2028 */
2029 static int
2030 random_ssample(i_fcolor *out, double x, double y, 
2031                struct fount_state *state) {
2032   i_fcolor *work = state->ssample_data;
2033   int i, ch;
2034   int maxsamples = state->parm;
2035   double rand_scale = 1.0 / RAND_MAX;
2036   int samp_count = 0;
2037   for (i = 0; i < maxsamples; ++i) {
2038     if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale, 
2039                     y - 0.5 + rand() * rand_scale, state)) {
2040       ++samp_count;
2041     }
2042   }
2043   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2044     out->channel[ch] = 0;
2045     for (i = 0; i < samp_count; ++i) {
2046       out->channel[ch] += work[i].channel[ch];
2047     }
2048     /* we divide by maxsamples rather than samp_count since if there's
2049        only one valid sample it should be mostly transparent */
2050     out->channel[ch] /= maxsamples;
2051   }
2052   return samp_count;
2053 }
2054
2055 /*
2056 =item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2057
2058 Super-sampling around the circumference of a circle.
2059
2060 I considered saving the sin()/cos() values and transforming step-size
2061 around the circle, but that's inaccurate, though it may not matter
2062 much.
2063
2064 =cut
2065  */
2066 static int
2067 circle_ssample(i_fcolor *out, double x, double y, 
2068                struct fount_state *state) {
2069   i_fcolor *work = state->ssample_data;
2070   int i, ch;
2071   int maxsamples = state->parm;
2072   double angle = 2 * PI / maxsamples;
2073   double radius = 0.3; /* semi-random */
2074   int samp_count = 0;
2075   for (i = 0; i < maxsamples; ++i) {
2076     if (fount_getat(work+samp_count, x + radius * cos(angle * i), 
2077                     y + radius * sin(angle * i), state)) {
2078       ++samp_count;
2079     }
2080   }
2081   for (ch = 0; ch < MAXCHANNELS; ++ch) {
2082     out->channel[ch] = 0;
2083     for (i = 0; i < samp_count; ++i) {
2084       out->channel[ch] += work[i].channel[ch];
2085     }
2086     /* we divide by maxsamples rather than samp_count since if there's
2087        only one valid sample it should be mostly transparent */
2088     out->channel[ch] /= maxsamples;
2089   }
2090   return samp_count;
2091 }
2092
2093 /*
2094 =item fount_r_none(v)
2095
2096 Implements no repeats.  Simply clamps the fill value.
2097
2098 =cut
2099 */
2100 static double
2101 fount_r_none(double v) {
2102   return v < 0 ? 0 : v > 1 ? 1 : v;
2103 }
2104
2105 /*
2106 =item fount_r_sawtooth(v)
2107
2108 Implements sawtooth repeats.  Clamps negative values and uses fmod()
2109 on others.
2110
2111 =cut
2112 */
2113 static double
2114 fount_r_sawtooth(double v) {
2115   return v < 0 ? 0 : fmod(v, 1.0);
2116 }
2117
2118 /*
2119 =item fount_r_triangle(v)
2120
2121 Implements triangle repeats.  Clamps negative values, uses fmod to get
2122 a range 0 through 2 and then adjusts values > 1.
2123
2124 =cut
2125 */
2126 static double
2127 fount_r_triangle(double v) {
2128   if (v < 0)
2129     return 0;
2130   else {
2131     v = fmod(v, 2.0);
2132     return v > 1.0 ? 2.0 - v : v;
2133   }
2134 }
2135
2136 /*
2137 =item fount_r_saw_both(v)
2138
2139 Implements sawtooth repeats in the both postive and negative directions.
2140
2141 Adjusts the value to be postive and then just uses fmod().
2142
2143 =cut
2144 */
2145 static double
2146 fount_r_saw_both(double v) {
2147   if (v < 0)
2148     v += 1+(int)(-v);
2149   return fmod(v, 1.0);
2150 }
2151
2152 /*
2153 =item fount_r_tri_both(v)
2154
2155 Implements triangle repeats in the both postive and negative directions.
2156
2157 Uses fmod on the absolute value, and then adjusts values > 1.
2158
2159 =cut
2160 */
2161 static double
2162 fount_r_tri_both(double v) {
2163   v = fmod(fabs(v), 2.0);
2164   return v > 1.0 ? 2.0 - v : v;
2165 }
2166
2167 /*
2168 =item fill_fountf(fill, x, y, width, channels, data)
2169
2170 The fill function for fountain fills.
2171
2172 =cut
2173 */
2174 static void
2175 fill_fountf(i_fill_t *fill, int x, int y, int width, int channels, 
2176             i_fcolor *data) {
2177   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2178   
2179   while (width--) {
2180     i_fcolor c;
2181     int got_one;
2182     
2183     if (f->state.ssfunc)
2184       got_one = f->state.ssfunc(&c, x, y, &f->state);
2185     else
2186       got_one = fount_getat(&c, x, y, &f->state);
2187     
2188     *data++ = c;
2189     
2190     ++x;
2191   }
2192 }
2193
2194 /*
2195 =item fount_fill_destroy(fill)
2196
2197 =cut
2198 */
2199 static void
2200 fount_fill_destroy(i_fill_t *fill) {
2201   i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2202   fount_finish_state(&f->state);
2203 }
2204
2205 /*
2206 =back
2207
2208 =head1 AUTHOR
2209
2210 Arnar M. Hrafnkelsson <addi@umich.edu>
2211
2212 Tony Cook <tony@develop-help.com> (i_fountain())
2213
2214 =head1 SEE ALSO
2215
2216 Imager(3)
2217
2218 =cut
2219 */