gaussian-blur.c 8.02 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/* This file is an image processing operation for GEGL
 * GEGL is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * GEGL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GEGL; if not, write to the
 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 * Boston, MA 02110-1301, USA.
 *
 * Copyright 2006 Dominik Ernst <dernst@gmx.de>
18 19 20
 * 
 * Recursive Gauss IIR Filter as described by Young / van Vliet
 * in "Signal Processing 44 (1995) 139 - 151"
21
 */
22

23
#ifdef GEGL_CHANT_SELF
24
 
25 26
gegl_chant_double (radius_x, -0.91675, 500.0, 4.0, "blur radius in horizontal direction.")
gegl_chant_double (radius_y, -0.91675, 500.0, 4.0, "blur radius in vertical direction.")
27 28 29

#else

30 31 32 33 34 35
#define GEGL_CHANT_FILTER
#define GEGL_CHANT_NAME            gaussian_blur
#define GEGL_CHANT_DESCRIPTION     "Performs an averaging of neighbouring pixels with the normal distribution as weighting."
#define GEGL_CHANT_SELF            "gaussian-blur.c"
#define GEGL_CHANT_CATEGORIES      "blur"
#define GEGL_CHANT_CLASS_CONSTRUCT
36 37 38
#include "gegl-chant.h"

#include <math.h>
39
#include <stdio.h>
40 41 42 43

static void
hor_blur (GeglBuffer *src,
          GeglBuffer *dst,
44 45
          gdouble     B,
          gdouble    *b);
46 47 48 49

static void
ver_blur (GeglBuffer *src,
          GeglBuffer *dst,
50 51 52 53 54 55 56
          gdouble     B,
          gdouble    *b);

static void
find_iir_constants (gfloat   radius, 
                    gdouble *B,
                    gdouble *b);
57 58 59 60 61 62 63 64 65 66 67 68 69 70




/* Actual image processing code
 ************************************************************************/
static gboolean
evaluate (GeglOperation *operation,
          const gchar   *output_prop)
{
  GeglOperationFilter *filter;
  ChantInstance *self;

  filter = GEGL_OPERATION_FILTER (operation);
71
  self   = GEGL_CHANT_INSTANCE (operation);
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113

  GeglBuffer *input  = filter->input;
  GeglBuffer *output;

  if(strcmp("output", output_prop))
    return FALSE;

    {
      GeglRect   *result = gegl_operation_result_rect (operation);
      GeglBuffer *temp_in;
      GeglBuffer *temp;

      if (result->w==0 ||
          result->h==0)
        {
          output = g_object_ref (input);
        }
      else
        {
          temp_in = g_object_new (GEGL_TYPE_BUFFER,
                                 "source", input,
                                 "x",      result->x,
                                 "y",      result->y,
                                 "width",  result->w,
                                 "height", result->h,
                                 NULL);
          temp   = g_object_new (GEGL_TYPE_BUFFER,
                                 "format", gegl_buffer_get_format (input),
                                 "x",      result->x,
                                 "y",      result->y,
                                 "width",  result->w,
                                 "height", result->h,
                                 NULL);

          output = g_object_new (GEGL_TYPE_BUFFER,
                                 "format", gegl_buffer_get_format (input),
                                 "x",      result->x,
                                 "y",      result->y,
                                 "width",  result->w,
                                 "height", result->h,
                                 NULL);

114 115
          {
            gdouble B, b[4];
116

117 118 119 120 121 122
            find_iir_constants (self->radius_x, &B, b);
            hor_blur (temp_in, temp,   B, b);
            
            find_iir_constants (self->radius_y, &B, b);
            ver_blur (temp,    output, B, b);
          }
123 124 125 126 127 128 129 130 131 132 133 134
          
          g_object_unref (temp);
          g_object_unref (temp_in);
        }

      if (filter->output)
        g_object_unref (filter->output);
      filter->output = output;
    }
  return  TRUE;
}

135 136 137 138
static void
find_iir_constants (gfloat   radius, 
                    gdouble *B,
                    gdouble *b)
139
{
140 141
  gdouble sigma0, q;
  
142
  sigma0 = (radius/3.0)+0.305584;
143 144 145 146 147 148 149 150 151 152

  if(sigma0 >= 2.5)
    q = 0.98711*sigma0 - 0.96330;
  else
    q = 3.97156 - 4.14554*sqrt(1-0.26891*sigma0);

  b[0] = 1.57825 + (2.44413*q) + (1.4281*q*q) + (0.422205*q*q*q);
  b[1] = (2.44413*q) + (2.85619*q*q) + (1.26661*q*q*q);
  b[2] = -((1.4281*q*q) + (1.26661*q*q*q));
  b[3] = 0.422205*q*q*q;
153

154 155 156 157 158 159 160 161 162 163 164 165 166 167
  *B = 1 - ( (b[1]+b[2]+b[3])/b[0] );
}

