?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
filters: Fix typo in VH Convolution.
[gfxprim.git] / libs / filters / GP_Linear.c
blob6f906b835a6866c2a8dc39ae76727342b6f17064
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-2012 Cyril Hrubis <metan@ucw.cz> *
20 * *
21 *****************************************************************************/
22
23 #include <math.h>
24
25 #include "core/GP_Context.h"
26 #include "core/GP_GetPutPixel.h"
27
28 #include "core/GP_Debug.h"
29
30 #include "GP_Linear.h"
31
32 static inline unsigned int gaussian_kernel_size(float sigma)
33 {
34 int center = 3 * sigma;
35
36 return 2 * center + 1;
37 }
38
39 static inline float gaussian_kernel_init(float sigma, float *kernel)
40 {
41 int i, center = 3 * sigma;
42 int N = 2 * center + 1;
43 float ret = 0;
44
45 double sigma2 = sigma * sigma;
46
47 for (i = 0; i < N; i++) {
48 double r = center - i;
49 kernel[i] = exp(-0.5 * (r * r) / sigma2);
50 ret += kernel[i];
51 }
52
53 return ret;
54 }
55
56 static int gaussian_callback_horiz(GP_ProgressCallback *self)
57 {
58 GP_ProgressCallback *callback = self->priv;
59
60 callback->percentage = self->percentage / 2;
61 return callback->callback(callback);
62 }
63
64 static int gaussian_callback_vert(GP_ProgressCallback *self)
65 {
66 GP_ProgressCallback *callback = self->priv;
67
68 callback->percentage = self->percentage / 2 + 50;
69 return callback->callback(callback);
70 }
71
72 int GP_FilterGaussianBlur_Raw(const GP_Context *src, GP_Context *dst,
73 float sigma_x, float sigma_y,
74 GP_ProgressCallback *callback)
75 {
76 unsigned int size_x = gaussian_kernel_size(sigma_x);
77 unsigned int size_y = gaussian_kernel_size(sigma_y);
78
79 GP_DEBUG(1, "Gaussian blur sigma_x=%2.3f sigma_y=%2.3f kernel %ix%i image %ux%u",
80 sigma_x, sigma_y, size_x, size_y, src->w, src->h);
81
82 GP_ProgressCallback *new_callback = NULL;
83
84 GP_ProgressCallback gaussian_callback = {
85 .callback = gaussian_callback_horiz,
86 .priv = callback
87 };
88
89 if (callback != NULL)
90 new_callback = &gaussian_callback;
91
92 /* compute kernel and apply in horizontal direction */
93 if (sigma_x > 0) {
94 float kernel_x[size_x];
95 float sum = gaussian_kernel_init(sigma_x, kernel_x);
96
97 if (GP_FilterHLinearConvolution_Raw(src, dst, kernel_x, size_x,
98 sum, new_callback))
99 return 1;
100 }
101
102 if (new_callback != NULL)
103 new_callback->callback = gaussian_callback_vert;
104
105 /* compute kernel and apply in vertical direction */
106 if (sigma_y > 0) {
107 float kernel_y[size_y];
108 float sum = gaussian_kernel_init(sigma_y, kernel_y);
109
110 if (GP_FilterVLinearConvolution_Raw(dst, dst, kernel_y, size_y,
111 sum, new_callback))
112 return 1;
113 }
114
115 GP_ProgressCallbackDone(callback);
116 return 0;
117 }
118
119 GP_Context *GP_FilterGaussianBlur(const GP_Context *src, GP_Context *dst,
120 float sigma_x, float sigma_y,
121 GP_ProgressCallback *callback)
122 {
123 /* TODO: templatetize */
124 if (src->pixel_type != GP_PIXEL_RGB888)
125 return NULL;
126
127 if (dst == NULL) {
128 dst = GP_ContextCopy(src, 0);
129
130 if (dst == NULL)
131 return NULL;
132 } else {
133 GP_ASSERT(src->pixel_type == dst->pixel_type,
134 "The src and dst pixel types must match");
135 GP_ASSERT(src->w <= dst->w && src->h <= dst->h,
136 "Destination is not big enough");
137 }
138
139 GP_FilterGaussianBlur_Raw(src, dst, sigma_x, sigma_y, callback);
140
141 return dst;
142 }
143
144 #define MUL 1024
145
146 #define CLAMP(val) do { \
147 if (val > 255) \
148 val = 255; \
149 if (val < 0) \
150 val = 0; \
151 } while (0)
152
153 int GP_FilterHLinearConvolution_Raw(const GP_Context *src, GP_Context *dst,
154 float kernel[], uint32_t kw, float kern_div,
155 GP_ProgressCallback *callback)
156 {
157 GP_Coord x, y;
158 uint32_t i;
159 int32_t ikernel[kw], ikern_div;
160 uint32_t size = dst->w + kw - 1;
161
162 GP_DEBUG(1, "Horizontal linear convolution kernel width %i image %ux%u",
163 kw, src->w, src->h);
164
165 for (i = 0; i < kw; i++)
166 ikernel[i] = kernel[i] * MUL + 0.5;
167
168 ikern_div = kern_div * MUL + 0.5;
169
170 /* do linear convolution */
171 for (y = 0; y < (GP_Coord)dst->h; y++) {
172 uint8_t R[size], G[size], B[size];
173
174 /* Fetch the whole row */
175 GP_Pixel pix = GP_GetPixel_Raw_24BPP(src, 0, y);
176
177 for (i = 0; i < kw/2; i++) {
178 R[i] = GP_Pixel_GET_R_RGB888(pix);
179 G[i] = GP_Pixel_GET_G_RGB888(pix);
180 B[i] = GP_Pixel_GET_B_RGB888(pix);
181 }
182
183 for (i = 0; i < src->w; i++) {
184 pix = GP_GetPixel_Raw_24BPP(src, i, y);
185
186 R[i+kw/2] = GP_Pixel_GET_R_RGB888(pix);
187 G[i+kw/2] = GP_Pixel_GET_G_RGB888(pix);
188 B[i+kw/2] = GP_Pixel_GET_B_RGB888(pix);
189 }
190
191 for (i = src->w + kw/2; i < size; i++) {
192 R[i] = GP_Pixel_GET_R_RGB888(pix);
193 G[i] = GP_Pixel_GET_G_RGB888(pix);
194 B[i] = GP_Pixel_GET_B_RGB888(pix);
195 }
196
197 for (x = 0; x < (GP_Coord)dst->w; x++) {
198 int32_t r = MUL/2, g = MUL/2, b = MUL/2;
199
200 /* count the pixel value from neighbours weighted by kernel */
201 for (i = 0; i < kw; i++) {
202 r += R[i + x] * ikernel[i];
203 g += G[i + x] * ikernel[i];
204 b += B[i + x] * ikernel[i];
205 }
206
207 /* divide the result */
208 r /= ikern_div;
209 g /= ikern_div;
210 b /= ikern_div;
211
212 /* and clamp just to be extra sure */
213 CLAMP(r);
214 CLAMP(g);
215 CLAMP(b);
216
217 GP_PutPixel_Raw_24BPP(dst, x, y,
218 GP_Pixel_CREATE_RGB888(r, g, b));
219 }
220
221 if (GP_ProgressCallbackReport(callback, y, dst->h, dst->w))
222 return 1;
223 }
224
225 GP_ProgressCallbackDone(callback);
226 return 0;
227 }
228
229 int GP_FilterVLinearConvolution_Raw(const GP_Context *src, GP_Context *dst,
230 float kernel[], uint32_t kh, float kern_div,
231 GP_ProgressCallback *callback)
232 {
233 GP_Coord x, y;
234 uint32_t i;
235 int32_t ikernel[kh], ikern_div;
236 uint32_t size = dst->h + kh - 1;
237
238 for (i = 0; i < kh; i++)
239 ikernel[i] = kernel[i] * MUL + 0.5;
240
241 GP_DEBUG(1, "Vertical linear convolution kernel width %i image %ux%u",
242 kh, src->w, src->h);
243
244 ikern_div = kern_div * MUL + 0.5;
245
246 /* do linear convolution */
247 for (x = 0; x < (GP_Coord)dst->w; x++) {
248 uint8_t R[size], G[size], B[size];
249
250 /* Fetch the whole column */
251 GP_Pixel pix = GP_GetPixel_Raw_24BPP(src, x, 0);
252
253 for (i = 0; i < kh/2; i++) {
254 R[i] = GP_Pixel_GET_R_RGB888(pix);
255 G[i] = GP_Pixel_GET_G_RGB888(pix);
256 B[i] = GP_Pixel_GET_B_RGB888(pix);
257 }
258
259 for (i = 0; i < src->h; i++) {
260 pix = GP_GetPixel_Raw_24BPP(src, x, i);
261
262 R[i+kh/2] = GP_Pixel_GET_R_RGB888(pix);
263 G[i+kh/2] = GP_Pixel_GET_G_RGB888(pix);
264 B[i+kh/2] = GP_Pixel_GET_B_RGB888(pix);
265 }
266
267 for (i = src->h + kh/2; i < size; i++) {
268 R[i] = GP_Pixel_GET_R_RGB888(pix);
269 G[i] = GP_Pixel_GET_G_RGB888(pix);
270 B[i] = GP_Pixel_GET_B_RGB888(pix);
271 }
272
273 for (y = 0; y < (GP_Coord)dst->h; y++) {
274 int32_t r = MUL/2, g = MUL/2, b = MUL/2;
275
276 /* count the pixel value from neighbours weighted by kernel */
277 for (i = 0; i < kh; i++) {
278 r += R[y + i] * ikernel[i];
279 g += G[y + i] * ikernel[i];
280 b += B[y + i] * ikernel[i];
281 }
282
283 /* divide the result */
284 r /= ikern_div;
285 g /= ikern_div;
286 b /= ikern_div;
287
288 /* and clamp just to be extra sure */
289 CLAMP(r);
290 CLAMP(g);
291 CLAMP(b);
292
293 GP_PutPixel_Raw_24BPP(dst, x, y,
294 GP_Pixel_CREATE_RGB888(r, g, b));
295 }
296
297 if (GP_ProgressCallbackReport(callback, x, dst->w, dst->h))
298 return 1;
299 }
300
301 GP_ProgressCallbackDone(callback);
302 return 0;
303 }
304
305 static int h_callback(GP_ProgressCallback *self)
306 {
307 GP_ProgressCallback *callback = self->priv;
308
309 callback->percentage = self->percentage / 2;
310 return callback->callback(callback);
311 }
312
313 static int v_callback(GP_ProgressCallback *self)
314 {
315 GP_ProgressCallback *callback = self->priv;
316
317 callback->percentage = self->percentage / 2 + 50;
318 return callback->callback(callback);
319 }
320
321 int GP_FilterVHLinearConvolution_Raw(const GP_Context *src, GP_Context *dst,
322 float hkernel[], uint32_t kw, float hkern_div,
323 float vkernel[], uint32_t kh, float vkern_div,
324 GP_ProgressCallback *callback)
325 {
326 GP_ProgressCallback *new_callback;
327
328 GP_ProgressCallback conv_callback = {
329 .callback = h_callback,
330 .priv = callback,
331 };
332
333 new_callback = callback ? &conv_callback : NULL;
334
335 if (GP_FilterVLinearConvolution_Raw(src, dst, hkernel, kw, hkern_div, new_callback))
336 return 1;
337
338 conv_callback.callback = v_callback;
339
340 if (GP_FilterHLinearConvolution_Raw(dst, dst, vkernel, kh, vkern_div, new_callback))
341 return 1;
342
343 GP_ProgressCallbackDone(callback);
344 return 0;
345 }
346
347 /*
348 * Linear convolution.
349 *
350 * Can be used in-place.
351 */
352 int GP_FilterLinearConvolution_Raw(const GP_Context *src, GP_Context *dst,
353 float kernel[], uint32_t kw, uint32_t kh,
354 float kern_div, GP_ProgressCallback *callback)
355 {
356 GP_Coord x, y;
357 uint32_t i, j;
358
359 GP_DEBUG(1, "Linear convolution kernel %ix%i image %ux%u",
360 kw, kh, src->w, src->h);
361
362 /* do linear convolution */
363 for (y = 0; y < (GP_Coord)dst->h; y++) {
364 GP_Pixel pix;
365 uint32_t R[kw][kh], G[kw][kh], B[kw][kh];
366
367 /* prefill the buffer on the start */
368 for (j = 0; j < kh; j++) {
369 for (i = 0; i < kw - 1; i++) {
370 int cx = i - kw/2;
371 int cy = y + j - kh/2;
372
373 if (cx < 0)
374 cx = 0;
375
376 if (cy < 0)
377 cy = 0;
378
379 if (cy >= (int)src->h)
380 cy = src->h - 1;
381
382 pix = GP_GetPixel_Raw_24BPP(src, cx, cy);
383
384 R[i][j] = GP_Pixel_GET_R_RGB888(pix);
385 G[i][j] = GP_Pixel_GET_G_RGB888(pix);
386 B[i][j] = GP_Pixel_GET_B_RGB888(pix);
387 }
388 }
389
390 int idx = kw - 1;
391
392 for (x = 0; x < (GP_Coord)dst->w; x++) {
393 float r = 0, g = 0, b = 0;
394
395 for (j = 0; j < kh; j++) {
396 int cy = y + j - kh/2;
397 int cx = x + kw/2;
398
399 if (cy < 0)
400 cy = 0;
401
402 if (cy >= (int)src->h)
403 cy = src->h - 1;
404
405 if (cx >= (int)src->w)
406 cx = src->w - 1;
407
408 pix = GP_GetPixel_Raw_24BPP(src, cx, cy);
409
410 R[idx][j] = GP_Pixel_GET_R_RGB888(pix);
411 G[idx][j] = GP_Pixel_GET_G_RGB888(pix);
412 B[idx][j] = GP_Pixel_GET_B_RGB888(pix);
413 }
414
415 /* count the pixel value from neighbours weighted by kernel */
416 for (i = 0; i < kw; i++) {
417 int k;
418
419 if ((int)i < idx + 1)
420 k = kw - idx - 1 + i;
421 else
422 k = i - idx - 1;
423
424 for (j = 0; j < kh; j++) {
425 r += R[i][j] * kernel[k + j * kw];
426 g += G[i][j] * kernel[k + j * kw];
427 b += B[i][j] * kernel[k + j * kw];
428 }
429 }
430
431 /* divide the result */
432 r /= kern_div;
433 g /= kern_div;
434 b /= kern_div;
435
436 /* and clamp just to be extra sure */
437 if (r > 255)
438 r = 255;
439 if (r < 0)
440 r = 0;
441 if (g > 255)
442 g = 255;
443 if (g < 0)
444 g = 0;
445 if (b > 255)
446 b = 255;
447 if (b < 0)
448 b = 0;
449
450 pix = GP_Pixel_CREATE_RGB888((uint32_t)r, (uint32_t)g, (uint32_t)b);
451
452 GP_PutPixel_Raw_24BPP(dst, x, y, pix);
453
454 idx++;
455
456 if (idx >= (int)kw)
457 idx = 0;
458 }
459
460 if (GP_ProgressCallbackReport(callback, y, dst->h, dst->w))
461 return 1;
462 }
463
464 GP_ProgressCallbackDone(callback);
465 return 0;
466 }
A simple 2D graphics library with emphasis on correctness and well-defined operation.