1 /*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 - This software is distributed in the hope that it will be
4 - useful, but with NO WARRANTY OF ANY KIND.
5 - No author or distributor accepts responsibility to anyone for the
6 - consequences of using this software, or for whether it serves any
7 - particular purpose or works at all, unless he or she says so in
8 - writing. Everyone is granted permission to copy, modify and
9 - redistribute this source code, for commercial or non-commercial
10 - purposes, with the following restrictions: (1) the origin of this
11 - source code must not be misrepresented; (2) modified versions must
12 - be plainly marked as such; and (3) this notice may not be removed
13 - or altered from any source or modified source distribution.
14 *====================================================================*/
15
16 /*
17 * convolve.c
18 *
19 * Top level grayscale or color block convolution
20 * PIX *pixBlockconv()
21 *
22 * Grayscale block convolution
23 * PIX *pixBlockconvGray()
24 *
25 * Accumulator for 1, 8 and 32 bpp convolution
26 * PIX *pixBlockconvAccum()
27 *
28 * Un-normalized grayscale block convolution
29 * PIX *pixBlockconvGrayUnnormalized()
30 *
31 * Tiled grayscale or color block convolution
32 * PIX *pixBlockconvTiled()
33 * PIX *pixBlockconvGrayTile()
34 *
35 * Convolution for average in specified window
36 * PIX *pixWindowedMean()
37 *
38 * Convolution for average square value in specified window
39 * PIX *pixWindowedMeanSquare()
40 * DPIX *pixMeanSquareAccum()
41 *
42 * Binary block sum and rank filter
43 * PIX *pixBlockrank()
44 * PIX *pixBlocksum()
45 *
46 * Woodfill transform
47 * PIX *pixWoodfillTransform()
48 *
49 * Generic convolution (with Pix)
50 * PIX *pixConvolve()
51 * PIX *pixConvolveSep()
52 *
53 * Generic convolution (with float arrays)
54 * FPIX *fpixConvolve()
55 * FPIX *fpixConvolveSep()
56 */
57
58 #include <stdio.h>
59 #include <stdlib.h>
60 #include "allheaders.h"
61
62
63 /*----------------------------------------------------------------------*
64 * Top-level grayscale or color block convolution *
65 *----------------------------------------------------------------------*/
66 /*!
67 * pixBlockconv()
68 *
69 * Input: pix (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
70 * wc, hc (half width/height of convolution kernel)
71 * Return: pixd, or null on error
72 *
73 * Notes:
74 * (1) The full width and height of the convolution kernel
75 * are (2 * wc + 1) and (2 * hc + 1)
76 * (2) Returns a copy if both wc and hc are 0
77 * (3) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1,
78 * where (w,h) are the dimensions of pixs.
79 */
80 PIX *
pixBlockconv(PIX * pix,l_int32 wc,l_int32 hc)81 pixBlockconv(PIX *pix,
82 l_int32 wc,
83 l_int32 hc)
84 {
85 l_int32 w, h, d;
86 PIX *pixs, *pixd, *pixr, *pixrc, *pixg, *pixgc, *pixb, *pixbc;
87
88 PROCNAME("pixBlockconv");
89
90 if (!pix)
91 return (PIX *)ERROR_PTR("pix not defined", procName, NULL);
92 if (wc < 0) wc = 0;
93 if (hc < 0) hc = 0;
94 pixGetDimensions(pix, &w, &h, &d);
95 if (w < 2 * wc + 1 || h < 2 * hc + 1) {
96 wc = L_MIN(wc, (w - 1) / 2);
97 hc = L_MIN(hc, (h - 1) / 2);
98 L_WARNING("kernel too large; reducing!", procName);
99 L_INFO_INT2("wc = %d, hc = %d", procName, wc, hc);
100 }
101 if (wc == 0 && hc == 0) /* no-op */
102 return pixCopy(NULL, pix);
103
104 /* Remove colormap if necessary */
105 if ((d == 2 || d == 4 || d == 8) && pixGetColormap(pix)) {
106 L_WARNING("pix has colormap; removing", procName);
107 pixs = pixRemoveColormap(pix, REMOVE_CMAP_BASED_ON_SRC);
108 d = pixGetDepth(pixs);
109 }
110 else
111 pixs = pixClone(pix);
112
113 if (d != 8 && d != 32) {
114 pixDestroy(&pixs);
115 return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, NULL);
116 }
117
118 if (d == 8)
119 pixd = pixBlockconvGray(pixs, NULL, wc, hc);
120 else { /* d == 32 */
121 pixr = pixGetRGBComponent(pixs, COLOR_RED);
122 pixrc = pixBlockconvGray(pixr, NULL, wc, hc);
123 pixDestroy(&pixr);
124 pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
125 pixgc = pixBlockconvGray(pixg, NULL, wc, hc);
126 pixDestroy(&pixg);
127 pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
128 pixbc = pixBlockconvGray(pixb, NULL, wc, hc);
129 pixDestroy(&pixb);
130 pixd = pixCreateRGBImage(pixrc, pixgc, pixbc);
131 pixDestroy(&pixrc);
132 pixDestroy(&pixgc);
133 pixDestroy(&pixbc);
134 }
135
136 pixDestroy(&pixs);
137 return pixd;
138 }
139
140
141 /*----------------------------------------------------------------------*
142 * Grayscale block convolution *
143 *----------------------------------------------------------------------*/
144 /*!
145 * pixBlockconvGray()
146 *
147 * Input: pix (8 bpp)
148 * accum pix (32 bpp; can be null)
149 * wc, hc (half width/height of convolution kernel)
150 * Return: pix (8 bpp), or null on error
151 *
152 * Notes:
153 * (1) If accum pix is null, make one and destroy it before
154 * returning; otherwise, just use the input accum pix.
155 * (2) The full width and height of the convolution kernel
156 * are (2 * wc + 1) and (2 * hc + 1).
157 * (3) Returns a copy if both wc and hc are 0.
158 * (4) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1,
159 * where (w,h) are the dimensions of pixs.
160 */
161 PIX *
pixBlockconvGray(PIX * pixs,PIX * pixacc,l_int32 wc,l_int32 hc)162 pixBlockconvGray(PIX *pixs,
163 PIX *pixacc,
164 l_int32 wc,
165 l_int32 hc)
166 {
167 l_int32 w, h, d, wpl, wpla;
168 l_uint32 *datad, *dataa;
169 PIX *pixd, *pixt;
170
171 PROCNAME("pixBlockconvGray");
172
173 if (!pixs)
174 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
175 pixGetDimensions(pixs, &w, &h, &d);
176 if (d != 8)
177 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
178 if (wc < 0) wc = 0;
179 if (hc < 0) hc = 0;
180 if (w < 2 * wc + 1 || h < 2 * hc + 1) {
181 wc = L_MIN(wc, (w - 1) / 2);
182 hc = L_MIN(hc, (h - 1) / 2);
183 L_WARNING("kernel too large; reducing!", procName);
184 L_INFO_INT2("wc = %d, hc = %d", procName, wc, hc);
185 }
186 if (wc == 0 && hc == 0) /* no-op */
187 return pixCopy(NULL, pixs);
188
189 if (pixacc) {
190 if (pixGetDepth(pixacc) == 32)
191 pixt = pixClone(pixacc);
192 else {
193 L_WARNING("pixacc not 32 bpp; making new one", procName);
194 if ((pixt = pixBlockconvAccum(pixs)) == NULL)
195 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
196 }
197 }
198 else {
199 if ((pixt = pixBlockconvAccum(pixs)) == NULL)
200 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
201 }
202
203 if ((pixd = pixCreateTemplate(pixs)) == NULL)
204 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
205
206 wpl = pixGetWpl(pixs);
207 wpla = pixGetWpl(pixt);
208 datad = pixGetData(pixd);
209 dataa = pixGetData(pixt);
210 blockconvLow(datad, w, h, wpl, dataa, wpla, wc, hc);
211
212 pixDestroy(&pixt);
213 return pixd;
214 }
215
216
217 /*----------------------------------------------------------------------*
218 * Accumulator for 1, 8 and 32 bpp convolution *
219 *----------------------------------------------------------------------*/
220 /*!
221 * pixBlockconvAccum()
222 *
223 * Input: pixs (1, 8 or 32 bpp)
224 * Return: accum pix (32 bpp), or null on error.
225 *
226 * Notes:
227 * (1) The general recursion relation is
228 * a(i,j) = v(i,j) + a(i-1, j) + a(i, j-1) - a(i-1, j-1)
229 * For the first line, this reduces to the special case
230 * a(i,j) = v(i,j) + a(i, j-1)
231 * For the first column, the special case is
232 * a(i,j) = v(i,j) + a(i-1, j)
233 */
234 PIX *
pixBlockconvAccum(PIX * pixs)235 pixBlockconvAccum(PIX *pixs)
236 {
237 l_int32 w, h, d, wpls, wpld;
238 l_uint32 *datas, *datad;
239 PIX *pixd;
240
241 PROCNAME("pixBlockconvAccum");
242
243 if (!pixs)
244 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
245
246 pixGetDimensions(pixs, &w, &h, &d);
247 if (d != 1 && d != 8 && d != 32)
248 return (PIX *)ERROR_PTR("pixs not 1, 8 or 32 bpp", procName, NULL);
249 if ((pixd = pixCreate(w, h, 32)) == NULL)
250 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
251
252 datas = pixGetData(pixs);
253 datad = pixGetData(pixd);
254 wpls = pixGetWpl(pixs);
255 wpld = pixGetWpl(pixd);
256 blockconvAccumLow(datad, w, h, wpld, datas, d, wpls);
257
258 return pixd;
259 }
260
261
262 /*----------------------------------------------------------------------*
263 * Un-normalized grayscale block convolution *
264 *----------------------------------------------------------------------*/
265 /*!
266 * pixBlockconvGrayUnnormalized()
267 *
268 * Input: pixs (8 bpp)
269 * wc, hc (half width/height of convolution kernel)
270 * Return: pix (32 bpp; containing the convolution without normalizing
271 * for the window size), or null on error
272 *
273 * Notes:
274 * (1) The full width and height of the convolution kernel
275 * are (2 * wc + 1) and (2 * hc + 1).
276 * (2) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1,
277 * where (w,h) are the dimensions of pixs.
278 * (3) Returns a copy if both wc and hc are 0.
279 * (3) Adds mirrored border to avoid treating the boundary pixels
280 * specially. Note that we add wc + 1 pixels to the left
281 * and wc to the right. The added width is 2 * wc + 1 pixels,
282 * and the particular choice simplifies the indexing in the loop.
283 * Likewise, add hc + 1 pixels to the top and hc to the bottom.
284 * (4) To get the normalized result, divide by the area of the
285 * convolution kernel: (2 * wc + 1) * (2 * hc + 1)
286 * Specifically, do this:
287 * pixc = pixBlockconvGrayUnnormalized(pixs, wc, hc);
288 * fract = 1. / ((2 * wc + 1) * (2 * hc + 1));
289 * pixMultConstantGray(pixc, fract);
290 * pixd = pixGetRGBComponent(pixc, L_ALPHA_CHANNEL);
291 * (5) Unlike pixBlockconvGray(), this always computes the accumulation
292 * pix because its size is tied to wc and hc.
293 * (6) Compare this implementation with pixBlockconvGray(), where
294 * most of the code in blockconvLow() is special casing for
295 * efficiently handling the boundary. Here, the use of
296 * mirrored borders and destination indexing makes the
297 * implementation very simple.
298 */
299 PIX *
pixBlockconvGrayUnnormalized(PIX * pixs,l_int32 wc,l_int32 hc)300 pixBlockconvGrayUnnormalized(PIX *pixs,
301 l_int32 wc,
302 l_int32 hc)
303 {
304 l_int32 i, j, w, h, d, wpla, wpld, jmax;
305 l_uint32 *linemina, *linemaxa, *lined, *dataa, *datad;
306 PIX *pixsb, *pixacc, *pixd;
307
308 PROCNAME("pixBlockconvGrayUnnormalized");
309
310 if (!pixs)
311 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
312 pixGetDimensions(pixs, &w, &h, &d);
313 if (d != 8)
314 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
315 if (wc < 0) wc = 0;
316 if (hc < 0) hc = 0;
317 if (w < 2 * wc + 1 || h < 2 * hc + 1) {
318 wc = L_MIN(wc, (w - 1) / 2);
319 hc = L_MIN(hc, (h - 1) / 2);
320 L_WARNING("kernel too large; reducing!", procName);
321 L_INFO_INT2("wc = %d, hc = %d", procName, wc, hc);
322 }
323 if (wc == 0 && hc == 0) /* no-op */
324 return pixCopy(NULL, pixs);
325
326 if ((pixsb = pixAddMirroredBorder(pixs, wc + 1, wc, hc + 1, hc)) == NULL)
327 return (PIX *)ERROR_PTR("pixsb not made", procName, NULL);
328 if ((pixacc = pixBlockconvAccum(pixsb)) == NULL)
329 return (PIX *)ERROR_PTR("pixacc not made", procName, NULL);
330 pixDestroy(&pixsb);
331 if ((pixd = pixCreate(w, h, 32)) == NULL)
332 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
333
334 wpla = pixGetWpl(pixacc);
335 wpld = pixGetWpl(pixd);
336 datad = pixGetData(pixd);
337 dataa = pixGetData(pixacc);
338 for (i = 0; i < h; i++) {
339 lined = datad + i * wpld;
340 linemina = dataa + i * wpla;
341 linemaxa = dataa + (i + 2 * hc + 1) * wpla;
342 for (j = 0; j < w; j++) {
343 jmax = j + 2 * wc + 1;
344 lined[j] = linemaxa[jmax] - linemaxa[j] -
345 linemina[jmax] + linemina[j];
346 }
347 }
348
349 pixDestroy(&pixacc);
350 return pixd;
351 }
352
353
354 /*----------------------------------------------------------------------*
355 * Tiled grayscale or color block convolution *
356 *----------------------------------------------------------------------*/
357 /*!
358 * pixBlockconvTiled()
359 *
360 * Input: pix (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
361 * wc, hc (half width/height of convolution kernel)
362 * nx, ny (subdivision into tiles)
363 * Return: pixd, or null on error
364 *
365 * Notes:
366 * (1) The full width and height of the convolution kernel
367 * are (2 * wc + 1) and (2 * hc + 1)
368 * (2) Returns a copy if both wc and hc are 0
369 * (3) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1,
370 * where (w,h) are the dimensions of pixs.
371 * (4) For nx == ny == 1, this defaults to pixBlockconv(), which
372 * is typically about twice as fast, and gives nearly
373 * identical results as pixBlockconvGrayTile().
374 * (5) If the tiles are too small, nx and/or ny are reduced
375 * a minimum amount so that the tiles are expanded to the
376 * smallest workable size in the problematic direction(s).
377 * (6) Why a tiled version? Three reasons:
378 * (a) Because the accumulator is a uint32, overflow can occur
379 * for an image with more than 16M pixels.
380 * (b) The accumulator array for 16M pixels is 64 MB; using
381 * tiles reduces the size of this array.
382 * (c) Each tile can be processed independently, in parallel,
383 * on a multicore processor.
384 */
385 PIX *
pixBlockconvTiled(PIX * pix,l_int32 wc,l_int32 hc,l_int32 nx,l_int32 ny)386 pixBlockconvTiled(PIX *pix,
387 l_int32 wc,
388 l_int32 hc,
389 l_int32 nx,
390 l_int32 ny)
391 {
392 l_int32 i, j, w, h, d, xrat, yrat;
393 PIX *pixs, *pixd, *pixc, *pixt;
394 PIX *pixr, *pixrc, *pixg, *pixgc, *pixb, *pixbc;
395 PIXTILING *pt;
396
397 PROCNAME("pixBlockconvTiled");
398
399 if (!pix)
400 return (PIX *)ERROR_PTR("pix not defined", procName, NULL);
401 if (wc < 0) wc = 0;
402 if (hc < 0) hc = 0;
403 pixGetDimensions(pix, &w, &h, &d);
404 if (w < 2 * wc + 3 || h < 2 * hc + 3) {
405 wc = L_MAX(0, L_MIN(wc, (w - 3) / 2));
406 hc = L_MAX(0, L_MIN(hc, (h - 3) / 2));
407 L_WARNING("kernel too large; reducing!", procName);
408 L_INFO_INT2("wc = %d, hc = %d", procName, wc, hc);
409 }
410 if (wc == 0 && hc == 0) /* no-op */
411 return pixCopy(NULL, pix);
412 if (nx <= 1 && ny <= 1)
413 return pixBlockconv(pix, wc, hc);
414
415 /* Test to see if the tiles are too small. The required
416 * condition is that the tile dimensions must be at least
417 * (wc + 2) x (hc + 2). */
418 xrat = w / nx;
419 yrat = h / ny;
420 if (xrat < wc + 2) {
421 nx = w / (wc + 2);
422 L_WARNING_INT("tile width too small; nx reduced to %d", procName, nx);
423 }
424 if (yrat < hc + 2) {
425 ny = h / (hc + 2);
426 L_WARNING_INT("tile height too small; ny reduced to %d", procName, ny);
427 }
428
429 /* Remove colormap if necessary */
430 if ((d == 2 || d == 4 || d == 8) && pixGetColormap(pix)) {
431 L_WARNING("pix has colormap; removing", procName);
432 pixs = pixRemoveColormap(pix, REMOVE_CMAP_BASED_ON_SRC);
433 d = pixGetDepth(pixs);
434 }
435 else
436 pixs = pixClone(pix);
437
438 if (d != 8 && d != 32) {
439 pixDestroy(&pixs);
440 return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, NULL);
441 }
442
443 /* Note that the overlaps in the width and height that
444 * are added to the tile are (wc + 2) and (hc + 2).
445 * These overlaps are removed by pixTilingPaintTile().
446 * They are larger than the extent of the filter because
447 * although the filter is symmetric with respect to its origin,
448 * the implementation is asymmetric -- see the implementation in
449 * pixBlockconvGrayTile(). */
450 pixd = pixCreateTemplateNoInit(pixs);
451 pt = pixTilingCreate(pixs, nx, ny, 0, 0, wc + 2, hc + 2);
452 for (i = 0; i < ny; i++) {
453 for (j = 0; j < nx; j++) {
454 pixt = pixTilingGetTile(pt, i, j);
455
456 /* Convolve over the tile */
457 if (d == 8)
458 pixc = pixBlockconvGrayTile(pixt, NULL, wc, hc);
459 else { /* d == 32 */
460 pixr = pixGetRGBComponent(pixt, COLOR_RED);
461 pixrc = pixBlockconvGrayTile(pixr, NULL, wc, hc);
462 pixDestroy(&pixr);
463 pixg = pixGetRGBComponent(pixt, COLOR_GREEN);
464 pixgc = pixBlockconvGrayTile(pixg, NULL, wc, hc);
465 pixDestroy(&pixg);
466 pixb = pixGetRGBComponent(pixt, COLOR_BLUE);
467 pixbc = pixBlockconvGrayTile(pixb, NULL, wc, hc);
468 pixDestroy(&pixb);
469 pixc = pixCreateRGBImage(pixrc, pixgc, pixbc);
470 pixDestroy(&pixrc);
471 pixDestroy(&pixgc);
472 pixDestroy(&pixbc);
473 }
474
475 pixTilingPaintTile(pixd, i, j, pixc, pt);
476 pixDestroy(&pixt);
477 pixDestroy(&pixc);
478 }
479 }
480
481 pixDestroy(&pixs);
482 pixTilingDestroy(&pt);
483 return pixd;
484 }
485
486
487 /*!
488 * pixBlockconvGrayTile()
489 *
490 * Input: pixs (8 bpp gray)
491 * pixacc (32 bpp accum pix)
492 * wc, hc (half width/height of convolution kernel)
493 * Return: pixd, or null on error
494 *
495 * Notes:
496 * (1) The full width and height of the convolution kernel
497 * are (2 * wc + 1) and (2 * hc + 1)
498 * (2) Assumes that the input pixs is padded with (wc + 1) pixels on
499 * left and right, and with (hc + 1) pixels on top and bottom.
500 * The returned pix has these stripped off; they are only used
501 * for computation.
502 * (3) Returns a copy if both wc and hc are 0
503 * (4) Require that w > 2 * wc + 1 and h > 2 * hc + 1,
504 * where (w,h) are the dimensions of pixs.
505 */
506 PIX *
pixBlockconvGrayTile(PIX * pixs,PIX * pixacc,l_int32 wc,l_int32 hc)507 pixBlockconvGrayTile(PIX *pixs,
508 PIX *pixacc,
509 l_int32 wc,
510 l_int32 hc)
511 {
512 l_int32 w, h, d, wd, hd, i, j, imin, imax, jmin, jmax, wplt, wpld;
513 l_float32 norm;
514 l_uint32 val;
515 l_uint32 *datat, *datad, *lined, *linemint, *linemaxt;
516 PIX *pixt, *pixd;
517
518 PROCNAME("pixBlockconvGrayTile");
519
520 if (!pixs)
521 return (PIX *)ERROR_PTR("pix not defined", procName, NULL);
522 pixGetDimensions(pixs, &w, &h, &d);
523 if (d != 8)
524 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
525 if (wc < 0) wc = 0;
526 if (hc < 0) hc = 0;
527 if (w < 2 * wc + 3 || h < 2 * hc + 3) {
528 wc = L_MAX(0, L_MIN(wc, (w - 3) / 2));
529 hc = L_MAX(0, L_MIN(hc, (h - 3) / 2));
530 L_WARNING("kernel too large; reducing!", procName);
531 L_INFO_INT2("wc = %d, hc = %d", procName, wc, hc);
532 }
533 if (wc == 0 && hc == 0)
534 return pixCopy(NULL, pixs);
535 wd = w - 2 * wc;
536 hd = h - 2 * hc;
537
538 if (pixacc) {
539 if (pixGetDepth(pixacc) == 32)
540 pixt = pixClone(pixacc);
541 else {
542 L_WARNING("pixacc not 32 bpp; making new one", procName);
543 if ((pixt = pixBlockconvAccum(pixs)) == NULL)
544 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
545 }
546 }
547 else {
548 if ((pixt = pixBlockconvAccum(pixs)) == NULL)
549 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
550 }
551
552 if ((pixd = pixCreateTemplate(pixs)) == NULL)
553 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
554 datat = pixGetData(pixt);
555 wplt = pixGetWpl(pixt);
556 datad = pixGetData(pixd);
557 wpld = pixGetWpl(pixd);
558 norm = 1. / (l_float32)((2 * wc + 1) * (2 * hc + 1));
559
560 /* Do the convolution over the subregion of size (wd - 2, hd - 2),
561 * which exactly corresponds to the size of the subregion that
562 * will be extracted by pixTilingPaintTile(). Note that the
563 * region in which points are computed is not symmetric about
564 * the center of the images; instead the computation in
565 * the accumulator image is shifted up and to the left by 1,
566 * relative to the center, because the 4 accumulator sampling
567 * points are taken at the LL corner of the filter and at 3 other
568 * points that are shifted -wc and -hc to the left and above. */
569 for (i = hc; i < hc + hd - 2; i++) {
570 imin = L_MAX(i - hc - 1, 0);
571 imax = L_MIN(i + hc, h - 1);
572 lined = datad + i * wpld;
573 linemint = datat + imin * wplt;
574 linemaxt = datat + imax * wplt;
575 for (j = wc; j < wc + wd - 2; j++) {
576 jmin = L_MAX(j - wc - 1, 0);
577 jmax = L_MIN(j + wc, w - 1);
578 val = linemaxt[jmax] - linemaxt[jmin]
579 + linemint[jmin] - linemint[jmax];
580 val = (l_uint8)(norm * val + 0.5);
581 SET_DATA_BYTE(lined, j, val);
582 }
583 }
584
585 pixDestroy(&pixt);
586 return pixd;
587 }
588
589
590 /*----------------------------------------------------------------------*
591 * Convolution for average in specified window *
592 *----------------------------------------------------------------------*/
593 /*!
594 * pixWindowedMean()
595 *
596 * Input: pixs (8 or 32 bpp grayscale)
597 * wc, hc (half width/height of convolution kernel)
598 * normflag (1 for normalization to get average in window;
599 * 0 for the sum in the window (un-normalized))
600 * Return: pixd (8 or 32 bpp, average over kernel window)
601 *
602 * Notes:
603 * (1) The input and output depths are the same.
604 * (2) A set of border pixels of width (wc + 1) on left and right,
605 * and of height (hc + 1) on top and bottom, is included in pixs.
606 * The output pixd (after convolution) has this border removed.
607 * (3) Typically, @normflag == 1. However, if you want the sum
608 * within the window, rather than a normalized convolution,
609 * use @normflag == 0.
610 * (4) This builds a block accumulator pix, uses it here, and
611 * destroys it.
612 */
613 PIX *
pixWindowedMean(PIX * pixs,l_int32 wc,l_int32 hc,l_int32 normflag)614 pixWindowedMean(PIX *pixs,
615 l_int32 wc,
616 l_int32 hc,
617 l_int32 normflag)
618 {
619 l_int32 i, j, w, h, d, wd, hd, wplc, wpld, wincr, hincr;
620 l_uint32 val;
621 l_uint32 *datac, *datad, *linec1, *linec2, *lined;
622 l_float32 norm;
623 PIX *pixc, *pixd;
624
625 PROCNAME("pixWindowedMean");
626
627 if (!pixs)
628 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
629 pixGetDimensions(pixs, &w, &h, &d);
630 if (d != 8 && d != 32)
631 return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
632 if (wc < 2 || hc < 2)
633 return (PIX *)ERROR_PTR("wc and hc not >= 2", procName, NULL);
634
635 /* Strip off wc + 1 border pixels from each side and
636 * hc + 1 border pixels from top and bottom. */
637 wd = w - 2 * (wc + 1);
638 hd = h - 2 * (hc + 1);
639 if (wd < 2 || hd < 2)
640 return (PIX *)ERROR_PTR("w or h too small for kernel", procName, NULL);
641 if ((pixd = pixCreate(wd, hd, d)) == NULL)
642 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
643
644 /* Make the accumulator pix */
645 if ((pixc = pixBlockconvAccum(pixs)) == NULL)
646 return (PIX *)ERROR_PTR("pixc not made", procName, NULL);
647 wplc = pixGetWpl(pixc);
648 wpld = pixGetWpl(pixd);
649 datad = pixGetData(pixd);
650 datac = pixGetData(pixc);
651
652 wincr = 2 * wc + 1;
653 hincr = 2 * hc + 1;
654 norm = 1.0; /* use this for sum-in-window */
655 if (normflag)
656 norm = 1.0 / (wincr * hincr);
657 for (i = 0; i < hd; i++) {
658 linec1 = datac + i * wplc;
659 linec2 = datac + (i + hincr) * wplc;
660 lined = datad + i * wpld;
661 for (j = 0; j < wd; j++) {
662 val = linec2[j + wincr] - linec2[j] - linec1[j + wincr] + linec1[j];
663 if (d == 8) {
664 val = (l_uint8)(norm * val);
665 SET_DATA_BYTE(lined, j, val);
666 } else { /* d == 32 */
667 val = (l_uint32)(norm * val);
668 lined[j] = val;
669 }
670 }
671 }
672
673 pixDestroy(&pixc);
674 return pixd;
675 }
676
677
678 /*----------------------------------------------------------------------*
679 * Convolution for average square value in specified window *
680 *----------------------------------------------------------------------*/
681 /*!
682 * pixWindowedMeanSquare()
683 *
684 * Input: pixs (8 bpp grayscale)
685 * size (halfwidth of convolution kernel)
686 * Return: pixd (32 bpp, average over window of size (2 * size + 1))
687 *
688 * Notes:
689 * (1) A set of border pixels of width (size + 1) is included
690 * in pixs. The output pixd (after convolution) has this
691 * border removed.
692 * (2) The advantage is that we are unaffected by the boundary, and
693 * it is not necessary to treat pixels within @size of the
694 * border differently. This is because processing for pixd
695 * only takes place for pixels in pixs for which the
696 * kernel is entirely contained in pixs.
697 * (3) Why do we have an added border of width (@size + 1), when
698 * we only need @size pixels to satisfy this condition?
699 * Answer: the accumulators are asymmetric, requiring an
700 * extra row and column of pixels at top and left to work
701 * accurately.
702 */
703 PIX *
pixWindowedMeanSquare(PIX * pixs,l_int32 size)704 pixWindowedMeanSquare(PIX *pixs,
705 l_int32 size)
706 {
707 l_int32 i, j, w, h, wd, hd, wpl, wpld, incr;
708 l_uint32 ival;
709 l_uint32 *datad, *lined;
710 l_float64 norm;
711 l_float64 val;
712 l_float64 *data, *line1, *line2;
713 DPIX *dpix;
714 PIX *pixd;
715
716 PROCNAME("pixWindowedMeanSquare");
717
718
719 if (!pixs || (pixGetDepth(pixs) != 8))
720 return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", procName, NULL);
721 pixGetDimensions(pixs, &w, &h, NULL);
722 if (size < 2)
723 return (PIX *)ERROR_PTR("size not >= 2", procName, NULL);
724
725 if ((dpix = pixMeanSquareAccum(pixs)) == NULL)
726 return (PIX *)ERROR_PTR("dpix not made", procName, NULL);
727 wpl = dpixGetWpl(dpix);
728 data = dpixGetData(dpix);
729
730 /* Strip off 2 * (size + 1) border pixels */
731 wd = w - 2 * (size + 1);
732 hd = h - 2 * (size + 1);
733 if ((pixd = pixCreate(wd, hd, 32)) == NULL)
734 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
735 wpld = pixGetWpl(pixd);
736 datad = pixGetData(pixd);
737
738 incr = 2 * size + 1;
739 norm = 1.0 / (incr * incr);
740 for (i = 0; i < hd; i++) {
741 line1 = data + i * wpl;
742 line2 = data + (i + incr) * wpl;
743 lined = datad + i * wpld;
744 for (j = 0; j < wd; j++) {
745 val = line2[j + incr] - line2[j] - line1[j + incr] + line1[j];
746 ival = (l_uint32)(norm * val);
747 lined[j] = ival;
748 }
749 }
750
751 dpixDestroy(&dpix);
752 return pixd;
753 }
754
755
756 /*!
757 * pixMeanSquareAccum()
758 *
759 * Input: pixs (8 bpp grayscale)
760 * Return: dpix (64 bit array), or null on error
761 *
762 * Notes:
763 * (1) Similar to pixBlockconvAccum(), this computes the
764 * sum of the squares of the pixel values in such a way
765 * that the value at (i,j) is the sum of all squares in
766 * the rectangle from the origin to (i,j).
767 * (2) The general recursion relation (v are squared pixel values) is
768 * a(i,j) = v(i,j) + a(i-1, j) + a(i, j-1) - a(i-1, j-1)
769 * For the first line, this reduces to the special case
770 * a(i,j) = v(i,j) + a(i, j-1)
771 * For the first column, the special case is
772 * a(i,j) = v(i,j) + a(i-1, j)
773 */
774 DPIX *
pixMeanSquareAccum(PIX * pixs)775 pixMeanSquareAccum(PIX *pixs)
776 {
777 l_int32 i, j, w, h, wpl, wpls, val;
778 l_uint32 *datas, *lines;
779 l_float64 *data, *line, *linep;
780 DPIX *dpix;
781
782 PROCNAME("pixMeanSquareAccum");
783
784
785 if (!pixs || (pixGetDepth(pixs) != 8))
786 return (DPIX *)ERROR_PTR("pixs undefined or not 8 bpp", procName, NULL);
787 pixGetDimensions(pixs, &w, &h, NULL);
788 if ((dpix = dpixCreate(w, h)) == NULL)
789 return (DPIX *)ERROR_PTR("dpix not made", procName, NULL);
790
791 datas = pixGetData(pixs);
792 wpls = pixGetWpl(pixs);
793 data = dpixGetData(dpix);
794 wpl = dpixGetWpl(dpix);
795
796 lines = datas;
797 line = data;
798 for (j = 0; j < w; j++) { /* first line */
799 val = GET_DATA_BYTE(lines, j);
800 if (j == 0)
801 line[0] = val * val;
802 else
803 line[j] = line[j - 1] + val * val;
804 }
805
806 /* Do the other lines */
807 for (i = 1; i < h; i++) {
808 lines = datas + i * wpls;
809 line = data + i * wpl; /* current dest line */
810 linep = line - wpl;; /* prev dest line */
811 for (j = 0; j < w; j++) {
812 val = GET_DATA_BYTE(lines, j);
813 if (j == 0)
814 line[0] = linep[0] + val * val;
815 else
816 line[j] = line[j - 1] + linep[j] - linep[j - 1] + val * val;
817 }
818 }
819
820 return dpix;
821 }
822
823
824
825 /*----------------------------------------------------------------------*
826 * Binary block sum/rank *
827 *----------------------------------------------------------------------*/
828 /*!
829 * pixBlockrank()
830 *
831 * Input: pixs (1 bpp)
832 * accum pix (<optional> 32 bpp)
833 * wc, hc (half width/height of block sum/rank kernel)
834 * rank (between 0.0 and 1.0; 0.5 is median filter)
835 * Return: pixd (1 bpp)
836 *
837 * Notes:
838 * (1) The full width and height of the convolution kernel
839 * are (2 * wc + 1) and (2 * hc + 1)
840 * (2) This returns a pixd where each pixel is a 1 if the
841 * neighborhood (2 * wc + 1) x (2 * hc + 1)) pixels
842 * contains the rank fraction of 1 pixels. Otherwise,
843 * the returned pixel is 0. Note that the special case
844 * of rank = 0.0 is always satisfied, so the returned
845 * pixd has all pixels with value 1.
846 * (3) If accum pix is null, make one, use it, and destroy it
847 * before returning; otherwise, just use the input accum pix
848 * (4) If both wc and hc are 0, returns a copy unless rank == 0.0,
849 * in which case this returns an all-ones image.
850 * (5) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1,
851 * where (w,h) are the dimensions of pixs.
852 */
853 PIX *
pixBlockrank(PIX * pixs,PIX * pixacc,l_int32 wc,l_int32 hc,l_float32 rank)854 pixBlockrank(PIX *pixs,
855 PIX *pixacc,
856 l_int32 wc,
857 l_int32 hc,
858 l_float32 rank)
859 {
860 l_int32 w, h, d, thresh;
861 PIX *pixt, *pixd;
862
863 PROCNAME("pixBlockrank");
864
865 if (!pixs)
866 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
867 pixGetDimensions(pixs, &w, &h, &d);
868 if (d != 1)
869 return (PIX *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
870 if (rank < 0.0 || rank > 1.0)
871 return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", procName, NULL);
872 if (rank == 0.0) {
873 pixd = pixCreateTemplate(pixs);
874 pixSetAll(pixd);
875 }
876 if (wc < 0) wc = 0;
877 if (hc < 0) hc = 0;
878 if (w < 2 * wc + 1 || h < 2 * hc + 1) {
879 wc = L_MIN(wc, (w - 1) / 2);
880 hc = L_MIN(hc, (h - 1) / 2);
881 L_WARNING("kernel too large; reducing!", procName);
882 L_INFO_INT2("wc = %d, hc = %d", procName, wc, hc);
883 }
884 if (wc == 0 && hc == 0)
885 return pixCopy(NULL, pixs);
886
887 if ((pixt = pixBlocksum(pixs, pixacc, wc, hc)) == NULL)
888 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
889
890 /* 1 bpp block rank filter output.
891 * Must invert because threshold gives 1 for values < thresh,
892 * but we need a 1 if the value is >= thresh. */
893 thresh = (l_int32)(255. * rank);
894 pixd = pixThresholdToBinary(pixt, thresh);
895 pixInvert(pixd, pixd);
896 pixDestroy(&pixt);
897 return pixd;
898 }
899
900
901 /*!
902 * pixBlocksum()
903 *
904 * Input: pixs (1 bpp)
905 * accum pix (<optional> 32 bpp)
906 * wc, hc (half width/height of block sum/rank kernel)
907 * Return: pixd (8 bpp)
908 *
909 * Notes:
910 * (1) If accum pix is null, make one and destroy it before
911 * returning; otherwise, just use the input accum pix
912 * (2) The full width and height of the convolution kernel
913 * are (2 * wc + 1) and (2 * hc + 1)
914 * (3) Use of wc = hc = 1, followed by pixInvert() on the
915 * 8 bpp result, gives a nice anti-aliased, and somewhat
916 * darkened, result on text.
917 * (4) Require that w >= 2 * wc + 1 and h >= 2 * hc + 1,
918 * where (w,h) are the dimensions of pixs.
919 * (5) Returns in each dest pixel the sum of all src pixels
920 * that are within a block of size of the kernel, centered
921 * on the dest pixel. This sum is the number of src ON
922 * pixels in the block at each location, normalized to 255
923 * for a block containing all ON pixels. For pixels near
924 * the boundary, where the block is not entirely contained
925 * within the image, we then multiply by a second normalization
926 * factor that is greater than one, so that all results
927 * are normalized by the number of participating pixels
928 * within the block.
929 */
930 PIX *
pixBlocksum(PIX * pixs,PIX * pixacc,l_int32 wc,l_int32 hc)931 pixBlocksum(PIX *pixs,
932 PIX *pixacc,
933 l_int32 wc,
934 l_int32 hc)
935 {
936 l_int32 w, h, d, wplt, wpld;
937 l_uint32 *datat, *datad;
938 PIX *pixt, *pixd;
939
940 PROCNAME("pixBlocksum");
941
942 if (!pixs)
943 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
944 pixGetDimensions(pixs, &w, &h, &d);
945 if (d != 1)
946 return (PIX *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
947 if (wc < 0) wc = 0;
948 if (hc < 0) hc = 0;
949 if (w < 2 * wc + 1 || h < 2 * hc + 1) {
950 wc = L_MIN(wc, (w - 1) / 2);
951 hc = L_MIN(hc, (h - 1) / 2);
952 L_WARNING("kernel too large; reducing!", procName);
953 L_INFO_INT2("wc = %d, hc = %d", procName, wc, hc);
954 }
955 if (wc == 0 && hc == 0)
956 return pixCopy(NULL, pixs);
957
958 if (pixacc) {
959 if (pixGetDepth(pixacc) != 32)
960 return (PIX *)ERROR_PTR("pixacc not 32 bpp", procName, NULL);
961 pixt = pixClone(pixacc);
962 }
963 else {
964 if ((pixt = pixBlockconvAccum(pixs)) == NULL)
965 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
966 }
967
968 /* 8 bpp block sum output */
969 if ((pixd = pixCreate(w, h, 8)) == NULL)
970 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
971 pixCopyResolution(pixd, pixs);
972
973 wpld = pixGetWpl(pixd);
974 wplt = pixGetWpl(pixt);
975 datad = pixGetData(pixd);
976 datat = pixGetData(pixt);
977 blocksumLow(datad, w, h, wpld, datat, wplt, wc, hc);
978
979 pixDestroy(&pixt);
980 return pixd;
981 }
982
983
984 /*----------------------------------------------------------------------*
985 * Woodfill transform *
986 *----------------------------------------------------------------------*/
987 /*!
988 * pixWoodfillTransform()
989 *
990 * Input: pixs (8 bpp)
991 * halfsize (of square over which neighbors are averaged)
992 * accum pix (<optional> 32 bpp)
993 * Return: pixd (1 bpp)
994 *
995 * Notes:
996 * (1) The Woodfill transform compares each pixel against
997 * the average of its neighbors (in a square of odd
998 * dimension centered on the pixel). If the pixel is
999 * greater than the average of its neighbors, the output
1000 * pixel value is 1; otherwise it is 0.
1001 * (2) This can be used as an encoding for an image that is
1002 * fairly robust against slow illumination changes, with
1003 * applications in image comparison and mosaicing.
1004 * (3) The size of the convolution kernel is (2 * halfsize + 1)
1005 * on a side. The halfsize parameter must be >= 1.
1006 * (4) If accum pix is null, make one, use it, and destroy it
1007 * before returning; otherwise, just use the input accum pix
1008 */
1009 PIX *
pixWoodfillTransform(PIX * pixs,l_int32 halfsize,PIX * pixacc)1010 pixWoodfillTransform(PIX *pixs,
1011 l_int32 halfsize,
1012 PIX *pixacc)
1013 {
1014 l_int32 i, j, w, h, wpls, wplv, wpld;
1015 l_int32 vals, valv;
1016 l_uint32 *datas, *datav, *datad, *lines, *linev, *lined;
1017 PIX *pixav, *pixd;
1018
1019 PROCNAME("pixWoodfillTransform");
1020
1021 if (!pixs)
1022 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1023 if (pixGetDepth(pixs) != 8)
1024 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
1025 if (halfsize < 1)
1026 return (PIX *)ERROR_PTR("halfsize must be >= 1", procName, NULL);
1027
1028 /* Get the average of each pixel with its neighbors */
1029 if ((pixav = pixBlockconvGray(pixs, pixacc, halfsize, halfsize))
1030 == NULL)
1031 return (PIX *)ERROR_PTR("pixav not made", procName, NULL);
1032
1033 /* Subtract the pixel from the average, and then compare
1034 * the pixel value with the remaining average */
1035 pixGetDimensions(pixs, &w, &h, NULL);
1036 if ((pixd = pixCreate(w, h, 1)) == NULL)
1037 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
1038 datas = pixGetData(pixs);
1039 datav = pixGetData(pixav);
1040 datad = pixGetData(pixd);
1041 wpls = pixGetWpl(pixs);
1042 wplv = pixGetWpl(pixav);
1043 wpld = pixGetWpl(pixd);
1044 for (i = 0; i < h; i++) {
1045 lines = datas + i * wpls;
1046 linev = datav + i * wplv;
1047 lined = datad + i * wpld;
1048 for (j = 0; j < w; j++) {
1049 vals = GET_DATA_BYTE(lines, j);
1050 valv = GET_DATA_BYTE(linev, j);
1051 if (vals > valv)
1052 SET_DATA_BIT(lined, j);
1053 }
1054 }
1055
1056 pixDestroy(&pixav);
1057 return pixd;
1058 }
1059
1060
1061 /*----------------------------------------------------------------------*
1062 * Generic convolution *
1063 *----------------------------------------------------------------------*/
1064 /*!
1065 * pixConvolve()
1066 *
1067 * Input: pixs (8, 16, 32 bpp; no colormap)
1068 * kernel
1069 * outdepth (of pixd: 8, 16 or 32)
1070 * normflag (1 to normalize kernel to unit sum; 0 otherwise)
1071 * Return: pixd (8, 16 or 32 bpp)
1072 *
1073 * Notes:
1074 * (1) This gives a convolution with an arbitrary kernel.
1075 * (2) The parameter @outdepth determines the depth of the result.
1076 * (3) If normflag == 1, the result is normalized by scaling
1077 * all kernel values for a unit sum. Do not normalize if
1078 * the kernel has null sum, such as a DoG.
1079 * (4) If the kernel is normalized to unit sum, the output values
1080 * can never exceed 255, so an output depth of 8 bpp is sufficient.
1081 * If the kernel is not normalized, it may be necessary to use
1082 * 16 or 32 bpp output to avoid overflow.
1083 * (5) The kernel values can be positive or negative, but the
1084 * result for the convolution can only be stored as a positive
1085 * number. Consequently, if it goes negative, the choices are
1086 * to clip to 0 or take the absolute value. We're choosing
1087 * the former for now. Another possibility would be to output
1088 * a second unsigned image for the negative values.
1089 * (6) This uses a mirrored border to avoid special casing on
1090 * the boundaries.
1091 * (7) The function is slow, running at about 12 machine cycles for
1092 * each pixel-op in the convolution. For example, with a 3 GHz
1093 * cpu, a 1 Mpixel grayscale image, and a kernel with
1094 * (sx * sy) = 25 elements, the convolution takes about 100 msec.
1095 */
1096 PIX *
pixConvolve(PIX * pixs,L_KERNEL * kel,l_int32 outdepth,l_int32 normflag)1097 pixConvolve(PIX *pixs,
1098 L_KERNEL *kel,
1099 l_int32 outdepth,
1100 l_int32 normflag)
1101 {
1102 l_int32 i, j, k, m, w, h, d, sx, sy, cx, cy, wplt, wpld;
1103 l_int32 val;
1104 l_uint32 *datat, *datad, *linet, *lined;
1105 l_float32 sum;
1106 L_KERNEL *keli, *keln;
1107 PIX *pixt, *pixd;
1108
1109 PROCNAME("pixConvolve");
1110
1111 if (!pixs)
1112 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1113 if (pixGetColormap(pixs))
1114 return (PIX *)ERROR_PTR("pixs has colormap", procName, NULL);
1115 pixGetDimensions(pixs, &w, &h, &d);
1116 if (d != 8 && d != 16 && d != 32)
1117 return (PIX *)ERROR_PTR("pixs not 8, 16, or 32 bpp", procName, NULL);
1118 if (!kel)
1119 return (PIX *)ERROR_PTR("kel not defined", procName, NULL);
1120
1121 keli = kernelInvert(kel);
1122 kernelGetParameters(keli, &sy, &sx, &cy, &cx);
1123 if (normflag)
1124 keln = kernelNormalize(keli, 1.0);
1125 else
1126 keln = kernelCopy(keli);
1127
1128 if ((pixt = pixAddMirroredBorder(pixs, cx, sx - cx, cy, sy - cy)) == NULL)
1129 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
1130
1131 pixd = pixCreate(w, h, outdepth);
1132 datat = pixGetData(pixt);
1133 datad = pixGetData(pixd);
1134 wplt = pixGetWpl(pixt);
1135 wpld = pixGetWpl(pixd);
1136 for (i = 0; i < h; i++) {
1137 lined = datad + i * wpld;
1138 for (j = 0; j < w; j++) {
1139 sum = 0.0;
1140 for (k = 0; k < sy; k++) {
1141 linet = datat + (i + k) * wplt;
1142 if (d == 8) {
1143 for (m = 0; m < sx; m++) {
1144 val = GET_DATA_BYTE(linet, j + m);
1145 sum += val * keln->data[k][m];
1146 }
1147 }
1148 else if (d == 16) {
1149 for (m = 0; m < sx; m++) {
1150 val = GET_DATA_TWO_BYTES(linet, j + m);
1151 sum += val * keln->data[k][m];
1152 }
1153 }
1154 else { /* d == 32 */
1155 for (m = 0; m < sx; m++) {
1156 val = *(linet + j + m);
1157 sum += val * keln->data[k][m];
1158 }
1159 }
1160 }
1161 #if 1
1162 if (sum < 0.0) sum = -sum; /* make it non-negative */
1163 #endif
1164 if (outdepth == 8)
1165 SET_DATA_BYTE(lined, j, (l_int32)(sum + 0.5));
1166 else if (outdepth == 16)
1167 SET_DATA_TWO_BYTES(lined, j, (l_int32)(sum + 0.5));
1168 else /* outdepth == 32 */
1169 *(lined + j) = (l_uint32)(sum + 0.5);
1170 }
1171 }
1172
1173 kernelDestroy(&keli);
1174 kernelDestroy(&keln);
1175 pixDestroy(&pixt);
1176 return pixd;
1177 }
1178
1179
1180 /*!
1181 * pixConvolveSep()
1182 *
1183 * Input: pixs (8 bpp)
1184 * kelx (x-dependent kernel)
1185 * kely (y-dependent kernel)
1186 * outdepth (of pixd: 8, 16 or 32)
1187 * normflag (1 to normalize kernel to unit sum; 0 otherwise)
1188 * Return: pixd (8, 16 or 32 bpp)
1189 *
1190 * Notes:
1191 * (1) This does a convolution with a separable kernel that is
1192 * is a sequence of convolutions in x and y. The two
1193 * one-dimensional kernel components must be input separately;
1194 * the full kernel is the product of these components.
1195 * The support for the full kernel is thus a rectangular region.
1196 * (2) The @outdepth and @normflag parameters are used as in
1197 * pixConvolve().
1198 * (3) If the kernel is normalized to unit sum, the output values
1199 * can never exceed 255, so an output depth of 8 bpp is sufficient.
1200 * If the kernel is not normalized, it may be necessary to use
1201 * 16 or 32 bpp output to avoid overflow.
1202 * (4) The kernel values can be positive or negative, but the
1203 * result for the convolution can only be stored as a positive
1204 * number. Consequently, if it goes negative, the choices are
1205 * to clip to 0 or take the absolute value. We're choosing
1206 * the former for now. Another possibility would be to output
1207 * a second unsigned image for the negative values.
1208 * (5) This uses mirrored borders to avoid special casing on
1209 * the boundaries.
1210 */
1211 PIX *
pixConvolveSep(PIX * pixs,L_KERNEL * kelx,L_KERNEL * kely,l_int32 outdepth,l_int32 normflag)1212 pixConvolveSep(PIX *pixs,
1213 L_KERNEL *kelx,
1214 L_KERNEL *kely,
1215 l_int32 outdepth,
1216 l_int32 normflag)
1217 {
1218 l_int32 w, h, d;
1219 L_KERNEL *kelxn, *kelyn;
1220 PIX *pixt, *pixd;
1221
1222 PROCNAME("pixConvolveSep");
1223
1224 if (!pixs)
1225 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1226 pixGetDimensions(pixs, &w, &h, &d);
1227 if (d != 8 && d != 16 && d != 32)
1228 return (PIX *)ERROR_PTR("pixs not 8, 16, or 32 bpp", procName, NULL);
1229 if (!kelx)
1230 return (PIX *)ERROR_PTR("kelx not defined", procName, NULL);
1231 if (!kely)
1232 return (PIX *)ERROR_PTR("kely not defined", procName, NULL);
1233
1234 if (normflag) {
1235 kelxn = kernelNormalize(kelx, 1000.0);
1236 kelyn = kernelNormalize(kely, 0.001);
1237 pixt = pixConvolve(pixs, kelxn, 32, 0);
1238 pixd = pixConvolve(pixt, kelyn, outdepth, 0);
1239 kernelDestroy(&kelxn);
1240 kernelDestroy(&kelyn);
1241 }
1242 else { /* don't normalize */
1243 pixt = pixConvolve(pixs, kelx, 32, 0);
1244 pixd = pixConvolve(pixt, kely, outdepth, 0);
1245 }
1246
1247 pixDestroy(&pixt);
1248 return pixd;
1249 }
1250
1251
1252 /*----------------------------------------------------------------------*
1253 * Generic convolution with float array *
1254 *----------------------------------------------------------------------*/
1255 /*!
1256 * fpixConvolve()
1257 *
1258 * Input: fpixs (32 bit float array)
1259 * kernel
1260 * normflag (1 to normalize kernel to unit sum; 0 otherwise)
1261 * Return: fpixd (32 bit float array)
1262 *
1263 * Notes:
1264 * (1) This gives a float convolution with an arbitrary kernel.
1265 * (2) If normflag == 1, the result is normalized by scaling
1266 * all kernel values for a unit sum. Do not normalize if
1267 * the kernel has null sum, such as a DoG.
1268 * (3) With the FPix, there are no issues about negative
1269 * array or kernel values. The convolution is performed
1270 * with single precision arithmetic.
1271 * (4) This uses a mirrored border to avoid special casing on
1272 * the boundaries.
1273 */
1274 FPIX *
fpixConvolve(FPIX * fpixs,L_KERNEL * kel,l_int32 normflag)1275 fpixConvolve(FPIX *fpixs,
1276 L_KERNEL *kel,
1277 l_int32 normflag)
1278 {
1279 l_int32 i, j, k, m, w, h, sx, sy, cx, cy, wplt, wpld;
1280 l_float32 val;
1281 l_float32 *datat, *datad, *linet, *lined;
1282 l_float32 sum;
1283 L_KERNEL *keli, *keln;
1284 FPIX *fpixt, *fpixd;
1285
1286 PROCNAME("fpixConvolve");
1287
1288 if (!fpixs)
1289 return (FPIX *)ERROR_PTR("fpixs not defined", procName, NULL);
1290 if (!kel)
1291 return (FPIX *)ERROR_PTR("kel not defined", procName, NULL);
1292
1293 keli = kernelInvert(kel);
1294 kernelGetParameters(keli, &sy, &sx, &cy, &cx);
1295 if (normflag)
1296 keln = kernelNormalize(keli, 1.0);
1297 else
1298 keln = kernelCopy(keli);
1299
1300 fpixGetDimensions(fpixs, &w, &h);
1301 fpixt = fpixAddMirroredBorder(fpixs, cx, sx - cx, cy, sy - cy);
1302 if (!fpixt)
1303 return (FPIX *)ERROR_PTR("fpixt not made", procName, NULL);
1304
1305 fpixd = fpixCreate(w, h);
1306 datat = fpixGetData(fpixt);
1307 datad = fpixGetData(fpixd);
1308 wplt = fpixGetWpl(fpixt);
1309 wpld = fpixGetWpl(fpixd);
1310 for (i = 0; i < h; i++) {
1311 lined = datad + i * wpld;
1312 for (j = 0; j < w; j++) {
1313 sum = 0.0;
1314 for (k = 0; k < sy; k++) {
1315 linet = datat + (i + k) * wplt;
1316 for (m = 0; m < sx; m++) {
1317 val = *(linet + j + m);
1318 sum += val * keln->data[k][m];
1319 }
1320 }
1321 *(lined + j) = sum;
1322 }
1323 }
1324
1325 kernelDestroy(&keli);
1326 kernelDestroy(&keln);
1327 fpixDestroy(&fpixt);
1328 return fpixd;
1329 }
1330
1331
1332 /*!
1333 * fpixConvolveSep()
1334 *
1335 * Input: fpixs (32 bit float array)
1336 * kelx (x-dependent kernel)
1337 * kely (y-dependent kernel)
1338 * normflag (1 to normalize kernel to unit sum; 0 otherwise)
1339 * Return: fpixd (32 bit float array)
1340 *
1341 * Notes:
1342 * (1) This does a convolution with a separable kernel that is
1343 * is a sequence of convolutions in x and y. The two
1344 * one-dimensional kernel components must be input separately;
1345 * the full kernel is the product of these components.
1346 * The support for the full kernel is thus a rectangular region.
1347 * (2) The normflag parameter is used as in fpixConvolve().
1348 * (3) This uses mirrored borders to avoid special casing on
1349 * the boundaries.
1350 */
1351 FPIX *
fpixConvolveSep(FPIX * fpixs,L_KERNEL * kelx,L_KERNEL * kely,l_int32 normflag)1352 fpixConvolveSep(FPIX *fpixs,
1353 L_KERNEL *kelx,
1354 L_KERNEL *kely,
1355 l_int32 normflag)
1356 {
1357 L_KERNEL *kelxn, *kelyn;
1358 FPIX *fpixt, *fpixd;
1359
1360 PROCNAME("fpixConvolveSep");
1361
1362 if (!fpixs)
1363 return (FPIX *)ERROR_PTR("pixs not defined", procName, NULL);
1364 if (!kelx)
1365 return (FPIX *)ERROR_PTR("kelx not defined", procName, NULL);
1366 if (!kely)
1367 return (FPIX *)ERROR_PTR("kely not defined", procName, NULL);
1368
1369 if (normflag) {
1370 kelxn = kernelNormalize(kelx, 1.0);
1371 kelyn = kernelNormalize(kely, 1.0);
1372 fpixt = fpixConvolve(fpixs, kelxn, 0);
1373 fpixd = fpixConvolve(fpixt, kelyn, 0);
1374 kernelDestroy(&kelxn);
1375 kernelDestroy(&kelyn);
1376 }
1377 else { /* don't normalize */
1378 fpixt = fpixConvolve(fpixs, kelx, 0);
1379 fpixd = fpixConvolve(fpixt, kely, 0);
1380 }
1381
1382 fpixDestroy(&fpixt);
1383 return fpixd;
1384 }
1385
1386
1387
1388