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