static inline void
blur_1D (gfloat  * buf,
         gint      offset,
         gint      delta_offset,
         gdouble   B,
         gdouble * b,
         gfloat  * w,
         gint      w_len)
{
  gint wcount, i;
  gdouble tmp;
168

169 170 171 172
  /* forward filter */
  wcount = 0;
  
  while (wcount < w_len)
173
    {
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
      tmp = 0;

      for (i=1; i<4; i++)
        {
          if (wcount-i >= 0)
            tmp += b[i]*w[wcount-i];
        }

      tmp /= b[0];
      tmp += B*buf[offset];
      w[wcount] = tmp;

      wcount++;
      offset += delta_offset;
    }

  /* backward filter */
  wcount = w_len - 1;
  offset -= delta_offset;

  while (wcount >= 0)
    {
      tmp = 0;

      for (i=1; i<4; i++)
        {
          if (wcount+i < w_len)
            tmp += b[i]*buf[offset+delta_offset*i];
        }

      tmp /= b[0];
      tmp += B*w[wcount];
      buf[offset] = tmp;

      offset -= delta_offset;
      wcount--;
210 211 212 213 214 215
    }
}

static void
hor_blur (GeglBuffer *src,
          GeglBuffer *dst,
216 217
          gdouble     B,
          gdouble    *b)
218
{
219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
  gint v;
  gint c;
  gint w_len;
  gfloat *buf;
  gfloat *w;

  buf = g_malloc0 (src->width * src->height * 4 * 4);
  w   = g_malloc0 (src->width * 4);

  gegl_buffer_get_fmt (src, buf, babl_format ("RaGaBaA float"));

  w_len = src->width;
  for (v=0; v<src->height; v++)
    {
      for (c = 0; c < 4; c++)
        { 
          blur_1D (buf,
                   v*src->width*4+c,
                   4,
                   B,
                   b,
                   w,
                   w_len);
        }
    }
244

245 246 247 248
  gegl_buffer_set_fmt (dst, buf, babl_format ("RaGaBaA float"));
  g_free (buf);
  g_free (w);
}
249 250 251 252

static void
ver_blur (GeglBuffer *src,
          GeglBuffer *dst,
253 254
          gdouble     B,
          gdouble    *b)
255
{
256 257 258 259 260
  gint v;
  gint c;
  gint w_len;
  gfloat *buf;
  gfloat *w;
261

262 263 264 265 266 267
  buf = g_malloc0 (src->width * src->height * 4 * 4);
  w   = g_malloc0 (src->height * 4);

  gegl_buffer_get_fmt (src, buf, babl_format ("RaGaBaA float"));

  w_len = src->height;
268
  
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
  for (v=0; v<src->width; v++)
    {
      for (c = 0; c < 4; c++)
        {
          blur_1D (buf,
                   v*4 + c,
                   src->width*4,
                   B,
                   b,
                   w,
                   w_len);
        }
    }

  gegl_buffer_set_fmt (dst, buf, babl_format ("RaGaBaA float"));
  g_free (buf);
  g_free (w);
286 287 288 289 290 291
}

static gboolean
calc_have_rect (GeglOperation *operation)
{
  GeglRect *in_rect = gegl_operation_get_have_rect (operation, "input");
292
  ChantInstance* self = GEGL_CHANT_INSTANCE(operation);
293 294 295
  gint radius_x       = ceil(self->radius_x+0.5);
  gint radius_y       = ceil(self->radius_y+0.5);
  if(!in_rect)
296 297 298
    return FALSE;

  gegl_operation_set_have_rect (operation, 
299 300
     in_rect->x-radius_x,   in_rect->y-radius_y,
     in_rect->w+radius_x*2, in_rect->h+radius_y*2);
301 302 303 304 305 306 307 308
  return TRUE;
}


static gboolean
calc_need_rect (GeglOperation *self)
{
  GeglRect  *need   = gegl_operation_need_rect (self);
309
  ChantInstance *blur = GEGL_CHANT_INSTANCE(self);
310 311
  gint radius_x = ceil(blur->radius_x+0.5);
  gint radius_y = ceil(blur->radius_y+0.5);
312 313

  gegl_operation_set_need_rect (self, "input",
314 315
     need->x-radius_x,   need->y-radius_y,
     need->w+radius_x*2, need->h+radius_y*2);
316 317 318 319 320 321 322 323 324 325 326
  return TRUE;
}

static void class_init (GeglOperationClass *operation_class)
{
  operation_class->calc_have_rect = calc_have_rect;
  operation_class->calc_need_rect = calc_need_rect;
}


#endif