merge in tiff re-work branch
[imager.git] / gaussian.im
1 #include "imager.h"
2 #include <math.h>
3
4 static double
5 gauss(int x, double std) {
6   return 1.0/(sqrt(2.0*PI)*std)*exp(-(double)(x)*(float)(x)/(2*std*std));
7 }
8
9 /* Counters are as follows
10  l:  lines
11  i:  columns
12  c:  filter coeffs
13  ch: channels
14  pc: coeff equalization
15 */
16
17
18
19 int
20 i_gaussian(i_img *im, double stddev) {
21   int i,l,c,ch;
22   double pc;
23   double *coeff;
24   double res[MAXCHANNELS];
25   i_img *timg;
26   int radius, diameter;
27
28   mm_log((1,"i_gaussian(im %p, stdev %.2f)\n",im,stddev));
29   i_clear_error();
30
31   if (stddev <= 0) {
32     i_push_error(0, "stddev must be positive");
33     return 0;
34   }
35   /* totally silly cutoff */
36   if (stddev > 1000) {
37     stddev = 1000;
38   }
39  
40   timg = i_sametype(im, im->xsize, im->ysize);
41
42   if (im->bits <= 8)
43     radius = ceil(2 * stddev);
44   else
45     radius = ceil(3 * stddev);
46
47   diameter = 1 + radius * 2;
48
49   coeff = mymalloc(sizeof(double) * diameter);
50
51   for(i=0;i <= radius;i++) 
52     coeff[radius + i]=coeff[radius - i]=gauss(i, stddev);
53   pc=0;
54   for(i=0; i < diameter; i++)
55     pc+=coeff[i];
56   for(i=0;i < diameter;i++)
57     coeff[i] /= pc;
58
59
60 #code im->bits <= 8
61   IM_COLOR rcolor;
62   for(l=0;l<im->ysize;l++) {
63     for(i=0;i<im->xsize;i++) {
64       pc=0.0;
65       for(ch=0;ch<im->channels;ch++) 
66         res[ch]=0; 
67       for(c = 0;c < diameter; c++)
68         if (IM_GPIX(im,i+c-radius,l,&rcolor)!=-1) {
69           for(ch=0;ch<im->channels;ch++) 
70             res[ch]+= rcolor.channel[ch] * coeff[c];
71           pc+=coeff[c];
72         }
73       for(ch=0;ch<im->channels;ch++) {
74         double value = res[ch] / pc;
75         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
76       }
77       IM_PPIX(timg,i,l,&rcolor);
78     }
79   }
80   
81   for(l=0;l<im->xsize;l++) {
82     for(i=0;i<im->ysize;i++) {
83       pc=0.0;
84       for(ch=0; ch<im->channels; ch++)
85         res[ch]=0; 
86       for(c=0; c < diameter; c++)
87         if (IM_GPIX(timg,l,i+c-radius,&rcolor)!=-1) {
88           for(ch=0;ch<im->channels;ch++) 
89             res[ch]+= rcolor.channel[ch] * coeff[c];
90           pc+=coeff[c];
91         }
92       for(ch=0;ch<im->channels;ch++) {
93         double value = res[ch]/pc;
94         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
95       }
96       IM_PPIX(im,l,i,&rcolor);
97     }
98   }
99 #/code
100   myfree(coeff);
101   i_img_destroy(timg);
102   
103   return 1;
104 }
105
106
107
108
109
110
111