libs/filters/GP_WeightedMedian.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
  30 #include "core/GP_Debug.h"
  31
  32 #include "GP_WeightedMedian.h"
  33
  34 #include <string.h>
  35
  36 static unsigned int sum_weights(GP_MedianWeights *weights)
  37 {
  38         unsigned int i;
  39         unsigned int sum = 0;
  40
  41         for (i = 0; i < weights->w * weights->h; i++)
  42                 sum += weights->weights[i];
  43
  44         return sum;
  45 }
  46
  47 static inline void hist_add(unsigned int *hist, unsigned int val,
  48                             unsigned int count)
  49 {
  50         hist[val] += count;
  51 }
  52
  53 static inline unsigned int hist_med(unsigned int *hist, unsigned int size,
  54                                     unsigned int threshold)
  55 {
  56         unsigned int i;
  57         unsigned int acc = 0;
  58
  59         for (i = 0; i < size; i++) {
  60                 acc += hist[i];
  61                 if (acc >= threshold)
  62                         return i;
  63         }
  64
  65         GP_BUG("Threshold not reached");
  66         return 0;
  67 }
  68
  69 static inline void hist_clear(unsigned int *hist, unsigned int size)
  70 {
  71         memset(hist, 0, sizeof(unsigned int) * size);
  72 }
  73
  74 static inline unsigned int get_weight(GP_MedianWeights *weights,
  75                                       unsigned int x, unsigned int y)
  76 {
  77         return weights->weights[y * weights->w + x];
  78 }
  79
  80 static int GP_FilterWeightedMedian_Raw(const GP_Pixmap *src,
  81                                        GP_Coord x_src, GP_Coord y_src,
  82                                        GP_Size w_src, GP_Size h_src,
  83                                        GP_Pixmap *dst,
  84                                        GP_Coord x_dst, GP_Coord y_dst,
  85                                        GP_MedianWeights *weights,
  86                                        GP_ProgressCallback *callback)
  87 {
  88         int x, y, sum = sum_weights(weights);
  89         unsigned int x1, y1;
  90
  91         if (src->pixel_type != GP_PIXEL_RGB888) {
  92                 errno = ENOSYS;
  93                 return -1;
  94         }
  95
  96         GP_DEBUG(1, "Weighted Median filter size %ux%u xmed=%u ymed=%u sum=%u",
  97                     w_src, h_src, weights->w, weights->h, sum);
  98
  99         unsigned int w = w_src +  weights->w;
 100         unsigned int size = w * weights->h;
 101
 102         GP_TempAllocCreate(temp, 3 * size * sizeof(unsigned int));
 103
 104         unsigned int *R = GP_TempAllocGet(temp, size * sizeof(unsigned int));
 105         unsigned int *G = GP_TempAllocGet(temp, size * sizeof(unsigned int));
 106         unsigned int *B = GP_TempAllocGet(temp, size * sizeof(unsigned int));
 107
 108         /* prefil the sampled array */
 109         for (x = 0; x < (int)w; x++) {
 110                 int xi = GP_CLAMP(x_src + x - (int)weights->w/2, 0, (int)src->w - 1);
 111
 112                 for (y = 0; y < (int)weights->h; y++) {
 113                         int yi = GP_CLAMP(y_src + y - (int)weights->h, 0, (int)src->h - 1);
 114
 115                         GP_Pixel pix = GP_GetPixel_Raw_24BPP(src, xi, yi);
 116
 117                         R[y * w + x] = GP_Pixel_GET_R_RGB888(pix);
 118                         G[y * w + x] = GP_Pixel_GET_G_RGB888(pix);
 119                         B[y * w + x] = GP_Pixel_GET_B_RGB888(pix);
 120                 }
 121         }
 122
 123         unsigned int hist_R[256];
 124         unsigned int hist_G[256];
 125         unsigned int hist_B[256];
 126
 127         hist_clear(hist_R, 256);
 128         hist_clear(hist_G, 256);
 129         hist_clear(hist_B, 256);
 130
 131         /* Apply the weighted median filter */
 132         for (y = 0; y < (int)h_src; y++) {
 133                 for (x = 0; x < (int)w_src; x++) {
 134                         /* compute weighted histogram and then median */
 135                         for (x1 = 0; x1 < weights->w; x1++) {
 136                                 for (y1 = 0; y1 < weights->h; y1++) {
 137                                         unsigned int weight = get_weight(weights, x1, y1);
 138                                         hist_add(hist_R, R[y1 * w + x + x1], weight);
 139                                         hist_add(hist_G, G[y1 * w + x + x1], weight);
 140                                         hist_add(hist_B, B[y1 * w + x + x1], weight);
 141                                 }
 142                         }
 143
 144                         unsigned int r = hist_med(hist_R, 256, sum/2);
 145                         unsigned int g = hist_med(hist_G, 256, sum/2);
 146                         unsigned int b = hist_med(hist_B, 256, sum/2);
 147
 148                         GP_PutPixel_Raw_24BPP(dst, x_dst + x, y_dst + y,
 149                                               GP_Pixel_CREATE_RGB888(r, g, b));
 150
 151                         hist_clear(hist_R, 256);
 152                         hist_clear(hist_G, 256);
 153                         hist_clear(hist_B, 256);
 154                 }
 155
 156                 for (x = 0; x < (int)w; x++) {
 157                         int xi = GP_CLAMP(x_src + x - (int)weights->w/2, 0, (int)src->w - 1);
 158
 159                         for (y1 = 0; y1 < weights->h; y1++) {
 160                                 int yi = GP_CLAMP(y_src + y + (int)y1 - (int)weights->h/2, 0, (int)src->h - 1);
 161
 162                                 GP_Pixel pix = GP_GetPixel_Raw_24BPP(src, xi, yi);
 163
 164                                 R[y1 * w + x] = GP_Pixel_GET_R_RGB888(pix);
 165                                 G[y1 * w + x] = GP_Pixel_GET_G_RGB888(pix);
 166                                 B[y1 * w + x] = GP_Pixel_GET_B_RGB888(pix);
 167                         }
 168                 }
 169
 170                 if (GP_ProgressCallbackReport(callback, y, h_src, w_src)) {
 171                         GP_TempAllocFree(temp);
 172                         return 1;
 173                 }
 174         }
 175
 176         GP_TempAllocFree(temp);
 177         GP_ProgressCallbackDone(callback);
 178
 179         return 0;
 180 }
 181
 182 int GP_FilterWeightedMedianEx(const GP_Pixmap *src,
 183                               GP_Coord x_src, GP_Coord y_src,
 184                               GP_Size w_src, GP_Size h_src,
 185                               GP_Pixmap *dst,
 186                               GP_Coord x_dst, GP_Coord y_dst,
 187                               GP_MedianWeights *weights,
 188                               GP_ProgressCallback *callback)
 189 {
 190         GP_CHECK(src->pixel_type == dst->pixel_type);
 191
 192         /* Check that destination is large enough */
 193         GP_CHECK(x_dst + (GP_Coord)w_src <= (GP_Coord)dst->w);
 194         GP_CHECK(y_dst + (GP_Coord)h_src <= (GP_Coord)dst->h);
 195
 196         //GP_CHECK(xmed >= 0 && ymed >= 0);
 197
 198         return GP_FilterWeightedMedian_Raw(src, x_src, y_src, w_src, h_src,
 199                                            dst, x_dst, y_dst, weights, callback);
 200 }
 201
 202 GP_Pixmap *GP_FilterWeightedMedianExAlloc(const GP_Pixmap *src,
 203                                            GP_Coord x_src, GP_Coord y_src,
 204                                            GP_Size w_src, GP_Size h_src,
 205                                            GP_MedianWeights *weights,
 206                                            GP_ProgressCallback *callback)
 207 {
 208         int ret;
 209
 210         //GP_CHECK(xmed >= 0 && ymed >= 0);
 211
 212         GP_Pixmap *dst = GP_PixmapAlloc(w_src, h_src, src->pixel_type);
 213
 214         if (dst == NULL)
 215                 return NULL;
 216
 217         ret = GP_FilterWeightedMedian_Raw(src, x_src, y_src, w_src, h_src,
 218                                           dst, 0, 0, weights, callback);
 219
 220         if (ret) {
 221                 GP_PixmapFree(dst);
 222                 return NULL;
 223         }
 224
 225         return dst;
 226 }