libs/filters/GP_Sigma.c

   1 /*****************************************************************************
   2  * This file is part of gfxprim library.                                     *
   3  *                                                                           *
   4  * Gfxprim is free software; you can redistribute it and/or                  *
   5  * modify it under the terms of the GNU Lesser General Public                *
   6  * License as published by the Free Software Foundation; either              *
   7  * version 2.1 of the License, or (at your option) any later version.        *
   8  *                                                                           *
   9  * Gfxprim is distributed in the hope that it will be useful,                *
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU         *
  12  * Lesser General Public License for more details.                           *
  13  *                                                                           *
  14  * You should have received a copy of the GNU Lesser General Public          *
  15  * License along with gfxprim; if not, write to the Free Software            *
  16  * Foundation, Inc., 51 Franklin Street, Fifth Floor,                        *
  17  * Boston, MA  02110-1301  USA                                               *
  18  *                                                                           *
  19  * Copyright (C) 2009-2013 Cyril Hrubis <metan@ucw.cz>                       *
  20  *                                                                           *
  21  *****************************************************************************/
  22
  23 #include <errno.h>
  24 #include <string.h>
  25
  26 #include <core/GP_Pixmap.h>
  27 #include <core/GP_GetPutPixel.h>
  28 #include <core/GP_TempAlloc.h>
  29 #include <core/GP_Clamp.h>
  30 #include <core/GP_Common.h>
  31 #include <core/GP_Debug.h>
  32 #include <filters/GP_Sigma.h>
  33
  34 static int gp_filter_sigma_raw(const gp_pixmap *src,
  35                                gp_coord x_src, gp_coord y_src,
  36                                gp_size w_src, gp_size h_src,
  37                                gp_pixmap *dst,
  38                                gp_coord x_dst, gp_coord y_dst,
  39                                int xrad, int yrad,
  40                                unsigned int min, float sigma,
  41                                gp_progress_cb *callback)
  42 {
  43         int x, y;
  44         unsigned int x1, y1;
  45
  46         if (src->pixel_type != GP_PIXEL_RGB888) {
  47                 errno = ENOSYS;
  48                 return -1;
  49         }
  50
  51         GP_DEBUG(1, "Sigma Mean filter size %ux%u xrad=%u yrad=%u sigma=%.2f",
  52                  w_src, h_src, xrad, yrad, sigma);
  53
  54         unsigned int R_sigma = 255 * sigma;
  55         unsigned int G_sigma = 255 * sigma;
  56         unsigned int B_sigma = 255 * sigma;
  57
  58         unsigned int xdiam = 2 * xrad + 1;
  59         unsigned int ydiam = 2 * yrad + 1;
  60
  61         unsigned int w = w_src + xdiam;
  62         unsigned int size = w * ydiam;
  63
  64         gp_temp_alloc_create(temp, 3 * size * sizeof(unsigned int));
  65
  66         unsigned int *R = gp_temp_alloc_get(temp, size * sizeof(unsigned int));
  67         unsigned int *G = gp_temp_alloc_get(temp, size * sizeof(unsigned int));
  68         unsigned int *B = gp_temp_alloc_get(temp, size * sizeof(unsigned int));
  69
  70         /* prefil the sampled array */
  71         for (x = 0; x < (int)w; x++) {
  72                 int xi = GP_CLAMP(x_src + x - xrad, 0, (int)src->w - 1);
  73
  74                 for (y = 0; y < (int)ydiam; y++) {
  75                         int yi = GP_CLAMP(y_src + y - yrad, 0, (int)src->h - 1);
  76
  77                         gp_pixel pix = gp_getpixel_raw_24BPP(src, xi, yi);
  78
  79                         R[y * w + x] = GP_PIXEL_GET_R_RGB888(pix);
  80                         G[y * w + x] = GP_PIXEL_GET_G_RGB888(pix);
  81                         B[y * w + x] = GP_PIXEL_GET_B_RGB888(pix);
  82                 }
  83         }
  84
  85         unsigned int R_sum;
  86         unsigned int G_sum;
  87         unsigned int B_sum;
  88
  89         unsigned int R_ssum;
  90         unsigned int G_ssum;
  91         unsigned int B_ssum;
  92
  93         unsigned int R_cnt;
  94         unsigned int G_cnt;
  95         unsigned int B_cnt;
  96
  97         unsigned int cnt = xdiam * ydiam - 1;
  98
  99         /* center pixel ypsilon in the buffer */
 100         unsigned int yc = yrad;
 101         /* last sampled ypsilon in the buffer */
 102         unsigned int yl = 0;
 103
 104         /* Apply the sigma mean filter */
 105         for (y = 0; y < (int)h_src; y++) {
 106                 for (x = 0; x < (int)w_src; x++) {
 107                         /* Get center pixel */
 108                         unsigned int R_center = R[yc * w + x + xrad];
 109                         unsigned int G_center = G[yc * w + x + xrad];
 110                         unsigned int B_center = B[yc * w + x + xrad];
 111
 112                         /* Reset sum counters */
 113                         R_sum = 0;
 114                         G_sum = 0;
 115                         B_sum = 0;
 116
 117                         R_ssum = 0;
 118                         G_ssum = 0;
 119                         B_ssum = 0;
 120
 121                         R_cnt = 0;
 122                         G_cnt = 0;
 123                         B_cnt = 0;
 124
 125                         for (x1 = 0; x1 < xdiam; x1++) {
 126                                 for (y1 = 0; y1 < ydiam; y1++) {
 127                                         unsigned int R_cur = R[y1 * w + x + x1];
 128                                         unsigned int G_cur = G[y1 * w + x + x1];
 129                                         unsigned int B_cur = B[y1 * w + x + x1];
 130
 131                                         R_sum += R_cur;
 132                                         G_sum += G_cur;
 133                                         B_sum += B_cur;
 134
 135                                         if (abs(R_cur - R_center) < R_sigma) {
 136                                                 R_ssum += R_cur;
 137                                                 R_cnt++;
 138                                         }
 139
 140                                         if (abs(G_cur - G_center) < G_sigma) {
 141                                                 G_ssum += G_cur;
 142                                                 G_cnt++;
 143                                         }
 144
 145                                         if (abs(B_cur - B_center) < B_sigma) {
 146                                                 B_ssum += B_cur;
 147                                                 B_cnt++;
 148                                         }
 149                                 }
 150                         }
 151
 152                         R_sum -= R_center;
 153                         G_sum -= G_center;
 154                         B_sum -= B_center;
 155
 156                         unsigned int r;
 157                         unsigned int g;
 158                         unsigned int b;
 159
 160                         if (R_cnt >= min)
 161                                 r = R_ssum / R_cnt;
 162                         else
 163                                 r = R_sum / cnt;
 164
 165                         if (G_cnt >= min)
 166                                 g = G_ssum / G_cnt;
 167                         else
 168                                 g = G_sum / cnt;
 169
 170                         if (B_cnt >= min)
 171                                 b = B_ssum / B_cnt;
 172                         else
 173                                 b = B_sum / cnt;
 174
 175                         gp_putpixel_raw_24BPP(dst, x_dst + x, y_dst + y,
 176                                               GP_PIXEL_CREATE_RGB888(r, g, b));
 177                 }
 178
 179                 int yi = GP_CLAMP(y_src + y + yrad + 1, 0, (int)src->h - 1);
 180
 181                 for (x = 0; x < (int)w; x++) {
 182                         int xi = GP_CLAMP(x_src + x - xrad, 0, (int)src->w - 1);
 183
 184                         gp_pixel pix = gp_getpixel_raw_24BPP(src, xi, yi);
 185
 186                         R[yl * w + x] = GP_PIXEL_GET_R_RGB888(pix);
 187                         G[yl * w + x] = GP_PIXEL_GET_G_RGB888(pix);
 188                         B[yl * w + x] = GP_PIXEL_GET_B_RGB888(pix);
 189                 }
 190
 191                 yc = (yc+1) % ydiam;
 192                 yl = (yl+1) % ydiam;
 193
 194                 if (gp_progress_cb_report(callback, y, h_src, w_src)) {
 195                         gp_temp_alloc_free(temp);
 196                         return 1;
 197                 }
 198         }
 199
 200         gp_temp_alloc_free(temp);
 201         gp_progress_cb_done(callback);
 202
 203         return 0;
 204 }
 205
 206 int gp_filter_sigma_ex(const gp_pixmap *src,
 207                        gp_coord x_src, gp_coord y_src,
 208                        gp_size w_src, gp_size h_src,
 209                        gp_pixmap *dst,
 210                        gp_coord x_dst, gp_coord y_dst,
 211                        int xrad, int yrad,
 212                        unsigned int min, float sigma,
 213                        gp_progress_cb *callback)
 214 {
 215         GP_CHECK(src->pixel_type == dst->pixel_type);
 216
 217         /* Check that destination is large enough */
 218         GP_CHECK(x_dst + (gp_coord)w_src <= (gp_coord)dst->w);
 219         GP_CHECK(y_dst + (gp_coord)h_src <= (gp_coord)dst->h);
 220
 221         GP_CHECK(xrad >= 0 && yrad >= 0);
 222
 223         return gp_filter_sigma_raw(src, x_src, y_src, w_src, h_src,
 224                                    dst, x_dst, y_dst, xrad, yrad, min, sigma,
 225                                    callback);
 226 }
 227
 228 gp_pixmap *gp_filter_sigma_ex_alloc(const gp_pixmap *src,
 229                                     gp_coord x_src, gp_coord y_src,
 230                                     gp_size w_src, gp_size h_src,
 231                                     int xrad, int yrad,
 232                                     unsigned int min, float sigma,
 233                                     gp_progress_cb *callback)
 234 {
 235         int ret;
 236
 237         GP_CHECK(xrad >= 0 && yrad >= 0);
 238
 239         gp_pixmap *dst = gp_pixmap_alloc(w_src, h_src, src->pixel_type);
 240
 241         if (dst == NULL)
 242                 return NULL;
 243
 244         ret = gp_filter_sigma_raw(src, x_src, y_src, w_src, h_src,
 245                                   dst, 0, 0, xrad, yrad, min, sigma, callback);
 246
 247         if (ret) {
 248                 gp_pixmap_free(dst);
 249                 return NULL;
 250         }
 251
 252         return dst;
 253 }