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)*(double)(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, c, ch;
  22   i_img_dim x, y;
  23   double pc;
  24   double *coeff;
  25   double res[MAXCHANNELS];
  26   i_img *timg;
  27   int radius, diameter;
  28
  29   mm_log((1,"i_gaussian(im %p, stdev %.2f)\n",im,stddev));
  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;
  63   for(y = 0; y < im->ysize; y++) {
  64     for(x = 0; x < im->xsize; x++) {
  65       pc=0.0;
  66       for(ch=0;ch<im->channels;ch++)
  67         res[ch]=0;
  68       for(c = 0;c < diameter; c++)
  69         if (IM_GPIX(im,x+c-radius,y,&rcolor)!=-1) {
  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;
  76         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
  77       }
  78       IM_PPIX(timg, x, y, &rcolor);
  79     }
  80   }
  81
  82   for(x = 0;x < im->xsize; x++) {
  83     for(y = 0; y < im->ysize; y++) {
  84       pc=0.0;
  85       for(ch=0; ch<im->channels; ch++)
  86         res[ch]=0;
  87       for(c=0; c < diameter; c++)
  88         if (IM_GPIX(timg, x, y+c-radius, &rcolor)!=-1) {
  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;
  95         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
  96       }
  97       IM_PPIX(im, x, y, &rcolor);
  98     }
  99   }
 100 #/code
 101   myfree(coeff);
 102   i_img_destroy(timg);
 103
 104   return 1;
 105 }
 106
 107
 108
 109
 110
 111
 112