Coverage Report

Created: 2025-06-13 06:48

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