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
  25 #include "core/GP_Pixmap.h"
  26 #include "core/GP_GetPutPixel.h"
  27 #include "core/GP_TempAlloc.h"
  28 #include "core/GP_Clamp.h"
  29 #include "core/GP_Common.h"
  30
  31 #include "core/GP_Debug.h"
  32
  33 #include "GP_Sigma.h"
  34
  35 #include <string.h>
  36
  37 static int GP_FilterSigma_Raw(const GP_Pixmap *src,
  38                               GP_Coord x_src, GP_Coord y_src,
  39                               GP_Size w_src, GP_Size h_src,
  40                               GP_Pixmap *dst,
  41                               GP_Coord x_dst, GP_Coord y_dst,
  42                               int xrad, int yrad,
  43                               unsigned int min, float sigma,
  44                               GP_ProgressCallback *callback)
  45 {
  46         int x, y;
  47         unsigned int x1, y1;
  48
  49         if (src->pixel_type != GP_PIXEL_RGB888) {
  50                 errno = ENOSYS;
  51                 return -1;
  52         }
  53
  54         GP_DEBUG(1, "Sigma Mean filter size %ux%u xrad=%u yrad=%u sigma=%.2f",
  55                  w_src, h_src, xrad, yrad, sigma);
  56
  57         unsigned int R_sigma = 255 * sigma;
  58         unsigned int G_sigma = 255 * sigma;
  59         unsigned int B_sigma = 255 * sigma;
  60
  61         unsigned int xdiam = 2 * xrad + 1;
  62         unsigned int ydiam = 2 * yrad + 1;
  63
  64         unsigned int w = w_src + xdiam;
  65         unsigned int size = w * ydiam;
  66
  67         GP_TempAllocCreate(temp, 3 * size * sizeof(unsigned int));
  68
  69         unsigned int *R = GP_TempAllocGet(temp, size * sizeof(unsigned int));
  70         unsigned int *G = GP_TempAllocGet(temp, size * sizeof(unsigned int));
  71         unsigned int *B = GP_TempAllocGet(temp, size * sizeof(unsigned int));
  72
  73         /* prefil the sampled array */
  74         for (x = 0; x < (int)w; x++) {
  75                 int xi = GP_CLAMP(x_src + x - xrad, 0, (int)src->w - 1);
  76
  77                 for (y = 0; y < (int)ydiam; y++) {
  78                         int yi = GP_CLAMP(y_src + y - yrad, 0, (int)src->h - 1);
  79
  80                         GP_Pixel pix = GP_GetPixel_Raw_24BPP(src, xi, yi);
  81
  82                         R[y * w + x] = GP_Pixel_GET_R_RGB888(pix);
  83                         G[y * w + x] = GP_Pixel_GET_G_RGB888(pix);
  84                         B[y * w + x] = GP_Pixel_GET_B_RGB888(pix);
  85                 }
  86         }
  87
  88         unsigned int R_sum;
  89         unsigned int G_sum;
  90         unsigned int B_sum;
  91
  92         unsigned int R_ssum;
  93         unsigned int G_ssum;
  94         unsigned int B_ssum;
  95
  96         unsigned int R_cnt;
  97         unsigned int G_cnt;
  98         unsigned int B_cnt;
  99
 100         unsigned int cnt = xdiam * ydiam - 1;
 101
 102         /* center pixel ypsilon in the buffer */
 103         unsigned int yc = yrad;
 104         /* last sampled ypsilon in the buffer */
 105         unsigned int yl = 0;
 106
 107         /* Apply the sigma mean filter */
 108         for (y = 0; y < (int)h_src; y++) {
 109                 for (x = 0; x < (int)w_src; x++) {
 110                         /* Get center pixel */
 111                         unsigned int R_center = R[yc * w + x + xrad];
 112                         unsigned int G_center = G[yc * w + x + xrad];
 113                         unsigned int B_center = B[yc * w + x + xrad];
 114
 115                         /* Reset sum counters */
 116                         R_sum = 0;
 117                         G_sum = 0;
 118                         B_sum = 0;
 119
 120                         R_ssum = 0;
 121                         G_ssum = 0;
 122                         B_ssum = 0;
 123
 124                         R_cnt = 0;
 125                         G_cnt = 0;
 126                         B_cnt = 0;
 127
 128                         for (x1 = 0; x1 < xdiam; x1++) {
 129                                 for (y1 = 0; y1 < ydiam; y1++) {
 130                                         unsigned int R_cur = R[y1 * w + x + x1];
 131                                         unsigned int G_cur = G[y1 * w + x + x1];
 132                                         unsigned int B_cur = B[y1 * w + x + x1];
 133
 134                                         R_sum += R_cur;
 135                                         G_sum += G_cur;
 136                                         B_sum += B_cur;
 137
 138                                         if (abs(R_cur - R_center) < R_sigma) {
 139                                                 R_ssum += R_cur;
 140                                                 R_cnt++;
 141                                         }
 142
 143                                         if (abs(G_cur - G_center) < G_sigma) {
 144                                                 G_ssum += G_cur;
 145                                                 G_cnt++;
 146                                         }
 147
 148                                         if (abs(B_cur - B_center) < B_sigma) {
 149                                                 B_ssum += B_cur;
 150                                                 B_cnt++;
 151                                         }
 152                                 }
 153                         }
 154
 155                         R_sum -= R_center;
 156                         G_sum -= G_center;
 157                         B_sum -= B_center;
 158
 159                         unsigned int r;
 160                         unsigned int g;
 161                         unsigned int b;
 162
 163                         if (R_cnt >= min)
 164                                 r = R_ssum / R_cnt;
 165                         else
 166                                 r = R_sum / cnt;
 167
 168                         if (G_cnt >= min)
 169                                 g = G_ssum / G_cnt;
 170                         else
 171                                 g = G_sum / cnt;
 172
 173                         if (B_cnt >= min)
 174                                 b = B_ssum / B_cnt;
 175                         else
 176                                 b = B_sum / cnt;
 177
 178                         GP_PutPixel_Raw_24BPP(dst, x_dst + x, y_dst + y,
 179                                               GP_Pixel_CREATE_RGB888(r, g, b));
 180                 }
 181
 182                 int yi = GP_CLAMP(y_src + y + yrad + 1, 0, (int)src->h - 1);
 183
 184                 for (x = 0; x < (int)w; x++) {
 185                         int xi = GP_CLAMP(x_src + x - xrad, 0, (int)src->w - 1);
 186
 187                         GP_Pixel pix = GP_GetPixel_Raw_24BPP(src, xi, yi);
 188
 189                         R[yl * w + x] = GP_Pixel_GET_R_RGB888(pix);
 190                         G[yl * w + x] = GP_Pixel_GET_G_RGB888(pix);
 191                         B[yl * w + x] = GP_Pixel_GET_B_RGB888(pix);
 192                 }
 193
 194                 yc = (yc+1) % ydiam;
 195                 yl = (yl+1) % ydiam;
 196
 197                 if (GP_ProgressCallbackReport(callback, y, h_src, w_src)) {
 198                         GP_TempAllocFree(temp);
 199                         return 1;
 200                 }
 201         }
 202
 203         GP_TempAllocFree(temp);
 204         GP_ProgressCallbackDone(callback);
 205
 206         return 0;
 207 }
 208
 209 int GP_FilterSigmaEx(const GP_Pixmap *src,
 210                      GP_Coord x_src, GP_Coord y_src,
 211                      GP_Size w_src, GP_Size h_src,
 212                      GP_Pixmap *dst,
 213                      GP_Coord x_dst, GP_Coord y_dst,
 214                      int xrad, int yrad,
 215                      unsigned int min, float sigma,
 216                      GP_ProgressCallback *callback)
 217 {
 218         GP_CHECK(src->pixel_type == dst->pixel_type);
 219
 220         /* Check that destination is large enough */
 221         GP_CHECK(x_dst + (GP_Coord)w_src <= (GP_Coord)dst->w);
 222         GP_CHECK(y_dst + (GP_Coord)h_src <= (GP_Coord)dst->h);
 223
 224         GP_CHECK(xrad >= 0 && yrad >= 0);
 225
 226         return GP_FilterSigma_Raw(src, x_src, y_src, w_src, h_src,
 227                                   dst, x_dst, y_dst, xrad, yrad, min, sigma,
 228                                   callback);
 229 }
 230
 231 GP_Pixmap *GP_FilterSigmaExAlloc(const GP_Pixmap *src,
 232                                   GP_Coord x_src, GP_Coord y_src,
 233                                   GP_Size w_src, GP_Size h_src,
 234                                   int xrad, int yrad,
 235                                   unsigned int min, float sigma,
 236                                   GP_ProgressCallback *callback)
 237 {
 238         int ret;
 239
 240         GP_CHECK(xrad >= 0 && yrad >= 0);
 241
 242         GP_Pixmap *dst = GP_PixmapAlloc(w_src, h_src, src->pixel_type);
 243
 244         if (dst == NULL)
 245                 return NULL;
 246
 247         ret = GP_FilterSigma_Raw(src, x_src, y_src, w_src, h_src,
 248                                  dst, 0, 0, xrad, yrad, min, sigma, callback);
 249
 250         if (ret) {
 251                 GP_PixmapFree(dst);
 252                 return NULL;
 253         }
 254
 255         return dst;
 256 }