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