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