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 pix4.c
29  * <pre>
30  *
31  *    This file has these operations:
32  *
33  *      (1) Pixel histograms
34  *      (2) Pixel row/column statistics
35  *      (3) Foreground/background estimation
36  *
37  *    Pixel histogram, rank val, averaging and min/max
38  *           NUMA       *pixGetGrayHistogram()
39  *           NUMA       *pixGetGrayHistogramMasked()
40  *           NUMA       *pixGetGrayHistogramInRect()
41  *           NUMAA      *pixGetGrayHistogramTiled()
42  *           l_int32     pixGetColorHistogram()
43  *           l_int32     pixGetColorHistogramMasked()
44  *           NUMA       *pixGetCmapHistogram()
45  *           NUMA       *pixGetCmapHistogramMasked()
46  *           NUMA       *pixGetCmapHistogramInRect()
47  *           l_int32     pixCountRGBColors()
48  *           L_AMAP     *pixGetColorAmapHistogram()
49  *           l_int32     amapGetCountForColor()
50  *           l_int32     pixGetRankValue()
51  *           l_int32     pixGetRankValueMaskedRGB()
52  *           l_int32     pixGetRankValueMasked()
53  *           l_int32     pixGetPixelAverage()
54  *           l_int32     pixGetPixelStats()
55  *           l_int32     pixGetAverageMaskedRGB()
56  *           l_int32     pixGetAverageMasked()
57  *           l_int32     pixGetAverageTiledRGB()
58  *           PIX        *pixGetAverageTiled()
59  *           NUMA       *pixRowStats()
60  *           NUMA       *pixColumnStats()
61  *           l_int32     pixGetRangeValues()
62  *           l_int32     pixGetExtremeValue()
63  *           l_int32     pixGetMaxValueInRect()
64  *           l_int32     pixGetBinnedComponentRange()
65  *           l_int32     pixGetRankColorArray()
66  *           l_int32     pixGetBinnedColor()
67  *           PIX        *pixDisplayColorArray()
68  *           PIX        *pixRankBinByStrip()
69  *
70  *    Pixelwise aligned statistics
71  *           PIX        *pixaGetAlignedStats()
72  *           l_int32     pixaExtractColumnFromEachPix()
73  *           l_int32     pixGetRowStats()
74  *           l_int32     pixGetColumnStats()
75  *           l_int32     pixSetPixelColumn()
76  *
77  *    Foreground/background estimation
78  *           l_int32     pixThresholdForFgBg()
79  *           l_int32     pixSplitDistributionFgBg()
80  * </pre>
81  */
82 
83 #include <string.h>
84 #include <math.h>
85 #include "allheaders.h"
86 
87 
88 /*------------------------------------------------------------------*
89  *                  Pixel histogram and averaging                   *
90  *------------------------------------------------------------------*/
91 /*!
92  * \brief   pixGetGrayHistogram()
93  *
94  * \param[in]    pixs 1, 2, 4, 8, 16 bpp; can be colormapped
95  * \param[in]    factor subsampling factor; integer >= 1
96  * \return  na histogram, or NULL on error
97  *
98  * <pre>
99  * Notes:
100  *      (1) If pixs has a colormap, it is converted to 8 bpp gray.
101  *          If you want a histogram of the colormap indices, use
102  *          pixGetCmapHistogram().
103  *      (2) If pixs does not have a colormap, the output histogram is
104  *          of size 2^d, where d is the depth of pixs.
105  *      (3) Set the subsampling factor > 1 to reduce the amount of computation.
106  * </pre>
107  */
108 NUMA *
pixGetGrayHistogram(PIX * pixs,l_int32 factor)109 pixGetGrayHistogram(PIX     *pixs,
110                     l_int32  factor)
111 {
112 l_int32     i, j, w, h, d, wpl, val, size, count;
113 l_uint32   *data, *line;
114 l_float32  *array;
115 NUMA       *na;
116 PIX        *pixg;
117 
118     PROCNAME("pixGetGrayHistogram");
119 
120     if (!pixs)
121         return (NUMA *)ERROR_PTR("pixs not defined", procName, NULL);
122     d = pixGetDepth(pixs);
123     if (d > 16)
124         return (NUMA *)ERROR_PTR("depth not in {1,2,4,8,16}", procName, NULL);
125     if (factor < 1)
126         return (NUMA *)ERROR_PTR("sampling must be >= 1", procName, NULL);
127 
128     if (pixGetColormap(pixs))
129         pixg = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
130     else
131         pixg = pixClone(pixs);
132 
133     pixGetDimensions(pixg, &w, &h, &d);
134     size = 1 << d;
135     if ((na = numaCreate(size)) == NULL) {
136         pixDestroy(&pixg);
137         return (NUMA *)ERROR_PTR("na not made", procName, NULL);
138     }
139     numaSetCount(na, size);  /* all initialized to 0.0 */
140     array = numaGetFArray(na, L_NOCOPY);
141 
142     if (d == 1) {  /* special case */
143         pixCountPixels(pixg, &count, NULL);
144         array[0] = w * h - count;
145         array[1] = count;
146         pixDestroy(&pixg);
147         return na;
148     }
149 
150     wpl = pixGetWpl(pixg);
151     data = pixGetData(pixg);
152     for (i = 0; i < h; i += factor) {
153         line = data + i * wpl;
154         if (d == 2) {
155             for (j = 0; j < w; j += factor) {
156                 val = GET_DATA_DIBIT(line, j);
157                 array[val] += 1.0;
158             }
159         } else if (d == 4) {
160             for (j = 0; j < w; j += factor) {
161                 val = GET_DATA_QBIT(line, j);
162                 array[val] += 1.0;
163             }
164         } else if (d == 8) {
165             for (j = 0; j < w; j += factor) {
166                 val = GET_DATA_BYTE(line, j);
167                 array[val] += 1.0;
168             }
169         } else {  /* d == 16 */
170             for (j = 0; j < w; j += factor) {
171                 val = GET_DATA_TWO_BYTES(line, j);
172                 array[val] += 1.0;
173             }
174         }
175     }
176 
177     pixDestroy(&pixg);
178     return na;
179 }
180 
181 
182 /*!
183  * \brief   pixGetGrayHistogramMasked()
184  *
185  * \param[in]    pixs 8 bpp, or colormapped
186  * \param[in]    pixm [optional] 1 bpp mask over which histogram is
187  *                    to be computed; use all pixels if null
188  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
189  *                    can be < 0; these values are ignored if pixm is null
190  * \param[in]    factor subsampling factor; integer >= 1
191  * \return  na histogram, or NULL on error
192  *
193  * <pre>
194  * Notes:
195  *      (1) If pixs is cmapped, it is converted to 8 bpp gray.
196  *          If you want a histogram of the colormap indices, use
197  *          pixGetCmapHistogramMasked().
198  *      (2) This always returns a 256-value histogram of pixel values.
199  *      (3) Set the subsampling factor > 1 to reduce the amount of computation.
200  *      (4) Clipping of pixm (if it exists) to pixs is done in the inner loop.
201  *      (5) Input x,y are ignored unless pixm exists.
202  * </pre>
203  */
204 NUMA *
pixGetGrayHistogramMasked(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor)205 pixGetGrayHistogramMasked(PIX        *pixs,
206                           PIX        *pixm,
207                           l_int32     x,
208                           l_int32     y,
209                           l_int32     factor)
210 {
211 l_int32     i, j, w, h, wm, hm, dm, wplg, wplm, val;
212 l_uint32   *datag, *datam, *lineg, *linem;
213 l_float32  *array;
214 NUMA       *na;
215 PIX        *pixg;
216 
217     PROCNAME("pixGetGrayHistogramMasked");
218 
219     if (!pixm)
220         return pixGetGrayHistogram(pixs, factor);
221     if (!pixs)
222         return (NUMA *)ERROR_PTR("pixs not defined", procName, NULL);
223     if (pixGetDepth(pixs) != 8 && !pixGetColormap(pixs))
224         return (NUMA *)ERROR_PTR("pixs neither 8 bpp nor colormapped",
225                                  procName, NULL);
226     pixGetDimensions(pixm, &wm, &hm, &dm);
227     if (dm != 1)
228         return (NUMA *)ERROR_PTR("pixm not 1 bpp", procName, NULL);
229     if (factor < 1)
230         return (NUMA *)ERROR_PTR("sampling must be >= 1", procName, NULL);
231 
232     if ((na = numaCreate(256)) == NULL)
233         return (NUMA *)ERROR_PTR("na not made", procName, NULL);
234     numaSetCount(na, 256);  /* all initialized to 0.0 */
235     array = numaGetFArray(na, L_NOCOPY);
236 
237     if (pixGetColormap(pixs))
238         pixg = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
239     else
240         pixg = pixClone(pixs);
241     pixGetDimensions(pixg, &w, &h, NULL);
242     datag = pixGetData(pixg);
243     wplg = pixGetWpl(pixg);
244     datam = pixGetData(pixm);
245     wplm = pixGetWpl(pixm);
246 
247         /* Generate the histogram */
248     for (i = 0; i < hm; i += factor) {
249         if (y + i < 0 || y + i >= h) continue;
250         lineg = datag + (y + i) * wplg;
251         linem = datam + i * wplm;
252         for (j = 0; j < wm; j += factor) {
253             if (x + j < 0 || x + j >= w) continue;
254             if (GET_DATA_BIT(linem, j)) {
255                 val = GET_DATA_BYTE(lineg, x + j);
256                 array[val] += 1.0;
257             }
258         }
259     }
260 
261     pixDestroy(&pixg);
262     return na;
263 }
264 
265 
266 /*!
267  * \brief   pixGetGrayHistogramInRect()
268  *
269  * \param[in]    pixs 8 bpp, or colormapped
270  * \param[in]    box [optional] over which histogram is to be computed;
271  *                    use full image if NULL
272  * \param[in]    factor subsampling factor; integer >= 1
273  * \return  na histogram, or NULL on error
274  *
275  * <pre>
276  * Notes:
277  *      (1) If pixs is cmapped, it is converted to 8 bpp gray.
278  *          If you want a histogram of the colormap indices, use
279  *          pixGetCmapHistogramInRect().
280  *      (2) This always returns a 256-value histogram of pixel values.
281  *      (3) Set the subsampling %factor > 1 to reduce the amount of computation.
282  * </pre>
283  */
284 NUMA *
pixGetGrayHistogramInRect(PIX * pixs,BOX * box,l_int32 factor)285 pixGetGrayHistogramInRect(PIX     *pixs,
286                           BOX     *box,
287                           l_int32  factor)
288 {
289 l_int32     i, j, bx, by, bw, bh, w, h, wplg, val;
290 l_uint32   *datag, *lineg;
291 l_float32  *array;
292 NUMA       *na;
293 PIX        *pixg;
294 
295     PROCNAME("pixGetGrayHistogramInRect");
296 
297     if (!box)
298         return pixGetGrayHistogram(pixs, factor);
299     if (!pixs)
300         return (NUMA *)ERROR_PTR("pixs not defined", procName, NULL);
301     if (pixGetDepth(pixs) != 8 && !pixGetColormap(pixs))
302         return (NUMA *)ERROR_PTR("pixs neither 8 bpp nor colormapped",
303                                  procName, NULL);
304     if (factor < 1)
305         return (NUMA *)ERROR_PTR("sampling must be >= 1", procName, NULL);
306 
307     if ((na = numaCreate(256)) == NULL)
308         return (NUMA *)ERROR_PTR("na not made", procName, NULL);
309     numaSetCount(na, 256);  /* all initialized to 0.0 */
310     array = numaGetFArray(na, L_NOCOPY);
311 
312     if (pixGetColormap(pixs))
313         pixg = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
314     else
315         pixg = pixClone(pixs);
316     pixGetDimensions(pixg, &w, &h, NULL);
317     datag = pixGetData(pixg);
318     wplg = pixGetWpl(pixg);
319     boxGetGeometry(box, &bx, &by, &bw, &bh);
320 
321         /* Generate the histogram */
322     for (i = 0; i < bh; i += factor) {
323         if (by + i < 0 || by + i >= h) continue;
324         lineg = datag + (by + i) * wplg;
325         for (j = 0; j < bw; j += factor) {
326             if (bx + j < 0 || bx + j >= w) continue;
327             val = GET_DATA_BYTE(lineg, bx + j);
328             array[val] += 1.0;
329         }
330     }
331 
332     pixDestroy(&pixg);
333     return na;
334 }
335 
336 
337 /*!
338  * \brief   pixGetGrayHistogramTiled()
339  *
340  * \param[in]    pixs any depth, colormap OK
341  * \param[in]    factor subsampling factor; integer >= 1
342  * \param[in]    nx, ny tiling; >= 1; typically small
343  * \return  naa set of histograms, or NULL on error
344  *
345  * <pre>
346  * Notes:
347  *      (1) If pixs is cmapped, it is converted to 8 bpp gray.
348  *      (2) This returns a set of 256-value histograms of pixel values.
349  *      (3) Set the subsampling factor > 1 to reduce the amount of computation.
350  * </pre>
351  */
352 NUMAA *
pixGetGrayHistogramTiled(PIX * pixs,l_int32 factor,l_int32 nx,l_int32 ny)353 pixGetGrayHistogramTiled(PIX     *pixs,
354                          l_int32  factor,
355                          l_int32  nx,
356                          l_int32  ny)
357 {
358 l_int32  i, n;
359 NUMA    *na;
360 NUMAA   *naa;
361 PIX     *pix1, *pix2;
362 PIXA    *pixa;
363 
364     PROCNAME("pixGetGrayHistogramTiled");
365 
366     if (!pixs)
367         return (NUMAA *)ERROR_PTR("pixs not defined", procName, NULL);
368     if (factor < 1)
369         return (NUMAA *)ERROR_PTR("sampling must be >= 1", procName, NULL);
370     if (nx < 1 || ny < 1)
371         return (NUMAA *)ERROR_PTR("nx and ny must both be > 0", procName, NULL);
372 
373     n = nx * ny;
374     if ((naa = numaaCreate(n)) == NULL)
375         return (NUMAA *)ERROR_PTR("naa not made", procName, NULL);
376 
377     pix1 = pixConvertTo8(pixs, FALSE);
378     pixa = pixaSplitPix(pix1, nx, ny, 0, 0);
379     for (i = 0; i < n; i++) {
380         pix2 = pixaGetPix(pixa, i, L_CLONE);
381         na = pixGetGrayHistogram(pix2, factor);
382         numaaAddNuma(naa, na, L_INSERT);
383         pixDestroy(&pix2);
384     }
385 
386     pixDestroy(&pix1);
387     pixaDestroy(&pixa);
388     return naa;
389 }
390 
391 
392 /*!
393  * \brief   pixGetColorHistogram()
394  *
395  * \param[in]    pixs rgb or colormapped
396  * \param[in]    factor subsampling factor; integer >= 1
397  * \param[out]   pnar red histogram
398  * \param[out]   pnag green histogram
399  * \param[out]   pnab blue histogram
400  * \return  0 if OK, 1 on error
401  *
402  * <pre>
403  * Notes:
404  *      (1) This generates a set of three 256 entry histograms,
405  *          one for each color component (r,g,b).
406  *      (2) Set the subsampling %factor > 1 to reduce the amount of computation.
407  * </pre>
408  */
409 l_int32
pixGetColorHistogram(PIX * pixs,l_int32 factor,NUMA ** pnar,NUMA ** pnag,NUMA ** pnab)410 pixGetColorHistogram(PIX     *pixs,
411                      l_int32  factor,
412                      NUMA   **pnar,
413                      NUMA   **pnag,
414                      NUMA   **pnab)
415 {
416 l_int32     i, j, w, h, d, wpl, index, rval, gval, bval;
417 l_uint32   *data, *line;
418 l_float32  *rarray, *garray, *barray;
419 NUMA       *nar, *nag, *nab;
420 PIXCMAP    *cmap;
421 
422     PROCNAME("pixGetColorHistogram");
423 
424     if (pnar) *pnar = NULL;
425     if (pnag) *pnag = NULL;
426     if (pnab) *pnab = NULL;
427     if (!pnar || !pnag || !pnab)
428         return ERROR_INT("&nar, &nag, &nab not all defined", procName, 1);
429     if (!pixs)
430         return ERROR_INT("pixs not defined", procName, 1);
431     pixGetDimensions(pixs, &w, &h, &d);
432     cmap = pixGetColormap(pixs);
433     if (cmap && (d != 2 && d != 4 && d != 8))
434         return ERROR_INT("colormap and not 2, 4, or 8 bpp", procName, 1);
435     if (!cmap && d != 32)
436         return ERROR_INT("no colormap and not rgb", procName, 1);
437     if (factor < 1)
438         return ERROR_INT("sampling factor must be >= 1", procName, 1);
439 
440         /* Set up the histogram arrays */
441     nar = numaCreate(256);
442     nag = numaCreate(256);
443     nab = numaCreate(256);
444     numaSetCount(nar, 256);
445     numaSetCount(nag, 256);
446     numaSetCount(nab, 256);
447     rarray = numaGetFArray(nar, L_NOCOPY);
448     garray = numaGetFArray(nag, L_NOCOPY);
449     barray = numaGetFArray(nab, L_NOCOPY);
450     *pnar = nar;
451     *pnag = nag;
452     *pnab = nab;
453 
454         /* Generate the color histograms */
455     data = pixGetData(pixs);
456     wpl = pixGetWpl(pixs);
457     if (cmap) {
458         for (i = 0; i < h; i += factor) {
459             line = data + i * wpl;
460             for (j = 0; j < w; j += factor) {
461                 if (d == 8)
462                     index = GET_DATA_BYTE(line, j);
463                 else if (d == 4)
464                     index = GET_DATA_QBIT(line, j);
465                 else   /* 2 bpp */
466                     index = GET_DATA_DIBIT(line, j);
467                 pixcmapGetColor(cmap, index, &rval, &gval, &bval);
468                 rarray[rval] += 1.0;
469                 garray[gval] += 1.0;
470                 barray[bval] += 1.0;
471             }
472         }
473     } else {  /* 32 bpp rgb */
474         for (i = 0; i < h; i += factor) {
475             line = data + i * wpl;
476             for (j = 0; j < w; j += factor) {
477                 extractRGBValues(line[j], &rval, &gval, &bval);
478                 rarray[rval] += 1.0;
479                 garray[gval] += 1.0;
480                 barray[bval] += 1.0;
481             }
482         }
483     }
484 
485     return 0;
486 }
487 
488 
489 /*!
490  * \brief   pixGetColorHistogramMasked()
491  *
492  * \param[in]    pixs 32 bpp rgb, or colormapped
493  * \param[in]    pixm [optional] 1 bpp mask over which histogram is
494  *                    to be computed; use all pixels if null
495  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
496  *                    can be < 0; these values are ignored if pixm is null
497  * \param[in]    factor subsampling factor; integer >= 1
498  * \param[out]   pnar red histogram
499  * \param[out]   pnag green histogram
500  * \param[out]   pnab blue histogram
501  * \return  0 if OK, 1 on error
502  *
503  * <pre>
504  * Notes:
505  *      (1) This generates a set of three 256 entry histograms,
506  *      (2) Set the subsampling %factor > 1 to reduce the amount of computation.
507  *      (3) Clipping of pixm (if it exists) to pixs is done in the inner loop.
508  *      (4) Input x,y are ignored unless pixm exists.
509  * </pre>
510  */
511 l_int32
pixGetColorHistogramMasked(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor,NUMA ** pnar,NUMA ** pnag,NUMA ** pnab)512 pixGetColorHistogramMasked(PIX        *pixs,
513                            PIX        *pixm,
514                            l_int32     x,
515                            l_int32     y,
516                            l_int32     factor,
517                            NUMA      **pnar,
518                            NUMA      **pnag,
519                            NUMA      **pnab)
520 {
521 l_int32     i, j, w, h, d, wm, hm, dm, wpls, wplm, index, rval, gval, bval;
522 l_uint32   *datas, *datam, *lines, *linem;
523 l_float32  *rarray, *garray, *barray;
524 NUMA       *nar, *nag, *nab;
525 PIXCMAP    *cmap;
526 
527     PROCNAME("pixGetColorHistogramMasked");
528 
529     if (!pixm)
530         return pixGetColorHistogram(pixs, factor, pnar, pnag, pnab);
531 
532     if (pnar) *pnar = NULL;
533     if (pnag) *pnag = NULL;
534     if (pnab) *pnab = NULL;
535     if (!pnar || !pnag || !pnab)
536         return ERROR_INT("&nar, &nag, &nab not all defined", procName, 1);
537     if (!pixs)
538         return ERROR_INT("pixs not defined", procName, 1);
539     pixGetDimensions(pixs, &w, &h, &d);
540     cmap = pixGetColormap(pixs);
541     if (cmap && (d != 2 && d != 4 && d != 8))
542         return ERROR_INT("colormap and not 2, 4, or 8 bpp", procName, 1);
543     if (!cmap && d != 32)
544         return ERROR_INT("no colormap and not rgb", procName, 1);
545     pixGetDimensions(pixm, &wm, &hm, &dm);
546     if (dm != 1)
547         return ERROR_INT("pixm not 1 bpp", procName, 1);
548     if (factor < 1)
549         return ERROR_INT("sampling factor must be >= 1", procName, 1);
550 
551         /* Set up the histogram arrays */
552     nar = numaCreate(256);
553     nag = numaCreate(256);
554     nab = numaCreate(256);
555     numaSetCount(nar, 256);
556     numaSetCount(nag, 256);
557     numaSetCount(nab, 256);
558     rarray = numaGetFArray(nar, L_NOCOPY);
559     garray = numaGetFArray(nag, L_NOCOPY);
560     barray = numaGetFArray(nab, L_NOCOPY);
561     *pnar = nar;
562     *pnag = nag;
563     *pnab = nab;
564 
565         /* Generate the color histograms */
566     datas = pixGetData(pixs);
567     wpls = pixGetWpl(pixs);
568     datam = pixGetData(pixm);
569     wplm = pixGetWpl(pixm);
570     if (cmap) {
571         for (i = 0; i < hm; i += factor) {
572             if (y + i < 0 || y + i >= h) continue;
573             lines = datas + (y + i) * wpls;
574             linem = datam + i * wplm;
575             for (j = 0; j < wm; j += factor) {
576                 if (x + j < 0 || x + j >= w) continue;
577                 if (GET_DATA_BIT(linem, j)) {
578                     if (d == 8)
579                         index = GET_DATA_BYTE(lines, x + j);
580                     else if (d == 4)
581                         index = GET_DATA_QBIT(lines, x + j);
582                     else   /* 2 bpp */
583                         index = GET_DATA_DIBIT(lines, x + j);
584                     pixcmapGetColor(cmap, index, &rval, &gval, &bval);
585                     rarray[rval] += 1.0;
586                     garray[gval] += 1.0;
587                     barray[bval] += 1.0;
588                 }
589             }
590         }
591     } else {  /* 32 bpp rgb */
592         for (i = 0; i < hm; i += factor) {
593             if (y + i < 0 || y + i >= h) continue;
594             lines = datas + (y + i) * wpls;
595             linem = datam + i * wplm;
596             for (j = 0; j < wm; j += factor) {
597                 if (x + j < 0 || x + j >= w) continue;
598                 if (GET_DATA_BIT(linem, j)) {
599                     extractRGBValues(lines[x + j], &rval, &gval, &bval);
600                     rarray[rval] += 1.0;
601                     garray[gval] += 1.0;
602                     barray[bval] += 1.0;
603                 }
604             }
605         }
606     }
607 
608     return 0;
609 }
610 
611 
612 /*!
613  * \brief   pixGetCmapHistogram()
614  *
615  * \param[in]    pixs colormapped: d = 2, 4 or 8
616  * \param[in]    factor subsampling factor; integer >= 1
617  * \return  na histogram of cmap indices, or NULL on error
618  *
619  * <pre>
620  * Notes:
621  *      (1) This generates a histogram of colormap pixel indices,
622  *          and is of size 2^d.
623  *      (2) Set the subsampling %factor > 1 to reduce the amount of computation.
624  * </pre>
625  */
626 NUMA *
pixGetCmapHistogram(PIX * pixs,l_int32 factor)627 pixGetCmapHistogram(PIX     *pixs,
628                     l_int32  factor)
629 {
630 l_int32     i, j, w, h, d, wpl, val, size;
631 l_uint32   *data, *line;
632 l_float32  *array;
633 NUMA       *na;
634 
635     PROCNAME("pixGetCmapHistogram");
636 
637     if (!pixs)
638         return (NUMA *)ERROR_PTR("pixs not defined", procName, NULL);
639     if (pixGetColormap(pixs) == NULL)
640         return (NUMA *)ERROR_PTR("pixs not cmapped", procName, NULL);
641     if (factor < 1)
642         return (NUMA *)ERROR_PTR("sampling must be >= 1", procName, NULL);
643     pixGetDimensions(pixs, &w, &h, &d);
644     if (d != 2 && d != 4 && d != 8)
645         return (NUMA *)ERROR_PTR("d not 2, 4 or 8", procName, NULL);
646 
647     size = 1 << d;
648     if ((na = numaCreate(size)) == NULL)
649         return (NUMA *)ERROR_PTR("na not made", procName, NULL);
650     numaSetCount(na, size);  /* all initialized to 0.0 */
651     array = numaGetFArray(na, L_NOCOPY);
652 
653     wpl = pixGetWpl(pixs);
654     data = pixGetData(pixs);
655     for (i = 0; i < h; i += factor) {
656         line = data + i * wpl;
657         for (j = 0; j < w; j += factor) {
658             if (d == 8)
659                 val = GET_DATA_BYTE(line, j);
660             else if (d == 4)
661                 val = GET_DATA_QBIT(line, j);
662             else  /* d == 2 */
663                 val = GET_DATA_DIBIT(line, j);
664             array[val] += 1.0;
665         }
666     }
667 
668     return na;
669 }
670 
671 
672 /*!
673  * \brief   pixGetCmapHistogramMasked()
674  *
675  * \param[in]    pixs colormapped: d = 2, 4 or 8
676  * \param[in]    pixm [optional] 1 bpp mask over which histogram is
677  *                    to be computed; use all pixels if null
678  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
679  *                    can be < 0; these values are ignored if pixm is null
680  * \param[in]    factor subsampling factor; integer >= 1
681  * \return  na histogram, or NULL on error
682  *
683  * <pre>
684  * Notes:
685  *      (1) This generates a histogram of colormap pixel indices,
686  *          and is of size 2^d.
687  *      (2) Set the subsampling %factor > 1 to reduce the amount of computation.
688  *      (3) Clipping of pixm to pixs is done in the inner loop.
689  * </pre>
690  */
691 NUMA *
pixGetCmapHistogramMasked(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor)692 pixGetCmapHistogramMasked(PIX     *pixs,
693                           PIX     *pixm,
694                           l_int32  x,
695                           l_int32  y,
696                           l_int32  factor)
697 {
698 l_int32     i, j, w, h, d, wm, hm, dm, wpls, wplm, val, size;
699 l_uint32   *datas, *datam, *lines, *linem;
700 l_float32  *array;
701 NUMA       *na;
702 
703     PROCNAME("pixGetCmapHistogramMasked");
704 
705     if (!pixm)
706         return pixGetCmapHistogram(pixs, factor);
707 
708     if (!pixs)
709         return (NUMA *)ERROR_PTR("pixs not defined", procName, NULL);
710     if (pixGetColormap(pixs) == NULL)
711         return (NUMA *)ERROR_PTR("pixs not cmapped", procName, NULL);
712     pixGetDimensions(pixm, &wm, &hm, &dm);
713     if (dm != 1)
714         return (NUMA *)ERROR_PTR("pixm not 1 bpp", procName, NULL);
715     if (factor < 1)
716         return (NUMA *)ERROR_PTR("sampling must be >= 1", procName, NULL);
717     pixGetDimensions(pixs, &w, &h, &d);
718     if (d != 2 && d != 4 && d != 8)
719         return (NUMA *)ERROR_PTR("d not 2, 4 or 8", procName, NULL);
720 
721     size = 1 << d;
722     if ((na = numaCreate(size)) == NULL)
723         return (NUMA *)ERROR_PTR("na not made", procName, NULL);
724     numaSetCount(na, size);  /* all initialized to 0.0 */
725     array = numaGetFArray(na, L_NOCOPY);
726 
727     datas = pixGetData(pixs);
728     wpls = pixGetWpl(pixs);
729     datam = pixGetData(pixm);
730     wplm = pixGetWpl(pixm);
731 
732     for (i = 0; i < hm; i += factor) {
733         if (y + i < 0 || y + i >= h) continue;
734         lines = datas + (y + i) * wpls;
735         linem = datam + i * wplm;
736         for (j = 0; j < wm; j += factor) {
737             if (x + j < 0 || x + j >= w) continue;
738             if (GET_DATA_BIT(linem, j)) {
739                 if (d == 8)
740                     val = GET_DATA_BYTE(lines, x + j);
741                 else if (d == 4)
742                     val = GET_DATA_QBIT(lines, x + j);
743                 else  /* d == 2 */
744                     val = GET_DATA_DIBIT(lines, x + j);
745                 array[val] += 1.0;
746             }
747         }
748     }
749 
750     return na;
751 }
752 
753 
754 /*!
755  * \brief   pixGetCmapHistogramInRect()
756  *
757  * \param[in]    pixs colormapped: d = 2, 4 or 8
758  * \param[in]    box [optional] over which histogram is to be computed;
759  *                    use full image if NULL
760  * \param[in]    factor subsampling factor; integer >= 1
761  * \return  na histogram, or NULL on error
762  *
763  * <pre>
764  * Notes:
765  *      (1) This generates a histogram of colormap pixel indices,
766  *          and is of size 2^d.
767  *      (2) Set the subsampling %factor > 1 to reduce the amount of computation.
768  *      (3) Clipping to the box is done in the inner loop.
769  * </pre>
770  */
771 NUMA *
pixGetCmapHistogramInRect(PIX * pixs,BOX * box,l_int32 factor)772 pixGetCmapHistogramInRect(PIX     *pixs,
773                           BOX     *box,
774                           l_int32  factor)
775 {
776 l_int32     i, j, bx, by, bw, bh, w, h, d, wpls, val, size;
777 l_uint32   *datas, *lines;
778 l_float32  *array;
779 NUMA       *na;
780 
781     PROCNAME("pixGetCmapHistogramInRect");
782 
783     if (!box)
784         return pixGetCmapHistogram(pixs, factor);
785     if (!pixs)
786         return (NUMA *)ERROR_PTR("pixs not defined", procName, NULL);
787     if (pixGetColormap(pixs) == NULL)
788         return (NUMA *)ERROR_PTR("pixs not cmapped", procName, NULL);
789     if (factor < 1)
790         return (NUMA *)ERROR_PTR("sampling must be >= 1", procName, NULL);
791     pixGetDimensions(pixs, &w, &h, &d);
792     if (d != 2 && d != 4 && d != 8)
793         return (NUMA *)ERROR_PTR("d not 2, 4 or 8", procName, NULL);
794 
795     size = 1 << d;
796     if ((na = numaCreate(size)) == NULL)
797         return (NUMA *)ERROR_PTR("na not made", procName, NULL);
798     numaSetCount(na, size);  /* all initialized to 0.0 */
799     array = numaGetFArray(na, L_NOCOPY);
800 
801     datas = pixGetData(pixs);
802     wpls = pixGetWpl(pixs);
803     boxGetGeometry(box, &bx, &by, &bw, &bh);
804 
805     for (i = 0; i < bh; i += factor) {
806         if (by + i < 0 || by + i >= h) continue;
807         lines = datas + (by + i) * wpls;
808         for (j = 0; j < bw; j += factor) {
809             if (bx + j < 0 || bx + j >= w) continue;
810             if (d == 8)
811                 val = GET_DATA_BYTE(lines, bx + j);
812             else if (d == 4)
813                 val = GET_DATA_QBIT(lines, bx + j);
814             else  /* d == 2 */
815                 val = GET_DATA_DIBIT(lines, bx + j);
816             array[val] += 1.0;
817         }
818     }
819 
820     return na;
821 }
822 
823 
824 /*!
825  * \brief   pixCountRGBColors()
826  *
827  * \param[in]    pixs    rgb or rgba
828  * \return  ncolors, or -1 on error
829  */
830 l_int32
pixCountRGBColors(PIX * pixs)831 pixCountRGBColors(PIX  *pixs)
832 {
833 l_int32  ncolors;
834 L_AMAP  *amap;
835 
836     PROCNAME("pixCountRGBColors");
837 
838     if (!pixs || pixGetDepth(pixs) != 32)
839         return ERROR_INT("pixs not defined or not 32 bpp", procName, -1);
840     amap = pixGetColorAmapHistogram(pixs, 1);
841     ncolors = l_amapSize(amap);
842     l_amapDestroy(&amap);
843     return ncolors;
844 }
845 
846 
847 /*!
848  * \brief   pixGetColorAmapHistogram()
849  *
850  * \param[in]    pixs    rgb or rgba
851  * \param[in]    factor  subsampling factor; integer >= 1
852  * \return  amap, or NULL on error
853  *
854  * <pre>
855  * Notes:
856  *      (1) This generates an ordered map from pixel value to histogram count.
857  *      (2) Use amapGetCountForColor() to use the map to look up a count.
858  * </pre>
859  */
860 L_AMAP  *
pixGetColorAmapHistogram(PIX * pixs,l_int32 factor)861 pixGetColorAmapHistogram(PIX     *pixs,
862                          l_int32  factor)
863 {
864 l_int32    i, j, w, h, wpl;
865 l_uint32  *data, *line;
866 L_AMAP    *amap;
867 RB_TYPE    key, value;
868 RB_TYPE   *pval;
869 
870     PROCNAME("pixGetColorAmapHistogram");
871 
872     if (!pixs)
873         return (L_AMAP *)ERROR_PTR("pixs not defined", procName, NULL);
874     if (pixGetDepth(pixs) != 32)
875         return (L_AMAP *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
876     pixGetDimensions(pixs, &w, &h, NULL);
877     data = pixGetData(pixs);
878     wpl = pixGetWpl(pixs);
879     amap = l_amapCreate(L_UINT_TYPE);
880     for (i = 0; i < h; i += factor) {
881         line = data + i * wpl;
882         for (j = 0; j < w; j += factor) {
883             key.utype = line[j];
884             pval = l_amapFind(amap, key);
885             if (!pval)
886                 value.itype = 1;
887             else
888                 value.itype = 1 + pval->itype;
889             l_amapInsert(amap, key, value);
890         }
891     }
892 
893     return amap;
894 }
895 
896 
897 /*!
898  * \brief   amapGetCountForColor()
899  *
900  * \param[in]    amap    map from pixel value to count
901  * \param[in]    val     rgb or rgba pixel value
902  * \return  count, or -1 on error
903  *
904  * <pre>
905  * Notes:
906  *      (1) The ordered map is made by pixGetColorAmapHistogram().
907  * </pre>
908  */
909 l_int32
amapGetCountForColor(L_AMAP * amap,l_uint32 val)910 amapGetCountForColor(L_AMAP   *amap,
911                      l_uint32  val)
912 {
913 RB_TYPE   key;
914 RB_TYPE  *pval;
915 
916     PROCNAME("amapGetCountForColor");
917 
918     if (!amap)
919         return ERROR_INT("amap not defined", procName, -1);
920 
921     key.utype = val;
922     pval = l_amapFind(amap, key);
923     return (pval) ? pval->itype : 0;
924 }
925 
926 
927 /*!
928  * \brief   pixGetRankValue()
929  *
930  * \param[in]    pixs 8 bpp, 32 bpp or colormapped
931  * \param[in]    factor subsampling factor; integer >= 1
932  * \param[in]    rank between 0.0 and 1.0; 1.0 is brightest, 0.0 is darkest
933  * \param[out]   pvalue pixel value corresponding to input rank
934  * \return  0 if OK, 1 on error
935  *
936  * <pre>
937  * Notes:
938  *      (1) Simple function to get rank values of an image.
939  *          For a color image, the median value (rank = 0.5) can be
940  *          used to linearly remap the colors based on the median
941  *          of a target image, using pixLinearMapToTargetColor().
942  * </pre>
943  */
944 l_int32
pixGetRankValue(PIX * pixs,l_int32 factor,l_float32 rank,l_uint32 * pvalue)945 pixGetRankValue(PIX       *pixs,
946                 l_int32    factor,
947                 l_float32  rank,
948                 l_uint32  *pvalue)
949 {
950 l_int32    d;
951 l_float32  val, rval, gval, bval;
952 PIX       *pixt;
953 PIXCMAP   *cmap;
954 
955     PROCNAME("pixGetRankValue");
956 
957     if (!pvalue)
958         return ERROR_INT("&value not defined", procName, 1);
959     *pvalue = 0;
960     if (!pixs)
961         return ERROR_INT("pixs not defined", procName, 1);
962     d = pixGetDepth(pixs);
963     cmap = pixGetColormap(pixs);
964     if (d != 8 && d != 32 && !cmap)
965         return ERROR_INT("pixs not 8 or 32 bpp, or cmapped", procName, 1);
966     if (cmap)
967         pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
968     else
969         pixt = pixClone(pixs);
970     d = pixGetDepth(pixt);
971 
972     if (d == 8) {
973         pixGetRankValueMasked(pixt, NULL, 0, 0, factor, rank, &val, NULL);
974         *pvalue = lept_roundftoi(val);
975     } else {
976         pixGetRankValueMaskedRGB(pixt, NULL, 0, 0, factor, rank,
977                                  &rval, &gval, &bval);
978         composeRGBPixel(lept_roundftoi(rval), lept_roundftoi(gval),
979                         lept_roundftoi(bval), pvalue);
980     }
981 
982     pixDestroy(&pixt);
983     return 0;
984 }
985 
986 
987 /*!
988  * \brief   pixGetRankValueMaskedRGB()
989  *
990  * \param[in]    pixs 32 bpp
991  * \param[in]    pixm [optional] 1 bpp mask over which rank val is to be taken;
992  *                    use all pixels if null
993  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
994  *                    can be < 0; these values are ignored if pixm is null
995  * \param[in]    factor subsampling factor; integer >= 1
996  * \param[in]    rank between 0.0 and 1.0; 1.0 is brightest, 0.0 is darkest
997  * \param[out]   prval [optional] red component val for input rank
998  * \param[out]   pgval [optional] green component val for input rank
999  * \param[out]   pbval [optional] blue component val for input rank
1000  * \return  0 if OK, 1 on error
1001  *
1002  * <pre>
1003  * Notes:
1004  *      (1) Computes the rank component values of pixels in pixs that
1005  *          are under the fg of the optional mask.  If the mask is null, it
1006  *          computes the average of the pixels in pixs.
1007  *      (2) Set the subsampling %factor > 1 to reduce the amount of
1008  *          computation.
1009  *      (4) Input x,y are ignored unless pixm exists.
1010  *      (5) The rank must be in [0.0 ... 1.0], where the brightest pixel
1011  *          has rank 1.0.  For the median pixel value, use 0.5.
1012  * </pre>
1013  */
1014 l_int32
pixGetRankValueMaskedRGB(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor,l_float32 rank,l_float32 * prval,l_float32 * pgval,l_float32 * pbval)1015 pixGetRankValueMaskedRGB(PIX        *pixs,
1016                          PIX        *pixm,
1017                          l_int32     x,
1018                          l_int32     y,
1019                          l_int32     factor,
1020                          l_float32   rank,
1021                          l_float32  *prval,
1022                          l_float32  *pgval,
1023                          l_float32  *pbval)
1024 {
1025 l_float32  scale;
1026 PIX       *pixmt, *pixt;
1027 
1028     PROCNAME("pixGetRankValueMaskedRGB");
1029 
1030     if (prval) *prval = 0.0;
1031     if (pgval) *pgval = 0.0;
1032     if (pbval) *pbval = 0.0;
1033     if (!prval && !pgval && !pbval)
1034         return ERROR_INT("no results requested", procName, 1);
1035     if (!pixs)
1036         return ERROR_INT("pixs not defined", procName, 1);
1037     if (pixGetDepth(pixs) != 32)
1038         return ERROR_INT("pixs not 32 bpp", procName, 1);
1039     if (pixm && pixGetDepth(pixm) != 1)
1040         return ERROR_INT("pixm not 1 bpp", procName, 1);
1041     if (factor < 1)
1042         return ERROR_INT("sampling factor must be >= 1", procName, 1);
1043     if (rank < 0.0 || rank > 1.0)
1044         return ERROR_INT("rank not in [0.0 ... 1.0]", procName, 1);
1045 
1046     pixmt = NULL;
1047     if (pixm) {
1048         scale = 1.0 / (l_float32)factor;
1049         pixmt = pixScale(pixm, scale, scale);
1050     }
1051     if (prval) {
1052         pixt = pixScaleRGBToGrayFast(pixs, factor, COLOR_RED);
1053         pixGetRankValueMasked(pixt, pixmt, x / factor, y / factor,
1054                               factor, rank, prval, NULL);
1055         pixDestroy(&pixt);
1056     }
1057     if (pgval) {
1058         pixt = pixScaleRGBToGrayFast(pixs, factor, COLOR_GREEN);
1059         pixGetRankValueMasked(pixt, pixmt, x / factor, y / factor,
1060                               factor, rank, pgval, NULL);
1061         pixDestroy(&pixt);
1062     }
1063     if (pbval) {
1064         pixt = pixScaleRGBToGrayFast(pixs, factor, COLOR_BLUE);
1065         pixGetRankValueMasked(pixt, pixmt, x / factor, y / factor,
1066                               factor, rank, pbval, NULL);
1067         pixDestroy(&pixt);
1068     }
1069     pixDestroy(&pixmt);
1070     return 0;
1071 }
1072 
1073 
1074 /*!
1075  * \brief   pixGetRankValueMasked()
1076  *
1077  * \param[in]    pixs 8 bpp, or colormapped
1078  * \param[in]    pixm [optional] 1 bpp mask over which rank val is to be taken;
1079  *                    use all pixels if null
1080  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
1081  *                    can be < 0; these values are ignored if pixm is null
1082  * \param[in]    factor subsampling factor; integer >= 1
1083  * \param[in]    rank between 0.0 and 1.0; 1.0 is brightest, 0.0 is darkest
1084  * \param[out]   pval pixel value corresponding to input rank
1085  * \param[out]   pna [optional] of histogram
1086  * \return  0 if OK, 1 on error
1087  *
1088  * <pre>
1089  * Notes:
1090  *      (1) Computes the rank value of pixels in pixs that are under
1091  *          the fg of the optional mask.  If the mask is null, it
1092  *          computes the average of the pixels in pixs.
1093  *      (2) Set the subsampling %factor > 1 to reduce the amount of
1094  *          computation.
1095  *      (3) Clipping of pixm (if it exists) to pixs is done in the inner loop.
1096  *      (4) Input x,y are ignored unless pixm exists.
1097  *      (5) The rank must be in [0.0 ... 1.0], where the brightest pixel
1098  *          has rank 1.0.  For the median pixel value, use 0.5.
1099  *      (6) The histogram can optionally be returned, so that other rank
1100  *          values can be extracted without recomputing the histogram.
1101  *          In that case, just use
1102  *              numaHistogramGetValFromRank(na, rank, &val);
1103  *          on the returned Numa for additional rank values.
1104  * </pre>
1105  */
1106 l_int32
pixGetRankValueMasked(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor,l_float32 rank,l_float32 * pval,NUMA ** pna)1107 pixGetRankValueMasked(PIX        *pixs,
1108                       PIX        *pixm,
1109                       l_int32     x,
1110                       l_int32     y,
1111                       l_int32     factor,
1112                       l_float32   rank,
1113                       l_float32  *pval,
1114                       NUMA      **pna)
1115 {
1116 NUMA  *na;
1117 
1118     PROCNAME("pixGetRankValueMasked");
1119 
1120     if (pna) *pna = NULL;
1121     if (!pval)
1122         return ERROR_INT("&val not defined", procName, 1);
1123     *pval = 0.0;
1124     if (!pixs)
1125         return ERROR_INT("pixs not defined", procName, 1);
1126     if (pixGetDepth(pixs) != 8 && !pixGetColormap(pixs))
1127         return ERROR_INT("pixs neither 8 bpp nor colormapped", procName, 1);
1128     if (pixm && pixGetDepth(pixm) != 1)
1129         return ERROR_INT("pixm not 1 bpp", procName, 1);
1130     if (factor < 1)
1131         return ERROR_INT("sampling factor must be >= 1", procName, 1);
1132     if (rank < 0.0 || rank > 1.0)
1133         return ERROR_INT("rank not in [0.0 ... 1.0]", procName, 1);
1134 
1135     if ((na = pixGetGrayHistogramMasked(pixs, pixm, x, y, factor)) == NULL)
1136         return ERROR_INT("na not made", procName, 1);
1137     numaHistogramGetValFromRank(na, rank, pval);
1138     if (pna)
1139         *pna = na;
1140     else
1141         numaDestroy(&na);
1142 
1143     return 0;
1144 }
1145 
1146 
1147 /*!
1148  * \brief   pixGetPixelAverage()
1149  *
1150  * \param[in]    pixs 8 or 32 bpp, or colormapped
1151  * \param[in]    pixm [optional] 1 bpp mask over which average is to be taken;
1152  *                    use all pixels if null
1153  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
1154  *                    can be < 0
1155  * \param[in]    factor subsampling factor; >= 1
1156  * \param[out]   pval  average pixel value
1157  * \return  0 if OK, 1 on error
1158  *
1159  * <pre>
1160  * Notes:
1161  *      (1) For rgb pix, this is a more direct computation of the
1162  *          average value of the pixels in %pixs that are under the
1163  *          mask %pixm. It is faster than pixGetPixelStats(), which
1164  *          calls pixGetAverageMaskedRGB() and has the overhead of
1165  *          generating a temporary pix of each of the three components;
1166  *          this can take most of the time if %factor > 1.
1167  *      (2) If %pixm is null, this gives the average value of all
1168  *          pixels in %pixs.  The returned value is an integer.
1169  *      (3) For color %pixs, the returned pixel value is in the standard
1170  *          uint32 RGBA packing.
1171  *      (4) Clipping of pixm (if it exists) to pixs is done in the inner loop.
1172  *      (5) Input x,y are ignored if %pixm does not exist.
1173  * </pre>
1174  */
1175 l_int32
pixGetPixelAverage(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor,l_uint32 * pval)1176 pixGetPixelAverage(PIX       *pixs,
1177                    PIX       *pixm,
1178                    l_int32    x,
1179                    l_int32    y,
1180                    l_int32    factor,
1181                    l_uint32  *pval)
1182 {
1183 l_int32    i, j, w, h, d, wm, hm, wpl1, wplm, val, rval, gval, bval, count;
1184 l_uint32  *data1, *datam, *line1, *linem;
1185 l_float64  sum, rsum, gsum, bsum;
1186 PIX       *pix1;
1187 
1188     PROCNAME("pixGetPixelAverage");
1189 
1190     if (!pval)
1191         return ERROR_INT("&val not defined", procName, 1);
1192     *pval = 0;
1193     if (!pixs)
1194         return ERROR_INT("pixs not defined", procName, 1);
1195     d = pixGetDepth(pixs);
1196     if (d != 32 && !pixGetColormap(pixs))
1197         return ERROR_INT("pixs not rgb or colormapped", procName, 1);
1198     if (pixm && pixGetDepth(pixm) != 1)
1199         return ERROR_INT("pixm not 1 bpp", procName, 1);
1200     if (factor < 1)
1201         return ERROR_INT("sampling factor must be >= 1", procName, 1);
1202 
1203     if (pixGetColormap(pixs))
1204         pix1 = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
1205     else
1206         pix1 = pixClone(pixs);
1207     pixGetDimensions(pix1, &w, &h, &d);
1208     if (d == 1) {
1209         pixDestroy(&pix1);
1210         return ERROR_INT("pix1 is just 1 bpp", procName, 1);
1211     }
1212     data1 = pixGetData(pix1);
1213     wpl1 = pixGetWpl(pix1);
1214 
1215     sum = rsum = gsum = bsum = 0.0;
1216     count = 0;
1217     if (!pixm) {
1218         for (i = 0; i < h; i += factor) {
1219             line1 = data1 + i * wpl1;
1220             for (j = 0; j < w; j += factor) {
1221                 if (d == 8) {
1222                     val = GET_DATA_BYTE(line1, j);
1223                     sum += val;
1224                 } else {  /* rgb */
1225                     extractRGBValues(*(line1 + j), &rval, &gval, &bval);
1226                     rsum += rval;
1227                     gsum += gval;
1228                     bsum += bval;
1229                 }
1230                 count++;
1231             }
1232         }
1233     } else {  /* masked */
1234         pixGetDimensions(pixm, &wm, &hm, NULL);
1235         datam = pixGetData(pixm);
1236         wplm = pixGetWpl(pixm);
1237         for (i = 0; i < hm; i += factor) {
1238             if (y + i < 0 || y + i >= h) continue;
1239             line1 = data1 + (y + i) * wpl1;
1240             linem = datam + i * wplm;
1241             for (j = 0; j < wm; j += factor) {
1242                 if (x + j < 0 || x + j >= w) continue;
1243                 if (GET_DATA_BIT(linem, j)) {
1244                     if (d == 8) {
1245                         val = GET_DATA_BYTE(line1, x + j);
1246                         sum += val;
1247                     } else {  /* rgb */
1248                         extractRGBValues(*(line1 + x + j), &rval, &gval, &bval);
1249                         rsum += rval;
1250                         gsum += gval;
1251                         bsum += bval;
1252                     }
1253                     count++;
1254                 }
1255             }
1256         }
1257     }
1258 
1259     pixDestroy(&pix1);
1260     if (count == 0)
1261         return ERROR_INT("no pixels sampled", procName, 1);
1262     if (d == 8) {
1263         *pval = (l_uint32)((l_float64)sum / (l_float64)count);
1264     } else {  /* d == 32 */
1265         rval = (l_uint32)((l_float64)rsum / (l_float64)count);
1266         gval = (l_uint32)((l_float64)gsum / (l_float64)count);
1267         bval = (l_uint32)((l_float64)bsum / (l_float64)count);
1268         composeRGBPixel(rval, gval, bval, pval);
1269     }
1270 
1271     return 0;
1272 }
1273 
1274 
1275 /*!
1276  * \brief   pixGetPixelStats()
1277  *
1278  * \param[in]    pixs 8 bpp, 32 bpp or colormapped
1279  * \param[in]    factor subsampling factor; integer >= 1
1280  * \param[in]    type L_MEAN_ABSVAL, L_ROOT_MEAN_SQUARE,
1281  *                    L_STANDARD_DEVIATION, L_VARIANCE
1282  * \param[out]   pvalue pixel value corresponding to input type
1283  * \return  0 if OK, 1 on error
1284  *
1285  * <pre>
1286  * Notes:
1287  *      (1) Simple function to get one of four statistical values of an image.
1288  *      (2) It does not take a mask: it uses the entire image.
1289  *      (3) To get the average pixel value of an RGB image, suggest using
1290  *          pixGetPixelAverage(), which is considerably faster.
1291  * </pre>
1292  */
1293 l_int32
pixGetPixelStats(PIX * pixs,l_int32 factor,l_int32 type,l_uint32 * pvalue)1294 pixGetPixelStats(PIX       *pixs,
1295                  l_int32    factor,
1296                  l_int32    type,
1297                  l_uint32  *pvalue)
1298 {
1299 l_int32    d;
1300 l_float32  val, rval, gval, bval;
1301 PIX       *pixt;
1302 PIXCMAP   *cmap;
1303 
1304     PROCNAME("pixGetPixelStats");
1305 
1306     if (!pvalue)
1307         return ERROR_INT("&value not defined", procName, 1);
1308     *pvalue = 0;
1309     if (!pixs)
1310         return ERROR_INT("pixs not defined", procName, 1);
1311     d = pixGetDepth(pixs);
1312     cmap = pixGetColormap(pixs);
1313     if (d != 8 && d != 32 && !cmap)
1314         return ERROR_INT("pixs not 8 or 32 bpp, or cmapped", procName, 1);
1315     if (cmap)
1316         pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
1317     else
1318         pixt = pixClone(pixs);
1319     d = pixGetDepth(pixt);
1320 
1321     if (d == 8) {
1322         pixGetAverageMasked(pixt, NULL, 0, 0, factor, type, &val);
1323         *pvalue = lept_roundftoi(val);
1324     } else {
1325         pixGetAverageMaskedRGB(pixt, NULL, 0, 0, factor, type,
1326                                &rval, &gval, &bval);
1327         composeRGBPixel(lept_roundftoi(rval), lept_roundftoi(gval),
1328                         lept_roundftoi(bval), pvalue);
1329     }
1330 
1331     pixDestroy(&pixt);
1332     return 0;
1333 }
1334 
1335 
1336 /*!
1337  * \brief   pixGetAverageMaskedRGB()
1338  *
1339  * \param[in]    pixs 32 bpp, or colormapped
1340  * \param[in]    pixm [optional] 1 bpp mask over which average is to be taken;
1341  *                    use all pixels if null
1342  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
1343  *                    can be < 0
1344  * \param[in]    factor subsampling factor; >= 1
1345  * \param[in]    type L_MEAN_ABSVAL, L_ROOT_MEAN_SQUARE,
1346  *                    L_STANDARD_DEVIATION, L_VARIANCE
1347  * \param[out]   prval [optional] measured red value of given 'type'
1348  * \param[out]   pgval [optional] measured green value of given 'type'
1349  * \param[out]   pbval [optional] measured blue value of given 'type'
1350  * \return  0 if OK, 1 on error
1351  *
1352  * <pre>
1353  * Notes:
1354  *      (1) For usage, see pixGetAverageMasked().
1355  *      (2) If there is a colormap, it is removed before the 8 bpp
1356  *          component images are extracted.
1357  *      (3) A better name for this would be: pixGetPixelStatsRGB()
1358  * </pre>
1359  */
1360 l_int32
pixGetAverageMaskedRGB(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor,l_int32 type,l_float32 * prval,l_float32 * pgval,l_float32 * pbval)1361 pixGetAverageMaskedRGB(PIX        *pixs,
1362                        PIX        *pixm,
1363                        l_int32     x,
1364                        l_int32     y,
1365                        l_int32     factor,
1366                        l_int32     type,
1367                        l_float32  *prval,
1368                        l_float32  *pgval,
1369                        l_float32  *pbval)
1370 {
1371 PIX      *pixt;
1372 PIXCMAP  *cmap;
1373 
1374     PROCNAME("pixGetAverageMaskedRGB");
1375 
1376     if (prval) *prval = 0.0;
1377     if (pgval) *pgval = 0.0;
1378     if (pbval) *pbval = 0.0;
1379     if (!prval && !pgval && !pbval)
1380         return ERROR_INT("no values requested", procName, 1);
1381     if (!pixs)
1382         return ERROR_INT("pixs not defined", procName, 1);
1383     cmap = pixGetColormap(pixs);
1384     if (pixGetDepth(pixs) != 32 && !cmap)
1385         return ERROR_INT("pixs neither 32 bpp nor colormapped", procName, 1);
1386     if (pixm && pixGetDepth(pixm) != 1)
1387         return ERROR_INT("pixm not 1 bpp", procName, 1);
1388     if (factor < 1)
1389         return ERROR_INT("sampling factor must be >= 1", procName, 1);
1390     if (type != L_MEAN_ABSVAL && type != L_ROOT_MEAN_SQUARE &&
1391         type != L_STANDARD_DEVIATION && type != L_VARIANCE)
1392         return ERROR_INT("invalid measure type", procName, 1);
1393 
1394     if (prval) {
1395         if (cmap)
1396             pixt = pixGetRGBComponentCmap(pixs, COLOR_RED);
1397         else
1398             pixt = pixGetRGBComponent(pixs, COLOR_RED);
1399         pixGetAverageMasked(pixt, pixm, x, y, factor, type, prval);
1400         pixDestroy(&pixt);
1401     }
1402     if (pgval) {
1403         if (cmap)
1404             pixt = pixGetRGBComponentCmap(pixs, COLOR_GREEN);
1405         else
1406             pixt = pixGetRGBComponent(pixs, COLOR_GREEN);
1407         pixGetAverageMasked(pixt, pixm, x, y, factor, type, pgval);
1408         pixDestroy(&pixt);
1409     }
1410     if (pbval) {
1411         if (cmap)
1412             pixt = pixGetRGBComponentCmap(pixs, COLOR_BLUE);
1413         else
1414             pixt = pixGetRGBComponent(pixs, COLOR_BLUE);
1415         pixGetAverageMasked(pixt, pixm, x, y, factor, type, pbval);
1416         pixDestroy(&pixt);
1417     }
1418 
1419     return 0;
1420 }
1421 
1422 
1423 /*!
1424  * \brief   pixGetAverageMasked()
1425  *
1426  * \param[in]    pixs 8 or 16 bpp, or colormapped
1427  * \param[in]    pixm [optional] 1 bpp mask over which average is to be taken;
1428  *                    use all pixels if null
1429  * \param[in]    x, y UL corner of pixm relative to the UL corner of pixs;
1430  *                    can be < 0
1431  * \param[in]    factor subsampling factor; >= 1
1432  * \param[in]    type L_MEAN_ABSVAL, L_ROOT_MEAN_SQUARE,
1433  *                    L_STANDARD_DEVIATION, L_VARIANCE
1434  * \param[out]   pval measured value of given 'type'
1435  * \return  0 if OK, 1 on error
1436  *
1437  * <pre>
1438  * Notes:
1439  *      (1) Use L_MEAN_ABSVAL to get the average value of pixels in pixs
1440  *          that are under the fg of the optional mask.  If the mask
1441  *          is null, it finds the average of the pixels in pixs.
1442  *      (2) Likewise, use L_ROOT_MEAN_SQUARE to get the rms value of
1443  *          pixels in pixs, either masked or not; L_STANDARD_DEVIATION
1444  *          to get the standard deviation from the mean of the pixels;
1445  *          L_VARIANCE to get the average squared difference from the
1446  *          expected value.  The variance is the square of the stdev.
1447  *          For the standard deviation, we use
1448  *              sqrt(<(<x> - x)>^2) = sqrt(<x^2> - <x>^2)
1449  *      (3) Set the subsampling %factor > 1 to reduce the amount of
1450  *          computation.
1451  *      (4) Clipping of pixm (if it exists) to pixs is done in the inner loop.
1452  *      (5) Input x,y are ignored unless pixm exists.
1453  *      (6) A better name for this would be: pixGetPixelStatsGray()
1454  * </pre>
1455  */
1456 l_int32
pixGetAverageMasked(PIX * pixs,PIX * pixm,l_int32 x,l_int32 y,l_int32 factor,l_int32 type,l_float32 * pval)1457 pixGetAverageMasked(PIX        *pixs,
1458                     PIX        *pixm,
1459                     l_int32     x,
1460                     l_int32     y,
1461                     l_int32     factor,
1462                     l_int32     type,
1463                     l_float32  *pval)
1464 {
1465 l_int32    i, j, w, h, d, wm, hm, wplg, wplm, val, count;
1466 l_uint32  *datag, *datam, *lineg, *linem;
1467 l_float64  sumave, summs, ave, meansq, var;
1468 PIX       *pixg;
1469 
1470     PROCNAME("pixGetAverageMasked");
1471 
1472     if (!pval)
1473         return ERROR_INT("&val not defined", procName, 1);
1474     *pval = 0.0;
1475     if (!pixs)
1476         return ERROR_INT("pixs not defined", procName, 1);
1477     d = pixGetDepth(pixs);
1478     if (d != 8 && d != 16 && !pixGetColormap(pixs))
1479         return ERROR_INT("pixs not 8 or 16 bpp or colormapped", procName, 1);
1480     if (pixm && pixGetDepth(pixm) != 1)
1481         return ERROR_INT("pixm not 1 bpp", procName, 1);
1482     if (factor < 1)
1483         return ERROR_INT("sampling factor must be >= 1", procName, 1);
1484     if (type != L_MEAN_ABSVAL && type != L_ROOT_MEAN_SQUARE &&
1485         type != L_STANDARD_DEVIATION && type != L_VARIANCE)
1486         return ERROR_INT("invalid measure type", procName, 1);
1487 
1488     if (pixGetColormap(pixs))
1489         pixg = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
1490     else
1491         pixg = pixClone(pixs);
1492     pixGetDimensions(pixg, &w, &h, &d);
1493     datag = pixGetData(pixg);
1494     wplg = pixGetWpl(pixg);
1495 
1496     sumave = summs = 0.0;
1497     count = 0;
1498     if (!pixm) {
1499         for (i = 0; i < h; i += factor) {
1500             lineg = datag + i * wplg;
1501             for (j = 0; j < w; j += factor) {
1502                 if (d == 8)
1503                     val = GET_DATA_BYTE(lineg, j);
1504                 else  /* d == 16 */
1505                     val = GET_DATA_TWO_BYTES(lineg, j);
1506                 if (type != L_ROOT_MEAN_SQUARE)
1507                     sumave += val;
1508                 if (type != L_MEAN_ABSVAL)
1509                     summs += (l_float64)(val) * val;
1510                 count++;
1511             }
1512         }
1513     } else {
1514         pixGetDimensions(pixm, &wm, &hm, NULL);
1515         datam = pixGetData(pixm);
1516         wplm = pixGetWpl(pixm);
1517         for (i = 0; i < hm; i += factor) {
1518             if (y + i < 0 || y + i >= h) continue;
1519             lineg = datag + (y + i) * wplg;
1520             linem = datam + i * wplm;
1521             for (j = 0; j < wm; j += factor) {
1522                 if (x + j < 0 || x + j >= w) continue;
1523                 if (GET_DATA_BIT(linem, j)) {
1524                     if (d == 8)
1525                         val = GET_DATA_BYTE(lineg, x + j);
1526                     else  /* d == 16 */
1527                         val = GET_DATA_TWO_BYTES(lineg, x + j);
1528                     if (type != L_ROOT_MEAN_SQUARE)
1529                         sumave += val;
1530                     if (type != L_MEAN_ABSVAL)
1531                         summs += (l_float64)(val) * val;
1532                     count++;
1533                 }
1534             }
1535         }
1536     }
1537 
1538     pixDestroy(&pixg);
1539     if (count == 0)
1540         return ERROR_INT("no pixels sampled", procName, 1);
1541     ave = sumave / (l_float64)count;
1542     meansq = summs / (l_float64)count;
1543     var = meansq - ave * ave;
1544     if (type == L_MEAN_ABSVAL)
1545         *pval = (l_float32)ave;
1546     else if (type == L_ROOT_MEAN_SQUARE)
1547         *pval = (l_float32)sqrt(meansq);
1548     else if (type == L_STANDARD_DEVIATION)
1549         *pval = (l_float32)sqrt(var);
1550     else  /* type == L_VARIANCE */
1551         *pval = (l_float32)var;
1552 
1553     return 0;
1554 }
1555 
1556 
1557 /*!
1558  * \brief   pixGetAverageTiledRGB()
1559  *
1560  * \param[in]    pixs 32 bpp, or colormapped
1561  * \param[in]    sx, sy tile size; must be at least 2 x 2
1562  * \param[in]    type L_MEAN_ABSVAL, L_ROOT_MEAN_SQUARE, L_STANDARD_DEVIATION
1563  * \param[out]   ppixr [optional] tiled 'average' of red component
1564  * \param[out]   ppixg [optional] tiled 'average' of green component
1565  * \param[out]   ppixb [optional] tiled 'average' of blue component
1566  * \return  0 if OK, 1 on error
1567  *
1568  * <pre>
1569  * Notes:
1570  *      (1) For usage, see pixGetAverageTiled().
1571  *      (2) If there is a colormap, it is removed before the 8 bpp
1572  *          component images are extracted.
1573  * </pre>
1574  */
1575 l_int32
pixGetAverageTiledRGB(PIX * pixs,l_int32 sx,l_int32 sy,l_int32 type,PIX ** ppixr,PIX ** ppixg,PIX ** ppixb)1576 pixGetAverageTiledRGB(PIX     *pixs,
1577                       l_int32  sx,
1578                       l_int32  sy,
1579                       l_int32  type,
1580                       PIX    **ppixr,
1581                       PIX    **ppixg,
1582                       PIX    **ppixb)
1583 {
1584 PIX      *pixt;
1585 PIXCMAP  *cmap;
1586 
1587     PROCNAME("pixGetAverageTiledRGB");
1588 
1589     if (ppixr) *ppixr = NULL;
1590     if (ppixg) *ppixg = NULL;
1591     if (ppixb) *ppixb = NULL;
1592     if (!ppixr && !ppixg && !ppixb)
1593         return ERROR_INT("no data requested", procName, 1);
1594     if (!pixs)
1595         return ERROR_INT("pixs not defined", procName, 1);
1596     cmap = pixGetColormap(pixs);
1597     if (pixGetDepth(pixs) != 32 && !cmap)
1598         return ERROR_INT("pixs neither 32 bpp nor colormapped", procName, 1);
1599     if (sx < 2 || sy < 2)
1600         return ERROR_INT("sx and sy not both > 1", procName, 1);
1601     if (type != L_MEAN_ABSVAL && type != L_ROOT_MEAN_SQUARE &&
1602         type != L_STANDARD_DEVIATION)
1603         return ERROR_INT("invalid measure type", procName, 1);
1604 
1605     if (ppixr) {
1606         if (cmap)
1607             pixt = pixGetRGBComponentCmap(pixs, COLOR_RED);
1608         else
1609             pixt = pixGetRGBComponent(pixs, COLOR_RED);
1610         *ppixr = pixGetAverageTiled(pixt, sx, sy, type);
1611         pixDestroy(&pixt);
1612     }
1613     if (ppixg) {
1614         if (cmap)
1615             pixt = pixGetRGBComponentCmap(pixs, COLOR_GREEN);
1616         else
1617             pixt = pixGetRGBComponent(pixs, COLOR_GREEN);
1618         *ppixg = pixGetAverageTiled(pixt, sx, sy, type);
1619         pixDestroy(&pixt);
1620     }
1621     if (ppixb) {
1622         if (cmap)
1623             pixt = pixGetRGBComponentCmap(pixs, COLOR_BLUE);
1624         else
1625             pixt = pixGetRGBComponent(pixs, COLOR_BLUE);
1626         *ppixb = pixGetAverageTiled(pixt, sx, sy, type);
1627         pixDestroy(&pixt);
1628     }
1629 
1630     return 0;
1631 }
1632 
1633 
1634 /*!
1635  * \brief   pixGetAverageTiled()
1636  *
1637  * \param[in]    pixs 8 bpp, or colormapped
1638  * \param[in]    sx, sy tile size; must be at least 2 x 2
1639  * \param[in]    type L_MEAN_ABSVAL, L_ROOT_MEAN_SQUARE, L_STANDARD_DEVIATION
1640  * \return  pixd average values in each tile, or NULL on error
1641  *
1642  * <pre>
1643  * Notes:
1644  *      (1) Only computes for tiles that are entirely contained in pixs.
1645  *      (2) Use L_MEAN_ABSVAL to get the average abs value within the tile;
1646  *          L_ROOT_MEAN_SQUARE to get the rms value within each tile;
1647  *          L_STANDARD_DEVIATION to get the standard dev. from the average
1648  *          within each tile.
1649  *      (3) If colormapped, converts to 8 bpp gray.
1650  * </pre>
1651  */
1652 PIX *
pixGetAverageTiled(PIX * pixs,l_int32 sx,l_int32 sy,l_int32 type)1653 pixGetAverageTiled(PIX     *pixs,
1654                    l_int32  sx,
1655                    l_int32  sy,
1656                    l_int32  type)
1657 {
1658 l_int32    i, j, k, m, w, h, wd, hd, d, pos, wplt, wpld, valt;
1659 l_uint32  *datat, *datad, *linet, *lined, *startt;
1660 l_float64  sumave, summs, ave, meansq, normfact;
1661 PIX       *pixt, *pixd;
1662 
1663     PROCNAME("pixGetAverageTiled");
1664 
1665     if (!pixs)
1666         return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
1667     pixGetDimensions(pixs, &w, &h, &d);
1668     if (d != 8 && !pixGetColormap(pixs))
1669         return (PIX *)ERROR_PTR("pixs not 8 bpp or cmapped", procName, NULL);
1670     if (sx < 2 || sy < 2)
1671         return (PIX *)ERROR_PTR("sx and sy not both > 1", procName, NULL);
1672     wd = w / sx;
1673     hd = h / sy;
1674     if (wd < 1 || hd < 1)
1675         return (PIX *)ERROR_PTR("wd or hd == 0", procName, NULL);
1676     if (type != L_MEAN_ABSVAL && type != L_ROOT_MEAN_SQUARE &&
1677         type != L_STANDARD_DEVIATION)
1678         return (PIX *)ERROR_PTR("invalid measure type", procName, NULL);
1679 
1680     pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
1681     pixd = pixCreate(wd, hd, 8);
1682     datat = pixGetData(pixt);
1683     wplt = pixGetWpl(pixt);
1684     datad = pixGetData(pixd);
1685     wpld = pixGetWpl(pixd);
1686     normfact = 1. / (l_float64)(sx * sy);
1687     for (i = 0; i < hd; i++) {
1688         lined = datad + i * wpld;
1689         linet = datat + i * sy * wplt;
1690         for (j = 0; j < wd; j++) {
1691             if (type == L_MEAN_ABSVAL || type == L_STANDARD_DEVIATION) {
1692                 sumave = 0.0;
1693                 for (k = 0; k < sy; k++) {
1694                     startt = linet + k * wplt;
1695                     for (m = 0; m < sx; m++) {
1696                         pos = j * sx + m;
1697                         valt = GET_DATA_BYTE(startt, pos);
1698                         sumave += valt;
1699                     }
1700                 }
1701                 ave = normfact * sumave;
1702             }
1703             if (type == L_ROOT_MEAN_SQUARE || type == L_STANDARD_DEVIATION) {
1704                 summs = 0.0;
1705                 for (k = 0; k < sy; k++) {
1706                     startt = linet + k * wplt;
1707                     for (m = 0; m < sx; m++) {
1708                         pos = j * sx + m;
1709                         valt = GET_DATA_BYTE(startt, pos);
1710                         summs += (l_float64)(valt) * valt;
1711                     }
1712                 }
1713                 meansq = normfact * summs;
1714             }
1715             if (type == L_MEAN_ABSVAL)
1716                 valt = (l_int32)(ave + 0.5);
1717             else if (type == L_ROOT_MEAN_SQUARE)
1718                 valt = (l_int32)(sqrt(meansq) + 0.5);
1719             else  /* type == L_STANDARD_DEVIATION */
1720                 valt = (l_int32)(sqrt(meansq - ave * ave) + 0.5);
1721             SET_DATA_BYTE(lined, j, valt);
1722         }
1723     }
1724 
1725     pixDestroy(&pixt);
1726     return pixd;
1727 }
1728 
1729 
1730 /*!
1731  * \brief   pixRowStats()
1732  *
1733  * \param[in]    pixs 8 bpp; not cmapped
1734  * \param[in]    box [optional] clipping box; can be null
1735  * \param[out]   pnamean [optional] numa of mean values
1736  * \param[out]   pnamedian [optional] numa of median values
1737  * \param[out]   pnamode [optional] numa of mode intensity values
1738  * \param[out]   pnamodecount [optional] numa of mode counts
1739  * \param[out]   pnavar [optional] numa of variance
1740  * \param[out]   pnarootvar [optional] numa of square root of variance
1741  * \return  na numa of requested statistic for each row, or NULL on error
1742  *
1743  * <pre>
1744  * Notes:
1745  *      (1) This computes numas that represent column vectors of statistics,
1746  *          with each of its values derived from the corresponding row of a Pix.
1747  *      (2) Use NULL on input to prevent computation of any of the 5 numas.
1748  *      (3) Other functions that compute pixel row statistics are:
1749  *             pixCountPixelsByRow()
1750  *             pixAverageByRow()
1751  *             pixVarianceByRow()
1752  *             pixGetRowStats()
1753  * </pre>
1754  */
1755 l_int32
pixRowStats(PIX * pixs,BOX * box,NUMA ** pnamean,NUMA ** pnamedian,NUMA ** pnamode,NUMA ** pnamodecount,NUMA ** pnavar,NUMA ** pnarootvar)1756 pixRowStats(PIX    *pixs,
1757             BOX    *box,
1758             NUMA  **pnamean,
1759             NUMA  **pnamedian,
1760             NUMA  **pnamode,
1761             NUMA  **pnamodecount,
1762             NUMA  **pnavar,
1763             NUMA  **pnarootvar)
1764 {
1765 l_int32     i, j, k, w, h, val, wpls, sum, sumsq, target, max, modeval;
1766 l_int32     xstart, xend, ystart, yend, bw, bh;
1767 l_int32    *histo;
1768 l_uint32   *lines, *datas;
1769 l_float32   norm;
1770 l_float32  *famean, *fameansq, *favar, *farootvar;
1771 l_float32  *famedian, *famode, *famodecount;
1772 
1773     PROCNAME("pixRowStats");
1774 
1775     if (pnamean) *pnamean = NULL;
1776     if (pnamedian) *pnamedian = NULL;
1777     if (pnamode) *pnamode = NULL;
1778     if (pnamodecount) *pnamodecount = NULL;
1779     if (pnavar) *pnavar = NULL;
1780     if (pnarootvar) *pnarootvar = NULL;
1781     if (!pixs || pixGetDepth(pixs) != 8)
1782         return ERROR_INT("pixs undefined or not 8 bpp", procName, 1);
1783     famean = fameansq = favar = farootvar = NULL;
1784     famedian = famode = famodecount = NULL;
1785 
1786     pixGetDimensions(pixs, &w, &h, NULL);
1787     if (boxClipToRectangleParams(box, w, h, &xstart, &ystart, &xend, &yend,
1788                                  &bw, &bh) == 1)
1789         return ERROR_INT("invalid clipping box", procName, 1);
1790 
1791         /* We need the mean for variance and root variance */
1792     datas = pixGetData(pixs);
1793     wpls = pixGetWpl(pixs);
1794     if (pnamean || pnavar || pnarootvar) {
1795         norm = 1. / (l_float32)bw;
1796         famean = (l_float32 *)LEPT_CALLOC(bh, sizeof(l_float32));
1797         fameansq = (l_float32 *)LEPT_CALLOC(bh, sizeof(l_float32));
1798         if (pnavar || pnarootvar) {
1799             favar = (l_float32 *)LEPT_CALLOC(bh, sizeof(l_float32));
1800             if (pnarootvar)
1801                 farootvar = (l_float32 *)LEPT_CALLOC(bh, sizeof(l_float32));
1802         }
1803         for (i = ystart; i < yend; i++) {
1804             sum = sumsq = 0;
1805             lines = datas + i * wpls;
1806             for (j = xstart; j < xend; j++) {
1807                 val = GET_DATA_BYTE(lines, j);
1808                 sum += val;
1809                 sumsq += val * val;
1810             }
1811             famean[i] = norm * sum;
1812             fameansq[i] = norm * sumsq;
1813             if (pnavar || pnarootvar) {
1814                 favar[i] = fameansq[i] - famean[i] * famean[i];
1815                 if (pnarootvar)
1816                     farootvar[i] = sqrtf(favar[i]);
1817             }
1818         }
1819         LEPT_FREE(fameansq);
1820         if (pnamean)
1821             *pnamean = numaCreateFromFArray(famean, bh, L_INSERT);
1822         else
1823             LEPT_FREE(famean);
1824         if (pnavar)
1825             *pnavar = numaCreateFromFArray(favar, bh, L_INSERT);
1826         else
1827             LEPT_FREE(favar);
1828         if (pnarootvar)
1829             *pnarootvar = numaCreateFromFArray(farootvar, bh, L_INSERT);
1830     }
1831 
1832         /* We need a histogram to find the median and/or mode values */
1833     if (pnamedian || pnamode || pnamodecount) {
1834         histo = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
1835         if (pnamedian) {
1836             *pnamedian = numaMakeConstant(0, bh);
1837             famedian = numaGetFArray(*pnamedian, L_NOCOPY);
1838         }
1839         if (pnamode) {
1840             *pnamode = numaMakeConstant(0, bh);
1841             famode = numaGetFArray(*pnamode, L_NOCOPY);
1842         }
1843         if (pnamodecount) {
1844             *pnamodecount = numaMakeConstant(0, bh);
1845             famodecount = numaGetFArray(*pnamodecount, L_NOCOPY);
1846         }
1847         for (i = ystart; i < yend; i++) {
1848             lines = datas + i * wpls;
1849             memset(histo, 0, 1024);
1850             for (j = xstart; j < xend; j++) {
1851                 val = GET_DATA_BYTE(lines, j);
1852                 histo[val]++;
1853             }
1854 
1855             if (pnamedian) {
1856                 sum = 0;
1857                 target = (bw + 1) / 2;
1858                 for (k = 0; k < 256; k++) {
1859                     sum += histo[k];
1860                     if (sum >= target) {
1861                         famedian[i] = k;
1862                         break;
1863                     }
1864                 }
1865             }
1866 
1867             if (pnamode || pnamodecount) {
1868                 max = 0;
1869                 modeval = 0;
1870                 for (k = 0; k < 256; k++) {
1871                     if (histo[k] > max) {
1872                         max = histo[k];
1873                         modeval = k;
1874                     }
1875                 }
1876                 if (pnamode)
1877                     famode[i] = modeval;
1878                 if (pnamodecount)
1879                     famodecount[i] = max;
1880             }
1881         }
1882         LEPT_FREE(histo);
1883     }
1884 
1885     return 0;
1886 }
1887 
1888 
1889 /*!
1890  * \brief   pixColumnStats()
1891  *
1892  * \param[in]    pixs 8 bpp; not cmapped
1893  * \param[in]    box [optional] clipping box; can be null
1894  * \param[out]   pnamean [optional] numa of mean values
1895  * \param[out]   pnamedian [optional] numa of median values
1896  * \param[out]   pnamode [optional] numa of mode intensity values
1897  * \param[out]   pnamodecount [optional] numa of mode counts
1898  * \param[out]   pnavar [optional] numa of variance
1899  * \param[out]   pnarootvar [optional] numa of square root of variance
1900  * \return  na numa of requested statistic for each column,
1901  *                  or NULL on error
1902  *
1903  * <pre>
1904  * Notes:
1905  *      (1) This computes numas that represent row vectors of statistics,
1906  *          with each of its values derived from the corresponding col of a Pix.
1907  *      (2) Use NULL on input to prevent computation of any of the 5 numas.
1908  *      (3) Other functions that compute pixel column statistics are:
1909  *             pixCountPixelsByColumn()
1910  *             pixAverageByColumn()
1911  *             pixVarianceByColumn()
1912  *             pixGetColumnStats()
1913  * </pre>
1914  */
1915 l_int32
pixColumnStats(PIX * pixs,BOX * box,NUMA ** pnamean,NUMA ** pnamedian,NUMA ** pnamode,NUMA ** pnamodecount,NUMA ** pnavar,NUMA ** pnarootvar)1916 pixColumnStats(PIX    *pixs,
1917                BOX    *box,
1918                NUMA  **pnamean,
1919                NUMA  **pnamedian,
1920                NUMA  **pnamode,
1921                NUMA  **pnamodecount,
1922                NUMA  **pnavar,
1923                NUMA  **pnarootvar)
1924 {
1925 l_int32     i, j, k, w, h, val, wpls, sum, sumsq, target, max, modeval;
1926 l_int32     xstart, xend, ystart, yend, bw, bh;
1927 l_int32    *histo;
1928 l_uint32   *lines, *datas;
1929 l_float32   norm;
1930 l_float32  *famean, *fameansq, *favar, *farootvar;
1931 l_float32  *famedian, *famode, *famodecount;
1932 
1933     PROCNAME("pixColumnStats");
1934 
1935     if (pnamean) *pnamean = NULL;
1936     if (pnamedian) *pnamedian = NULL;
1937     if (pnamode) *pnamode = NULL;
1938     if (pnamodecount) *pnamodecount = NULL;
1939     if (pnavar) *pnavar = NULL;
1940     if (pnarootvar) *pnarootvar = NULL;
1941     if (!pixs || pixGetDepth(pixs) != 8)
1942         return ERROR_INT("pixs undefined or not 8 bpp", procName, 1);
1943     famean = fameansq = favar = farootvar = NULL;
1944     famedian = famode = famodecount = NULL;
1945 
1946     pixGetDimensions(pixs, &w, &h, NULL);
1947     if (boxClipToRectangleParams(box, w, h, &xstart, &ystart, &xend, &yend,
1948                                  &bw, &bh) == 1)
1949         return ERROR_INT("invalid clipping box", procName, 1);
1950 
1951         /* We need the mean for variance and root variance */
1952     datas = pixGetData(pixs);
1953     wpls = pixGetWpl(pixs);
1954     if (pnamean || pnavar || pnarootvar) {
1955         norm = 1. / (l_float32)bh;
1956         famean = (l_float32 *)LEPT_CALLOC(bw, sizeof(l_float32));
1957         fameansq = (l_float32 *)LEPT_CALLOC(bw, sizeof(l_float32));
1958         if (pnavar || pnarootvar) {
1959             favar = (l_float32 *)LEPT_CALLOC(bw, sizeof(l_float32));
1960             if (pnarootvar)
1961                 farootvar = (l_float32 *)LEPT_CALLOC(bw, sizeof(l_float32));
1962         }
1963         for (j = xstart; j < xend; j++) {
1964             sum = sumsq = 0;
1965             for (i = ystart, lines = datas; i < yend; lines += wpls, i++) {
1966                 val = GET_DATA_BYTE(lines, j);
1967                 sum += val;
1968                 sumsq += val * val;
1969             }
1970             famean[j] = norm * sum;
1971             fameansq[j] = norm * sumsq;
1972             if (pnavar || pnarootvar) {
1973                 favar[j] = fameansq[j] - famean[j] * famean[j];
1974                 if (pnarootvar)
1975                     farootvar[j] = sqrtf(favar[j]);
1976             }
1977         }
1978         LEPT_FREE(fameansq);
1979         if (pnamean)
1980             *pnamean = numaCreateFromFArray(famean, bw, L_INSERT);
1981         else
1982             LEPT_FREE(famean);
1983         if (pnavar)
1984             *pnavar = numaCreateFromFArray(favar, bw, L_INSERT);
1985         else
1986             LEPT_FREE(favar);
1987         if (pnarootvar)
1988             *pnarootvar = numaCreateFromFArray(farootvar, bw, L_INSERT);
1989     }
1990 
1991         /* We need a histogram to find the median and/or mode values */
1992     if (pnamedian || pnamode || pnamodecount) {
1993         histo = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
1994         if (pnamedian) {
1995             *pnamedian = numaMakeConstant(0, bw);
1996             famedian = numaGetFArray(*pnamedian, L_NOCOPY);
1997         }
1998         if (pnamode) {
1999             *pnamode = numaMakeConstant(0, bw);
2000             famode = numaGetFArray(*pnamode, L_NOCOPY);
2001         }
2002         if (pnamodecount) {
2003             *pnamodecount = numaMakeConstant(0, bw);
2004             famodecount = numaGetFArray(*pnamodecount, L_NOCOPY);
2005         }
2006         for (j = xstart; j < xend; j++) {
2007             memset(histo, 0, 1024);
2008             for (i = ystart, lines = datas; i < yend; lines += wpls, i++) {
2009                 val = GET_DATA_BYTE(lines, j);
2010                 histo[val]++;
2011             }
2012 
2013             if (pnamedian) {
2014                 sum = 0;
2015                 target = (bh + 1) / 2;
2016                 for (k = 0; k < 256; k++) {
2017                     sum += histo[k];
2018                     if (sum >= target) {
2019                         famedian[j] = k;
2020                         break;
2021                     }
2022                 }
2023             }
2024 
2025             if (pnamode || pnamodecount) {
2026                 max = 0;
2027                 modeval = 0;
2028                 for (k = 0; k < 256; k++) {
2029                     if (histo[k] > max) {
2030                         max = histo[k];
2031                         modeval = k;
2032                     }
2033                 }
2034                 if (pnamode)
2035                     famode[j] = modeval;
2036                 if (pnamodecount)
2037                     famodecount[j] = max;
2038             }
2039         }
2040         LEPT_FREE(histo);
2041     }
2042 
2043     return 0;
2044 }
2045 
2046 
2047 /*!
2048  * \brief   pixGetRangeValues()
2049  *
2050  * \param[in]    pixs 8 bpp grayscale, 32 bpp rgb, or colormapped
2051  * \param[in]    factor subsampling factor; >= 1; ignored if colormapped
2052  * \param[in]    color L_SELECT_RED, L_SELECT_GREEN or L_SELECT_BLUE
2053  * \param[out]   pminval [optional] minimum value of component
2054  * \param[out]   pmaxval [optional] maximum value of component
2055  * \return  0 if OK, 1 on error
2056  *
2057  * <pre>
2058  * Notes:
2059  *      (1) If pixs is 8 bpp grayscale, the color selection type is ignored.
2060  * </pre>
2061  */
2062 l_int32
pixGetRangeValues(PIX * pixs,l_int32 factor,l_int32 color,l_int32 * pminval,l_int32 * pmaxval)2063 pixGetRangeValues(PIX      *pixs,
2064                   l_int32   factor,
2065                   l_int32   color,
2066                   l_int32  *pminval,
2067                   l_int32  *pmaxval)
2068 {
2069 l_int32   d;
2070 PIXCMAP  *cmap;
2071 
2072     PROCNAME("pixGetRangeValues");
2073 
2074     if (pminval) *pminval = 0;
2075     if (pmaxval) *pmaxval = 0;
2076     if (!pminval && !pmaxval)
2077         return ERROR_INT("no result requested", procName, 1);
2078     if (!pixs)
2079         return ERROR_INT("pixs not defined", procName, 1);
2080 
2081     cmap = pixGetColormap(pixs);
2082     if (cmap)
2083         return pixcmapGetRangeValues(cmap, color, pminval, pmaxval,
2084                                      NULL, NULL);
2085 
2086     if (factor < 1)
2087         return ERROR_INT("sampling factor must be >= 1", procName, 1);
2088     d = pixGetDepth(pixs);
2089     if (d != 8 && d != 32)
2090         return ERROR_INT("pixs not 8 or 32 bpp", procName, 1);
2091 
2092     if (d == 8) {
2093         pixGetExtremeValue(pixs, factor, L_SELECT_MIN,
2094                            NULL, NULL, NULL, pminval);
2095         pixGetExtremeValue(pixs, factor, L_SELECT_MAX,
2096                            NULL, NULL, NULL, pmaxval);
2097     } else if (color == L_SELECT_RED) {
2098         pixGetExtremeValue(pixs, factor, L_SELECT_MIN,
2099                            pminval, NULL, NULL, NULL);
2100         pixGetExtremeValue(pixs, factor, L_SELECT_MAX,
2101                            pmaxval, NULL, NULL, NULL);
2102     } else if (color == L_SELECT_GREEN) {
2103         pixGetExtremeValue(pixs, factor, L_SELECT_MIN,
2104                            NULL, pminval, NULL, NULL);
2105         pixGetExtremeValue(pixs, factor, L_SELECT_MAX,
2106                            NULL, pmaxval, NULL, NULL);
2107     } else if (color == L_SELECT_BLUE) {
2108         pixGetExtremeValue(pixs, factor, L_SELECT_MIN,
2109                            NULL, NULL, pminval, NULL);
2110         pixGetExtremeValue(pixs, factor, L_SELECT_MAX,
2111                            NULL, NULL, pmaxval, NULL);
2112     } else {
2113         return ERROR_INT("invalid color", procName, 1);
2114     }
2115 
2116     return 0;
2117 }
2118 
2119 
2120 /*!
2121  * \brief   pixGetExtremeValue()
2122  *
2123  * \param[in]    pixs 8 bpp grayscale, 32 bpp rgb, or colormapped
2124  * \param[in]    factor subsampling factor; >= 1; ignored if colormapped
2125  * \param[in]    type L_SELECT_MIN or L_SELECT_MAX
2126  * \param[out]   prval [optional] red component
2127  * \param[out]   pgval [optional] green component
2128  * \param[out]   pbval [optional] blue component
2129  * \param[out]   pgrayval [optional] min or max gray value
2130  * \return  0 if OK, 1 on error
2131  *
2132  * <pre>
2133  * Notes:
2134  *      (1) If pixs is grayscale, the result is returned in &grayval.
2135  *          Otherwise, if there is a colormap or d == 32,
2136  *          each requested color component is returned.  At least
2137  *          one color component (address) must be input.
2138  * </pre>
2139  */
2140 l_int32
pixGetExtremeValue(PIX * pixs,l_int32 factor,l_int32 type,l_int32 * prval,l_int32 * pgval,l_int32 * pbval,l_int32 * pgrayval)2141 pixGetExtremeValue(PIX      *pixs,
2142                    l_int32   factor,
2143                    l_int32   type,
2144                    l_int32  *prval,
2145                    l_int32  *pgval,
2146                    l_int32  *pbval,
2147                    l_int32  *pgrayval)
2148 {
2149 l_int32    i, j, w, h, d, wpl;
2150 l_int32    val, extval, rval, gval, bval, extrval, extgval, extbval;
2151 l_uint32   pixel;
2152 l_uint32  *data, *line;
2153 PIXCMAP   *cmap;
2154 
2155     PROCNAME("pixGetExtremeValue");
2156 
2157     if (prval) *prval = -1;
2158     if (pgval) *pgval = -1;
2159     if (pbval) *pbval = -1;
2160     if (pgrayval) *pgrayval = -1;
2161     if (!pixs)
2162         return ERROR_INT("pixs not defined", procName, 1);
2163     if (type != L_SELECT_MIN && type != L_SELECT_MAX)
2164         return ERROR_INT("invalid type", procName, 1);
2165 
2166     cmap = pixGetColormap(pixs);
2167     if (cmap) {
2168         if (type == L_SELECT_MIN) {
2169             if (prval) pixcmapGetRangeValues(cmap, L_SELECT_RED, prval, NULL,
2170                                              NULL, NULL);
2171             if (pgval) pixcmapGetRangeValues(cmap, L_SELECT_GREEN, pgval, NULL,
2172                                              NULL, NULL);
2173             if (pbval) pixcmapGetRangeValues(cmap, L_SELECT_BLUE, pbval, NULL,
2174                                              NULL, NULL);
2175         } else {  /* type == L_SELECT_MAX */
2176             if (prval) pixcmapGetRangeValues(cmap, L_SELECT_RED, NULL, prval,
2177                                              NULL, NULL);
2178             if (pgval) pixcmapGetRangeValues(cmap, L_SELECT_GREEN, NULL, pgval,
2179                                              NULL, NULL);
2180             if (pbval) pixcmapGetRangeValues(cmap, L_SELECT_BLUE, NULL, pbval,
2181                                              NULL, NULL);
2182         }
2183         return 0;
2184     }
2185 
2186     pixGetDimensions(pixs, &w, &h, &d);
2187     if (factor < 1)
2188         return ERROR_INT("sampling factor must be >= 1", procName, 1);
2189     if (d != 8 && d != 32)
2190         return ERROR_INT("pixs not 8 or 32 bpp", procName, 1);
2191     if (d == 8 && !pgrayval)
2192         return ERROR_INT("can't return result in grayval", procName, 1);
2193     if (d == 32 && !prval && !pgval && !pbval)
2194         return ERROR_INT("can't return result in r/g/b-val", procName, 1);
2195 
2196     data = pixGetData(pixs);
2197     wpl = pixGetWpl(pixs);
2198     if (d == 8) {
2199         if (type == L_SELECT_MIN)
2200             extval = 100000;
2201         else  /* get max */
2202             extval = -1;
2203 
2204         for (i = 0; i < h; i += factor) {
2205             line = data + i * wpl;
2206             for (j = 0; j < w; j += factor) {
2207                 val = GET_DATA_BYTE(line, j);
2208                 if ((type == L_SELECT_MIN && val < extval) ||
2209                     (type == L_SELECT_MAX && val > extval))
2210                     extval = val;
2211             }
2212         }
2213         *pgrayval = extval;
2214         return 0;
2215     }
2216 
2217         /* 32 bpp rgb */
2218     if (type == L_SELECT_MIN) {
2219         extrval = 100000;
2220         extgval = 100000;
2221         extbval = 100000;
2222     } else {
2223         extrval = -1;
2224         extgval = -1;
2225         extbval = -1;
2226     }
2227     for (i = 0; i < h; i += factor) {
2228         line = data + i * wpl;
2229         for (j = 0; j < w; j += factor) {
2230             pixel = line[j];
2231             if (prval) {
2232                 rval = (pixel >> L_RED_SHIFT) & 0xff;
2233                 if ((type == L_SELECT_MIN && rval < extrval) ||
2234                     (type == L_SELECT_MAX && rval > extrval))
2235                     extrval = rval;
2236             }
2237             if (pgval) {
2238                 gval = (pixel >> L_GREEN_SHIFT) & 0xff;
2239                 if ((type == L_SELECT_MIN && gval < extgval) ||
2240                     (type == L_SELECT_MAX && gval > extgval))
2241                     extgval = gval;
2242             }
2243             if (pbval) {
2244                 bval = (pixel >> L_BLUE_SHIFT) & 0xff;
2245                 if ((type == L_SELECT_MIN && bval < extbval) ||
2246                     (type == L_SELECT_MAX && bval > extbval))
2247                     extbval = bval;
2248             }
2249         }
2250     }
2251     if (prval) *prval = extrval;
2252     if (pgval) *pgval = extgval;
2253     if (pbval) *pbval = extbval;
2254     return 0;
2255 }
2256 
2257 
2258 /*!
2259  * \brief   pixGetMaxValueInRect()
2260  *
2261  * \param[in]    pixs 8, 16 or 32 bpp grayscale; no color space components
2262  * \param[in]    box [optional] region; set box = NULL to use entire pixs
2263  * \param[out]   pmaxval [optional] max value in region
2264  * \param[out]   pxmax [optional] x location of max value
2265  * \param[out]   pymax [optional] y location of max value
2266  * \return  0 if OK, 1 on error
2267  *
2268  * <pre>
2269  * Notes:
2270  *      (1) This can be used to find the maximum and its location
2271  *          in a 2-dimensional histogram, where the x and y directions
2272  *          represent two color components (e.g., saturation and hue).
2273  *      (2) Note that here a 32 bpp pixs has pixel values that are simply
2274  *          numbers.  They are not 8 bpp components in a colorspace.
2275  * </pre>
2276  */
2277 l_int32
pixGetMaxValueInRect(PIX * pixs,BOX * box,l_uint32 * pmaxval,l_int32 * pxmax,l_int32 * pymax)2278 pixGetMaxValueInRect(PIX       *pixs,
2279                      BOX       *box,
2280                      l_uint32  *pmaxval,
2281                      l_int32   *pxmax,
2282                      l_int32   *pymax)
2283 {
2284 l_int32    i, j, w, h, d, wpl, bw, bh;
2285 l_int32    xstart, ystart, xend, yend, xmax, ymax;
2286 l_uint32   val, maxval;
2287 l_uint32  *data, *line;
2288 
2289     PROCNAME("pixGetMaxValueInRect");
2290 
2291     if (pmaxval) *pmaxval = 0;
2292     if (pxmax) *pxmax = 0;
2293     if (pymax) *pymax = 0;
2294     if (!pmaxval && !pxmax && !pymax)
2295         return ERROR_INT("no data requested", procName, 1);
2296     if (!pixs)
2297         return ERROR_INT("pixs not defined", procName, 1);
2298     if (pixGetColormap(pixs) != NULL)
2299         return ERROR_INT("pixs has colormap", procName, 1);
2300     pixGetDimensions(pixs, &w, &h, &d);
2301     if (d != 8 && d != 16 && d != 32)
2302         return ERROR_INT("pixs not 8, 16 or 32 bpp", procName, 1);
2303 
2304     xstart = ystart = 0;
2305     xend = w - 1;
2306     yend = h - 1;
2307     if (box) {
2308         boxGetGeometry(box, &xstart, &ystart, &bw, &bh);
2309         xend = xstart + bw - 1;
2310         yend = ystart + bh - 1;
2311     }
2312 
2313     data = pixGetData(pixs);
2314     wpl = pixGetWpl(pixs);
2315     maxval = 0;
2316     xmax = ymax = 0;
2317     for (i = ystart; i <= yend; i++) {
2318         line = data + i * wpl;
2319         for (j = xstart; j <= xend; j++) {
2320             if (d == 8)
2321                 val = GET_DATA_BYTE(line, j);
2322             else if (d == 16)
2323                 val = GET_DATA_TWO_BYTES(line, j);
2324             else  /* d == 32 */
2325                 val = line[j];
2326             if (val > maxval) {
2327                 maxval = val;
2328                 xmax = j;
2329                 ymax = i;
2330             }
2331         }
2332     }
2333     if (maxval == 0) {  /* no counts; pick the center of the rectangle */
2334         xmax = (xstart + xend) / 2;
2335         ymax = (ystart + yend) / 2;
2336     }
2337 
2338     if (pmaxval) *pmaxval = maxval;
2339     if (pxmax) *pxmax = xmax;
2340     if (pymax) *pymax = ymax;
2341     return 0;
2342 }
2343 
2344 
2345 /*!
2346  * \brief   pixGetBinnedComponentRange()
2347  *
2348  * \param[in]    pixs 32 bpp rgb
2349  * \param[in]    nbins number of equal population bins; must be > 1
2350  * \param[in]    factor subsampling factor; >= 1
2351  * \param[in]    color L_SELECT_RED, L_SELECT_GREEN or L_SELECT_BLUE
2352  * \param[out]   pminval [optional] minimum value of component
2353  * \param[out]   pmaxval [optional] maximum value of component
2354  * \param[out]   pcarray [optional] color array of bins
2355  * \param[in]    fontsize [optional] 0 for no debug; for debug, valid set
2356  *                        is {4,6,8,10,12,14,16,18,20}.
2357  * \return  0 if OK, 1 on error
2358  *
2359  * <pre>
2360  * Notes:
2361  *      (1) This returns the min and max average values of the
2362  *          selected color component in the set of rank bins,
2363  *          where the ranking is done using the specified component.
2364  * </pre>
2365  */
2366 l_int32
pixGetBinnedComponentRange(PIX * pixs,l_int32 nbins,l_int32 factor,l_int32 color,l_int32 * pminval,l_int32 * pmaxval,l_uint32 ** pcarray,l_int32 fontsize)2367 pixGetBinnedComponentRange(PIX        *pixs,
2368                            l_int32     nbins,
2369                            l_int32     factor,
2370                            l_int32     color,
2371                            l_int32    *pminval,
2372                            l_int32    *pmaxval,
2373                            l_uint32  **pcarray,
2374                            l_int32     fontsize)
2375 {
2376 l_int32    i, minval, maxval, rval, gval, bval;
2377 l_uint32  *carray;
2378 PIX       *pixt;
2379 
2380     PROCNAME("pixGetBinnedComponentRange");
2381 
2382     if (pminval) *pminval = 0;
2383     if (pmaxval) *pmaxval = 0;
2384     if (pcarray) *pcarray = NULL;
2385     if (!pminval && !pmaxval)
2386         return ERROR_INT("no result requested", procName, 1);
2387     if (!pixs || pixGetDepth(pixs) != 32)
2388         return ERROR_INT("pixs not defined or not 32 bpp", procName, 1);
2389     if (factor < 1)
2390         return ERROR_INT("sampling factor must be >= 1", procName, 1);
2391     if (color != L_SELECT_RED && color != L_SELECT_GREEN &&
2392         color != L_SELECT_BLUE)
2393         return ERROR_INT("invalid color", procName, 1);
2394     if (fontsize < 0 || fontsize > 20 || fontsize & 1 || fontsize == 2)
2395         return ERROR_INT("invalid fontsize", procName, 1);
2396 
2397     pixGetRankColorArray(pixs, nbins, color, factor, &carray, 0, 0);
2398     if (fontsize > 0) {
2399         for (i = 0; i < nbins; i++)
2400             L_INFO("c[%d] = %x\n", procName, i, carray[i]);
2401         pixt = pixDisplayColorArray(carray, nbins, 200, 5, fontsize);
2402         pixDisplay(pixt, 100, 100);
2403         pixDestroy(&pixt);
2404     }
2405 
2406     extractRGBValues(carray[0], &rval, &gval, &bval);
2407     minval = rval;
2408     if (color == L_SELECT_GREEN)
2409         minval = gval;
2410     else if (color == L_SELECT_BLUE)
2411         minval = bval;
2412     extractRGBValues(carray[nbins - 1], &rval, &gval, &bval);
2413     maxval = rval;
2414     if (color == L_SELECT_GREEN)
2415         maxval = gval;
2416     else if (color == L_SELECT_BLUE)
2417         maxval = bval;
2418 
2419     if (pminval) *pminval = minval;
2420     if (pmaxval) *pmaxval = maxval;
2421     if (pcarray)
2422         *pcarray = carray;
2423     else
2424         LEPT_FREE(carray);
2425     return 0;
2426 }
2427 
2428 
2429 /*!
2430  * \brief   pixGetRankColorArray()
2431  *
2432  * \param[in]    pixs      32 bpp or cmapped
2433  * \param[in]    nbins     number of equal population bins; must be > 1
2434  * \param[in]    type      color selection flag
2435  * \param[in]    factor    subsampling factor; integer >= 1
2436  * \param[out]   pcarray   array of colors, ranked by intensity
2437  * \param[in]    debugflag 1 to display color squares and plots of color
2438  *                         components; 2 to write them as png to file
2439  * \param[in]    fontsize  [optional] 0 for no debug; for debug, valid set
2440  *                         is {4,6,8,10,12,14,16,18,20}.  Ignored if
2441  *                         debugflag == 0.  fontsize == 6 is typical.
2442  * \return  0 if OK, 1 on error
2443  *
2444  * <pre>
2445  * Notes:
2446  *      (1) The color selection flag is one of: L_SELECT_RED, L_SELECT_GREEN,
2447  *          L_SELECT_BLUE, L_SELECT_MIN, L_SELECT_MAX, L_SELECT_AVERAGE,
2448  *          L_SELECT_HUE, L_SELECT_SATURATION.
2449  *      (2) Then it finds the histogram of the selected color type in each
2450  *          RGB pixel.  For each of the %nbins sets of pixels,
2451  *          ordered by this color type value, find the average RGB color,
2452  *          and return this as a "rank color" array.  The output array
2453  *          has %nbins colors.
2454  *      (3) Set the subsampling factor > 1 to reduce the amount of
2455  *          computation.  Typically you want at least 10,000 pixels
2456  *          for reasonable statistics.
2457  *      (4) The rank color as a function of rank can then be found from
2458  *             rankint = (l_int32)(rank * (nbins - 1) + 0.5);
2459  *             extractRGBValues(array[rankint], &rval, &gval, &bval);
2460  *          where the rank is in [0.0 ... 1.0].
2461  *          This function is meant to be simple and approximate.
2462  *      (5) Compare this with pixGetBinnedColor(), which generates equal
2463  *          width intensity bins and finds the average color in each bin.
2464  * </pre>
2465  */
2466 l_int32
pixGetRankColorArray(PIX * pixs,l_int32 nbins,l_int32 type,l_int32 factor,l_uint32 ** pcarray,l_int32 debugflag,l_int32 fontsize)2467 pixGetRankColorArray(PIX        *pixs,
2468                      l_int32     nbins,
2469                      l_int32     type,
2470                      l_int32     factor,
2471                      l_uint32  **pcarray,
2472                      l_int32     debugflag,
2473                      l_int32     fontsize)
2474 {
2475 l_int32    ret;
2476 l_uint32  *array;
2477 NUMA      *na, *nan, *narbin;
2478 PIX       *pixt, *pixc, *pixg, *pixd;
2479 PIXCMAP   *cmap;
2480 
2481     PROCNAME("pixGetRankColorArray");
2482 
2483     if (!pcarray)
2484         return ERROR_INT("&carray not defined", procName, 1);
2485     *pcarray = NULL;
2486     if (factor < 1)
2487         return ERROR_INT("sampling factor must be >= 1", procName, 1);
2488     if (nbins < 2)
2489         return ERROR_INT("nbins must be at least 2", procName, 1);
2490     if (!pixs)
2491         return ERROR_INT("pixs not defined", procName, 1);
2492     cmap = pixGetColormap(pixs);
2493     if (pixGetDepth(pixs) != 32 && !cmap)
2494         return ERROR_INT("pixs neither 32 bpp nor cmapped", procName, 1);
2495     if (type != L_SELECT_RED && type != L_SELECT_GREEN &&
2496         type != L_SELECT_BLUE && type != L_SELECT_MIN &&
2497         type != L_SELECT_MAX && type != L_SELECT_AVERAGE &&
2498         type != L_SELECT_HUE && type != L_SELECT_SATURATION)
2499         return ERROR_INT("invalid type", procName, 1);
2500     if (debugflag > 0) {
2501         if (fontsize < 0 || fontsize > 20 || fontsize & 1 || fontsize == 2)
2502             return ERROR_INT("invalid fontsize", procName, 1);
2503     }
2504 
2505         /* Downscale by factor and remove colormap if it exists */
2506     pixt = pixScaleByIntSampling(pixs, factor);
2507     if (cmap)
2508         pixc = pixRemoveColormap(pixt, REMOVE_CMAP_TO_FULL_COLOR);
2509     else
2510         pixc = pixClone(pixt);
2511     pixDestroy(&pixt);
2512 
2513         /* Get normalized histogram of the selected component */
2514     if (type == L_SELECT_RED)
2515         pixg = pixGetRGBComponent(pixc, COLOR_RED);
2516     else if (type == L_SELECT_GREEN)
2517         pixg = pixGetRGBComponent(pixc, COLOR_GREEN);
2518     else if (type == L_SELECT_BLUE)
2519         pixg = pixGetRGBComponent(pixc, COLOR_BLUE);
2520     else if (type == L_SELECT_MIN)
2521         pixg = pixConvertRGBToGrayMinMax(pixc, L_CHOOSE_MIN);
2522     else if (type == L_SELECT_MAX)
2523         pixg = pixConvertRGBToGrayMinMax(pixc, L_CHOOSE_MAX);
2524     else if (type == L_SELECT_AVERAGE)
2525         pixg = pixConvertRGBToGray(pixc, 0.34, 0.33, 0.33);
2526     else if (type == L_SELECT_HUE)
2527         pixg = pixConvertRGBToHue(pixc);
2528     else  /* L_SELECT_SATURATION */
2529         pixg = pixConvertRGBToSaturation(pixc);
2530     if ((na = pixGetGrayHistogram(pixg, 1)) == NULL) {
2531         pixDestroy(&pixc);
2532         pixDestroy(&pixg);
2533         return ERROR_INT("na not made", procName, 1);
2534     }
2535     nan = numaNormalizeHistogram(na, 1.0);
2536 
2537         /* Get the following arrays:
2538          * (1) nar: cumulative normalized histogram (rank vs intensity value).
2539          *     With 256 intensity values, we have 257 rank values.
2540          * (2) nai: "average" intensity as function of rank bin, for
2541          *     %nbins equally spaced in rank between 0.0 and 1.0.
2542          * (3) narbin: bin number of discretized rank as a function of
2543          *     intensity.  This is the 'inverse' of nai.
2544          * (4) nabb: intensity value of the right bin boundary, for each
2545          *     of the %nbins discretized rank bins. */
2546     if (!debugflag) {
2547         numaDiscretizeRankAndIntensity(nan, nbins, &narbin, NULL, NULL, NULL);
2548     } else {
2549         NUMA  *nai, *nar, *nabb;
2550         numaDiscretizeRankAndIntensity(nan, nbins, &narbin, &nai, &nar, &nabb);
2551         lept_mkdir("lept/regout");
2552         gplotSimple1(nan, GPLOT_PNG, "/tmp/lept/regout/rtnan",
2553                      "Normalized Histogram");
2554         gplotSimple1(nar, GPLOT_PNG, "/tmp/lept/regout/rtnar",
2555                      "Cumulative Histogram");
2556         gplotSimple1(nai, GPLOT_PNG, "/tmp/lept/regout/rtnai",
2557                      "Intensity vs. rank bin");
2558         gplotSimple1(narbin, GPLOT_PNG, "/tmp/lept/regout/rtnarbin",
2559                      "LUT: rank bin vs. Intensity");
2560         gplotSimple1(nabb, GPLOT_PNG, "/tmp/lept/regout/rtnabb",
2561                      "Intensity of right edge vs. rank bin");
2562         numaDestroy(&nai);
2563         numaDestroy(&nar);
2564         numaDestroy(&nabb);
2565     }
2566 
2567         /* Get the average color in each bin for pixels whose grayscale
2568          * values fall in the bin range.  %narbin is the LUT that
2569          * determines the bin number from the grayscale version of
2570          * the image.  Because this mapping may not be unique,
2571          * some bins may not be represented in the LUT. In use, to get fair
2572          * allocation into all the bins, bin population is monitored
2573          * as pixels are accumulated, and when bins fill up,
2574          * pixels are required to overflow into succeeding bins. */
2575     pixGetBinnedColor(pixc, pixg, 1, nbins, narbin, pcarray, debugflag);
2576     ret = 0;
2577     if ((array = *pcarray) == NULL) {
2578         L_ERROR("color array not returned\n", procName);
2579         ret = 1;
2580         debugflag = 0;  /* make sure to skip the following */
2581     }
2582     if (debugflag) {
2583         pixd = pixDisplayColorArray(array, nbins, 200, 5, fontsize);
2584         if (debugflag == 1)
2585             pixDisplayWithTitle(pixd, 0, 500, "binned colors", 1);
2586         else  /* debugflag == 2 */
2587             pixWriteDebug("/tmp/lept/regout/rankhisto.png", pixd, IFF_PNG);
2588         pixDestroy(&pixd);
2589     }
2590 
2591     pixDestroy(&pixc);
2592     pixDestroy(&pixg);
2593     numaDestroy(&na);
2594     numaDestroy(&nan);
2595     numaDestroy(&narbin);
2596     return ret;
2597 }
2598 
2599 
2600 /*!
2601  * \brief   pixGetBinnedColor()
2602  *
2603  * \param[in]    pixs 32 bpp
2604  * \param[in]    pixg 8 bpp grayscale version of pixs
2605  * \param[in]    factor sampling factor along pixel counting direction
2606  * \param[in]    nbins number of intensity bins
2607  * \param[in]    nalut LUT for mapping from intensity to bin number
2608  * \param[out]   pcarray array of average color values in each bin
2609  * \param[in]    debugflag 1 to display output debug plots of color
2610  *                         components; 2 to write them as png to file
2611  * \return  0 if OK; 1 on error
2612  *
2613  * <pre>
2614  * Notes:
2615  *      (1) This takes a color image, a grayscale (intensity) version,
2616  *          a LUT from intensity to bin number, and the number of bins.
2617  *          It computes the average color for pixels whose intensity
2618  *          is in each bin.  This is returned as an array of l_uint32
2619  *          colors in our standard RGBA ordering.
2620  *      (2) This function generates equal width intensity bins and
2621  *          finds the average color in each bin.  Compare this with
2622  *          pixGetRankColorArray(), which rank orders the pixels
2623  *          by the value of the selected component in each pixel,
2624  *          sets up bins with equal population (not intensity width!),
2625  *          and gets the average color in each bin.
2626  * </pre>
2627  */
2628 l_int32
pixGetBinnedColor(PIX * pixs,PIX * pixg,l_int32 factor,l_int32 nbins,NUMA * nalut,l_uint32 ** pcarray,l_int32 debugflag)2629 pixGetBinnedColor(PIX        *pixs,
2630                   PIX        *pixg,
2631                   l_int32     factor,
2632                   l_int32     nbins,
2633                   NUMA       *nalut,
2634                   l_uint32  **pcarray,
2635                   l_int32     debugflag)
2636 {
2637 l_int32     i, j, w, h, wpls, wplg, grayval, bin, rval, gval, bval, success;
2638 l_int32     npts, avepts, maxpts;
2639 l_uint32   *datas, *datag, *lines, *lineg, *carray;
2640 l_float64   norm;
2641 l_float64  *rarray, *garray, *barray, *narray;
2642 
2643     PROCNAME("pixGetBinnedColor");
2644 
2645     if (!pcarray)
2646         return ERROR_INT("&carray not defined", procName, 1);
2647     *pcarray = NULL;
2648     if (!pixs)
2649         return ERROR_INT("pixs not defined", procName, 1);
2650     if (!pixg)
2651         return ERROR_INT("pixg not defined", procName, 1);
2652     if (!nalut)
2653         return ERROR_INT("nalut not defined", procName, 1);
2654     if (factor < 1) {
2655         L_WARNING("sampling factor less than 1; setting to 1\n", procName);
2656         factor = 1;
2657     }
2658 
2659         /* Find the color for each rank bin.  Note that we can have
2660          * multiple bins filled with pixels having the same gray value.
2661          * Therefore, because in general the mapping from gray value
2662          * to bin number is not unique, if a bin fills up (actually,
2663          * we allow it to slightly overfill), we roll the excess
2664          * over to the next bin, etc. */
2665     pixGetDimensions(pixs, &w, &h, NULL);
2666     npts = (w + factor - 1) * (h + factor - 1) / (factor * factor);
2667     avepts = (npts + nbins - 1) / nbins;  /* average number of pts in a bin */
2668     maxpts = (l_int32)((1.0 + 0.5 / (l_float32)nbins) * avepts);
2669     datas = pixGetData(pixs);
2670     wpls = pixGetWpl(pixs);
2671     datag = pixGetData(pixg);
2672     wplg = pixGetWpl(pixg);
2673     rarray = (l_float64 *)LEPT_CALLOC(nbins, sizeof(l_float64));
2674     garray = (l_float64 *)LEPT_CALLOC(nbins, sizeof(l_float64));
2675     barray = (l_float64 *)LEPT_CALLOC(nbins, sizeof(l_float64));
2676     narray = (l_float64 *)LEPT_CALLOC(nbins, sizeof(l_float64));
2677     for (i = 0; i < h; i += factor) {
2678         lines = datas + i * wpls;
2679         lineg = datag + i * wplg;
2680         for (j = 0; j < w; j += factor) {
2681             grayval = GET_DATA_BYTE(lineg, j);
2682             numaGetIValue(nalut, grayval, &bin);
2683             extractRGBValues(lines[j], &rval, &gval, &bval);
2684             while (narray[bin] >= maxpts && bin < nbins - 1)
2685                 bin++;
2686             rarray[bin] += rval;
2687             garray[bin] += gval;
2688             barray[bin] += bval;
2689             narray[bin] += 1.0;  /* count samples in each bin */
2690         }
2691     }
2692 
2693     for (i = 0; i < nbins; i++) {
2694         norm = 1. / narray[i];
2695         rarray[i] *= norm;
2696         garray[i] *= norm;
2697         barray[i] *= norm;
2698 /*        fprintf(stderr, "narray[%d] = %f\n", i, narray[i]);  */
2699     }
2700 
2701     if (debugflag) {
2702         NUMA *nared, *nagreen, *nablue;
2703         nared = numaCreate(nbins);
2704         nagreen = numaCreate(nbins);
2705         nablue = numaCreate(nbins);
2706         for (i = 0; i < nbins; i++) {
2707             numaAddNumber(nared, rarray[i]);
2708             numaAddNumber(nagreen, garray[i]);
2709             numaAddNumber(nablue, barray[i]);
2710         }
2711         lept_mkdir("lept/regout");
2712         gplotSimple1(nared, GPLOT_PNG, "/tmp/lept/regout/rtnared",
2713                      "Average red val vs. rank bin");
2714         gplotSimple1(nagreen, GPLOT_PNG, "/tmp/lept/regout/rtnagreen",
2715                      "Average green val vs. rank bin");
2716         gplotSimple1(nablue, GPLOT_PNG, "/tmp/lept/regout/rtnablue",
2717                      "Average blue val vs. rank bin");
2718         numaDestroy(&nared);
2719         numaDestroy(&nagreen);
2720         numaDestroy(&nablue);
2721     }
2722 
2723         /* Save colors for all bins  in a single array */
2724     success = TRUE;
2725     if ((carray = (l_uint32 *)LEPT_CALLOC(nbins, sizeof(l_uint32))) == NULL) {
2726         success = FALSE;
2727         L_ERROR("carray not made\n", procName);
2728         goto cleanup_arrays;
2729     }
2730     *pcarray = carray;
2731     for (i = 0; i < nbins; i++) {
2732         rval = (l_int32)(rarray[i] + 0.5);
2733         gval = (l_int32)(garray[i] + 0.5);
2734         bval = (l_int32)(barray[i] + 0.5);
2735         composeRGBPixel(rval, gval, bval, carray + i);
2736     }
2737 
2738 cleanup_arrays:
2739     LEPT_FREE(rarray);
2740     LEPT_FREE(garray);
2741     LEPT_FREE(barray);
2742     LEPT_FREE(narray);
2743     return (success) ? 0 : 1;
2744 }
2745 
2746 
2747 /*!
2748  * \brief   pixDisplayColorArray()
2749  *
2750  * \param[in]    carray array of colors: 0xrrggbb00
2751  * \param[in]    ncolors size of array
2752  * \param[in]    side size of each color square; suggest 200
2753  * \param[in]    ncols number of columns in output color matrix
2754  * \param[in]    fontsize to label each square with text.  Valid set is
2755  *                        {4,6,8,10,12,14,16,18,20}.  Use 0 to disable.
2756  * \return  pixd color array, or NULL on error
2757  */
2758 PIX *
pixDisplayColorArray(l_uint32 * carray,l_int32 ncolors,l_int32 side,l_int32 ncols,l_int32 fontsize)2759 pixDisplayColorArray(l_uint32  *carray,
2760                      l_int32    ncolors,
2761                      l_int32    side,
2762                      l_int32    ncols,
2763                      l_int32    fontsize)
2764 {
2765 char     textstr[256];
2766 l_int32  i, rval, gval, bval;
2767 L_BMF   *bmf;
2768 PIX     *pixt, *pixd;
2769 PIXA    *pixa;
2770 
2771     PROCNAME("pixDisplayColorArray");
2772 
2773     if (!carray)
2774         return (PIX *)ERROR_PTR("carray not defined", procName, NULL);
2775     if (fontsize < 0 || fontsize > 20 || fontsize & 1 || fontsize == 2)
2776         return (PIX *)ERROR_PTR("invalid fontsize", procName, NULL);
2777 
2778     bmf = (fontsize == 0) ? NULL : bmfCreate(NULL, fontsize);
2779     pixa = pixaCreate(ncolors);
2780     for (i = 0; i < ncolors; i++) {
2781         pixt = pixCreate(side, side, 32);
2782         pixSetAllArbitrary(pixt, carray[i]);
2783         if (bmf) {
2784             extractRGBValues(carray[i], &rval, &gval, &bval);
2785             snprintf(textstr, sizeof(textstr),
2786                      "%d: (%d %d %d)", i, rval, gval, bval);
2787             pixSaveTiledWithText(pixt, pixa, side, (i % ncols == 0) ? 1 : 0,
2788                                  20, 2, bmf, textstr, 0xff000000, L_ADD_BELOW);
2789         } else {
2790             pixSaveTiled(pixt, pixa, 1.0, (i % ncols == 0) ? 1 : 0, 20, 32);
2791         }
2792         pixDestroy(&pixt);
2793     }
2794     pixd = pixaDisplay(pixa, 0, 0);
2795 
2796     pixaDestroy(&pixa);
2797     bmfDestroy(&bmf);
2798     return pixd;
2799 }
2800 
2801 
2802 /*!
2803  * \brief   pixRankBinByStrip()
2804  *
2805  * \param[in]    pixs 32 bpp or cmapped
2806  * \param[in]    direction L_SCAN_HORIZONTAL or L_SCAN_VERTICAL
2807  * \param[in]    size of strips in scan direction
2808  * \param[in]    nbins number of equal population bins; must be > 1
2809  * \param[in]    type color selection flag
2810  * \return  pixd result, or NULL on error
2811  *
2812  * <pre>
2813  * Notes:
2814  *      (1) This generates a pix where each column represents a strip of
2815  *          the input image.  If %direction == L_SCAN_HORIZONTAL, the
2816  *          input impage is tiled into vertical strips of width %size,
2817  *          where %size is a compromise between getting better spatial
2818  *          columnwise resolution (small %size) and getting better
2819  *          columnwise statistical information (larger %size).  Likewise
2820  *          with rows of the image if %direction == L_SCAN_VERTICAL.
2821  *      (2) For L_HORIZONTAL_SCAN, the output pix contains rank binned
2822  *          median colors in each column that correspond to a vertical
2823  *          strip of width %size in the input image.
2824  *      (3) The color selection flag is one of: L_SELECT_RED, L_SELECT_GREEN,
2825  *          L_SELECT_BLUE, L_SELECT_MIN, L_SELECT_MAX, L_SELECT_AVERAGE.
2826  *          It determines how the rank ordering is done.
2827  *      (4) Typical input values might be %size = 5, %nbins = 10.
2828  * </pre>
2829  */
2830 PIX *
pixRankBinByStrip(PIX * pixs,l_int32 direction,l_int32 size,l_int32 nbins,l_int32 type)2831 pixRankBinByStrip(PIX     *pixs,
2832                   l_int32  direction,
2833                   l_int32  size,
2834                   l_int32  nbins,
2835                   l_int32  type)
2836 {
2837 l_int32    i, j, w, h, nstrips;
2838 l_uint32  *array;
2839 BOXA      *boxa;
2840 PIX       *pix1, *pix2, *pixd;
2841 PIXA      *pixa;
2842 PIXCMAP   *cmap;
2843 
2844     PROCNAME("pixRankBinByStrip");
2845 
2846     if (!pixs)
2847         return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
2848     cmap = pixGetColormap(pixs);
2849     if (pixGetDepth(pixs) != 32 && !cmap)
2850         return (PIX *)ERROR_PTR("pixs neither 32 bpp nor cmapped",
2851                                 procName, NULL);
2852     if (direction != L_SCAN_HORIZONTAL && direction != L_SCAN_VERTICAL)
2853         return (PIX *)ERROR_PTR("invalid direction", procName, NULL);
2854     if (size < 1)
2855         return (PIX *)ERROR_PTR("size < 1", procName, NULL);
2856     if (nbins < 2)
2857         return (PIX *)ERROR_PTR("nbins must be at least 2", procName, NULL);
2858     if (type != L_SELECT_RED && type != L_SELECT_GREEN &&
2859         type != L_SELECT_BLUE && type != L_SELECT_MIN &&
2860         type != L_SELECT_MAX && type != L_SELECT_AVERAGE)
2861         return (PIX *)ERROR_PTR("invalid type", procName, NULL);
2862 
2863         /* Downscale by factor and remove colormap if it exists */
2864     if (cmap)
2865         pix1 = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);
2866     else
2867         pix1 = pixClone(pixs);
2868     pixGetDimensions(pixs, &w, &h, NULL);
2869 
2870     pixd = NULL;
2871     boxa = makeMosaicStrips(w, h, direction, size);
2872     pixa = pixClipRectangles(pix1, boxa);
2873     nstrips = pixaGetCount(pixa);
2874     if (direction == L_SCAN_HORIZONTAL) {
2875         pixd = pixCreate(nstrips, nbins, 32);
2876         for (i = 0; i < nstrips; i++) {
2877             pix2 = pixaGetPix(pixa, i, L_CLONE);
2878             pixGetRankColorArray(pix2, nbins, type, 1, &array, 0, 0);
2879             for (j = 0; j < nbins; j++)
2880                 pixSetPixel(pixd, i, j, array[j]);
2881             LEPT_FREE(array);
2882             pixDestroy(&pix2);
2883         }
2884     } else {  /* L_SCAN_VERTICAL */
2885         pixd = pixCreate(nbins, nstrips, 32);
2886         for (i = 0; i < nstrips; i++) {
2887             pix2 = pixaGetPix(pixa, i, L_CLONE);
2888             pixGetRankColorArray(pix2, nbins, type, 1, &array, 0, 0);
2889             for (j = 0; j < nbins; j++)
2890                 pixSetPixel(pixd, j, i, array[j]);
2891             LEPT_FREE(array);
2892             pixDestroy(&pix2);
2893         }
2894     }
2895     pixDestroy(&pix1);
2896     boxaDestroy(&boxa);
2897     pixaDestroy(&pixa);
2898     return pixd;
2899 }
2900 
2901 
2902 
2903 /*-------------------------------------------------------------*
2904  *                 Pixelwise aligned statistics                *
2905  *-------------------------------------------------------------*/
2906 /*!
2907  * \brief   pixaGetAlignedStats()
2908  *
2909  * \param[in]    pixa of identically sized, 8 bpp pix; not cmapped
2910  * \param[in]    type L_MEAN_ABSVAL, L_MEDIAN_VAL, L_MODE_VAL, L_MODE_COUNT
2911  * \param[in]    nbins of histogram for median and mode; ignored for mean
2912  * \param[in]    thresh on histogram for mode val; ignored for all other types
2913  * \return  pix with pixelwise aligned stats, or NULL on error.
2914  *
2915  * <pre>
2916  * Notes:
2917  *      (1) Each pixel in the returned pix represents an average
2918  *          (or median, or mode) over the corresponding pixels in each
2919  *          pix in the pixa.
2920  *      (2) The %thresh parameter works with L_MODE_VAL only, and
2921  *          sets a minimum occupancy of the mode bin.
2922  *          If the occupancy of the mode bin is less than %thresh, the
2923  *          mode value is returned as 0.  To always return the actual
2924  *          mode value, set %thresh = 0.  See pixGetRowStats().
2925  * </pre>
2926  */
2927 PIX *
pixaGetAlignedStats(PIXA * pixa,l_int32 type,l_int32 nbins,l_int32 thresh)2928 pixaGetAlignedStats(PIXA     *pixa,
2929                     l_int32   type,
2930                     l_int32   nbins,
2931                     l_int32   thresh)
2932 {
2933 l_int32     j, n, w, h, d;
2934 l_float32  *colvect;
2935 PIX        *pixt, *pixd;
2936 
2937     PROCNAME("pixaGetAlignedStats");
2938 
2939     if (!pixa)
2940         return (PIX *)ERROR_PTR("pixa not defined", procName, NULL);
2941     if (type != L_MEAN_ABSVAL && type != L_MEDIAN_VAL &&
2942         type != L_MODE_VAL && type != L_MODE_COUNT)
2943         return (PIX *)ERROR_PTR("invalid type", procName, NULL);
2944     n = pixaGetCount(pixa);
2945     if (n == 0)
2946         return (PIX *)ERROR_PTR("no pix in pixa", procName, NULL);
2947     pixaGetPixDimensions(pixa, 0, &w, &h, &d);
2948     if (d != 8)
2949         return (PIX *)ERROR_PTR("pix not 8 bpp", procName, NULL);
2950 
2951     pixd = pixCreate(w, h, 8);
2952     pixt = pixCreate(n, h, 8);
2953     colvect = (l_float32 *)LEPT_CALLOC(h, sizeof(l_float32));
2954     for (j = 0; j < w; j++) {
2955         pixaExtractColumnFromEachPix(pixa, j, pixt);
2956         pixGetRowStats(pixt, type, nbins, thresh, colvect);
2957         pixSetPixelColumn(pixd, j, colvect);
2958     }
2959 
2960     LEPT_FREE(colvect);
2961     pixDestroy(&pixt);
2962     return pixd;
2963 }
2964 
2965 
2966 /*!
2967  * \brief   pixaExtractColumnFromEachPix()
2968  *
2969  * \param[in]    pixa of identically sized, 8 bpp; not cmapped
2970  * \param[in]    col column index
2971  * \param[in]    pixd pix into which each column is inserted
2972  * \return  0 if OK, 1 on error
2973  */
2974 l_int32
pixaExtractColumnFromEachPix(PIXA * pixa,l_int32 col,PIX * pixd)2975 pixaExtractColumnFromEachPix(PIXA    *pixa,
2976                              l_int32  col,
2977                              PIX     *pixd)
2978 {
2979 l_int32    i, k, n, w, h, ht, val, wplt, wpld;
2980 l_uint32  *datad, *datat;
2981 PIX       *pixt;
2982 
2983     PROCNAME("pixaExtractColumnFromEachPix");
2984 
2985     if (!pixa)
2986         return ERROR_INT("pixa not defined", procName, 1);
2987     if (!pixd || pixGetDepth(pixd) != 8)
2988         return ERROR_INT("pixd not defined or not 8 bpp", procName, 1);
2989     n = pixaGetCount(pixa);
2990     pixGetDimensions(pixd, &w, &h, NULL);
2991     if (n != w)
2992         return ERROR_INT("pix width != n", procName, 1);
2993     pixt = pixaGetPix(pixa, 0, L_CLONE);
2994     wplt = pixGetWpl(pixt);
2995     pixGetDimensions(pixt, NULL, &ht, NULL);
2996     pixDestroy(&pixt);
2997     if (h != ht)
2998         return ERROR_INT("pixd height != column height", procName, 1);
2999 
3000     datad = pixGetData(pixd);
3001     wpld = pixGetWpl(pixd);
3002     for (k = 0; k < n; k++) {
3003         pixt = pixaGetPix(pixa, k, L_CLONE);
3004         datat = pixGetData(pixt);
3005         for (i = 0; i < h; i++) {
3006             val = GET_DATA_BYTE(datat, col);
3007             SET_DATA_BYTE(datad + i * wpld, k, val);
3008             datat += wplt;
3009         }
3010         pixDestroy(&pixt);
3011     }
3012 
3013     return 0;
3014 }
3015 
3016 
3017 /*!
3018  * \brief   pixGetRowStats()
3019  *
3020  * \param[in]    pixs 8 bpp; not cmapped
3021  * \param[in]    type L_MEAN_ABSVAL, L_MEDIAN_VAL, L_MODE_VAL, L_MODE_COUNT
3022  * \param[in]    nbins of histogram for median and mode; ignored for mean
3023  * \param[in]    thresh on histogram for mode; ignored for mean and median
3024  * \param[in]    colvect vector of results gathered across the rows of pixs
3025  * \return  0 if OK, 1 on error
3026  *
3027  * <pre>
3028  * Notes:
3029  *      (1) This computes a column vector of statistics using each
3030  *          row of a Pix.  The result is put in %colvect.
3031  *      (2) The %thresh parameter works with L_MODE_VAL only, and
3032  *          sets a minimum occupancy of the mode bin.
3033  *          If the occupancy of the mode bin is less than %thresh, the
3034  *          mode value is returned as 0.  To always return the actual
3035  *          mode value, set %thresh = 0.
3036  *      (3) What is the meaning of this %thresh parameter?
3037  *          For each row, the total count in the histogram is w, the
3038  *          image width.  So %thresh, relative to w, gives a measure
3039  *          of the ratio of the bin width to the width of the distribution.
3040  *          The larger %thresh, the narrower the distribution must be
3041  *          for the mode value to be returned (instead of returning 0).
3042  *      (4) If the Pix consists of a set of corresponding columns,
3043  *          one for each Pix in a Pixa, the width of the Pix is the
3044  *          number of Pix in the Pixa and the column vector can
3045  *          be stored as a column in a Pix of the same size as
3046  *          each Pix in the Pixa.
3047  * </pre>
3048  */
3049 l_int32
pixGetRowStats(PIX * pixs,l_int32 type,l_int32 nbins,l_int32 thresh,l_float32 * colvect)3050 pixGetRowStats(PIX        *pixs,
3051                l_int32     type,
3052                l_int32     nbins,
3053                l_int32     thresh,
3054                l_float32  *colvect)
3055 {
3056 l_int32    i, j, k, w, h, val, wpls, sum, target, max, modeval;
3057 l_int32   *histo, *gray2bin, *bin2gray;
3058 l_uint32  *lines, *datas;
3059 
3060     PROCNAME("pixGetRowStats");
3061 
3062     if (!pixs || pixGetDepth(pixs) != 8)
3063         return ERROR_INT("pixs not defined or not 8 bpp", procName, 1);
3064     if (!colvect)
3065         return ERROR_INT("colvect not defined", procName, 1);
3066     if (type != L_MEAN_ABSVAL && type != L_MEDIAN_VAL &&
3067         type != L_MODE_VAL && type != L_MODE_COUNT)
3068         return ERROR_INT("invalid type", procName, 1);
3069     if (type != L_MEAN_ABSVAL && (nbins < 1 || nbins > 256))
3070         return ERROR_INT("invalid nbins", procName, 1);
3071     pixGetDimensions(pixs, &w, &h, NULL);
3072 
3073     datas = pixGetData(pixs);
3074     wpls = pixGetWpl(pixs);
3075     if (type == L_MEAN_ABSVAL) {
3076         for (i = 0; i < h; i++) {
3077             sum = 0;
3078             lines = datas + i * wpls;
3079             for (j = 0; j < w; j++)
3080                 sum += GET_DATA_BYTE(lines, j);
3081             colvect[i] = (l_float32)sum / (l_float32)w;
3082         }
3083         return 0;
3084     }
3085 
3086         /* We need a histogram; binwidth ~ 256 / nbins */
3087     histo = (l_int32 *)LEPT_CALLOC(nbins, sizeof(l_int32));
3088     gray2bin = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
3089     bin2gray = (l_int32 *)LEPT_CALLOC(nbins, sizeof(l_int32));
3090     for (i = 0; i < 256; i++)  /* gray value --> histo bin */
3091         gray2bin[i] = (i * nbins) / 256;
3092     for (i = 0; i < nbins; i++)  /* histo bin --> gray value */
3093         bin2gray[i] = (i * 256 + 128) / nbins;
3094 
3095     for (i = 0; i < h; i++) {
3096         lines = datas + i * wpls;
3097         for (k = 0; k < nbins; k++)
3098             histo[k] = 0;
3099         for (j = 0; j < w; j++) {
3100             val = GET_DATA_BYTE(lines, j);
3101             histo[gray2bin[val]]++;
3102         }
3103 
3104         if (type == L_MEDIAN_VAL) {
3105             sum = 0;
3106             target = (w + 1) / 2;
3107             for (k = 0; k < nbins; k++) {
3108                 sum += histo[k];
3109                 if (sum >= target) {
3110                     colvect[i] = bin2gray[k];
3111                     break;
3112                 }
3113             }
3114         } else if (type == L_MODE_VAL) {
3115             max = 0;
3116             modeval = 0;
3117             for (k = 0; k < nbins; k++) {
3118                 if (histo[k] > max) {
3119                     max = histo[k];
3120                     modeval = k;
3121                 }
3122             }
3123             if (max < thresh)
3124                 colvect[i] = 0;
3125             else
3126                 colvect[i] = bin2gray[modeval];
3127         } else {  /* type == L_MODE_COUNT */
3128             max = 0;
3129             for (k = 0; k < nbins; k++) {
3130                 if (histo[k] > max)
3131                     max = histo[k];
3132             }
3133             colvect[i] = max;
3134         }
3135     }
3136 
3137     LEPT_FREE(histo);
3138     LEPT_FREE(gray2bin);
3139     LEPT_FREE(bin2gray);
3140     return 0;
3141 }
3142 
3143 
3144 /*!
3145  * \brief   pixGetColumnStats()
3146  *
3147  * \param[in]    pixs 8 bpp; not cmapped
3148  * \param[in]    type L_MEAN_ABSVAL, L_MEDIAN_VAL, L_MODE_VAL, L_MODE_COUNT
3149  * \param[in]    nbins of histogram for median and mode; ignored for mean
3150  * \param[in]    thresh on histogram for mode val; ignored for all other types
3151  * \param[in]    rowvect vector of results gathered down the columns of pixs
3152  * \return  0 if OK, 1 on error
3153  *
3154  * <pre>
3155  * Notes:
3156  *      (1) This computes a row vector of statistics using each
3157  *          column of a Pix.  The result is put in %rowvect.
3158  *      (2) The %thresh parameter works with L_MODE_VAL only, and
3159  *          sets a minimum occupancy of the mode bin.
3160  *          If the occupancy of the mode bin is less than %thresh, the
3161  *          mode value is returned as 0.  To always return the actual
3162  *          mode value, set %thresh = 0.
3163  *      (3) What is the meaning of this %thresh parameter?
3164  *          For each column, the total count in the histogram is h, the
3165  *          image height.  So %thresh, relative to h, gives a measure
3166  *          of the ratio of the bin width to the width of the distribution.
3167  *          The larger %thresh, the narrower the distribution must be
3168  *          for the mode value to be returned (instead of returning 0).
3169  * </pre>
3170  */
3171 l_int32
pixGetColumnStats(PIX * pixs,l_int32 type,l_int32 nbins,l_int32 thresh,l_float32 * rowvect)3172 pixGetColumnStats(PIX        *pixs,
3173                   l_int32     type,
3174                   l_int32     nbins,
3175                   l_int32     thresh,
3176                   l_float32  *rowvect)
3177 {
3178 l_int32    i, j, k, w, h, val, wpls, sum, target, max, modeval;
3179 l_int32   *histo, *gray2bin, *bin2gray;
3180 l_uint32  *datas;
3181 
3182     PROCNAME("pixGetColumnStats");
3183 
3184     if (!pixs || pixGetDepth(pixs) != 8)
3185         return ERROR_INT("pixs not defined or not 8 bpp", procName, 1);
3186     if (!rowvect)
3187         return ERROR_INT("rowvect not defined", procName, 1);
3188     if (type != L_MEAN_ABSVAL && type != L_MEDIAN_VAL &&
3189         type != L_MODE_VAL && type != L_MODE_COUNT)
3190         return ERROR_INT("invalid type", procName, 1);
3191     if (type != L_MEAN_ABSVAL && (nbins < 1 || nbins > 256))
3192         return ERROR_INT("invalid nbins", procName, 1);
3193     pixGetDimensions(pixs, &w, &h, NULL);
3194 
3195     datas = pixGetData(pixs);
3196     wpls = pixGetWpl(pixs);
3197     if (type == L_MEAN_ABSVAL) {
3198         for (j = 0; j < w; j++) {
3199             sum = 0;
3200             for (i = 0; i < h; i++)
3201                 sum += GET_DATA_BYTE(datas + i * wpls, j);
3202             rowvect[j] = (l_float32)sum / (l_float32)h;
3203         }
3204         return 0;
3205     }
3206 
3207         /* We need a histogram; binwidth ~ 256 / nbins */
3208     histo = (l_int32 *)LEPT_CALLOC(nbins, sizeof(l_int32));
3209     gray2bin = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
3210     bin2gray = (l_int32 *)LEPT_CALLOC(nbins, sizeof(l_int32));
3211     for (i = 0; i < 256; i++)  /* gray value --> histo bin */
3212         gray2bin[i] = (i * nbins) / 256;
3213     for (i = 0; i < nbins; i++)  /* histo bin --> gray value */
3214         bin2gray[i] = (i * 256 + 128) / nbins;
3215 
3216     for (j = 0; j < w; j++) {
3217         for (i = 0; i < h; i++) {
3218             val = GET_DATA_BYTE(datas + i * wpls, j);
3219             histo[gray2bin[val]]++;
3220         }
3221 
3222         if (type == L_MEDIAN_VAL) {
3223             sum = 0;
3224             target = (h + 1) / 2;
3225             for (k = 0; k < nbins; k++) {
3226                 sum += histo[k];
3227                 if (sum >= target) {
3228                     rowvect[j] = bin2gray[k];
3229                     break;
3230                 }
3231             }
3232         } else if (type == L_MODE_VAL) {
3233             max = 0;
3234             modeval = 0;
3235             for (k = 0; k < nbins; k++) {
3236                 if (histo[k] > max) {
3237                     max = histo[k];
3238                     modeval = k;
3239                 }
3240             }
3241             if (max < thresh)
3242                 rowvect[j] = 0;
3243             else
3244                 rowvect[j] = bin2gray[modeval];
3245         } else {  /* type == L_MODE_COUNT */
3246             max = 0;
3247             for (k = 0; k < nbins; k++) {
3248                 if (histo[k] > max)
3249                     max = histo[k];
3250             }
3251             rowvect[j] = max;
3252         }
3253         for (k = 0; k < nbins; k++)
3254             histo[k] = 0;
3255     }
3256 
3257     LEPT_FREE(histo);
3258     LEPT_FREE(gray2bin);
3259     LEPT_FREE(bin2gray);
3260     return 0;
3261 }
3262 
3263 
3264 /*!
3265  * \brief   pixSetPixelColumn()
3266  *
3267  * \param[in]    pix 8 bpp; not cmapped
3268  * \param[in]    col column index
3269  * \param[in]    colvect vector of floats
3270  * \return  0 if OK, 1 on error
3271  */
3272 l_int32
pixSetPixelColumn(PIX * pix,l_int32 col,l_float32 * colvect)3273 pixSetPixelColumn(PIX        *pix,
3274                   l_int32     col,
3275                   l_float32  *colvect)
3276 {
3277 l_int32    i, w, h, wpl;
3278 l_uint32  *data;
3279 
3280     PROCNAME("pixSetCPixelColumn");
3281 
3282     if (!pix || pixGetDepth(pix) != 8)
3283         return ERROR_INT("pix not defined or not 8 bpp", procName, 1);
3284     if (!colvect)
3285         return ERROR_INT("colvect not defined", procName, 1);
3286     pixGetDimensions(pix, &w, &h, NULL);
3287     if (col < 0 || col > w)
3288         return ERROR_INT("invalid col", procName, 1);
3289 
3290     data = pixGetData(pix);
3291     wpl = pixGetWpl(pix);
3292     for (i = 0; i < h; i++)
3293         SET_DATA_BYTE(data + i * wpl, col, (l_int32)colvect[i]);
3294 
3295     return 0;
3296 }
3297 
3298 
3299 /*-------------------------------------------------------------*
3300  *              Foreground/background estimation               *
3301  *-------------------------------------------------------------*/
3302 /*!
3303  * \brief   pixThresholdForFgBg()
3304  *
3305  * \param[in]    pixs any depth; cmapped ok
3306  * \param[in]    factor subsampling factor; integer >= 1
3307  * \param[in]    thresh threshold for generating foreground mask
3308  * \param[out]   pfgval [optional] average foreground value
3309  * \param[out]   pbgval [optional] average background value
3310  * \return  0 if OK, 1 on error
3311  */
3312 l_int32
pixThresholdForFgBg(PIX * pixs,l_int32 factor,l_int32 thresh,l_int32 * pfgval,l_int32 * pbgval)3313 pixThresholdForFgBg(PIX      *pixs,
3314                     l_int32   factor,
3315                     l_int32   thresh,
3316                     l_int32  *pfgval,
3317                     l_int32  *pbgval)
3318 {
3319 l_float32  fval;
3320 PIX       *pixg, *pixm;
3321 
3322     PROCNAME("pixThresholdForFgBg");
3323 
3324     if (pfgval) *pfgval = 0;
3325     if (pbgval) *pbgval = 0;
3326     if (!pfgval && !pbgval)
3327         return ERROR_INT("no data requested", procName, 1);
3328     if (!pixs)
3329         return ERROR_INT("pixs not defined", procName, 1);
3330 
3331         /* Generate a subsampled 8 bpp version and a mask over the fg */
3332     pixg = pixConvertTo8BySampling(pixs, factor, 0);
3333     pixm = pixThresholdToBinary(pixg, thresh);
3334 
3335     if (pfgval) {
3336         pixGetAverageMasked(pixg, pixm, 0, 0, 1, L_MEAN_ABSVAL, &fval);
3337         *pfgval = (l_int32)(fval + 0.5);
3338     }
3339 
3340     if (pbgval) {
3341         pixInvert(pixm, pixm);
3342         pixGetAverageMasked(pixg, pixm, 0, 0, 1, L_MEAN_ABSVAL, &fval);
3343         *pbgval = (l_int32)(fval + 0.5);
3344     }
3345 
3346     pixDestroy(&pixg);
3347     pixDestroy(&pixm);
3348     return 0;
3349 }
3350 
3351 
3352 /*!
3353  * \brief   pixSplitDistributionFgBg()
3354  *
3355  * \param[in]    pixs any depth; cmapped ok
3356  * \param[in]    scorefract fraction of the max score, used to determine
3357  *                          the range over which the histogram min is searched
3358  * \param[in]    factor subsampling factor; integer >= 1
3359  * \param[out]   pthresh [optional] best threshold for separating
3360  * \param[out]   pfgval [optional] average foreground value
3361  * \param[out]   pbgval [optional] average background value
3362  * \param[out]   ppixdb [optional] plot of distribution and split point
3363  * \return  0 if OK, 1 on error
3364  *
3365  * <pre>
3366  * Notes:
3367  *      (1) See numaSplitDistribution() for details on the underlying
3368  *          method of choosing a threshold.
3369  * </pre>
3370  */
3371 l_int32
pixSplitDistributionFgBg(PIX * pixs,l_float32 scorefract,l_int32 factor,l_int32 * pthresh,l_int32 * pfgval,l_int32 * pbgval,PIX ** ppixdb)3372 pixSplitDistributionFgBg(PIX       *pixs,
3373                          l_float32  scorefract,
3374                          l_int32    factor,
3375                          l_int32   *pthresh,
3376                          l_int32   *pfgval,
3377                          l_int32   *pbgval,
3378                          PIX      **ppixdb)
3379 {
3380 char       buf[256];
3381 l_int32    thresh;
3382 l_float32  avefg, avebg, maxnum;
3383 GPLOT     *gplot;
3384 NUMA      *na, *nascore, *nax, *nay;
3385 PIX       *pixg;
3386 
3387     PROCNAME("pixSplitDistributionFgBg");
3388 
3389     if (pthresh) *pthresh = 0;
3390     if (pfgval) *pfgval = 0;
3391     if (pbgval) *pbgval = 0;
3392     if (ppixdb) *ppixdb = NULL;
3393     if (!pthresh && !pfgval && !pbgval)
3394         return ERROR_INT("no data requested", procName, 1);
3395     if (!pixs)
3396         return ERROR_INT("pixs not defined", procName, 1);
3397 
3398         /* Generate a subsampled 8 bpp version */
3399     pixg = pixConvertTo8BySampling(pixs, factor, 0);
3400 
3401         /* Make the fg/bg estimates */
3402     na = pixGetGrayHistogram(pixg, 1);
3403     if (ppixdb) {
3404         numaSplitDistribution(na, scorefract, &thresh, &avefg, &avebg,
3405                               NULL, NULL, &nascore);
3406         numaDestroy(&nascore);
3407     } else {
3408         numaSplitDistribution(na, scorefract, &thresh, &avefg, &avebg,
3409                               NULL, NULL, NULL);
3410     }
3411 
3412     if (pthresh) *pthresh = thresh;
3413     if (pfgval) *pfgval = (l_int32)(avefg + 0.5);
3414     if (pbgval) *pbgval = (l_int32)(avebg + 0.5);
3415 
3416     if (ppixdb) {
3417         lept_mkdir("lept/redout");
3418         gplot = gplotCreate("/tmp/lept/redout/histplot", GPLOT_PNG, "Histogram",
3419                             "Grayscale value", "Number of pixels");
3420         gplotAddPlot(gplot, NULL, na, GPLOT_LINES, NULL);
3421         nax = numaMakeConstant(thresh, 2);
3422         numaGetMax(na, &maxnum, NULL);
3423         nay = numaMakeConstant(0, 2);
3424         numaReplaceNumber(nay, 1, (l_int32)(0.5 * maxnum));
3425         snprintf(buf, sizeof(buf), "score fract = %3.1f", scorefract);
3426         gplotAddPlot(gplot, nax, nay, GPLOT_LINES, buf);
3427         gplotMakeOutput(gplot);
3428         gplotDestroy(&gplot);
3429         numaDestroy(&nax);
3430         numaDestroy(&nay);
3431         *ppixdb = pixRead("/tmp/lept/redout/histplot.png");
3432     }
3433 
3434     pixDestroy(&pixg);
3435     numaDestroy(&na);
3436     return 0;
3437 }
3438