/src/leptonica/src/convolve.c
Line | Count | Source (jump to first uncovered line) |
1 | | /*====================================================================* |
2 | | - Copyright (C) 2001 Leptonica. All rights reserved. |
3 | | - |
4 | | - Redistribution and use in source and binary forms, with or without |
5 | | - modification, are permitted provided that the following conditions |
6 | | - are met: |
7 | | - 1. Redistributions of source code must retain the above copyright |
8 | | - notice, this list of conditions and the following disclaimer. |
9 | | - 2. Redistributions in binary form must reproduce the above |
10 | | - copyright notice, this list of conditions and the following |
11 | | - disclaimer in the documentation and/or other materials |
12 | | - provided with the distribution. |
13 | | - |
14 | | - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
15 | | - ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
16 | | - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
17 | | - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ANY |
18 | | - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
19 | | - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
20 | | - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
21 | | - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY |
22 | | - OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
23 | | - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
24 | | - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
25 | | *====================================================================*/ |
26 | | |
27 | | /*! |
28 | | * \file convolve.c |
29 | | * <pre> |
30 | | * |
31 | | * Top level grayscale or color block convolution |
32 | | * PIX *pixBlockconv() |
33 | | * |
34 | | * Grayscale block convolution |
35 | | * PIX *pixBlockconvGray() |
36 | | * static void blockconvLow() |
37 | | * |
38 | | * Accumulator for 1, 8 and 32 bpp convolution |
39 | | * PIX *pixBlockconvAccum() |
40 | | * static void blockconvAccumLow() |
41 | | * |
42 | | * Un-normalized grayscale block convolution |
43 | | * PIX *pixBlockconvGrayUnnormalized() |
44 | | * |
45 | | * Tiled grayscale or color block convolution |
46 | | * PIX *pixBlockconvTiled() |
47 | | * PIX *pixBlockconvGrayTile() |
48 | | * |
49 | | * Convolution for mean, mean square, variance and rms deviation |
50 | | * in specified window |
51 | | * l_int32 pixWindowedStats() |
52 | | * PIX *pixWindowedMean() |
53 | | * PIX *pixWindowedMeanSquare() |
54 | | * l_int32 pixWindowedVariance() |
55 | | * DPIX *pixMeanSquareAccum() |
56 | | * |
57 | | * Binary block sum and rank filter |
58 | | * PIX *pixBlockrank() |
59 | | * PIX *pixBlocksum() |
60 | | * static void blocksumLow() |
61 | | * |
62 | | * Census transform |
63 | | * PIX *pixCensusTransform() |
64 | | * |
65 | | * Generic convolution (with Pix) |
66 | | * PIX *pixConvolve() |
67 | | * PIX *pixConvolveSep() |
68 | | * PIX *pixConvolveRGB() |
69 | | * PIX *pixConvolveRGBSep() |
70 | | * |
71 | | * Generic convolution (with float arrays) |
72 | | * FPIX *fpixConvolve() |
73 | | * FPIX *fpixConvolveSep() |
74 | | * |
75 | | * Convolution with bias (for non-negative output) |
76 | | * PIX *pixConvolveWithBias() |
77 | | * |
78 | | * Set parameter for convolution subsampling |
79 | | * void l_setConvolveSampling() |
80 | | * |
81 | | * Additive gaussian noise |
82 | | * PIX *pixAddGaussNoise() |
83 | | * l_float32 gaussDistribSampling() |
84 | | * </pre> |
85 | | */ |
86 | | |
87 | | #ifdef HAVE_CONFIG_H |
88 | | #include <config_auto.h> |
89 | | #endif /* HAVE_CONFIG_H */ |
90 | | |
91 | | #include <math.h> |
92 | | #include "allheaders.h" |
93 | | |
94 | | /* These globals determine the subsampling factors for |
95 | | * generic convolution of pix and fpix. Declare extern to use. |
96 | | * To change the values, use l_setConvolveSampling(). */ |
97 | | LEPT_DLL l_int32 ConvolveSamplingFactX = 1; |
98 | | LEPT_DLL l_int32 ConvolveSamplingFactY = 1; |
99 | | |
100 | | /* Low-level static functions */ |
101 | | static void blockconvLow(l_uint32 *data, l_int32 w, l_int32 h, l_int32 wpl, |
102 | | l_uint32 *dataa, l_int32 wpla, l_int32 wc, |
103 | | l_int32 hc); |
104 | | static void blockconvAccumLow(l_uint32 *datad, l_int32 w, l_int32 h, |
105 | | l_int32 wpld, l_uint32 *datas, l_int32 d, |
106 | | l_int32 wpls); |
107 | | static void blocksumLow(l_uint32 *datad, l_int32 w, l_int32 h, l_int32 wpl, |
108 | | l_uint32 *dataa, l_int32 wpla, l_int32 wc, l_int32 hc); |
109 | | |
110 | | |
111 | | /*----------------------------------------------------------------------* |
112 | | * Top-level grayscale or color block convolution * |
113 | | *----------------------------------------------------------------------*/ |
114 | | /*! |
115 | | * \brief pixBlockconv() |
116 | | * |
117 | | * \param[in] pix 8 or 32 bpp; or 2, 4 or 8 bpp with colormap |
118 | | * \param[in] wc, hc half width/height of convolution kernel |
119 | | * \return pixd, or NULL on error |
120 | | * |
121 | | * <pre> |
122 | | * Notes: |
123 | | * (1) The full width and height of the convolution kernel |
124 | | * are (2 * wc + 1) and (2 * hc + 1) |
125 | | * (2) Returns a copy if either wc or hc are 0 |
126 | | * (3) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1, |
127 | | * where (w,h) are the dimensions of pixs. Attempt to |
128 | | * reduce the kernel size if necessary. |
129 | | * </pre> |
130 | | */ |
131 | | PIX * |
132 | | pixBlockconv(PIX *pix, |
133 | | l_int32 wc, |
134 | | l_int32 hc) |
135 | 0 | { |
136 | 0 | l_int32 w, h, d; |
137 | 0 | PIX *pixs, *pixd, *pixr, *pixrc, *pixg, *pixgc, *pixb, *pixbc; |
138 | |
|
139 | 0 | if (!pix) |
140 | 0 | return (PIX *)ERROR_PTR("pix not defined", __func__, NULL); |
141 | 0 | if (wc <= 0 || hc <= 0) |
142 | 0 | return pixCopy(NULL, pix); |
143 | 0 | pixGetDimensions(pix, &w, &h, &d); |
144 | 0 | if (w < 2 * wc + 1 || h < 2 * hc + 1) { |
145 | 0 | L_WARNING("kernel too large: wc = %d, hc = %d, w = %d, h = %d; " |
146 | 0 | "reducing!\n", __func__, wc, hc, w, h); |
147 | 0 | wc = L_MIN(wc, (w - 1) / 2); |
148 | 0 | hc = L_MIN(hc, (h - 1) / 2); |
149 | 0 | } |
150 | 0 | if (wc == 0 || hc == 0) /* no-op */ |
151 | 0 | return pixCopy(NULL, pix); |
152 | | |
153 | | /* Remove colormap if necessary */ |
154 | 0 | if ((d == 2 || d == 4 || d == 8) && pixGetColormap(pix)) { |
155 | 0 | L_WARNING("pix has colormap; removing\n", __func__); |
156 | 0 | pixs = pixRemoveColormap(pix, REMOVE_CMAP_BASED_ON_SRC); |
157 | 0 | d = pixGetDepth(pixs); |
158 | 0 | } else { |
159 | 0 | pixs = pixClone(pix); |
160 | 0 | } |
161 | |
|
162 | 0 | if (d != 8 && d != 32) { |
163 | 0 | pixDestroy(&pixs); |
164 | 0 | return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, NULL); |
165 | 0 | } |
166 | | |
167 | 0 | if (d == 8) { |
168 | 0 | pixd = pixBlockconvGray(pixs, NULL, wc, hc); |
169 | 0 | } else { /* d == 32 */ |
170 | 0 | pixr = pixGetRGBComponent(pixs, COLOR_RED); |
171 | 0 | pixrc = pixBlockconvGray(pixr, NULL, wc, hc); |
172 | 0 | pixDestroy(&pixr); |
173 | 0 | pixg = pixGetRGBComponent(pixs, COLOR_GREEN); |
174 | 0 | pixgc = pixBlockconvGray(pixg, NULL, wc, hc); |
175 | 0 | pixDestroy(&pixg); |
176 | 0 | pixb = pixGetRGBComponent(pixs, COLOR_BLUE); |
177 | 0 | pixbc = pixBlockconvGray(pixb, NULL, wc, hc); |
178 | 0 | pixDestroy(&pixb); |
179 | 0 | pixd = pixCreateRGBImage(pixrc, pixgc, pixbc); |
180 | 0 | pixDestroy(&pixrc); |
181 | 0 | pixDestroy(&pixgc); |
182 | 0 | pixDestroy(&pixbc); |
183 | 0 | } |
184 | |
|
185 | 0 | pixDestroy(&pixs); |
186 | 0 | return pixd; |
187 | 0 | } |
188 | | |
189 | | |
190 | | /*----------------------------------------------------------------------* |
191 | | * Grayscale block convolution * |
192 | | *----------------------------------------------------------------------*/ |
193 | | /*! |
194 | | * \brief pixBlockconvGray() |
195 | | * |
196 | | * \param[in] pixs 8 bpp |
197 | | * \param[in] pixacc pix 32 bpp; can be null |
198 | | * \param[in] wc, hc half width/height of convolution kernel |
199 | | * \return pix 8 bpp, or NULL on error |
200 | | * |
201 | | * <pre> |
202 | | * Notes: |
203 | | * (1) If accum pix is null, make one and destroy it before |
204 | | * returning; otherwise, just use the input accum pix. |
205 | | * (2) The full width and height of the convolution kernel |
206 | | * are (2 * wc + 1) and (2 * hc + 1). |
207 | | * (3) Returns a copy if either wc or hc are 0 |
208 | | * (4) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1, |
209 | | * where (w,h) are the dimensions of pixs. Attempt to |
210 | | * reduce the kernel size if necessary. |
211 | | * </pre> |
212 | | */ |
213 | | PIX * |
214 | | pixBlockconvGray(PIX *pixs, |
215 | | PIX *pixacc, |
216 | | l_int32 wc, |
217 | | l_int32 hc) |
218 | 0 | { |
219 | 0 | l_int32 w, h, d, wpl, wpla; |
220 | 0 | l_uint32 *datad, *dataa; |
221 | 0 | PIX *pixd, *pixt; |
222 | |
|
223 | 0 | if (!pixs) |
224 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
225 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
226 | 0 | if (d != 8) |
227 | 0 | return (PIX *)ERROR_PTR("pixs not 8 bpp", __func__, NULL); |
228 | 0 | if (wc <= 0 || hc <= 0) /* no-op */ |
229 | 0 | return pixCopy(NULL, pixs); |
230 | 0 | if (w < 2 * wc + 1 || h < 2 * hc + 1) { |
231 | 0 | L_WARNING("kernel too large: wc = %d, hc = %d, w = %d, h = %d; " |
232 | 0 | "reducing!\n", __func__, wc, hc, w, h); |
233 | 0 | wc = L_MIN(wc, (w - 1) / 2); |
234 | 0 | hc = L_MIN(hc, (h - 1) / 2); |
235 | 0 | } |
236 | 0 | if (wc == 0 || hc == 0) |
237 | 0 | return pixCopy(NULL, pixs); |
238 | | |
239 | 0 | if (pixacc) { |
240 | 0 | if (pixGetDepth(pixacc) == 32) { |
241 | 0 | pixt = pixClone(pixacc); |
242 | 0 | } else { |
243 | 0 | L_WARNING("pixacc not 32 bpp; making new one\n", __func__); |
244 | 0 | if ((pixt = pixBlockconvAccum(pixs)) == NULL) |
245 | 0 | return (PIX *)ERROR_PTR("pixt not made", __func__, NULL); |
246 | 0 | } |
247 | 0 | } else { |
248 | 0 | if ((pixt = pixBlockconvAccum(pixs)) == NULL) |
249 | 0 | return (PIX *)ERROR_PTR("pixt not made", __func__, NULL); |
250 | 0 | } |
251 | | |
252 | 0 | if ((pixd = pixCreateTemplate(pixs)) == NULL) { |
253 | 0 | pixDestroy(&pixt); |
254 | 0 | return (PIX *)ERROR_PTR("pixd not made", __func__, NULL); |
255 | 0 | } |
256 | | |
257 | 0 | pixSetPadBits(pixt, 0); |
258 | 0 | wpl = pixGetWpl(pixd); |
259 | 0 | wpla = pixGetWpl(pixt); |
260 | 0 | datad = pixGetData(pixd); |
261 | 0 | dataa = pixGetData(pixt); |
262 | 0 | blockconvLow(datad, w, h, wpl, dataa, wpla, wc, hc); |
263 | |
|
264 | 0 | pixDestroy(&pixt); |
265 | 0 | return pixd; |
266 | 0 | } |
267 | | |
268 | | |
269 | | /*! |
270 | | * \brief blockconvLow() |
271 | | * |
272 | | * \param[in] data data of input image, to be convolved |
273 | | * \param[in] w, h, wpl |
274 | | * \param[in] dataa data of 32 bpp accumulator |
275 | | * \param[in] wpla accumulator |
276 | | * \param[in] wc convolution "half-width" |
277 | | * \param[in] hc convolution "half-height" |
278 | | * \return void |
279 | | * |
280 | | * <pre> |
281 | | * Notes: |
282 | | * (1) The full width and height of the convolution kernel |
283 | | * are (2 * wc + 1) and (2 * hc + 1). |
284 | | * (2) The lack of symmetry between the handling of the |
285 | | * first (hc + 1) lines and the last (hc) lines, |
286 | | * and similarly with the columns, is due to fact that |
287 | | * for the pixel at (x,y), the accumulator values are |
288 | | * taken at (x + wc, y + hc), (x - wc - 1, y + hc), |
289 | | * (x + wc, y - hc - 1) and (x - wc - 1, y - hc - 1). |
290 | | * (3) We compute sums, normalized as if there were no reduced |
291 | | * area at the boundary. This under-estimates the value |
292 | | * of the boundary pixels, so we multiply them by another |
293 | | * normalization factor that is greater than 1. |
294 | | * (4) This second normalization is done first for the first |
295 | | * hc + 1 lines; then for the last hc lines; and finally |
296 | | * for the first wc + 1 and last wc columns in the intermediate |
297 | | * lines. |
298 | | * (5) The caller should verify that wc < w and hc < h. |
299 | | * Failing either condition, illegal reads and writes can occur. |
300 | | * (6) Implementation note: to get the same results in the interior |
301 | | * between this function and pixConvolve(), it is necessary to |
302 | | * add 0.5 for roundoff in the main loop that runs over all pixels. |
303 | | * However, if we do that and have white (255) pixels near the |
304 | | * image boundary, some overflow occurs for pixels very close |
305 | | * to the boundary. We can't fix this by subtracting from the |
306 | | * normalized values for the boundary pixels, because this results |
307 | | * in underflow if the boundary pixels are black (0). Empirically, |
308 | | * adding 0.25 (instead of 0.5) before truncating in the main |
309 | | * loop will not cause overflow, but this gives some |
310 | | * off-by-1-level errors in interior pixel values. So we add |
311 | | * 0.5 for roundoff in the main loop, and for pixels within a |
312 | | * half filter width of the boundary, use a L_MIN of the |
313 | | * computed value and 255 to avoid overflow during normalization. |
314 | | * </pre> |
315 | | */ |
316 | | static void |
317 | | blockconvLow(l_uint32 *data, |
318 | | l_int32 w, |
319 | | l_int32 h, |
320 | | l_int32 wpl, |
321 | | l_uint32 *dataa, |
322 | | l_int32 wpla, |
323 | | l_int32 wc, |
324 | | l_int32 hc) |
325 | 0 | { |
326 | 0 | l_int32 i, j, imax, imin, jmax, jmin; |
327 | 0 | l_int32 wn, hn, fwc, fhc, wmwc, hmhc; |
328 | 0 | l_float32 norm, normh, normw; |
329 | 0 | l_uint32 val; |
330 | 0 | l_uint32 *linemina, *linemaxa, *line; |
331 | |
|
332 | 0 | wmwc = w - wc; |
333 | 0 | hmhc = h - hc; |
334 | 0 | if (wmwc <= 0 || hmhc <= 0) { |
335 | 0 | L_ERROR("wc >= w || hc >= h\n", __func__); |
336 | 0 | return; |
337 | 0 | } |
338 | 0 | fwc = 2 * wc + 1; |
339 | 0 | fhc = 2 * hc + 1; |
340 | 0 | norm = 1.0 / ((l_float32)(fwc) * fhc); |
341 | | |
342 | | /*------------------------------------------------------------* |
343 | | * Compute, using b.c. only to set limits on the accum image * |
344 | | *------------------------------------------------------------*/ |
345 | 0 | for (i = 0; i < h; i++) { |
346 | 0 | imin = L_MAX(i - 1 - hc, 0); |
347 | 0 | imax = L_MIN(i + hc, h - 1); |
348 | 0 | line = data + wpl * i; |
349 | 0 | linemina = dataa + wpla * imin; |
350 | 0 | linemaxa = dataa + wpla * imax; |
351 | 0 | for (j = 0; j < w; j++) { |
352 | 0 | jmin = L_MAX(j - 1 - wc, 0); |
353 | 0 | jmax = L_MIN(j + wc, w - 1); |
354 | 0 | val = linemaxa[jmax] - linemaxa[jmin] |
355 | 0 | + linemina[jmin] - linemina[jmax]; |
356 | 0 | val = (l_uint8)(norm * val + 0.5); /* see comment above */ |
357 | 0 | SET_DATA_BYTE(line, j, val); |
358 | 0 | } |
359 | 0 | } |
360 | | |
361 | | /*------------------------------------------------------------* |
362 | | * Fix normalization for boundary pixels * |
363 | | *------------------------------------------------------------*/ |
364 | 0 | for (i = 0; i <= hc; i++) { /* first hc + 1 lines */ |
365 | 0 | hn = L_MAX(1, hc + i); |
366 | 0 | normh = (l_float32)fhc / (l_float32)hn; /* >= 1 */ |
367 | 0 | line = data + wpl * i; |
368 | 0 | for (j = 0; j <= wc; j++) { |
369 | 0 | wn = L_MAX(1, wc + j); |
370 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* >= 1 */ |
371 | 0 | val = GET_DATA_BYTE(line, j); |
372 | 0 | val = (l_uint8)L_MIN(val * normh * normw, 255); |
373 | 0 | SET_DATA_BYTE(line, j, val); |
374 | 0 | } |
375 | 0 | for (j = wc + 1; j < wmwc; j++) { |
376 | 0 | val = GET_DATA_BYTE(line, j); |
377 | 0 | val = (l_uint8)L_MIN(val * normh, 255); |
378 | 0 | SET_DATA_BYTE(line, j, val); |
379 | 0 | } |
380 | 0 | for (j = wmwc; j < w; j++) { |
381 | 0 | wn = wc + w - j; |
382 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
383 | 0 | val = GET_DATA_BYTE(line, j); |
384 | 0 | val = (l_uint8)L_MIN(val * normh * normw, 255); |
385 | 0 | SET_DATA_BYTE(line, j, val); |
386 | 0 | } |
387 | 0 | } |
388 | |
|
389 | 0 | for (i = hmhc; i < h; i++) { /* last hc lines */ |
390 | 0 | hn = hc + h - i; |
391 | 0 | normh = (l_float32)fhc / (l_float32)hn; /* > 1 */ |
392 | 0 | line = data + wpl * i; |
393 | 0 | for (j = 0; j <= wc; j++) { |
394 | 0 | wn = wc + j; |
395 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
396 | 0 | val = GET_DATA_BYTE(line, j); |
397 | 0 | val = (l_uint8)L_MIN(val * normh * normw, 255); |
398 | 0 | SET_DATA_BYTE(line, j, val); |
399 | 0 | } |
400 | 0 | for (j = wc + 1; j < wmwc; j++) { |
401 | 0 | val = GET_DATA_BYTE(line, j); |
402 | 0 | val = (l_uint8)L_MIN(val * normh, 255); |
403 | 0 | SET_DATA_BYTE(line, j, val); |
404 | 0 | } |
405 | 0 | for (j = wmwc; j < w; j++) { |
406 | 0 | wn = wc + w - j; |
407 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
408 | 0 | val = GET_DATA_BYTE(line, j); |
409 | 0 | val = (l_uint8)L_MIN(val * normh * normw, 255); |
410 | 0 | SET_DATA_BYTE(line, j, val); |
411 | 0 | } |
412 | 0 | } |
413 | |
|
414 | 0 | for (i = hc + 1; i < hmhc; i++) { /* intermediate lines */ |
415 | 0 | line = data + wpl * i; |
416 | 0 | for (j = 0; j <= wc; j++) { /* first wc + 1 columns */ |
417 | 0 | wn = wc + j; |
418 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
419 | 0 | val = GET_DATA_BYTE(line, j); |
420 | 0 | val = (l_uint8)L_MIN(val * normw, 255); |
421 | 0 | SET_DATA_BYTE(line, j, val); |
422 | 0 | } |
423 | 0 | for (j = wmwc; j < w; j++) { /* last wc columns */ |
424 | 0 | wn = wc + w - j; |
425 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
426 | 0 | val = GET_DATA_BYTE(line, j); |
427 | 0 | val = (l_uint8)L_MIN(val * normw, 255); |
428 | 0 | SET_DATA_BYTE(line, j, val); |
429 | 0 | } |
430 | 0 | } |
431 | 0 | } |
432 | | |
433 | | |
434 | | /*----------------------------------------------------------------------* |
435 | | * Accumulator for 1, 8 and 32 bpp convolution * |
436 | | *----------------------------------------------------------------------*/ |
437 | | /*! |
438 | | * \brief pixBlockconvAccum() |
439 | | * |
440 | | * \param[in] pixs 1, 8 or 32 bpp |
441 | | * \return accum pix 32 bpp, or NULL on error. |
442 | | * |
443 | | * <pre> |
444 | | * Notes: |
445 | | * (1) The general recursion relation is |
446 | | * a(i,j) = v(i,j) + a(i-1, j) + a(i, j-1) - a(i-1, j-1) |
447 | | * For the first line, this reduces to the special case |
448 | | * a(i,j) = v(i,j) + a(i, j-1) |
449 | | * For the first column, the special case is |
450 | | * a(i,j) = v(i,j) + a(i-1, j) |
451 | | * </pre> |
452 | | */ |
453 | | PIX * |
454 | | pixBlockconvAccum(PIX *pixs) |
455 | 0 | { |
456 | 0 | l_int32 w, h, d, wpls, wpld; |
457 | 0 | l_uint32 *datas, *datad; |
458 | 0 | PIX *pixd; |
459 | |
|
460 | 0 | if (!pixs) |
461 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
462 | | |
463 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
464 | 0 | if (d != 1 && d != 8 && d != 32) |
465 | 0 | return (PIX *)ERROR_PTR("pixs not 1, 8 or 32 bpp", __func__, NULL); |
466 | 0 | if ((pixd = pixCreate(w, h, 32)) == NULL) |
467 | 0 | return (PIX *)ERROR_PTR("pixd not made", __func__, NULL); |
468 | | |
469 | 0 | datas = pixGetData(pixs); |
470 | 0 | datad = pixGetData(pixd); |
471 | 0 | wpls = pixGetWpl(pixs); |
472 | 0 | wpld = pixGetWpl(pixd); |
473 | 0 | blockconvAccumLow(datad, w, h, wpld, datas, d, wpls); |
474 | |
|
475 | 0 | return pixd; |
476 | 0 | } |
477 | | |
478 | | |
479 | | /* |
480 | | * \brief blockconvAccumLow() |
481 | | * |
482 | | * \param[in] datad 32 bpp dest |
483 | | * \param[in] w, h, wpld of 32 bpp dest |
484 | | * \param[in] datas 1, 8 or 32 bpp src |
485 | | * \param[in] d bpp of src |
486 | | * \param[in] wpls of src |
487 | | * \return void |
488 | | * |
489 | | * <pre> |
490 | | * Notes: |
491 | | * (1) The general recursion relation is |
492 | | * a(i,j) = v(i,j) + a(i-1, j) + a(i, j-1) - a(i-1, j-1) |
493 | | * For the first line, this reduces to the special case |
494 | | * a(0,j) = v(0,j) + a(0, j-1), j > 0 |
495 | | * For the first column, the special case is |
496 | | * a(i,0) = v(i,0) + a(i-1, 0), i > 0 |
497 | | * </pre> |
498 | | */ |
499 | | static void |
500 | | blockconvAccumLow(l_uint32 *datad, |
501 | | l_int32 w, |
502 | | l_int32 h, |
503 | | l_int32 wpld, |
504 | | l_uint32 *datas, |
505 | | l_int32 d, |
506 | | l_int32 wpls) |
507 | 0 | { |
508 | 0 | l_uint8 val; |
509 | 0 | l_int32 i, j; |
510 | 0 | l_uint32 val32; |
511 | 0 | l_uint32 *lines, *lined, *linedp; |
512 | |
|
513 | 0 | lines = datas; |
514 | 0 | lined = datad; |
515 | |
|
516 | 0 | if (d == 1) { |
517 | | /* Do the first line */ |
518 | 0 | for (j = 0; j < w; j++) { |
519 | 0 | val = GET_DATA_BIT(lines, j); |
520 | 0 | if (j == 0) |
521 | 0 | lined[0] = val; |
522 | 0 | else |
523 | 0 | lined[j] = lined[j - 1] + val; |
524 | 0 | } |
525 | | |
526 | | /* Do the other lines */ |
527 | 0 | for (i = 1; i < h; i++) { |
528 | 0 | lines = datas + i * wpls; |
529 | 0 | lined = datad + i * wpld; /* curr dest line */ |
530 | 0 | linedp = lined - wpld; /* prev dest line */ |
531 | 0 | for (j = 0; j < w; j++) { |
532 | 0 | val = GET_DATA_BIT(lines, j); |
533 | 0 | if (j == 0) |
534 | 0 | lined[0] = val + linedp[0]; |
535 | 0 | else |
536 | 0 | lined[j] = val + lined[j - 1] + linedp[j] - linedp[j - 1]; |
537 | 0 | } |
538 | 0 | } |
539 | 0 | } else if (d == 8) { |
540 | | /* Do the first line */ |
541 | 0 | for (j = 0; j < w; j++) { |
542 | 0 | val = GET_DATA_BYTE(lines, j); |
543 | 0 | if (j == 0) |
544 | 0 | lined[0] = val; |
545 | 0 | else |
546 | 0 | lined[j] = lined[j - 1] + val; |
547 | 0 | } |
548 | | |
549 | | /* Do the other lines */ |
550 | 0 | for (i = 1; i < h; i++) { |
551 | 0 | lines = datas + i * wpls; |
552 | 0 | lined = datad + i * wpld; /* curr dest line */ |
553 | 0 | linedp = lined - wpld; /* prev dest line */ |
554 | 0 | for (j = 0; j < w; j++) { |
555 | 0 | val = GET_DATA_BYTE(lines, j); |
556 | 0 | if (j == 0) |
557 | 0 | lined[0] = val + linedp[0]; |
558 | 0 | else |
559 | 0 | lined[j] = val + lined[j - 1] + linedp[j] - linedp[j - 1]; |
560 | 0 | } |
561 | 0 | } |
562 | 0 | } else if (d == 32) { |
563 | | /* Do the first line */ |
564 | 0 | for (j = 0; j < w; j++) { |
565 | 0 | val32 = lines[j]; |
566 | 0 | if (j == 0) |
567 | 0 | lined[0] = val32; |
568 | 0 | else |
569 | 0 | lined[j] = lined[j - 1] + val32; |
570 | 0 | } |
571 | | |
572 | | /* Do the other lines */ |
573 | 0 | for (i = 1; i < h; i++) { |
574 | 0 | lines = datas + i * wpls; |
575 | 0 | lined = datad + i * wpld; /* curr dest line */ |
576 | 0 | linedp = lined - wpld; /* prev dest line */ |
577 | 0 | for (j = 0; j < w; j++) { |
578 | 0 | val32 = lines[j]; |
579 | 0 | if (j == 0) |
580 | 0 | lined[0] = val32 + linedp[0]; |
581 | 0 | else |
582 | 0 | lined[j] = val32 + lined[j - 1] + linedp[j] - linedp[j - 1]; |
583 | 0 | } |
584 | 0 | } |
585 | 0 | } else { |
586 | 0 | L_ERROR("depth not 1, 8 or 32 bpp\n", __func__); |
587 | 0 | } |
588 | 0 | } |
589 | | |
590 | | |
591 | | /*----------------------------------------------------------------------* |
592 | | * Un-normalized grayscale block convolution * |
593 | | *----------------------------------------------------------------------*/ |
594 | | /*! |
595 | | * \brief pixBlockconvGrayUnnormalized() |
596 | | * |
597 | | * \param[in] pixs 8 bpp |
598 | | * \param[in] wc, hc half width/height of convolution kernel |
599 | | * \return pix 32 bpp; containing the convolution without normalizing |
600 | | * for the window size, or NULL on error |
601 | | * |
602 | | * <pre> |
603 | | * Notes: |
604 | | * (1) The full width and height of the convolution kernel |
605 | | * are (2 * wc + 1) and (2 * hc + 1). |
606 | | * (2) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1, |
607 | | * where (w,h) are the dimensions of pixs. Attempt to |
608 | | * reduce the kernel size if necessary. |
609 | | * (3) Returns a copy if either wc or hc are 0. |
610 | | * (3) Adds mirrored border to avoid treating the boundary pixels |
611 | | * specially. Note that we add wc + 1 pixels to the left |
612 | | * and wc to the right. The added width is 2 * wc + 1 pixels, |
613 | | * and the particular choice simplifies the indexing in the loop. |
614 | | * Likewise, add hc + 1 pixels to the top and hc to the bottom. |
615 | | * (4) To get the normalized result, divide by the area of the |
616 | | * convolution kernel: (2 * wc + 1) * (2 * hc + 1) |
617 | | * Specifically, do this: |
618 | | * pixc = pixBlockconvGrayUnnormalized(pixs, wc, hc); |
619 | | * fract = 1. / ((2 * wc + 1) * (2 * hc + 1)); |
620 | | * pixMultConstantGray(pixc, fract); |
621 | | * pixd = pixGetRGBComponent(pixc, L_ALPHA_CHANNEL); |
622 | | * (5) Unlike pixBlockconvGray(), this always computes the accumulation |
623 | | * pix because its size is tied to wc and hc. |
624 | | * (6) Compare this implementation with pixBlockconvGray(), where |
625 | | * most of the code in blockconvLow() is special casing for |
626 | | * efficiently handling the boundary. Here, the use of |
627 | | * mirrored borders and destination indexing makes the |
628 | | * implementation very simple. |
629 | | * </pre> |
630 | | */ |
631 | | PIX * |
632 | | pixBlockconvGrayUnnormalized(PIX *pixs, |
633 | | l_int32 wc, |
634 | | l_int32 hc) |
635 | 0 | { |
636 | 0 | l_int32 i, j, w, h, d, wpla, wpld, jmax; |
637 | 0 | l_uint32 *linemina, *linemaxa, *lined, *dataa, *datad; |
638 | 0 | PIX *pixsb, *pixacc, *pixd; |
639 | |
|
640 | 0 | if (!pixs) |
641 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
642 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
643 | 0 | if (d != 8) |
644 | 0 | return (PIX *)ERROR_PTR("pixs not 8 bpp", __func__, NULL); |
645 | 0 | if (wc <= 0 || hc <= 0) /* no-op */ |
646 | 0 | return pixCopy(NULL, pixs); |
647 | 0 | if (w < 2 * wc + 1 || h < 2 * hc + 1) { |
648 | 0 | L_WARNING("kernel too large: wc = %d, hc = %d, w = %d, h = %d; " |
649 | 0 | "reducing!\n", __func__, wc, hc, w, h); |
650 | 0 | wc = L_MIN(wc, (w - 1) / 2); |
651 | 0 | hc = L_MIN(hc, (h - 1) / 2); |
652 | 0 | } |
653 | 0 | if (wc == 0 || hc == 0) |
654 | 0 | return pixCopy(NULL, pixs); |
655 | | |
656 | 0 | if ((pixsb = pixAddMirroredBorder(pixs, wc + 1, wc, hc + 1, hc)) == NULL) |
657 | 0 | return (PIX *)ERROR_PTR("pixsb not made", __func__, NULL); |
658 | 0 | pixacc = pixBlockconvAccum(pixsb); |
659 | 0 | pixDestroy(&pixsb); |
660 | 0 | if (!pixacc) |
661 | 0 | return (PIX *)ERROR_PTR("pixacc not made", __func__, NULL); |
662 | 0 | if ((pixd = pixCreate(w, h, 32)) == NULL) { |
663 | 0 | pixDestroy(&pixacc); |
664 | 0 | return (PIX *)ERROR_PTR("pixd not made", __func__, NULL); |
665 | 0 | } |
666 | | |
667 | 0 | wpla = pixGetWpl(pixacc); |
668 | 0 | wpld = pixGetWpl(pixd); |
669 | 0 | datad = pixGetData(pixd); |
670 | 0 | dataa = pixGetData(pixacc); |
671 | 0 | for (i = 0; i < h; i++) { |
672 | 0 | lined = datad + i * wpld; |
673 | 0 | linemina = dataa + i * wpla; |
674 | 0 | linemaxa = dataa + (i + 2 * hc + 1) * wpla; |
675 | 0 | for (j = 0; j < w; j++) { |
676 | 0 | jmax = j + 2 * wc + 1; |
677 | 0 | lined[j] = linemaxa[jmax] - linemaxa[j] - |
678 | 0 | linemina[jmax] + linemina[j]; |
679 | 0 | } |
680 | 0 | } |
681 | |
|
682 | 0 | pixDestroy(&pixacc); |
683 | 0 | return pixd; |
684 | 0 | } |
685 | | |
686 | | |
687 | | /*----------------------------------------------------------------------* |
688 | | * Tiled grayscale or color block convolution * |
689 | | *----------------------------------------------------------------------*/ |
690 | | /*! |
691 | | * \brief pixBlockconvTiled() |
692 | | * |
693 | | * \param[in] pix 8 or 32 bpp; or 2, 4 or 8 bpp with colormap |
694 | | * \param[in] wc, hc half width/height of convolution kernel |
695 | | * \param[in] nx, ny subdivision into tiles |
696 | | * \return pixd, or NULL on error |
697 | | * |
698 | | * <pre> |
699 | | * Notes: |
700 | | * (1) The full width and height of the convolution kernel |
701 | | * are (2 * wc + 1) and (2 * hc + 1) |
702 | | * (2) Returns a copy if either wc or hc are 0. |
703 | | * (3) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1, |
704 | | * where (w,h) are the dimensions of pixs. Attempt to |
705 | | * reduce the kernel size if necessary. |
706 | | * (4) For nx == ny == 1, this defaults to pixBlockconv(), which |
707 | | * is typically about twice as fast, and gives nearly |
708 | | * identical results as pixBlockconvGrayTile(). |
709 | | * (5) If the tiles are too small, nx and/or ny are reduced |
710 | | * a minimum amount so that the tiles are expanded to the |
711 | | * smallest workable size in the problematic direction(s). |
712 | | * (6) Why a tiled version? Three reasons: |
713 | | * (a) Because the accumulator is a uint32, overflow can occur |
714 | | * for an image with more than 16M pixels. |
715 | | * (b) The accumulator array for 16M pixels is 64 MB; using |
716 | | * tiles reduces the size of this array. |
717 | | * (c) Each tile can be processed independently, in parallel, |
718 | | * on a multicore processor. |
719 | | * </pre> |
720 | | */ |
721 | | PIX * |
722 | | pixBlockconvTiled(PIX *pix, |
723 | | l_int32 wc, |
724 | | l_int32 hc, |
725 | | l_int32 nx, |
726 | | l_int32 ny) |
727 | 0 | { |
728 | 0 | l_int32 i, j, w, h, d, xrat, yrat; |
729 | 0 | PIX *pixs, *pixd, *pixc, *pixt; |
730 | 0 | PIX *pixr, *pixrc, *pixg, *pixgc, *pixb, *pixbc; |
731 | 0 | PIXTILING *pt; |
732 | |
|
733 | 0 | if (!pix) |
734 | 0 | return (PIX *)ERROR_PTR("pix not defined", __func__, NULL); |
735 | 0 | if (wc <= 0 || hc <= 0) /* no-op */ |
736 | 0 | return pixCopy(NULL, pix); |
737 | 0 | if (nx <= 1 && ny <= 1) |
738 | 0 | return pixBlockconv(pix, wc, hc); |
739 | 0 | pixGetDimensions(pix, &w, &h, &d); |
740 | 0 | if (w < 2 * wc + 3 || h < 2 * hc + 3) { |
741 | 0 | L_WARNING("kernel too large: wc = %d, hc = %d, w = %d, h = %d; " |
742 | 0 | "reducing!\n", __func__, wc, hc, w, h); |
743 | 0 | wc = L_MIN(wc, (w - 1) / 2); |
744 | 0 | hc = L_MIN(hc, (h - 1) / 2); |
745 | 0 | } |
746 | 0 | if (wc == 0 || hc == 0) |
747 | 0 | return pixCopy(NULL, pix); |
748 | | |
749 | | /* Test to see if the tiles are too small. The required |
750 | | * condition is that the tile dimensions must be at least |
751 | | * (wc + 2) x (hc + 2). */ |
752 | 0 | xrat = w / nx; |
753 | 0 | yrat = h / ny; |
754 | 0 | if (xrat < wc + 2) { |
755 | 0 | nx = w / (wc + 2); |
756 | 0 | L_WARNING("tile width too small; nx reduced to %d\n", __func__, nx); |
757 | 0 | } |
758 | 0 | if (yrat < hc + 2) { |
759 | 0 | ny = h / (hc + 2); |
760 | 0 | L_WARNING("tile height too small; ny reduced to %d\n", __func__, ny); |
761 | 0 | } |
762 | | |
763 | | /* Remove colormap if necessary */ |
764 | 0 | if ((d == 2 || d == 4 || d == 8) && pixGetColormap(pix)) { |
765 | 0 | L_WARNING("pix has colormap; removing\n", __func__); |
766 | 0 | pixs = pixRemoveColormap(pix, REMOVE_CMAP_BASED_ON_SRC); |
767 | 0 | d = pixGetDepth(pixs); |
768 | 0 | } else { |
769 | 0 | pixs = pixClone(pix); |
770 | 0 | } |
771 | |
|
772 | 0 | if (d != 8 && d != 32) { |
773 | 0 | pixDestroy(&pixs); |
774 | 0 | return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", __func__, NULL); |
775 | 0 | } |
776 | | |
777 | | /* Note that the overlaps in the width and height that |
778 | | * are added to the tile are (wc + 2) and (hc + 2). |
779 | | * These overlaps are removed by pixTilingPaintTile(). |
780 | | * They are larger than the extent of the filter because |
781 | | * although the filter is symmetric with respect to its origin, |
782 | | * the implementation is asymmetric -- see the implementation in |
783 | | * pixBlockconvGrayTile(). */ |
784 | 0 | if ((pixd = pixCreateTemplate(pixs)) == NULL) { |
785 | 0 | pixDestroy(&pixs); |
786 | 0 | return (PIX *)ERROR_PTR("pixd not made", __func__, NULL); |
787 | 0 | } |
788 | 0 | pt = pixTilingCreate(pixs, nx, ny, 0, 0, wc + 2, hc + 2); |
789 | 0 | for (i = 0; i < ny; i++) { |
790 | 0 | for (j = 0; j < nx; j++) { |
791 | 0 | pixt = pixTilingGetTile(pt, i, j); |
792 | | |
793 | | /* Convolve over the tile */ |
794 | 0 | if (d == 8) { |
795 | 0 | pixc = pixBlockconvGrayTile(pixt, NULL, wc, hc); |
796 | 0 | } else { /* d == 32 */ |
797 | 0 | pixr = pixGetRGBComponent(pixt, COLOR_RED); |
798 | 0 | pixrc = pixBlockconvGrayTile(pixr, NULL, wc, hc); |
799 | 0 | pixDestroy(&pixr); |
800 | 0 | pixg = pixGetRGBComponent(pixt, COLOR_GREEN); |
801 | 0 | pixgc = pixBlockconvGrayTile(pixg, NULL, wc, hc); |
802 | 0 | pixDestroy(&pixg); |
803 | 0 | pixb = pixGetRGBComponent(pixt, COLOR_BLUE); |
804 | 0 | pixbc = pixBlockconvGrayTile(pixb, NULL, wc, hc); |
805 | 0 | pixDestroy(&pixb); |
806 | 0 | pixc = pixCreateRGBImage(pixrc, pixgc, pixbc); |
807 | 0 | pixDestroy(&pixrc); |
808 | 0 | pixDestroy(&pixgc); |
809 | 0 | pixDestroy(&pixbc); |
810 | 0 | } |
811 | |
|
812 | 0 | pixTilingPaintTile(pixd, i, j, pixc, pt); |
813 | 0 | pixDestroy(&pixt); |
814 | 0 | pixDestroy(&pixc); |
815 | 0 | } |
816 | 0 | } |
817 | |
|
818 | 0 | pixDestroy(&pixs); |
819 | 0 | pixTilingDestroy(&pt); |
820 | 0 | return pixd; |
821 | 0 | } |
822 | | |
823 | | |
824 | | /*! |
825 | | * \brief pixBlockconvGrayTile() |
826 | | * |
827 | | * \param[in] pixs 8 bpp gray |
828 | | * \param[in] pixacc 32 bpp accum pix |
829 | | * \param[in] wc, hc half width/height of convolution kernel |
830 | | * \return pixd, or NULL on error |
831 | | * |
832 | | * <pre> |
833 | | * Notes: |
834 | | * (1) The full width and height of the convolution kernel |
835 | | * are (2 * wc + 1) and (2 * hc + 1) |
836 | | * (2) Assumes that the input pixs is padded with (wc + 1) pixels on |
837 | | * left and right, and with (hc + 1) pixels on top and bottom. |
838 | | * The returned pix has these stripped off; they are only used |
839 | | * for computation. |
840 | | * (3) Returns a copy if either wc or hc are 0. |
841 | | * (4) Require that w > 2 * wc + 3 and h > 2 * hc + 3, |
842 | | * where (w,h) are the dimensions of pixs. Attempt to |
843 | | * reduce the kernel size if necessary. |
844 | | * </pre> |
845 | | */ |
846 | | PIX * |
847 | | pixBlockconvGrayTile(PIX *pixs, |
848 | | PIX *pixacc, |
849 | | l_int32 wc, |
850 | | l_int32 hc) |
851 | 0 | { |
852 | 0 | l_int32 w, h, d, wd, hd, i, j, imin, imax, jmin, jmax, wplt, wpld; |
853 | 0 | l_float32 norm; |
854 | 0 | l_uint32 val; |
855 | 0 | l_uint32 *datat, *datad, *lined, *linemint, *linemaxt; |
856 | 0 | PIX *pixt, *pixd; |
857 | |
|
858 | 0 | if (!pixs) |
859 | 0 | return (PIX *)ERROR_PTR("pix not defined", __func__, NULL); |
860 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
861 | 0 | if (d != 8) |
862 | 0 | return (PIX *)ERROR_PTR("pixs not 8 bpp", __func__, NULL); |
863 | 0 | if (wc <= 0 || hc <= 0) /* no-op */ |
864 | 0 | return pixCopy(NULL, pixs); |
865 | 0 | if (w < 2 * wc + 3 || h < 2 * hc + 3) { |
866 | 0 | L_WARNING("kernel too large: wc = %d, hc = %d, w = %d, h = %d; " |
867 | 0 | "reducing!\n", __func__, wc, hc, w, h); |
868 | 0 | wc = L_MIN(wc, (w - 1) / 2); |
869 | 0 | hc = L_MIN(hc, (h - 1) / 2); |
870 | 0 | } |
871 | 0 | if (wc == 0 || hc == 0) |
872 | 0 | return pixCopy(NULL, pixs); |
873 | 0 | wd = w - 2 * wc; |
874 | 0 | hd = h - 2 * hc; |
875 | |
|
876 | 0 | if (pixacc) { |
877 | 0 | if (pixGetDepth(pixacc) == 32) { |
878 | 0 | pixt = pixClone(pixacc); |
879 | 0 | } else { |
880 | 0 | L_WARNING("pixacc not 32 bpp; making new one\n", __func__); |
881 | 0 | if ((pixt = pixBlockconvAccum(pixs)) == NULL) |
882 | 0 | return (PIX *)ERROR_PTR("pixt not made", __func__, NULL); |
883 | 0 | } |
884 | 0 | } else { |
885 | 0 | if ((pixt = pixBlockconvAccum(pixs)) == NULL) |
886 | 0 | return (PIX *)ERROR_PTR("pixt not made", __func__, NULL); |
887 | 0 | } |
888 | | |
889 | 0 | if ((pixd = pixCreateTemplate(pixs)) == NULL) { |
890 | 0 | pixDestroy(&pixt); |
891 | 0 | return (PIX *)ERROR_PTR("pixd not made", __func__, NULL); |
892 | 0 | } |
893 | 0 | datat = pixGetData(pixt); |
894 | 0 | wplt = pixGetWpl(pixt); |
895 | 0 | datad = pixGetData(pixd); |
896 | 0 | wpld = pixGetWpl(pixd); |
897 | 0 | norm = 1. / (l_float32)((2 * wc + 1) * (2 * hc + 1)); |
898 | | |
899 | | /* Do the convolution over the subregion of size (wd - 2, hd - 2), |
900 | | * which exactly corresponds to the size of the subregion that |
901 | | * will be extracted by pixTilingPaintTile(). Note that the |
902 | | * region in which points are computed is not symmetric about |
903 | | * the center of the images; instead the computation in |
904 | | * the accumulator image is shifted up and to the left by 1, |
905 | | * relative to the center, because the 4 accumulator sampling |
906 | | * points are taken at the LL corner of the filter and at 3 other |
907 | | * points that are shifted -wc and -hc to the left and above. */ |
908 | 0 | for (i = hc; i < hc + hd - 2; i++) { |
909 | 0 | imin = L_MAX(i - hc - 1, 0); |
910 | 0 | imax = L_MIN(i + hc, h - 1); |
911 | 0 | lined = datad + i * wpld; |
912 | 0 | linemint = datat + imin * wplt; |
913 | 0 | linemaxt = datat + imax * wplt; |
914 | 0 | for (j = wc; j < wc + wd - 2; j++) { |
915 | 0 | jmin = L_MAX(j - wc - 1, 0); |
916 | 0 | jmax = L_MIN(j + wc, w - 1); |
917 | 0 | val = linemaxt[jmax] - linemaxt[jmin] |
918 | 0 | + linemint[jmin] - linemint[jmax]; |
919 | 0 | val = (l_uint8)(norm * val + 0.5); |
920 | 0 | SET_DATA_BYTE(lined, j, val); |
921 | 0 | } |
922 | 0 | } |
923 | |
|
924 | 0 | pixDestroy(&pixt); |
925 | 0 | return pixd; |
926 | 0 | } |
927 | | |
928 | | |
929 | | /*----------------------------------------------------------------------* |
930 | | * Convolution for mean, mean square, variance and rms deviation * |
931 | | *----------------------------------------------------------------------*/ |
932 | | /*! |
933 | | * \brief pixWindowedStats() |
934 | | * |
935 | | * \param[in] pixs 8 bpp grayscale |
936 | | * \param[in] wc, hc half width/height of convolution kernel |
937 | | * \param[in] hasborder use 1 if it already has (wc + 1 border pixels |
938 | | * on left and right, and hc + 1 on top and bottom; |
939 | | * use 0 to add kernel-dependent border) |
940 | | * \param[out] ppixm [optional] 8 bpp mean value in window |
941 | | * \param[out] ppixms [optional] 32 bpp mean square value in window |
942 | | * \param[out] pfpixv [optional] float variance in window |
943 | | * \param[out] pfpixrv [optional] float rms deviation from the mean |
944 | | * \return 0 if OK, 1 on error |
945 | | * |
946 | | * <pre> |
947 | | * Notes: |
948 | | * (1) This is a high-level convenience function for calculating |
949 | | * any or all of these derived images. |
950 | | * (2) If %hasborder = 0, a border is added and the result is |
951 | | * computed over all pixels in pixs. Otherwise, no border is |
952 | | * added and the border pixels are removed from the output images. |
953 | | * (3) These statistical measures over the pixels in the |
954 | | * rectangular window are: |
955 | | * ~ average value: <p> (pixm) |
956 | | * ~ average squared value: <p*p> (pixms) |
957 | | * ~ variance: <(p - <p>)*(p - <p>)> = <p*p> - <p>*<p> (pixv) |
958 | | * ~ square-root of variance: (pixrv) |
959 | | * where the brackets < .. > indicate that the average value is |
960 | | * to be taken over the window. |
961 | | * (4) Note that the variance is just the mean square difference from |
962 | | * the mean value; and the square root of the variance is the |
963 | | * root mean square difference from the mean, sometimes also |
964 | | * called the 'standard deviation'. |
965 | | * (5) The added border, along with the use of an accumulator array, |
966 | | * allows computation without special treatment of pixels near |
967 | | * the image boundary, and runs in a time that is independent |
968 | | * of the size of the convolution kernel. |
969 | | * </pre> |
970 | | */ |
971 | | l_ok |
972 | | pixWindowedStats(PIX *pixs, |
973 | | l_int32 wc, |
974 | | l_int32 hc, |
975 | | l_int32 hasborder, |
976 | | PIX **ppixm, |
977 | | PIX **ppixms, |
978 | | FPIX **pfpixv, |
979 | | FPIX **pfpixrv) |
980 | 0 | { |
981 | 0 | PIX *pixb, *pixm, *pixms; |
982 | |
|
983 | 0 | if (!ppixm && !ppixms && !pfpixv && !pfpixrv) |
984 | 0 | return ERROR_INT("no output requested", __func__, 1); |
985 | 0 | if (ppixm) *ppixm = NULL; |
986 | 0 | if (ppixms) *ppixms = NULL; |
987 | 0 | if (pfpixv) *pfpixv = NULL; |
988 | 0 | if (pfpixrv) *pfpixrv = NULL; |
989 | 0 | if (!pixs || pixGetDepth(pixs) != 8) |
990 | 0 | return ERROR_INT("pixs not defined or not 8 bpp", __func__, 1); |
991 | 0 | if (wc < 2 || hc < 2) |
992 | 0 | return ERROR_INT("wc and hc not >= 2", __func__, 1); |
993 | | |
994 | | /* Add border if requested */ |
995 | 0 | if (!hasborder) |
996 | 0 | pixb = pixAddBorderGeneral(pixs, wc + 1, wc + 1, hc + 1, hc + 1, 0); |
997 | 0 | else |
998 | 0 | pixb = pixClone(pixs); |
999 | |
|
1000 | 0 | if (!pfpixv && !pfpixrv) { |
1001 | 0 | if (ppixm) *ppixm = pixWindowedMean(pixb, wc, hc, 1, 1); |
1002 | 0 | if (ppixms) *ppixms = pixWindowedMeanSquare(pixb, wc, hc, 1); |
1003 | 0 | pixDestroy(&pixb); |
1004 | 0 | return 0; |
1005 | 0 | } |
1006 | | |
1007 | 0 | pixm = pixWindowedMean(pixb, wc, hc, 1, 1); |
1008 | 0 | pixms = pixWindowedMeanSquare(pixb, wc, hc, 1); |
1009 | 0 | pixWindowedVariance(pixm, pixms, pfpixv, pfpixrv); |
1010 | 0 | if (ppixm) |
1011 | 0 | *ppixm = pixm; |
1012 | 0 | else |
1013 | 0 | pixDestroy(&pixm); |
1014 | 0 | if (ppixms) |
1015 | 0 | *ppixms = pixms; |
1016 | 0 | else |
1017 | 0 | pixDestroy(&pixms); |
1018 | 0 | pixDestroy(&pixb); |
1019 | 0 | return 0; |
1020 | 0 | } |
1021 | | |
1022 | | |
1023 | | /*! |
1024 | | * \brief pixWindowedMean() |
1025 | | * |
1026 | | * \param[in] pixs 8 or 32 bpp grayscale |
1027 | | * \param[in] wc, hc half width/height of convolution kernel |
1028 | | * \param[in] hasborder use 1 if it already has (wc + 1 border pixels |
1029 | | * on left and right, and hc + 1 on top and bottom; |
1030 | | * use 0 to add kernel-dependent border) |
1031 | | * \param[in] normflag 1 for normalization to get average in window; |
1032 | | * 0 for the sum in the window (un-normalized) |
1033 | | * \return pixd 8 or 32 bpp, average over kernel window |
1034 | | * |
1035 | | * <pre> |
1036 | | * Notes: |
1037 | | * (1) The input and output depths are the same. |
1038 | | * (2) A set of border pixels of width (wc + 1) on left and right, |
1039 | | * and of height (hc + 1) on top and bottom, must be on the |
1040 | | * pix before the accumulator is found. The output pixd |
1041 | | * (after convolution) has this border removed. |
1042 | | * If %hasborder = 0, the required border is added. |
1043 | | * (3) Typically, %normflag == 1. However, if you want the sum |
1044 | | * within the window, rather than a normalized convolution, |
1045 | | * use %normflag == 0. |
1046 | | * (4) This builds a block accumulator pix, uses it here, and |
1047 | | * destroys it. |
1048 | | * (5) The added border, along with the use of an accumulator array, |
1049 | | * allows computation without special treatment of pixels near |
1050 | | * the image boundary, and runs in a time that is independent |
1051 | | * of the size of the convolution kernel. |
1052 | | * </pre> |
1053 | | */ |
1054 | | PIX * |
1055 | | pixWindowedMean(PIX *pixs, |
1056 | | l_int32 wc, |
1057 | | l_int32 hc, |
1058 | | l_int32 hasborder, |
1059 | | l_int32 normflag) |
1060 | 0 | { |
1061 | 0 | l_int32 i, j, w, h, d, wd, hd, wplc, wpld, wincr, hincr; |
1062 | 0 | l_uint32 val; |
1063 | 0 | l_uint32 *datac, *datad, *linec1, *linec2, *lined; |
1064 | 0 | l_float32 norm; |
1065 | 0 | PIX *pixb, *pixc, *pixd; |
1066 | |
|
1067 | 0 | if (!pixs) |
1068 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
1069 | 0 | d = pixGetDepth(pixs); |
1070 | 0 | if (d != 8 && d != 32) |
1071 | 0 | return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL); |
1072 | 0 | if (wc < 2 || hc < 2) |
1073 | 0 | return (PIX *)ERROR_PTR("wc and hc not >= 2", __func__, NULL); |
1074 | | |
1075 | 0 | pixb = pixc = pixd = NULL; |
1076 | | |
1077 | | /* Add border if requested */ |
1078 | 0 | if (!hasborder) |
1079 | 0 | pixb = pixAddBorderGeneral(pixs, wc + 1, wc + 1, hc + 1, hc + 1, 0); |
1080 | 0 | else |
1081 | 0 | pixb = pixClone(pixs); |
1082 | | |
1083 | | /* Make the accumulator pix from pixb */ |
1084 | 0 | if ((pixc = pixBlockconvAccum(pixb)) == NULL) { |
1085 | 0 | L_ERROR("pixc not made\n", __func__); |
1086 | 0 | goto cleanup; |
1087 | 0 | } |
1088 | 0 | wplc = pixGetWpl(pixc); |
1089 | 0 | datac = pixGetData(pixc); |
1090 | | |
1091 | | /* The output has wc + 1 border pixels stripped from each side |
1092 | | * of pixb, and hc + 1 border pixels stripped from top and bottom. */ |
1093 | 0 | pixGetDimensions(pixb, &w, &h, NULL); |
1094 | 0 | wd = w - 2 * (wc + 1); |
1095 | 0 | hd = h - 2 * (hc + 1); |
1096 | 0 | if (wd < 2 || hd < 2) { |
1097 | 0 | L_ERROR("w or h is too small for the kernel\n", __func__); |
1098 | 0 | goto cleanup; |
1099 | 0 | } |
1100 | 0 | if ((pixd = pixCreate(wd, hd, d)) == NULL) { |
1101 | 0 | L_ERROR("pixd not made\n", __func__); |
1102 | 0 | goto cleanup; |
1103 | 0 | } |
1104 | 0 | wpld = pixGetWpl(pixd); |
1105 | 0 | datad = pixGetData(pixd); |
1106 | |
|
1107 | 0 | wincr = 2 * wc + 1; |
1108 | 0 | hincr = 2 * hc + 1; |
1109 | 0 | norm = 1.0; /* use this for sum-in-window */ |
1110 | 0 | if (normflag) |
1111 | 0 | norm = 1.0 / ((l_float32)(wincr) * hincr); |
1112 | 0 | for (i = 0; i < hd; i++) { |
1113 | 0 | linec1 = datac + i * wplc; |
1114 | 0 | linec2 = datac + (i + hincr) * wplc; |
1115 | 0 | lined = datad + i * wpld; |
1116 | 0 | for (j = 0; j < wd; j++) { |
1117 | 0 | val = linec2[j + wincr] - linec2[j] - linec1[j + wincr] + linec1[j]; |
1118 | 0 | if (d == 8) { |
1119 | 0 | val = (l_uint8)(norm * val); |
1120 | 0 | SET_DATA_BYTE(lined, j, val); |
1121 | 0 | } else { /* d == 32 */ |
1122 | 0 | val = (l_uint32)(norm * val); |
1123 | 0 | lined[j] = val; |
1124 | 0 | } |
1125 | 0 | } |
1126 | 0 | } |
1127 | |
|
1128 | 0 | cleanup: |
1129 | 0 | pixDestroy(&pixb); |
1130 | 0 | pixDestroy(&pixc); |
1131 | 0 | return pixd; |
1132 | 0 | } |
1133 | | |
1134 | | |
1135 | | /*! |
1136 | | * \brief pixWindowedMeanSquare() |
1137 | | * |
1138 | | * \param[in] pixs 8 bpp grayscale |
1139 | | * \param[in] wc, hc half width/height of convolution kernel |
1140 | | * \param[in] hasborder use 1 if it already has (wc + 1 border pixels |
1141 | | * on left and right, and hc + 1 on top and bottom; |
1142 | | * use 0 to add kernel-dependent border) |
1143 | | * \return pixd 32 bpp, average over rectangular window of |
1144 | | * width = 2 * wc + 1 and height = 2 * hc + 1 |
1145 | | * |
1146 | | * <pre> |
1147 | | * Notes: |
1148 | | * (1) A set of border pixels of width (wc + 1) on left and right, |
1149 | | * and of height (hc + 1) on top and bottom, must be on the |
1150 | | * pix before the accumulator is found. The output pixd |
1151 | | * (after convolution) has this border removed. |
1152 | | * If %hasborder = 0, the required border is added. |
1153 | | * (2) The advantage is that we are unaffected by the boundary, and |
1154 | | * it is not necessary to treat pixels within %wc and %hc of the |
1155 | | * border differently. This is because processing for pixd |
1156 | | * only takes place for pixels in pixs for which the |
1157 | | * kernel is entirely contained in pixs. |
1158 | | * (3) Why do we have an added border of width (%wc + 1) and |
1159 | | * height (%hc + 1), when we only need %wc and %hc pixels |
1160 | | * to satisfy this condition? Answer: the accumulators |
1161 | | * are asymmetric, requiring an extra row and column of |
1162 | | * pixels at top and left to work accurately. |
1163 | | * (4) The added border, along with the use of an accumulator array, |
1164 | | * allows computation without special treatment of pixels near |
1165 | | * the image boundary, and runs in a time that is independent |
1166 | | * of the size of the convolution kernel. |
1167 | | * </pre> |
1168 | | */ |
1169 | | PIX * |
1170 | | pixWindowedMeanSquare(PIX *pixs, |
1171 | | l_int32 wc, |
1172 | | l_int32 hc, |
1173 | | l_int32 hasborder) |
1174 | 0 | { |
1175 | 0 | l_int32 i, j, w, h, wd, hd, wpl, wpld, wincr, hincr; |
1176 | 0 | l_uint32 ival; |
1177 | 0 | l_uint32 *datad, *lined; |
1178 | 0 | l_float64 norm; |
1179 | 0 | l_float64 val; |
1180 | 0 | l_float64 *data, *line1, *line2; |
1181 | 0 | DPIX *dpix; |
1182 | 0 | PIX *pixb, *pixd; |
1183 | |
|
1184 | 0 | if (!pixs || (pixGetDepth(pixs) != 8)) |
1185 | 0 | return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", __func__, NULL); |
1186 | 0 | if (wc < 2 || hc < 2) |
1187 | 0 | return (PIX *)ERROR_PTR("wc and hc not >= 2", __func__, NULL); |
1188 | | |
1189 | 0 | pixd = NULL; |
1190 | | |
1191 | | /* Add border if requested */ |
1192 | 0 | if (!hasborder) |
1193 | 0 | pixb = pixAddBorderGeneral(pixs, wc + 1, wc + 1, hc + 1, hc + 1, 0); |
1194 | 0 | else |
1195 | 0 | pixb = pixClone(pixs); |
1196 | |
|
1197 | 0 | if ((dpix = pixMeanSquareAccum(pixb)) == NULL) { |
1198 | 0 | L_ERROR("dpix not made\n", __func__); |
1199 | 0 | goto cleanup; |
1200 | 0 | } |
1201 | 0 | wpl = dpixGetWpl(dpix); |
1202 | 0 | data = dpixGetData(dpix); |
1203 | | |
1204 | | /* The output has wc + 1 border pixels stripped from each side |
1205 | | * of pixb, and hc + 1 border pixels stripped from top and bottom. */ |
1206 | 0 | pixGetDimensions(pixb, &w, &h, NULL); |
1207 | 0 | wd = w - 2 * (wc + 1); |
1208 | 0 | hd = h - 2 * (hc + 1); |
1209 | 0 | if (wd < 2 || hd < 2) { |
1210 | 0 | L_ERROR("w or h too small for kernel\n", __func__); |
1211 | 0 | goto cleanup; |
1212 | 0 | } |
1213 | 0 | if ((pixd = pixCreate(wd, hd, 32)) == NULL) { |
1214 | 0 | L_ERROR("pixd not made\n", __func__); |
1215 | 0 | goto cleanup; |
1216 | 0 | } |
1217 | 0 | wpld = pixGetWpl(pixd); |
1218 | 0 | datad = pixGetData(pixd); |
1219 | |
|
1220 | 0 | wincr = 2 * wc + 1; |
1221 | 0 | hincr = 2 * hc + 1; |
1222 | 0 | norm = 1.0 / ((l_float32)(wincr) * hincr); |
1223 | 0 | for (i = 0; i < hd; i++) { |
1224 | 0 | line1 = data + i * wpl; |
1225 | 0 | line2 = data + (i + hincr) * wpl; |
1226 | 0 | lined = datad + i * wpld; |
1227 | 0 | for (j = 0; j < wd; j++) { |
1228 | 0 | val = line2[j + wincr] - line2[j] - line1[j + wincr] + line1[j]; |
1229 | 0 | ival = (l_uint32)(norm * val + 0.5); /* to round up */ |
1230 | 0 | lined[j] = ival; |
1231 | 0 | } |
1232 | 0 | } |
1233 | |
|
1234 | 0 | cleanup: |
1235 | 0 | dpixDestroy(&dpix); |
1236 | 0 | pixDestroy(&pixb); |
1237 | 0 | return pixd; |
1238 | 0 | } |
1239 | | |
1240 | | |
1241 | | /*! |
1242 | | * \brief pixWindowedVariance() |
1243 | | * |
1244 | | * \param[in] pixm mean over window; 8 or 32 bpp grayscale |
1245 | | * \param[in] pixms mean square over window; 32 bpp |
1246 | | * \param[out] pfpixv [optional] float variance -- the ms deviation |
1247 | | * from the mean |
1248 | | * \param[out] pfpixrv [optional] float rms deviation from the mean |
1249 | | * \return 0 if OK, 1 on error |
1250 | | * |
1251 | | * <pre> |
1252 | | * Notes: |
1253 | | * (1) The mean and mean square values are precomputed, using |
1254 | | * pixWindowedMean() and pixWindowedMeanSquare(). |
1255 | | * (2) Either or both of the variance and square-root of variance |
1256 | | * are returned as an fpix, where the variance is the |
1257 | | * average over the window of the mean square difference of |
1258 | | * the pixel value from the mean: |
1259 | | * <(p - <p>)*(p - <p>)> = <p*p> - <p>*<p> |
1260 | | * (3) To visualize the results: |
1261 | | * ~ for both, use fpixDisplayMaxDynamicRange(). |
1262 | | * ~ for rms deviation, simply convert the output fpix to pix, |
1263 | | * </pre> |
1264 | | */ |
1265 | | l_ok |
1266 | | pixWindowedVariance(PIX *pixm, |
1267 | | PIX *pixms, |
1268 | | FPIX **pfpixv, |
1269 | | FPIX **pfpixrv) |
1270 | 0 | { |
1271 | 0 | l_int32 i, j, w, h, ws, hs, ds, wplm, wplms, wplv, wplrv, valm, valms; |
1272 | 0 | l_float32 var; |
1273 | 0 | l_uint32 *linem, *linems, *datam, *datams; |
1274 | 0 | l_float32 *linev = NULL, *linerv = NULL, *datav = NULL, *datarv = NULL; |
1275 | 0 | FPIX *fpixv, *fpixrv; /* variance and square root of variance */ |
1276 | |
|
1277 | 0 | if (!pfpixv && !pfpixrv) |
1278 | 0 | return ERROR_INT("no output requested", __func__, 1); |
1279 | 0 | if (pfpixv) *pfpixv = NULL; |
1280 | 0 | if (pfpixrv) *pfpixrv = NULL; |
1281 | 0 | if (!pixm || pixGetDepth(pixm) != 8) |
1282 | 0 | return ERROR_INT("pixm undefined or not 8 bpp", __func__, 1); |
1283 | 0 | if (!pixms || pixGetDepth(pixms) != 32) |
1284 | 0 | return ERROR_INT("pixms undefined or not 32 bpp", __func__, 1); |
1285 | 0 | pixGetDimensions(pixm, &w, &h, NULL); |
1286 | 0 | pixGetDimensions(pixms, &ws, &hs, &ds); |
1287 | 0 | if (w != ws || h != hs) |
1288 | 0 | return ERROR_INT("pixm and pixms sizes differ", __func__, 1); |
1289 | | |
1290 | 0 | if (pfpixv) { |
1291 | 0 | fpixv = fpixCreate(w, h); |
1292 | 0 | *pfpixv = fpixv; |
1293 | 0 | wplv = fpixGetWpl(fpixv); |
1294 | 0 | datav = fpixGetData(fpixv); |
1295 | 0 | } |
1296 | 0 | if (pfpixrv) { |
1297 | 0 | fpixrv = fpixCreate(w, h); |
1298 | 0 | *pfpixrv = fpixrv; |
1299 | 0 | wplrv = fpixGetWpl(fpixrv); |
1300 | 0 | datarv = fpixGetData(fpixrv); |
1301 | 0 | } |
1302 | |
|
1303 | 0 | wplm = pixGetWpl(pixm); |
1304 | 0 | wplms = pixGetWpl(pixms); |
1305 | 0 | datam = pixGetData(pixm); |
1306 | 0 | datams = pixGetData(pixms); |
1307 | 0 | for (i = 0; i < h; i++) { |
1308 | 0 | linem = datam + i * wplm; |
1309 | 0 | linems = datams + i * wplms; |
1310 | 0 | if (pfpixv) |
1311 | 0 | linev = datav + i * wplv; |
1312 | 0 | if (pfpixrv) |
1313 | 0 | linerv = datarv + i * wplrv; |
1314 | 0 | for (j = 0; j < w; j++) { |
1315 | 0 | valm = GET_DATA_BYTE(linem, j); |
1316 | 0 | if (ds == 8) |
1317 | 0 | valms = GET_DATA_BYTE(linems, j); |
1318 | 0 | else /* ds == 32 */ |
1319 | 0 | valms = (l_int32)linems[j]; |
1320 | 0 | var = (l_float32)valms - (l_float32)valm * valm; |
1321 | 0 | if (pfpixv) |
1322 | 0 | linev[j] = var; |
1323 | 0 | if (pfpixrv) |
1324 | 0 | linerv[j] = (l_float32)sqrt(var); |
1325 | 0 | } |
1326 | 0 | } |
1327 | |
|
1328 | 0 | return 0; |
1329 | 0 | } |
1330 | | |
1331 | | |
1332 | | /*! |
1333 | | * \brief pixMeanSquareAccum() |
1334 | | * |
1335 | | * \param[in] pixs 8 bpp grayscale |
1336 | | * \return dpix 64 bit array, or NULL on error |
1337 | | * |
1338 | | * <pre> |
1339 | | * Notes: |
1340 | | * (1) Similar to pixBlockconvAccum(), this computes the |
1341 | | * sum of the squares of the pixel values in such a way |
1342 | | * that the value at (i,j) is the sum of all squares in |
1343 | | * the rectangle from the origin to (i,j). |
1344 | | * (2) The general recursion relation (v are squared pixel values) is |
1345 | | * a(i,j) = v(i,j) + a(i-1, j) + a(i, j-1) - a(i-1, j-1) |
1346 | | * For the first line, this reduces to the special case |
1347 | | * a(i,j) = v(i,j) + a(i, j-1) |
1348 | | * For the first column, the special case is |
1349 | | * a(i,j) = v(i,j) + a(i-1, j) |
1350 | | * </pre> |
1351 | | */ |
1352 | | DPIX * |
1353 | | pixMeanSquareAccum(PIX *pixs) |
1354 | 0 | { |
1355 | 0 | l_int32 i, j, w, h, wpl, wpls, val; |
1356 | 0 | l_uint32 *datas, *lines; |
1357 | 0 | l_float64 *data, *line, *linep; |
1358 | 0 | DPIX *dpix; |
1359 | |
|
1360 | 0 | if (!pixs || (pixGetDepth(pixs) != 8)) |
1361 | 0 | return (DPIX *)ERROR_PTR("pixs undefined or not 8 bpp", __func__, NULL); |
1362 | 0 | pixGetDimensions(pixs, &w, &h, NULL); |
1363 | 0 | if ((dpix = dpixCreate(w, h)) == NULL) |
1364 | 0 | return (DPIX *)ERROR_PTR("dpix not made", __func__, NULL); |
1365 | | |
1366 | 0 | datas = pixGetData(pixs); |
1367 | 0 | wpls = pixGetWpl(pixs); |
1368 | 0 | data = dpixGetData(dpix); |
1369 | 0 | wpl = dpixGetWpl(dpix); |
1370 | |
|
1371 | 0 | lines = datas; |
1372 | 0 | line = data; |
1373 | 0 | for (j = 0; j < w; j++) { /* first line */ |
1374 | 0 | val = GET_DATA_BYTE(lines, j); |
1375 | 0 | if (j == 0) |
1376 | 0 | line[0] = (l_float64)(val) * val; |
1377 | 0 | else |
1378 | 0 | line[j] = line[j - 1] + (l_float64)(val) * val; |
1379 | 0 | } |
1380 | | |
1381 | | /* Do the other lines */ |
1382 | 0 | for (i = 1; i < h; i++) { |
1383 | 0 | lines = datas + i * wpls; |
1384 | 0 | line = data + i * wpl; /* current dest line */ |
1385 | 0 | linep = line - wpl;; /* prev dest line */ |
1386 | 0 | for (j = 0; j < w; j++) { |
1387 | 0 | val = GET_DATA_BYTE(lines, j); |
1388 | 0 | if (j == 0) |
1389 | 0 | line[0] = linep[0] + (l_float64)(val) * val; |
1390 | 0 | else |
1391 | 0 | line[j] = line[j - 1] + linep[j] - linep[j - 1] |
1392 | 0 | + (l_float64)(val) * val; |
1393 | 0 | } |
1394 | 0 | } |
1395 | |
|
1396 | 0 | return dpix; |
1397 | 0 | } |
1398 | | |
1399 | | |
1400 | | /*----------------------------------------------------------------------* |
1401 | | * Binary block sum/rank * |
1402 | | *----------------------------------------------------------------------*/ |
1403 | | /*! |
1404 | | * \brief pixBlockrank() |
1405 | | * |
1406 | | * \param[in] pixs 1 bpp |
1407 | | * \param[in] pixacc pix [optional] 32 bpp |
1408 | | * \param[in] wc, hc half width/height of block sum/rank kernel |
1409 | | * \param[in] rank between 0.0 and 1.0; 0.5 is median filter |
1410 | | * \return pixd 1 bpp |
1411 | | * |
1412 | | * <pre> |
1413 | | * Notes: |
1414 | | * (1) The full width and height of the convolution kernel |
1415 | | * are (2 * wc + 1) and (2 * hc + 1) |
1416 | | * (2) This returns a pixd where each pixel is a 1 if the |
1417 | | * neighborhood (2 * wc + 1) x (2 * hc + 1)) pixels |
1418 | | * contains the rank fraction of 1 pixels. Otherwise, |
1419 | | * the returned pixel is 0. Note that the special case |
1420 | | * of rank = 0.0 is always satisfied, so the returned |
1421 | | * pixd has all pixels with value 1. |
1422 | | * (3) If accum pix is null, make one, use it, and destroy it |
1423 | | * before returning; otherwise, just use the input accum pix |
1424 | | * (4) If both wc and hc are 0, returns a copy unless rank == 0.0, |
1425 | | * in which case this returns an all-ones image. |
1426 | | * (5) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1, |
1427 | | * where (w,h) are the dimensions of pixs. Attempt to |
1428 | | * reduce the kernel size if necessary. |
1429 | | * </pre> |
1430 | | */ |
1431 | | PIX * |
1432 | | pixBlockrank(PIX *pixs, |
1433 | | PIX *pixacc, |
1434 | | l_int32 wc, |
1435 | | l_int32 hc, |
1436 | | l_float32 rank) |
1437 | 0 | { |
1438 | 0 | l_int32 w, h, d, thresh; |
1439 | 0 | PIX *pixt, *pixd; |
1440 | |
|
1441 | 0 | if (!pixs) |
1442 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
1443 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
1444 | 0 | if (d != 1) |
1445 | 0 | return (PIX *)ERROR_PTR("pixs not 1 bpp", __func__, NULL); |
1446 | 0 | if (rank < 0.0 || rank > 1.0) |
1447 | 0 | return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", __func__, NULL); |
1448 | | |
1449 | 0 | if (rank == 0.0) { |
1450 | 0 | pixd = pixCreateTemplate(pixs); |
1451 | 0 | pixSetAll(pixd); |
1452 | 0 | return pixd; |
1453 | 0 | } |
1454 | | |
1455 | 0 | if (wc <= 0 || hc <= 0) |
1456 | 0 | return pixCopy(NULL, pixs); |
1457 | 0 | if (w < 2 * wc + 1 || h < 2 * hc + 1) { |
1458 | 0 | L_WARNING("kernel too large: wc = %d, hc = %d, w = %d, h = %d; " |
1459 | 0 | "reducing!\n", __func__, wc, hc, w, h); |
1460 | 0 | wc = L_MIN(wc, (w - 1) / 2); |
1461 | 0 | hc = L_MIN(hc, (h - 1) / 2); |
1462 | 0 | } |
1463 | 0 | if (wc == 0 || hc == 0) |
1464 | 0 | return pixCopy(NULL, pixs); |
1465 | | |
1466 | 0 | if ((pixt = pixBlocksum(pixs, pixacc, wc, hc)) == NULL) |
1467 | 0 | return (PIX *)ERROR_PTR("pixt not made", __func__, NULL); |
1468 | | |
1469 | | /* 1 bpp block rank filter output. |
1470 | | * Must invert because threshold gives 1 for values < thresh, |
1471 | | * but we need a 1 if the value is >= thresh. */ |
1472 | 0 | thresh = (l_int32)(255. * rank); |
1473 | 0 | pixd = pixThresholdToBinary(pixt, thresh); |
1474 | 0 | pixInvert(pixd, pixd); |
1475 | 0 | pixDestroy(&pixt); |
1476 | 0 | return pixd; |
1477 | 0 | } |
1478 | | |
1479 | | |
1480 | | /*! |
1481 | | * \brief pixBlocksum() |
1482 | | * |
1483 | | * \param[in] pixs 1 bpp |
1484 | | * \param[in] pixacc pix [optional] 32 bpp |
1485 | | * \param[in] wc, hc half width/height of block sum/rank kernel |
1486 | | * \return pixd 8 bpp |
1487 | | * |
1488 | | * <pre> |
1489 | | * Notes: |
1490 | | * (1) If accum pix is null, make one and destroy it before |
1491 | | * returning; otherwise, just use the input accum pix |
1492 | | * (2) The full width and height of the convolution kernel |
1493 | | * are (2 * wc + 1) and (2 * hc + 1) |
1494 | | * (3) Use of wc = hc = 1, followed by pixInvert() on the |
1495 | | * 8 bpp result, gives a nice anti-aliased, and somewhat |
1496 | | * darkened, result on text. |
1497 | | * (4) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1, |
1498 | | * where (w,h) are the dimensions of pixs. Attempt to |
1499 | | * reduce the kernel size if necessary. |
1500 | | * (5) Returns in each dest pixel the sum of all src pixels |
1501 | | * that are within a block of size of the kernel, centered |
1502 | | * on the dest pixel. This sum is the number of src ON |
1503 | | * pixels in the block at each location, normalized to 255 |
1504 | | * for a block containing all ON pixels. For pixels near |
1505 | | * the boundary, where the block is not entirely contained |
1506 | | * within the image, we then multiply by a second normalization |
1507 | | * factor that is greater than one, so that all results |
1508 | | * are normalized by the number of participating pixels |
1509 | | * within the block. |
1510 | | * </pre> |
1511 | | */ |
1512 | | PIX * |
1513 | | pixBlocksum(PIX *pixs, |
1514 | | PIX *pixacc, |
1515 | | l_int32 wc, |
1516 | | l_int32 hc) |
1517 | 0 | { |
1518 | 0 | l_int32 w, h, d, wplt, wpld; |
1519 | 0 | l_uint32 *datat, *datad; |
1520 | 0 | PIX *pixt, *pixd; |
1521 | |
|
1522 | 0 | if (!pixs) |
1523 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
1524 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
1525 | 0 | if (d != 1) |
1526 | 0 | return (PIX *)ERROR_PTR("pixs not 1 bpp", __func__, NULL); |
1527 | 0 | if (wc <= 0 || hc <= 0) |
1528 | 0 | return pixCopy(NULL, pixs); |
1529 | 0 | if (w < 2 * wc + 1 || h < 2 * hc + 1) { |
1530 | 0 | L_WARNING("kernel too large: wc = %d, hc = %d, w = %d, h = %d; " |
1531 | 0 | "reducing!\n", __func__, wc, hc, w, h); |
1532 | 0 | wc = L_MIN(wc, (w - 1) / 2); |
1533 | 0 | hc = L_MIN(hc, (h - 1) / 2); |
1534 | 0 | } |
1535 | 0 | if (wc == 0 || hc == 0) |
1536 | 0 | return pixCopy(NULL, pixs); |
1537 | | |
1538 | 0 | if (pixacc) { |
1539 | 0 | if (pixGetDepth(pixacc) != 32) |
1540 | 0 | return (PIX *)ERROR_PTR("pixacc not 32 bpp", __func__, NULL); |
1541 | 0 | pixt = pixClone(pixacc); |
1542 | 0 | } else { |
1543 | 0 | if ((pixt = pixBlockconvAccum(pixs)) == NULL) |
1544 | 0 | return (PIX *)ERROR_PTR("pixt not made", __func__, NULL); |
1545 | 0 | } |
1546 | | |
1547 | | /* 8 bpp block sum output */ |
1548 | 0 | if ((pixd = pixCreate(w, h, 8)) == NULL) { |
1549 | 0 | pixDestroy(&pixt); |
1550 | 0 | return (PIX *)ERROR_PTR("pixd not made", __func__, NULL); |
1551 | 0 | } |
1552 | 0 | pixCopyResolution(pixd, pixs); |
1553 | |
|
1554 | 0 | wpld = pixGetWpl(pixd); |
1555 | 0 | wplt = pixGetWpl(pixt); |
1556 | 0 | datad = pixGetData(pixd); |
1557 | 0 | datat = pixGetData(pixt); |
1558 | 0 | blocksumLow(datad, w, h, wpld, datat, wplt, wc, hc); |
1559 | |
|
1560 | 0 | pixDestroy(&pixt); |
1561 | 0 | return pixd; |
1562 | 0 | } |
1563 | | |
1564 | | |
1565 | | /*! |
1566 | | * \brief blocksumLow() |
1567 | | * |
1568 | | * \param[in] datad of 8 bpp dest |
1569 | | * \param[in] w, h, wpl of 8 bpp dest |
1570 | | * \param[in] dataa of 32 bpp accum |
1571 | | * \param[in] wpla of 32 bpp accum |
1572 | | * \param[in] wc, hc convolution "half-width" and "half-height" |
1573 | | * \return void |
1574 | | * |
1575 | | * <pre> |
1576 | | * Notes: |
1577 | | * (1) The full width and height of the convolution kernel |
1578 | | * are (2 * wc + 1) and (2 * hc + 1). |
1579 | | * (2) The lack of symmetry between the handling of the |
1580 | | * first (hc + 1) lines and the last (hc) lines, |
1581 | | * and similarly with the columns, is due to fact that |
1582 | | * for the pixel at (x,y), the accumulator values are |
1583 | | * taken at (x + wc, y + hc), (x - wc - 1, y + hc), |
1584 | | * (x + wc, y - hc - 1) and (x - wc - 1, y - hc - 1). |
1585 | | * (3) Compute sums of ON pixels within the block filter size, |
1586 | | * normalized between 0 and 255, as if there were no reduced |
1587 | | * area at the boundary. This under-estimates the value |
1588 | | * of the boundary pixels, so we multiply them by another |
1589 | | * normalization factor that is greater than 1. |
1590 | | * (4) This second normalization is done first for the first |
1591 | | * hc + 1 lines; then for the last hc lines; and finally |
1592 | | * for the first wc + 1 and last wc columns in the intermediate |
1593 | | * lines. |
1594 | | * (5) Required constraints are: wc < w and hc < h. |
1595 | | * </pre> |
1596 | | */ |
1597 | | static void |
1598 | | blocksumLow(l_uint32 *datad, |
1599 | | l_int32 w, |
1600 | | l_int32 h, |
1601 | | l_int32 wpl, |
1602 | | l_uint32 *dataa, |
1603 | | l_int32 wpla, |
1604 | | l_int32 wc, |
1605 | | l_int32 hc) |
1606 | 0 | { |
1607 | 0 | l_int32 i, j, imax, imin, jmax, jmin; |
1608 | 0 | l_int32 wn, hn, fwc, fhc, wmwc, hmhc; |
1609 | 0 | l_float32 norm, normh, normw; |
1610 | 0 | l_uint32 val; |
1611 | 0 | l_uint32 *linemina, *linemaxa, *lined; |
1612 | |
|
1613 | 0 | wmwc = w - wc; |
1614 | 0 | hmhc = h - hc; |
1615 | 0 | if (wmwc <= 0 || hmhc <= 0) { |
1616 | 0 | L_ERROR("wc >= w || hc >=h\n", __func__); |
1617 | 0 | return; |
1618 | 0 | } |
1619 | 0 | fwc = 2 * wc + 1; |
1620 | 0 | fhc = 2 * hc + 1; |
1621 | 0 | norm = 255. / ((l_float32)(fwc) * fhc); |
1622 | | |
1623 | | /*------------------------------------------------------------* |
1624 | | * Compute, using b.c. only to set limits on the accum image * |
1625 | | *------------------------------------------------------------*/ |
1626 | 0 | for (i = 0; i < h; i++) { |
1627 | 0 | imin = L_MAX(i - 1 - hc, 0); |
1628 | 0 | imax = L_MIN(i + hc, h - 1); |
1629 | 0 | lined = datad + wpl * i; |
1630 | 0 | linemina = dataa + wpla * imin; |
1631 | 0 | linemaxa = dataa + wpla * imax; |
1632 | 0 | for (j = 0; j < w; j++) { |
1633 | 0 | jmin = L_MAX(j - 1 - wc, 0); |
1634 | 0 | jmax = L_MIN(j + wc, w - 1); |
1635 | 0 | val = linemaxa[jmax] - linemaxa[jmin] |
1636 | 0 | - linemina[jmax] + linemina[jmin]; |
1637 | 0 | val = (l_uint8)(norm * val); |
1638 | 0 | SET_DATA_BYTE(lined, j, val); |
1639 | 0 | } |
1640 | 0 | } |
1641 | | |
1642 | | /*------------------------------------------------------------* |
1643 | | * Fix normalization for boundary pixels * |
1644 | | *------------------------------------------------------------*/ |
1645 | 0 | for (i = 0; i <= hc; i++) { /* first hc + 1 lines */ |
1646 | 0 | hn = hc + i; |
1647 | 0 | normh = (l_float32)fhc / (l_float32)hn; /* > 1 */ |
1648 | 0 | lined = datad + wpl * i; |
1649 | 0 | for (j = 0; j <= wc; j++) { |
1650 | 0 | wn = wc + j; |
1651 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
1652 | 0 | val = GET_DATA_BYTE(lined, j); |
1653 | 0 | val = (l_uint8)(val * normh * normw); |
1654 | 0 | SET_DATA_BYTE(lined, j, val); |
1655 | 0 | } |
1656 | 0 | for (j = wc + 1; j < wmwc; j++) { |
1657 | 0 | val = GET_DATA_BYTE(lined, j); |
1658 | 0 | val = (l_uint8)(val * normh); |
1659 | 0 | SET_DATA_BYTE(lined, j, val); |
1660 | 0 | } |
1661 | 0 | for (j = wmwc; j < w; j++) { |
1662 | 0 | wn = wc + w - j; |
1663 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
1664 | 0 | val = GET_DATA_BYTE(lined, j); |
1665 | 0 | val = (l_uint8)(val * normh * normw); |
1666 | 0 | SET_DATA_BYTE(lined, j, val); |
1667 | 0 | } |
1668 | 0 | } |
1669 | |
|
1670 | 0 | for (i = hmhc; i < h; i++) { /* last hc lines */ |
1671 | 0 | hn = hc + h - i; |
1672 | 0 | normh = (l_float32)fhc / (l_float32)hn; /* > 1 */ |
1673 | 0 | lined = datad + wpl * i; |
1674 | 0 | for (j = 0; j <= wc; j++) { |
1675 | 0 | wn = wc + j; |
1676 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
1677 | 0 | val = GET_DATA_BYTE(lined, j); |
1678 | 0 | val = (l_uint8)(val * normh * normw); |
1679 | 0 | SET_DATA_BYTE(lined, j, val); |
1680 | 0 | } |
1681 | 0 | for (j = wc + 1; j < wmwc; j++) { |
1682 | 0 | val = GET_DATA_BYTE(lined, j); |
1683 | 0 | val = (l_uint8)(val * normh); |
1684 | 0 | SET_DATA_BYTE(lined, j, val); |
1685 | 0 | } |
1686 | 0 | for (j = wmwc; j < w; j++) { |
1687 | 0 | wn = wc + w - j; |
1688 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
1689 | 0 | val = GET_DATA_BYTE(lined, j); |
1690 | 0 | val = (l_uint8)(val * normh * normw); |
1691 | 0 | SET_DATA_BYTE(lined, j, val); |
1692 | 0 | } |
1693 | 0 | } |
1694 | |
|
1695 | 0 | for (i = hc + 1; i < hmhc; i++) { /* intermediate lines */ |
1696 | 0 | lined = datad + wpl * i; |
1697 | 0 | for (j = 0; j <= wc; j++) { /* first wc + 1 columns */ |
1698 | 0 | wn = wc + j; |
1699 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
1700 | 0 | val = GET_DATA_BYTE(lined, j); |
1701 | 0 | val = (l_uint8)(val * normw); |
1702 | 0 | SET_DATA_BYTE(lined, j, val); |
1703 | 0 | } |
1704 | 0 | for (j = wmwc; j < w; j++) { /* last wc columns */ |
1705 | 0 | wn = wc + w - j; |
1706 | 0 | normw = (l_float32)fwc / (l_float32)wn; /* > 1 */ |
1707 | 0 | val = GET_DATA_BYTE(lined, j); |
1708 | 0 | val = (l_uint8)(val * normw); |
1709 | 0 | SET_DATA_BYTE(lined, j, val); |
1710 | 0 | } |
1711 | 0 | } |
1712 | 0 | } |
1713 | | |
1714 | | |
1715 | | /*----------------------------------------------------------------------* |
1716 | | * Census transform * |
1717 | | *----------------------------------------------------------------------*/ |
1718 | | /*! |
1719 | | * \brief pixCensusTransform() |
1720 | | * |
1721 | | * \param[in] pixs 8 bpp |
1722 | | * \param[in] halfsize of square over which neighbors are averaged |
1723 | | * \param[in] pixacc [optional] 32 bpp pix |
1724 | | * \return pixd 1 bpp |
1725 | | * |
1726 | | * <pre> |
1727 | | * Notes: |
1728 | | * (1) The Census transform was invented by Ramin Zabih and John Woodfill |
1729 | | * ("Non-parametric local transforms for computing visual |
1730 | | * correspondence", Third European Conference on Computer Vision, |
1731 | | * Stockholm, Sweden, May 1994); see publications at |
1732 | | * http://www.cs.cornell.edu/~rdz/index.htm |
1733 | | * This compares each pixel against the average of its neighbors, |
1734 | | * in a square of odd dimension centered on the pixel. |
1735 | | * If the pixel is greater than the average of its neighbors, |
1736 | | * the output pixel value is 1; otherwise it is 0. |
1737 | | * (2) This can be used as an encoding for an image that is |
1738 | | * fairly robust against slow illumination changes, with |
1739 | | * applications in image comparison and mosaicing. |
1740 | | * (3) The size of the convolution kernel is (2 * halfsize + 1) |
1741 | | * on a side. The halfsize parameter must be >= 1. |
1742 | | * (4) If accum pix is null, make one, use it, and destroy it |
1743 | | * before returning; otherwise, just use the input accum pix |
1744 | | * </pre> |
1745 | | */ |
1746 | | PIX * |
1747 | | pixCensusTransform(PIX *pixs, |
1748 | | l_int32 halfsize, |
1749 | | PIX *pixacc) |
1750 | 0 | { |
1751 | 0 | l_int32 i, j, w, h, wpls, wplv, wpld; |
1752 | 0 | l_int32 vals, valv; |
1753 | 0 | l_uint32 *datas, *datav, *datad, *lines, *linev, *lined; |
1754 | 0 | PIX *pixav, *pixd; |
1755 | |
|
1756 | 0 | if (!pixs) |
1757 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
1758 | 0 | if (pixGetDepth(pixs) != 8) |
1759 | 0 | return (PIX *)ERROR_PTR("pixs not 8 bpp", __func__, NULL); |
1760 | 0 | if (halfsize < 1) |
1761 | 0 | return (PIX *)ERROR_PTR("halfsize must be >= 1", __func__, NULL); |
1762 | | |
1763 | | /* Get the average of each pixel with its neighbors */ |
1764 | 0 | if ((pixav = pixBlockconvGray(pixs, pixacc, halfsize, halfsize)) |
1765 | 0 | == NULL) |
1766 | 0 | return (PIX *)ERROR_PTR("pixav not made", __func__, NULL); |
1767 | | |
1768 | | /* Subtract the pixel from the average, and then compare |
1769 | | * the pixel value with the remaining average */ |
1770 | 0 | pixGetDimensions(pixs, &w, &h, NULL); |
1771 | 0 | if ((pixd = pixCreate(w, h, 1)) == NULL) { |
1772 | 0 | pixDestroy(&pixav); |
1773 | 0 | return (PIX *)ERROR_PTR("pixd not made", __func__, NULL); |
1774 | 0 | } |
1775 | 0 | datas = pixGetData(pixs); |
1776 | 0 | datav = pixGetData(pixav); |
1777 | 0 | datad = pixGetData(pixd); |
1778 | 0 | wpls = pixGetWpl(pixs); |
1779 | 0 | wplv = pixGetWpl(pixav); |
1780 | 0 | wpld = pixGetWpl(pixd); |
1781 | 0 | for (i = 0; i < h; i++) { |
1782 | 0 | lines = datas + i * wpls; |
1783 | 0 | linev = datav + i * wplv; |
1784 | 0 | lined = datad + i * wpld; |
1785 | 0 | for (j = 0; j < w; j++) { |
1786 | 0 | vals = GET_DATA_BYTE(lines, j); |
1787 | 0 | valv = GET_DATA_BYTE(linev, j); |
1788 | 0 | if (vals > valv) |
1789 | 0 | SET_DATA_BIT(lined, j); |
1790 | 0 | } |
1791 | 0 | } |
1792 | |
|
1793 | 0 | pixDestroy(&pixav); |
1794 | 0 | return pixd; |
1795 | 0 | } |
1796 | | |
1797 | | |
1798 | | /*----------------------------------------------------------------------* |
1799 | | * Generic convolution * |
1800 | | *----------------------------------------------------------------------*/ |
1801 | | /*! |
1802 | | * \brief pixConvolve() |
1803 | | * |
1804 | | * \param[in] pixs 8, 16, 32 bpp; no colormap |
1805 | | * \param[in] kel kernel |
1806 | | * \param[in] outdepth of pixd: 8, 16 or 32 |
1807 | | * \param[in] normflag 1 to normalize kernel to unit sum; 0 otherwise |
1808 | | * \return pixd 8, 16 or 32 bpp |
1809 | | * |
1810 | | * <pre> |
1811 | | * Notes: |
1812 | | * (1) This gives a convolution with an arbitrary kernel. |
1813 | | * (2) The input pixs must have only one sample/pixel. |
1814 | | * To do a convolution on an RGB image, use pixConvolveRGB(). |
1815 | | * (3) The parameter %outdepth determines the depth of the result. |
1816 | | * If the kernel is normalized to unit sum, the output values |
1817 | | * can never exceed 255, so an output depth of 8 bpp is sufficient. |
1818 | | * If the kernel is not normalized, it may be necessary to use |
1819 | | * 16 or 32 bpp output to avoid overflow. |
1820 | | * (4) If normflag == 1, the result is normalized by scaling all |
1821 | | * kernel values for a unit sum. If the sum of kernel values |
1822 | | * is very close to zero, the kernel can not be normalized and |
1823 | | * the convolution will not be performed. A warning is issued. |
1824 | | * (5) The kernel values can be positive or negative, but the |
1825 | | * result for the convolution can only be stored as a positive |
1826 | | * number. Consequently, if it goes negative, the choices are |
1827 | | * to clip to 0 or take the absolute value. We're choosing |
1828 | | * to take the absolute value. (Another possibility would be |
1829 | | * to output a second unsigned image for the negative values.) |
1830 | | * If you want to get a clipped result, or to keep the negative |
1831 | | * values in the result, use fpixConvolve(), with the |
1832 | | * converters in fpix2.c between pix and fpix. |
1833 | | * (6) This uses a mirrored border to avoid special casing on |
1834 | | * the boundaries. |
1835 | | * (7) To get a subsampled output, call l_setConvolveSampling(). |
1836 | | * The time to make a subsampled output is reduced by the |
1837 | | * product of the sampling factors. |
1838 | | * (8) The function is slow, running at about 12 machine cycles for |
1839 | | * each pixel-op in the convolution. For example, with a 3 GHz |
1840 | | * cpu, a 1 Mpixel grayscale image, and a kernel with |
1841 | | * (sx * sy) = 25 elements, the convolution takes about 100 msec. |
1842 | | * </pre> |
1843 | | */ |
1844 | | PIX * |
1845 | | pixConvolve(PIX *pixs, |
1846 | | L_KERNEL *kel, |
1847 | | l_int32 outdepth, |
1848 | | l_int32 normflag) |
1849 | 0 | { |
1850 | 0 | l_int32 i, j, id, jd, k, m, w, h, d, wd, hd, sx, sy, cx, cy, wplt, wpld; |
1851 | 0 | l_int32 val; |
1852 | 0 | l_uint32 *datat, *datad, *linet, *lined; |
1853 | 0 | l_float32 sum; |
1854 | 0 | L_KERNEL *keli, *keln; |
1855 | 0 | PIX *pixt, *pixd; |
1856 | |
|
1857 | 0 | if (!pixs) |
1858 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
1859 | 0 | if (pixGetColormap(pixs)) |
1860 | 0 | return (PIX *)ERROR_PTR("pixs has colormap", __func__, NULL); |
1861 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
1862 | 0 | if (d != 8 && d != 16 && d != 32) |
1863 | 0 | return (PIX *)ERROR_PTR("pixs not 8, 16, or 32 bpp", __func__, NULL); |
1864 | 0 | if (!kel) |
1865 | 0 | return (PIX *)ERROR_PTR("kel not defined", __func__, NULL); |
1866 | | |
1867 | 0 | pixd = NULL; |
1868 | |
|
1869 | 0 | keli = kernelInvert(kel); |
1870 | 0 | kernelGetParameters(keli, &sy, &sx, &cy, &cx); |
1871 | 0 | if (normflag) |
1872 | 0 | keln = kernelNormalize(keli, 1.0); |
1873 | 0 | else |
1874 | 0 | keln = kernelCopy(keli); |
1875 | |
|
1876 | 0 | if ((pixt = pixAddMirroredBorder(pixs, cx, sx - cx, cy, sy - cy)) == NULL) { |
1877 | 0 | L_ERROR("pixt not made\n", __func__); |
1878 | 0 | goto cleanup; |
1879 | 0 | } |
1880 | | |
1881 | 0 | wd = (w + ConvolveSamplingFactX - 1) / ConvolveSamplingFactX; |
1882 | 0 | hd = (h + ConvolveSamplingFactY - 1) / ConvolveSamplingFactY; |
1883 | 0 | pixd = pixCreate(wd, hd, outdepth); |
1884 | 0 | datat = pixGetData(pixt); |
1885 | 0 | datad = pixGetData(pixd); |
1886 | 0 | wplt = pixGetWpl(pixt); |
1887 | 0 | wpld = pixGetWpl(pixd); |
1888 | 0 | for (i = 0, id = 0; id < hd; i += ConvolveSamplingFactY, id++) { |
1889 | 0 | lined = datad + id * wpld; |
1890 | 0 | for (j = 0, jd = 0; jd < wd; j += ConvolveSamplingFactX, jd++) { |
1891 | 0 | sum = 0.0; |
1892 | 0 | for (k = 0; k < sy; k++) { |
1893 | 0 | linet = datat + (i + k) * wplt; |
1894 | 0 | if (d == 8) { |
1895 | 0 | for (m = 0; m < sx; m++) { |
1896 | 0 | val = GET_DATA_BYTE(linet, j + m); |
1897 | 0 | sum += val * keln->data[k][m]; |
1898 | 0 | } |
1899 | 0 | } else if (d == 16) { |
1900 | 0 | for (m = 0; m < sx; m++) { |
1901 | 0 | val = GET_DATA_TWO_BYTES(linet, j + m); |
1902 | 0 | sum += val * keln->data[k][m]; |
1903 | 0 | } |
1904 | 0 | } else { /* d == 32 */ |
1905 | 0 | for (m = 0; m < sx; m++) { |
1906 | 0 | val = *(linet + j + m); |
1907 | 0 | sum += val * keln->data[k][m]; |
1908 | 0 | } |
1909 | 0 | } |
1910 | 0 | } |
1911 | 0 | if (sum < 0.0) sum = -sum; /* make it non-negative */ |
1912 | 0 | if (outdepth == 8) |
1913 | 0 | SET_DATA_BYTE(lined, jd, (l_int32)(sum + 0.5)); |
1914 | 0 | else if (outdepth == 16) |
1915 | 0 | SET_DATA_TWO_BYTES(lined, jd, (l_int32)(sum + 0.5)); |
1916 | 0 | else /* outdepth == 32 */ |
1917 | 0 | *(lined + jd) = (l_uint32)(sum + 0.5); |
1918 | 0 | } |
1919 | 0 | } |
1920 | |
|
1921 | 0 | cleanup: |
1922 | 0 | kernelDestroy(&keli); |
1923 | 0 | kernelDestroy(&keln); |
1924 | 0 | pixDestroy(&pixt); |
1925 | 0 | return pixd; |
1926 | 0 | } |
1927 | | |
1928 | | |
1929 | | /*! |
1930 | | * \brief pixConvolveSep() |
1931 | | * |
1932 | | * \param[in] pixs 8, 16, 32 bpp; no colormap |
1933 | | * \param[in] kelx x-dependent kernel |
1934 | | * \param[in] kely y-dependent kernel |
1935 | | * \param[in] outdepth of pixd: 8, 16 or 32 |
1936 | | * \param[in] normflag 1 to normalize kernel to unit sum; 0 otherwise |
1937 | | * \return pixd 8, 16 or 32 bpp |
1938 | | * |
1939 | | * <pre> |
1940 | | * Notes: |
1941 | | * (1) This does a convolution with a separable kernel that is |
1942 | | * is a sequence of convolutions in x and y. The two |
1943 | | * one-dimensional kernel components must be input separately; |
1944 | | * the full kernel is the product of these components. |
1945 | | * The support for the full kernel is thus a rectangular region. |
1946 | | * (2) The input pixs must have only one sample/pixel. |
1947 | | * To do a convolution on an RGB image, use pixConvolveSepRGB(). |
1948 | | * (3) The parameter %outdepth determines the depth of the result. |
1949 | | * If the kernel is normalized to unit sum, the output values |
1950 | | * can never exceed 255, so an output depth of 8 bpp is sufficient. |
1951 | | * If the kernel is not normalized, it may be necessary to use |
1952 | | * 16 or 32 bpp output to avoid overflow. |
1953 | | * (2) The %normflag parameter is used as in pixConvolve(). |
1954 | | * (4) The kernel values can be positive or negative, but the |
1955 | | * result for the convolution can only be stored as a positive |
1956 | | * number. Consequently, if it goes negative, the choices are |
1957 | | * to clip to 0 or take the absolute value. We're choosing |
1958 | | * the former for now. Another possibility would be to output |
1959 | | * a second unsigned image for the negative values. |
1960 | | * (5) Warning: if you use l_setConvolveSampling() to get a |
1961 | | * subsampled output, and the sampling factor is larger than |
1962 | | * the kernel half-width, it is faster to use the non-separable |
1963 | | * version pixConvolve(). This is because the first convolution |
1964 | | * here must be done on every raster line, regardless of the |
1965 | | * vertical sampling factor. If the sampling factor is smaller |
1966 | | * than kernel half-width, it's faster to use the separable |
1967 | | * convolution. |
1968 | | * (6) This uses mirrored borders to avoid special casing on |
1969 | | * the boundaries. |
1970 | | * </pre> |
1971 | | */ |
1972 | | PIX * |
1973 | | pixConvolveSep(PIX *pixs, |
1974 | | L_KERNEL *kelx, |
1975 | | L_KERNEL *kely, |
1976 | | l_int32 outdepth, |
1977 | | l_int32 normflag) |
1978 | 0 | { |
1979 | 0 | l_int32 d, xfact, yfact; |
1980 | 0 | L_KERNEL *kelxn, *kelyn; |
1981 | 0 | PIX *pixt, *pixd; |
1982 | |
|
1983 | 0 | if (!pixs) |
1984 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
1985 | 0 | d = pixGetDepth(pixs); |
1986 | 0 | if (d != 8 && d != 16 && d != 32) |
1987 | 0 | return (PIX *)ERROR_PTR("pixs not 8, 16, or 32 bpp", __func__, NULL); |
1988 | 0 | if (!kelx) |
1989 | 0 | return (PIX *)ERROR_PTR("kelx not defined", __func__, NULL); |
1990 | 0 | if (!kely) |
1991 | 0 | return (PIX *)ERROR_PTR("kely not defined", __func__, NULL); |
1992 | | |
1993 | 0 | xfact = ConvolveSamplingFactX; |
1994 | 0 | yfact = ConvolveSamplingFactY; |
1995 | 0 | if (normflag) { |
1996 | 0 | kelxn = kernelNormalize(kelx, 1000.0f); |
1997 | 0 | kelyn = kernelNormalize(kely, 0.001f); |
1998 | 0 | l_setConvolveSampling(xfact, 1); |
1999 | 0 | pixt = pixConvolve(pixs, kelxn, 32, 0); |
2000 | 0 | l_setConvolveSampling(1, yfact); |
2001 | 0 | pixd = pixConvolve(pixt, kelyn, outdepth, 0); |
2002 | 0 | l_setConvolveSampling(xfact, yfact); /* restore */ |
2003 | 0 | kernelDestroy(&kelxn); |
2004 | 0 | kernelDestroy(&kelyn); |
2005 | 0 | } else { /* don't normalize */ |
2006 | 0 | l_setConvolveSampling(xfact, 1); |
2007 | 0 | pixt = pixConvolve(pixs, kelx, 32, 0); |
2008 | 0 | l_setConvolveSampling(1, yfact); |
2009 | 0 | pixd = pixConvolve(pixt, kely, outdepth, 0); |
2010 | 0 | l_setConvolveSampling(xfact, yfact); |
2011 | 0 | } |
2012 | |
|
2013 | 0 | pixDestroy(&pixt); |
2014 | 0 | return pixd; |
2015 | 0 | } |
2016 | | |
2017 | | |
2018 | | /*! |
2019 | | * \brief pixConvolveRGB() |
2020 | | * |
2021 | | * \param[in] pixs 32 bpp rgb |
2022 | | * \param[in] kel kernel |
2023 | | * \return pixd 32 bpp rgb |
2024 | | * |
2025 | | * <pre> |
2026 | | * Notes: |
2027 | | * (1) This gives a convolution on an RGB image using an |
2028 | | * arbitrary kernel (which we normalize to keep each |
2029 | | * component within the range [0 ... 255]. |
2030 | | * (2) The input pixs must be RGB. |
2031 | | * (3) The kernel values can be positive or negative, but the |
2032 | | * result for the convolution can only be stored as a positive |
2033 | | * number. Consequently, if it goes negative, we clip the |
2034 | | * result to 0. |
2035 | | * (4) To get a subsampled output, call l_setConvolveSampling(). |
2036 | | * The time to make a subsampled output is reduced by the |
2037 | | * product of the sampling factors. |
2038 | | * (5) This uses a mirrored border to avoid special casing on |
2039 | | * the boundaries. |
2040 | | * </pre> |
2041 | | */ |
2042 | | PIX * |
2043 | | pixConvolveRGB(PIX *pixs, |
2044 | | L_KERNEL *kel) |
2045 | 0 | { |
2046 | 0 | PIX *pixt, *pixr, *pixg, *pixb, *pixd; |
2047 | |
|
2048 | 0 | if (!pixs) |
2049 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
2050 | 0 | if (pixGetDepth(pixs) != 32) |
2051 | 0 | return (PIX *)ERROR_PTR("pixs is not 32 bpp", __func__, NULL); |
2052 | 0 | if (!kel) |
2053 | 0 | return (PIX *)ERROR_PTR("kel not defined", __func__, NULL); |
2054 | | |
2055 | 0 | pixt = pixGetRGBComponent(pixs, COLOR_RED); |
2056 | 0 | pixr = pixConvolve(pixt, kel, 8, 1); |
2057 | 0 | pixDestroy(&pixt); |
2058 | 0 | pixt = pixGetRGBComponent(pixs, COLOR_GREEN); |
2059 | 0 | pixg = pixConvolve(pixt, kel, 8, 1); |
2060 | 0 | pixDestroy(&pixt); |
2061 | 0 | pixt = pixGetRGBComponent(pixs, COLOR_BLUE); |
2062 | 0 | pixb = pixConvolve(pixt, kel, 8, 1); |
2063 | 0 | pixDestroy(&pixt); |
2064 | 0 | pixd = pixCreateRGBImage(pixr, pixg, pixb); |
2065 | |
|
2066 | 0 | pixDestroy(&pixr); |
2067 | 0 | pixDestroy(&pixg); |
2068 | 0 | pixDestroy(&pixb); |
2069 | 0 | return pixd; |
2070 | 0 | } |
2071 | | |
2072 | | |
2073 | | /*! |
2074 | | * \brief pixConvolveRGBSep() |
2075 | | * |
2076 | | * \param[in] pixs 32 bpp rgb |
2077 | | * \param[in] kelx x-dependent kernel |
2078 | | * \param[in] kely y-dependent kernel |
2079 | | * \return pixd 32 bpp rgb |
2080 | | * |
2081 | | * <pre> |
2082 | | * Notes: |
2083 | | * (1) This does a convolution on an RGB image using a separable |
2084 | | * kernel that is a sequence of convolutions in x and y. The two |
2085 | | * one-dimensional kernel components must be input separately; |
2086 | | * the full kernel is the product of these components. |
2087 | | * The support for the full kernel is thus a rectangular region. |
2088 | | * (2) The kernel values can be positive or negative, but the |
2089 | | * result for the convolution can only be stored as a positive |
2090 | | * number. Consequently, if it goes negative, we clip the |
2091 | | * result to 0. |
2092 | | * (3) To get a subsampled output, call l_setConvolveSampling(). |
2093 | | * The time to make a subsampled output is reduced by the |
2094 | | * product of the sampling factors. |
2095 | | * (4) This uses a mirrored border to avoid special casing on |
2096 | | * the boundaries. |
2097 | | * </pre> |
2098 | | */ |
2099 | | PIX * |
2100 | | pixConvolveRGBSep(PIX *pixs, |
2101 | | L_KERNEL *kelx, |
2102 | | L_KERNEL *kely) |
2103 | 0 | { |
2104 | 0 | PIX *pixt, *pixr, *pixg, *pixb, *pixd; |
2105 | |
|
2106 | 0 | if (!pixs) |
2107 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
2108 | 0 | if (pixGetDepth(pixs) != 32) |
2109 | 0 | return (PIX *)ERROR_PTR("pixs is not 32 bpp", __func__, NULL); |
2110 | 0 | if (!kelx || !kely) |
2111 | 0 | return (PIX *)ERROR_PTR("kelx, kely not both defined", __func__, NULL); |
2112 | | |
2113 | 0 | pixt = pixGetRGBComponent(pixs, COLOR_RED); |
2114 | 0 | pixr = pixConvolveSep(pixt, kelx, kely, 8, 1); |
2115 | 0 | pixDestroy(&pixt); |
2116 | 0 | pixt = pixGetRGBComponent(pixs, COLOR_GREEN); |
2117 | 0 | pixg = pixConvolveSep(pixt, kelx, kely, 8, 1); |
2118 | 0 | pixDestroy(&pixt); |
2119 | 0 | pixt = pixGetRGBComponent(pixs, COLOR_BLUE); |
2120 | 0 | pixb = pixConvolveSep(pixt, kelx, kely, 8, 1); |
2121 | 0 | pixDestroy(&pixt); |
2122 | 0 | pixd = pixCreateRGBImage(pixr, pixg, pixb); |
2123 | |
|
2124 | 0 | pixDestroy(&pixr); |
2125 | 0 | pixDestroy(&pixg); |
2126 | 0 | pixDestroy(&pixb); |
2127 | 0 | return pixd; |
2128 | 0 | } |
2129 | | |
2130 | | |
2131 | | /*----------------------------------------------------------------------* |
2132 | | * Generic convolution with float array * |
2133 | | *----------------------------------------------------------------------*/ |
2134 | | /*! |
2135 | | * \brief fpixConvolve() |
2136 | | * |
2137 | | * \param[in] fpixs 32 bit float array |
2138 | | * \param[in] kel kernel |
2139 | | * \param[in] normflag 1 to normalize kernel to unit sum; 0 otherwise |
2140 | | * \return fpixd 32 bit float array |
2141 | | * |
2142 | | * <pre> |
2143 | | * Notes: |
2144 | | * (1) This gives a float convolution with an arbitrary kernel. |
2145 | | * (2) If normflag == 1, the result is normalized by scaling all |
2146 | | * kernel values for a unit sum. If the sum of kernel values |
2147 | | * is very close to zero, the kernel can not be normalized and |
2148 | | * the convolution will not be performed. A warning is issued. |
2149 | | * (3) With the FPix, there are no issues about negative |
2150 | | * array or kernel values. The convolution is performed |
2151 | | * with single precision arithmetic. |
2152 | | * (4) To get a subsampled output, call l_setConvolveSampling(). |
2153 | | * The time to make a subsampled output is reduced by the |
2154 | | * product of the sampling factors. |
2155 | | * (5) This uses a mirrored border to avoid special casing on |
2156 | | * the boundaries. |
2157 | | * </pre> |
2158 | | */ |
2159 | | FPIX * |
2160 | | fpixConvolve(FPIX *fpixs, |
2161 | | L_KERNEL *kel, |
2162 | | l_int32 normflag) |
2163 | 0 | { |
2164 | 0 | l_int32 i, j, id, jd, k, m, w, h, wd, hd, sx, sy, cx, cy, wplt, wpld; |
2165 | 0 | l_float32 val; |
2166 | 0 | l_float32 *datat, *datad, *linet, *lined; |
2167 | 0 | l_float32 sum; |
2168 | 0 | L_KERNEL *keli, *keln; |
2169 | 0 | FPIX *fpixt, *fpixd; |
2170 | |
|
2171 | 0 | if (!fpixs) |
2172 | 0 | return (FPIX *)ERROR_PTR("fpixs not defined", __func__, NULL); |
2173 | 0 | if (!kel) |
2174 | 0 | return (FPIX *)ERROR_PTR("kel not defined", __func__, NULL); |
2175 | | |
2176 | 0 | fpixd = NULL; |
2177 | |
|
2178 | 0 | keli = kernelInvert(kel); |
2179 | 0 | kernelGetParameters(keli, &sy, &sx, &cy, &cx); |
2180 | 0 | if (normflag) |
2181 | 0 | keln = kernelNormalize(keli, 1.0); |
2182 | 0 | else |
2183 | 0 | keln = kernelCopy(keli); |
2184 | |
|
2185 | 0 | fpixGetDimensions(fpixs, &w, &h); |
2186 | 0 | fpixt = fpixAddMirroredBorder(fpixs, cx, sx - cx, cy, sy - cy); |
2187 | 0 | if (!fpixt) { |
2188 | 0 | L_ERROR("fpixt not made\n", __func__); |
2189 | 0 | goto cleanup; |
2190 | 0 | } |
2191 | | |
2192 | 0 | wd = (w + ConvolveSamplingFactX - 1) / ConvolveSamplingFactX; |
2193 | 0 | hd = (h + ConvolveSamplingFactY - 1) / ConvolveSamplingFactY; |
2194 | 0 | fpixd = fpixCreate(wd, hd); |
2195 | 0 | datat = fpixGetData(fpixt); |
2196 | 0 | datad = fpixGetData(fpixd); |
2197 | 0 | wplt = fpixGetWpl(fpixt); |
2198 | 0 | wpld = fpixGetWpl(fpixd); |
2199 | 0 | for (i = 0, id = 0; id < hd; i += ConvolveSamplingFactY, id++) { |
2200 | 0 | lined = datad + id * wpld; |
2201 | 0 | for (j = 0, jd = 0; jd < wd; j += ConvolveSamplingFactX, jd++) { |
2202 | 0 | sum = 0.0; |
2203 | 0 | for (k = 0; k < sy; k++) { |
2204 | 0 | linet = datat + (i + k) * wplt; |
2205 | 0 | for (m = 0; m < sx; m++) { |
2206 | 0 | val = *(linet + j + m); |
2207 | 0 | sum += val * keln->data[k][m]; |
2208 | 0 | } |
2209 | 0 | } |
2210 | 0 | *(lined + jd) = sum; |
2211 | 0 | } |
2212 | 0 | } |
2213 | |
|
2214 | 0 | cleanup: |
2215 | 0 | kernelDestroy(&keli); |
2216 | 0 | kernelDestroy(&keln); |
2217 | 0 | fpixDestroy(&fpixt); |
2218 | 0 | return fpixd; |
2219 | 0 | } |
2220 | | |
2221 | | |
2222 | | /*! |
2223 | | * \brief fpixConvolveSep() |
2224 | | * |
2225 | | * \param[in] fpixs 32 bit float array |
2226 | | * \param[in] kelx x-dependent kernel |
2227 | | * \param[in] kely y-dependent kernel |
2228 | | * \param[in] normflag 1 to normalize kernel to unit sum; 0 otherwise |
2229 | | * \return fpixd 32 bit float array |
2230 | | * |
2231 | | * <pre> |
2232 | | * Notes: |
2233 | | * (1) This does a convolution with a separable kernel that is |
2234 | | * is a sequence of convolutions in x and y. The two |
2235 | | * one-dimensional kernel components must be input separately; |
2236 | | * the full kernel is the product of these components. |
2237 | | * The support for the full kernel is thus a rectangular region. |
2238 | | * (2) The normflag parameter is used as in fpixConvolve(). |
2239 | | * (3) Warning: if you use l_setConvolveSampling() to get a |
2240 | | * subsampled output, and the sampling factor is larger than |
2241 | | * the kernel half-width, it is faster to use the non-separable |
2242 | | * version pixConvolve(). This is because the first convolution |
2243 | | * here must be done on every raster line, regardless of the |
2244 | | * vertical sampling factor. If the sampling factor is smaller |
2245 | | * than kernel half-width, it's faster to use the separable |
2246 | | * convolution. |
2247 | | * (4) This uses mirrored borders to avoid special casing on |
2248 | | * the boundaries. |
2249 | | * </pre> |
2250 | | */ |
2251 | | FPIX * |
2252 | | fpixConvolveSep(FPIX *fpixs, |
2253 | | L_KERNEL *kelx, |
2254 | | L_KERNEL *kely, |
2255 | | l_int32 normflag) |
2256 | 0 | { |
2257 | 0 | l_int32 xfact, yfact; |
2258 | 0 | L_KERNEL *kelxn, *kelyn; |
2259 | 0 | FPIX *fpixt, *fpixd; |
2260 | |
|
2261 | 0 | if (!fpixs) |
2262 | 0 | return (FPIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
2263 | 0 | if (!kelx) |
2264 | 0 | return (FPIX *)ERROR_PTR("kelx not defined", __func__, NULL); |
2265 | 0 | if (!kely) |
2266 | 0 | return (FPIX *)ERROR_PTR("kely not defined", __func__, NULL); |
2267 | | |
2268 | 0 | xfact = ConvolveSamplingFactX; |
2269 | 0 | yfact = ConvolveSamplingFactY; |
2270 | 0 | if (normflag) { |
2271 | 0 | kelxn = kernelNormalize(kelx, 1.0); |
2272 | 0 | kelyn = kernelNormalize(kely, 1.0); |
2273 | 0 | l_setConvolveSampling(xfact, 1); |
2274 | 0 | fpixt = fpixConvolve(fpixs, kelxn, 0); |
2275 | 0 | l_setConvolveSampling(1, yfact); |
2276 | 0 | fpixd = fpixConvolve(fpixt, kelyn, 0); |
2277 | 0 | l_setConvolveSampling(xfact, yfact); /* restore */ |
2278 | 0 | kernelDestroy(&kelxn); |
2279 | 0 | kernelDestroy(&kelyn); |
2280 | 0 | } else { /* don't normalize */ |
2281 | 0 | l_setConvolveSampling(xfact, 1); |
2282 | 0 | fpixt = fpixConvolve(fpixs, kelx, 0); |
2283 | 0 | l_setConvolveSampling(1, yfact); |
2284 | 0 | fpixd = fpixConvolve(fpixt, kely, 0); |
2285 | 0 | l_setConvolveSampling(xfact, yfact); |
2286 | 0 | } |
2287 | |
|
2288 | 0 | fpixDestroy(&fpixt); |
2289 | 0 | return fpixd; |
2290 | 0 | } |
2291 | | |
2292 | | |
2293 | | /*------------------------------------------------------------------------* |
2294 | | * Convolution with bias (for non-negative output) * |
2295 | | *------------------------------------------------------------------------*/ |
2296 | | /*! |
2297 | | * \brief pixConvolveWithBias() |
2298 | | * |
2299 | | * \param[in] pixs 8 bpp; no colormap |
2300 | | * \param[in] kel1 |
2301 | | * \param[in] kel2 can be null; use if separable |
2302 | | * \param[in] force8 if 1, force output to 8 bpp; otherwise, determine |
2303 | | * output depth by the dynamic range of pixel values |
2304 | | * \param[out] pbias applied bias |
2305 | | * \return pixd 8 or 16 bpp |
2306 | | * |
2307 | | * <pre> |
2308 | | * Notes: |
2309 | | * (1) This does a convolution with either a single kernel or |
2310 | | * a pair of separable kernels, and automatically applies whatever |
2311 | | * bias (shift) is required so that the resulting pixel values |
2312 | | * are non-negative. |
2313 | | * (2) The kernel is always normalized. If there are no negative |
2314 | | * values in the kernel, a standard normalized convolution is |
2315 | | * performed, with 8 bpp output. If the sum of kernel values is |
2316 | | * very close to zero, the kernel can not be normalized and |
2317 | | * the convolution will not be performed. An error message results. |
2318 | | * (3) If there are negative values in the kernel, the pix is |
2319 | | * converted to an fpix, the convolution is done on the fpix, and |
2320 | | * a bias (shift) may need to be applied. |
2321 | | * (4) If force8 == TRUE and the range of values after the convolution |
2322 | | * is > 255, the output values will be scaled to fit in [0 ... 255]. |
2323 | | * If force8 == FALSE, the output will be either 8 or 16 bpp, |
2324 | | * to accommodate the dynamic range of output values without scaling. |
2325 | | * </pre> |
2326 | | */ |
2327 | | PIX * |
2328 | | pixConvolveWithBias(PIX *pixs, |
2329 | | L_KERNEL *kel1, |
2330 | | L_KERNEL *kel2, |
2331 | | l_int32 force8, |
2332 | | l_int32 *pbias) |
2333 | 0 | { |
2334 | 0 | l_int32 outdepth; |
2335 | 0 | l_float32 min1, min2, min, minval, maxval, range; |
2336 | 0 | FPIX *fpix1, *fpix2; |
2337 | 0 | PIX *pixd; |
2338 | |
|
2339 | 0 | if (!pbias) |
2340 | 0 | return (PIX *)ERROR_PTR("&bias not defined", __func__, NULL); |
2341 | 0 | *pbias = 0; |
2342 | 0 | if (!pixs || pixGetDepth(pixs) != 8) |
2343 | 0 | return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", __func__, NULL); |
2344 | 0 | if (pixGetColormap(pixs)) |
2345 | 0 | return (PIX *)ERROR_PTR("pixs has colormap", __func__, NULL); |
2346 | 0 | if (!kel1) |
2347 | 0 | return (PIX *)ERROR_PTR("kel1 not defined", __func__, NULL); |
2348 | | |
2349 | | /* Determine if negative values can be produced in the convolution */ |
2350 | 0 | kernelGetMinMax(kel1, &min1, NULL); |
2351 | 0 | min2 = 0.0; |
2352 | 0 | if (kel2) |
2353 | 0 | kernelGetMinMax(kel2, &min2, NULL); |
2354 | 0 | min = L_MIN(min1, min2); |
2355 | |
|
2356 | 0 | if (min >= 0.0) { |
2357 | 0 | if (!kel2) |
2358 | 0 | return pixConvolve(pixs, kel1, 8, 1); |
2359 | 0 | else |
2360 | 0 | return pixConvolveSep(pixs, kel1, kel2, 8, 1); |
2361 | 0 | } |
2362 | | |
2363 | | /* Bias may need to be applied; convert to fpix and convolve */ |
2364 | 0 | fpix1 = pixConvertToFPix(pixs, 1); |
2365 | 0 | if (!kel2) |
2366 | 0 | fpix2 = fpixConvolve(fpix1, kel1, 1); |
2367 | 0 | else |
2368 | 0 | fpix2 = fpixConvolveSep(fpix1, kel1, kel2, 1); |
2369 | 0 | fpixDestroy(&fpix1); |
2370 | | |
2371 | | /* Determine the bias and the dynamic range. |
2372 | | * If the dynamic range is <= 255, just shift the values by the |
2373 | | * bias, if any. |
2374 | | * If the dynamic range is > 255, there are two cases: |
2375 | | * (1) the output depth is not forced to 8 bpp |
2376 | | * ==> apply the bias without scaling; outdepth = 16 |
2377 | | * (2) the output depth is forced to 8 |
2378 | | * ==> linearly map the pixel values to [0 ... 255]. */ |
2379 | 0 | fpixGetMin(fpix2, &minval, NULL, NULL); |
2380 | 0 | fpixGetMax(fpix2, &maxval, NULL, NULL); |
2381 | 0 | range = maxval - minval; |
2382 | 0 | *pbias = (minval < 0.0) ? -minval : 0.0; |
2383 | 0 | fpixAddMultConstant(fpix2, *pbias, 1.0); /* shift: min val ==> 0 */ |
2384 | 0 | if (range <= 255 || !force8) { /* no scaling of output values */ |
2385 | 0 | outdepth = (range > 255) ? 16 : 8; |
2386 | 0 | } else { /* scale output values to fit in 8 bpp */ |
2387 | 0 | fpixAddMultConstant(fpix2, 0.0, (255.0 / range)); |
2388 | 0 | outdepth = 8; |
2389 | 0 | } |
2390 | | |
2391 | | /* Convert back to pix; it won't do any clipping */ |
2392 | 0 | pixd = fpixConvertToPix(fpix2, outdepth, L_CLIP_TO_ZERO, 0); |
2393 | 0 | fpixDestroy(&fpix2); |
2394 | |
|
2395 | 0 | return pixd; |
2396 | 0 | } |
2397 | | |
2398 | | |
2399 | | /*------------------------------------------------------------------------* |
2400 | | * Set parameter for convolution subsampling * |
2401 | | *------------------------------------------------------------------------*/ |
2402 | | /*! |
2403 | | * \brief l_setConvolveSampling() |
2404 | | |
2405 | | * |
2406 | | * \param[in] xfact, yfact integer >= 1 |
2407 | | * \return void |
2408 | | * |
2409 | | * <pre> |
2410 | | * Notes: |
2411 | | * (1) This sets the x and y output subsampling factors for generic pix |
2412 | | * and fpix convolution. The default values are 1 (no subsampling). |
2413 | | * </pre> |
2414 | | */ |
2415 | | void |
2416 | | l_setConvolveSampling(l_int32 xfact, |
2417 | | l_int32 yfact) |
2418 | 0 | { |
2419 | 0 | if (xfact < 1) xfact = 1; |
2420 | 0 | if (yfact < 1) yfact = 1; |
2421 | 0 | ConvolveSamplingFactX = xfact; |
2422 | 0 | ConvolveSamplingFactY = yfact; |
2423 | 0 | } |
2424 | | |
2425 | | |
2426 | | /*------------------------------------------------------------------------* |
2427 | | * Additive gaussian noise * |
2428 | | *------------------------------------------------------------------------*/ |
2429 | | /*! |
2430 | | * \brief pixAddGaussianNoise() |
2431 | | * |
2432 | | * \param[in] pixs 8 bpp gray or 32 bpp rgb; no colormap |
2433 | | * \param[in] stdev of noise |
2434 | | * \return pixd 8 or 32 bpp, or NULL on error |
2435 | | * |
2436 | | * <pre> |
2437 | | * Notes: |
2438 | | * (1) This adds noise to each pixel, taken from a normal |
2439 | | * distribution with zero mean and specified standard deviation. |
2440 | | * </pre> |
2441 | | */ |
2442 | | PIX * |
2443 | | pixAddGaussianNoise(PIX *pixs, |
2444 | | l_float32 stdev) |
2445 | 0 | { |
2446 | 0 | l_int32 i, j, w, h, d, wpls, wpld, val, rval, gval, bval; |
2447 | 0 | l_uint32 pixel; |
2448 | 0 | l_uint32 *datas, *datad, *lines, *lined; |
2449 | 0 | PIX *pixd; |
2450 | |
|
2451 | 0 | if (!pixs) |
2452 | 0 | return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL); |
2453 | 0 | if (pixGetColormap(pixs)) |
2454 | 0 | return (PIX *)ERROR_PTR("pixs has colormap", __func__, NULL); |
2455 | 0 | pixGetDimensions(pixs, &w, &h, &d); |
2456 | 0 | if (d != 8 && d != 32) |
2457 | 0 | return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL); |
2458 | | |
2459 | 0 | pixd = pixCreateTemplate(pixs); |
2460 | 0 | datas = pixGetData(pixs); |
2461 | 0 | datad = pixGetData(pixd); |
2462 | 0 | wpls = pixGetWpl(pixs); |
2463 | 0 | wpld = pixGetWpl(pixd); |
2464 | 0 | for (i = 0; i < h; i++) { |
2465 | 0 | lines = datas + i * wpls; |
2466 | 0 | lined = datad + i * wpld; |
2467 | 0 | for (j = 0; j < w; j++) { |
2468 | 0 | if (d == 8) { |
2469 | 0 | val = GET_DATA_BYTE(lines, j); |
2470 | 0 | val += (l_int32)(stdev * gaussDistribSampling() + 0.5); |
2471 | 0 | val = L_MIN(255, L_MAX(0, val)); |
2472 | 0 | SET_DATA_BYTE(lined, j, val); |
2473 | 0 | } else { /* d = 32 */ |
2474 | 0 | pixel = *(lines + j); |
2475 | 0 | extractRGBValues(pixel, &rval, &gval, &bval); |
2476 | 0 | rval += (l_int32)(stdev * gaussDistribSampling() + 0.5); |
2477 | 0 | rval = L_MIN(255, L_MAX(0, rval)); |
2478 | 0 | gval += (l_int32)(stdev * gaussDistribSampling() + 0.5); |
2479 | 0 | gval = L_MIN(255, L_MAX(0, gval)); |
2480 | 0 | bval += (l_int32)(stdev * gaussDistribSampling() + 0.5); |
2481 | 0 | bval = L_MIN(255, L_MAX(0, bval)); |
2482 | 0 | composeRGBPixel(rval, gval, bval, lined + j); |
2483 | 0 | } |
2484 | 0 | } |
2485 | 0 | } |
2486 | 0 | return pixd; |
2487 | 0 | } |
2488 | | |
2489 | | |
2490 | | /*! |
2491 | | * \brief gaussDistribSampling() |
2492 | | * |
2493 | | * \return gaussian distributed variable with zero mean and unit stdev |
2494 | | * |
2495 | | * <pre> |
2496 | | * Notes: |
2497 | | * (1) For an explanation of the Box-Muller method for generating |
2498 | | * a normally distributed random variable with zero mean and |
2499 | | * unit standard deviation, see Numerical Recipes in C, |
2500 | | * 2nd edition, p. 288ff. |
2501 | | * (2) This can be called sequentially to get samples that can be |
2502 | | * used for adding noise to each pixel of an image, for example. |
2503 | | * </pre> |
2504 | | */ |
2505 | | l_float32 |
2506 | | gaussDistribSampling(void) |
2507 | 0 | { |
2508 | 0 | static l_int32 select = 0; /* flips between 0 and 1 on successive calls */ |
2509 | 0 | static l_float32 saveval; |
2510 | 0 | l_float32 frand, xval, yval, rsq, factor; |
2511 | |
|
2512 | 0 | if (select == 0) { |
2513 | 0 | while (1) { /* choose a point in a 2x2 square, centered at origin */ |
2514 | 0 | frand = (l_float32)rand() / (l_float32)RAND_MAX; |
2515 | 0 | xval = 2.0 * frand - 1.0; |
2516 | 0 | frand = (l_float32)rand() / (l_float32)RAND_MAX; |
2517 | 0 | yval = 2.0 * frand - 1.0; |
2518 | 0 | rsq = xval * xval + yval * yval; |
2519 | 0 | if (rsq > 0.0 && rsq < 1.0) /* point is inside the unit circle */ |
2520 | 0 | break; |
2521 | 0 | } |
2522 | 0 | factor = sqrt(-2.0 * log(rsq) / rsq); |
2523 | 0 | saveval = xval * factor; |
2524 | 0 | select = 1; |
2525 | 0 | return yval * factor; |
2526 | 0 | } |
2527 | 0 | else { |
2528 | 0 | select = 0; |
2529 | 0 | return saveval; |
2530 | 0 | } |
2531 | 0 | } |