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