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 /*!
29 * \file morphapp.c
30 * <pre>
31 *
32 * These are some useful and/or interesting composite
33 * image processing operations, of the type that are often
34 * useful in applications. Most are morphological in
35 * nature.
36 *
37 * Extraction of boundary pixels
38 * PIX *pixExtractBoundary()
39 *
40 * Selective morph sequence operation under mask
41 * PIX *pixMorphSequenceMasked()
42 *
43 * Selective morph sequence operation on each component
44 * PIX *pixMorphSequenceByComponent()
45 * PIXA *pixaMorphSequenceByComponent()
46 *
47 * Selective morph sequence operation on each region
48 * PIX *pixMorphSequenceByRegion()
49 * PIXA *pixaMorphSequenceByRegion()
50 *
51 * Union and intersection of parallel composite operations
52 * PIX *pixUnionOfMorphOps()
53 * PIX *pixIntersectionOfMorphOps()
54 *
55 * Selective connected component filling
56 * PIX *pixSelectiveConnCompFill()
57 *
58 * Removal of matched patterns
59 * PIX *pixRemoveMatchedPattern()
60 *
61 * Display of matched patterns
62 * PIX *pixDisplayMatchedPattern()
63 *
64 * Extension of pixa by iterative erosion or dilation (and by scaling)
65 * PIXA *pixaExtendByMorph()
66 * PIXA *pixaExtendByScaling()
67 *
68 * Iterative morphological seed filling (don't use for real work)
69 * PIX *pixSeedfillMorph()
70 *
71 * Granulometry on binary images
72 * NUMA *pixRunHistogramMorph()
73 *
74 * Composite operations on grayscale images
75 * PIX *pixTophat()
76 * PIX *pixHDome()
77 * PIX *pixFastTophat()
78 * PIX *pixMorphGradient()
79 *
80 * Centroid of component
81 * PTA *pixaCentroids()
82 * l_int32 pixCentroid()
83 * </pre>
84 */
85
86 #include "allheaders.h"
87
88 #define SWAP(x, y) {temp = (x); (x) = (y); (y) = temp;}
89
90
91 /*-----------------------------------------------------------------*
92 * Extraction of boundary pixels *
93 *-----------------------------------------------------------------*/
94 /*!
95 * \brief pixExtractBoundary()
96 *
97 * \param[in] pixs 1 bpp
98 * \param[in] type 0 for background pixels; 1 for foreground pixels
99 * \return pixd, or NULL on error
100 *
101 * <pre>
102 * Notes:
103 * (1) Extracts the fg or bg boundary pixels for each component.
104 * Components are assumed to end at the boundary of pixs.
105 * </pre>
106 */
107 PIX *
pixExtractBoundary(PIX * pixs,l_int32 type)108 pixExtractBoundary(PIX *pixs,
109 l_int32 type)
110 {
111 PIX *pixd;
112
113 PROCNAME("pixExtractBoundary");
114
115 if (!pixs)
116 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
117
118 if (type == 0)
119 pixd = pixDilateBrick(NULL, pixs, 3, 3);
120 else
121 pixd = pixErodeBrick(NULL, pixs, 3, 3);
122 pixXor(pixd, pixd, pixs);
123 return pixd;
124 }
125
126
127 /*-----------------------------------------------------------------*
128 * Selective morph sequence operation under mask *
129 *-----------------------------------------------------------------*/
130 /*!
131 * \brief pixMorphSequenceMasked()
132 *
133 * \param[in] pixs 1 bpp
134 * \param[in] pixm [optional] 1 bpp mask
135 * \param[in] sequence string specifying sequence of operations
136 * \param[in] dispsep horizontal separation in pixels between
137 * successive displays; use zero to suppress display
138 * \return pixd, or NULL on error
139 *
140 * <pre>
141 * Notes:
142 * (1) This applies the morph sequence to the image, but only allows
143 * changes in pixs for pixels under the background of pixm.
144 * (5) If pixm is NULL, this is just pixMorphSequence().
145 * </pre>
146 */
147 PIX *
pixMorphSequenceMasked(PIX * pixs,PIX * pixm,const char * sequence,l_int32 dispsep)148 pixMorphSequenceMasked(PIX *pixs,
149 PIX *pixm,
150 const char *sequence,
151 l_int32 dispsep)
152 {
153 PIX *pixd;
154
155 PROCNAME("pixMorphSequenceMasked");
156
157 if (!pixs)
158 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
159 if (!sequence)
160 return (PIX *)ERROR_PTR("sequence not defined", procName, NULL);
161
162 pixd = pixMorphSequence(pixs, sequence, dispsep);
163 pixCombineMasked(pixd, pixs, pixm); /* restore src pixels under mask fg */
164 return pixd;
165 }
166
167
168 /*-----------------------------------------------------------------*
169 * Morph sequence operation on each component *
170 *-----------------------------------------------------------------*/
171 /*!
172 * \brief pixMorphSequenceByComponent()
173 *
174 * \param[in] pixs 1 bpp
175 * \param[in] sequence string specifying sequence
176 * \param[in] connectivity 4 or 8
177 * \param[in] minw minimum width to consider; use 0 or 1 for any width
178 * \param[in] minh minimum height to consider; use 0 or 1 for any height
179 * \param[out] pboxa [optional] return boxa of c.c. in pixs
180 * \return pixd, or NULL on error
181 *
182 * <pre>
183 * Notes:
184 * (1) See pixMorphSequence() for composing operation sequences.
185 * (2) This operates separately on each c.c. in the input pix.
186 * (3) The dilation does NOT increase the c.c. size; it is clipped
187 * to the size of the original c.c. This is necessary to
188 * keep the c.c. independent after the operation.
189 * (4) You can specify that the width and/or height must equal
190 * or exceed a minimum size for the operation to take place.
191 * (5) Use NULL for boxa to avoid returning the boxa.
192 * </pre>
193 */
194 PIX *
pixMorphSequenceByComponent(PIX * pixs,const char * sequence,l_int32 connectivity,l_int32 minw,l_int32 minh,BOXA ** pboxa)195 pixMorphSequenceByComponent(PIX *pixs,
196 const char *sequence,
197 l_int32 connectivity,
198 l_int32 minw,
199 l_int32 minh,
200 BOXA **pboxa)
201 {
202 l_int32 n, i, x, y, w, h;
203 BOXA *boxa;
204 PIX *pix, *pixd;
205 PIXA *pixas, *pixad;
206
207 PROCNAME("pixMorphSequenceByComponent");
208
209 if (!pixs)
210 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
211 if (!sequence)
212 return (PIX *)ERROR_PTR("sequence not defined", procName, NULL);
213
214 if (minw <= 0) minw = 1;
215 if (minh <= 0) minh = 1;
216
217 /* Get the c.c. */
218 if ((boxa = pixConnComp(pixs, &pixas, connectivity)) == NULL)
219 return (PIX *)ERROR_PTR("boxa not made", procName, NULL);
220
221 /* Operate on each c.c. independently */
222 pixad = pixaMorphSequenceByComponent(pixas, sequence, minw, minh);
223 pixaDestroy(&pixas);
224 boxaDestroy(&boxa);
225 if (!pixad)
226 return (PIX *)ERROR_PTR("pixad not made", procName, NULL);
227
228 /* Display the result out into pixd */
229 pixd = pixCreateTemplate(pixs);
230 n = pixaGetCount(pixad);
231 for (i = 0; i < n; i++) {
232 pixaGetBoxGeometry(pixad, i, &x, &y, &w, &h);
233 pix = pixaGetPix(pixad, i, L_CLONE);
234 pixRasterop(pixd, x, y, w, h, PIX_PAINT, pix, 0, 0);
235 pixDestroy(&pix);
236 }
237
238 if (pboxa)
239 *pboxa = pixaGetBoxa(pixad, L_CLONE);
240 pixaDestroy(&pixad);
241 return pixd;
242 }
243
244
245 /*!
246 * \brief pixaMorphSequenceByComponent()
247 *
248 * \param[in] pixas of 1 bpp pix
249 * \param[in] sequence string specifying sequence
250 * \param[in] minw minimum width to consider; use 0 or 1 for any width
251 * \param[in] minh minimum height to consider; use 0 or 1 for any height
252 * \return pixad, or NULL on error
253 *
254 * <pre>
255 * Notes:
256 * (1) See pixMorphSequence() for composing operation sequences.
257 * (2) This operates separately on each c.c. in the input pixa.
258 * (3) You can specify that the width and/or height must equal
259 * or exceed a minimum size for the operation to take place.
260 * (4) The input pixa should have a boxa giving the locations
261 * of the pix components.
262 * </pre>
263 */
264 PIXA *
pixaMorphSequenceByComponent(PIXA * pixas,const char * sequence,l_int32 minw,l_int32 minh)265 pixaMorphSequenceByComponent(PIXA *pixas,
266 const char *sequence,
267 l_int32 minw,
268 l_int32 minh)
269 {
270 l_int32 n, i, w, h, d;
271 BOX *box;
272 PIX *pix1, *pix2;
273 PIXA *pixad;
274
275 PROCNAME("pixaMorphSequenceByComponent");
276
277 if (!pixas)
278 return (PIXA *)ERROR_PTR("pixas not defined", procName, NULL);
279 if ((n = pixaGetCount(pixas)) == 0)
280 return (PIXA *)ERROR_PTR("no pix in pixas", procName, NULL);
281 if (n != pixaGetBoxaCount(pixas))
282 L_WARNING("boxa size != n\n", procName);
283 pixaGetPixDimensions(pixas, 0, NULL, NULL, &d);
284 if (d != 1)
285 return (PIXA *)ERROR_PTR("depth not 1 bpp", procName, NULL);
286
287 if (!sequence)
288 return (PIXA *)ERROR_PTR("sequence not defined", procName, NULL);
289 if (minw <= 0) minw = 1;
290 if (minh <= 0) minh = 1;
291
292 if ((pixad = pixaCreate(n)) == NULL)
293 return (PIXA *)ERROR_PTR("pixad not made", procName, NULL);
294 for (i = 0; i < n; i++) {
295 pixaGetPixDimensions(pixas, i, &w, &h, NULL);
296 if (w >= minw && h >= minh) {
297 if ((pix1 = pixaGetPix(pixas, i, L_CLONE)) == NULL) {
298 pixaDestroy(&pixad);
299 return (PIXA *)ERROR_PTR("pix1 not found", procName, NULL);
300 }
301 if ((pix2 = pixMorphCompSequence(pix1, sequence, 0)) == NULL) {
302 pixaDestroy(&pixad);
303 return (PIXA *)ERROR_PTR("pix2 not made", procName, NULL);
304 }
305 pixaAddPix(pixad, pix2, L_INSERT);
306 box = pixaGetBox(pixas, i, L_COPY);
307 pixaAddBox(pixad, box, L_INSERT);
308 pixDestroy(&pix1);
309 }
310 }
311
312 return pixad;
313 }
314
315
316 /*-----------------------------------------------------------------*
317 * Morph sequence operation on each region *
318 *-----------------------------------------------------------------*/
319 /*!
320 * \brief pixMorphSequenceByRegion()
321 *
322 * \param[in] pixs 1 bpp
323 * \param[in] pixm mask specifying regions
324 * \param[in] sequence string specifying sequence
325 * \param[in] connectivity 4 or 8, used on mask
326 * \param[in] minw minimum width to consider; use 0 or 1 for any width
327 * \param[in] minh minimum height to consider; use 0 or 1 for any height
328 * \param[out] pboxa [optional] return boxa of c.c. in pixm
329 * \return pixd, or NULL on error
330 *
331 * <pre>
332 * Notes:
333 * (1) See pixMorphCompSequence() for composing operation sequences.
334 * (2) This operates separately on the region in pixs corresponding
335 * to each c.c. in the mask pixm. It differs from
336 * pixMorphSequenceByComponent() in that the latter does not have
337 * a pixm (mask), but instead operates independently on each
338 * component in pixs.
339 * (3) Dilation will NOT increase the region size; the result
340 * is clipped to the size of the mask region. This is necessary
341 * to make regions independent after the operation.
342 * (4) You can specify that the width and/or height of a region must
343 * equal or exceed a minimum size for the operation to take place.
344 * (5) Use NULL for %pboxa to avoid returning the boxa.
345 * </pre>
346 */
347 PIX *
pixMorphSequenceByRegion(PIX * pixs,PIX * pixm,const char * sequence,l_int32 connectivity,l_int32 minw,l_int32 minh,BOXA ** pboxa)348 pixMorphSequenceByRegion(PIX *pixs,
349 PIX *pixm,
350 const char *sequence,
351 l_int32 connectivity,
352 l_int32 minw,
353 l_int32 minh,
354 BOXA **pboxa)
355 {
356 l_int32 n, i, x, y, w, h;
357 BOXA *boxa;
358 PIX *pix, *pixd;
359 PIXA *pixam, *pixad;
360
361 PROCNAME("pixMorphSequenceByRegion");
362
363 if (pboxa) *pboxa = NULL;
364 if (!pixs)
365 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
366 if (!pixm)
367 return (PIX *)ERROR_PTR("pixm not defined", procName, NULL);
368 if (pixGetDepth(pixs) != 1 || pixGetDepth(pixm) != 1)
369 return (PIX *)ERROR_PTR("pixs and pixm not both 1 bpp", procName, NULL);
370 if (!sequence)
371 return (PIX *)ERROR_PTR("sequence not defined", procName, NULL);
372
373 if (minw <= 0) minw = 1;
374 if (minh <= 0) minh = 1;
375
376 /* Get the c.c. of the mask */
377 if ((boxa = pixConnComp(pixm, &pixam, connectivity)) == NULL)
378 return (PIX *)ERROR_PTR("boxa not made", procName, NULL);
379
380 /* Operate on each region in pixs independently */
381 pixad = pixaMorphSequenceByRegion(pixs, pixam, sequence, minw, minh);
382 pixaDestroy(&pixam);
383 boxaDestroy(&boxa);
384 if (!pixad)
385 return (PIX *)ERROR_PTR("pixad not made", procName, NULL);
386
387 /* Display the result out into pixd */
388 pixd = pixCreateTemplate(pixs);
389 n = pixaGetCount(pixad);
390 for (i = 0; i < n; i++) {
391 pixaGetBoxGeometry(pixad, i, &x, &y, &w, &h);
392 pix = pixaGetPix(pixad, i, L_CLONE);
393 pixRasterop(pixd, x, y, w, h, PIX_PAINT, pix, 0, 0);
394 pixDestroy(&pix);
395 }
396
397 if (pboxa)
398 *pboxa = pixaGetBoxa(pixad, L_CLONE);
399 pixaDestroy(&pixad);
400 return pixd;
401 }
402
403
404 /*!
405 * \brief pixaMorphSequenceByRegion()
406 *
407 * \param[in] pixs 1 bpp
408 * \param[in] pixam of 1 bpp mask elements
409 * \param[in] sequence string specifying sequence
410 * \param[in] minw minimum width to consider; use 0 or 1 for any width
411 * \param[in] minh minimum height to consider; use 0 or 1 for any height
412 * \return pixad, or NULL on error
413 *
414 * <pre>
415 * Notes:
416 * (1) See pixMorphSequence() for composing operation sequences.
417 * (2) This operates separately on each region in the input pixs
418 * defined by the components in pixam.
419 * (3) You can specify that the width and/or height of a mask
420 * component must equal or exceed a minimum size for the
421 * operation to take place.
422 * (4) The input pixam should have a boxa giving the locations
423 * of the regions in pixs.
424 * </pre>
425 */
426 PIXA *
pixaMorphSequenceByRegion(PIX * pixs,PIXA * pixam,const char * sequence,l_int32 minw,l_int32 minh)427 pixaMorphSequenceByRegion(PIX *pixs,
428 PIXA *pixam,
429 const char *sequence,
430 l_int32 minw,
431 l_int32 minh)
432 {
433 l_int32 n, i, w, h, same, maxd, fullpa, fullba;
434 BOX *box;
435 PIX *pix1, *pix2, *pix3;
436 PIXA *pixad;
437
438 PROCNAME("pixaMorphSequenceByRegion");
439
440 if (!pixs)
441 return (PIXA *)ERROR_PTR("pixs not defined", procName, NULL);
442 if (pixGetDepth(pixs) != 1)
443 return (PIXA *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
444 if (!sequence)
445 return (PIXA *)ERROR_PTR("sequence not defined", procName, NULL);
446 if (!pixam)
447 return (PIXA *)ERROR_PTR("pixam not defined", procName, NULL);
448 pixaVerifyDepth(pixam, &same, &maxd);
449 if (maxd != 1)
450 return (PIXA *)ERROR_PTR("mask depth not 1 bpp", procName, NULL);
451 pixaIsFull(pixam, &fullpa, &fullba);
452 if (!fullpa || !fullba)
453 return (PIXA *)ERROR_PTR("missing comps in pixam", procName, NULL);
454 n = pixaGetCount(pixam);
455 if (minw <= 0) minw = 1;
456 if (minh <= 0) minh = 1;
457
458 if ((pixad = pixaCreate(n)) == NULL)
459 return (PIXA *)ERROR_PTR("pixad not made", procName, NULL);
460
461 /* Use the rectangle to remove the appropriate part of pixs;
462 * then AND with the mask component to get the actual fg
463 * of pixs that is under the mask component. */
464 for (i = 0; i < n; i++) {
465 pixaGetPixDimensions(pixam, i, &w, &h, NULL);
466 if (w >= minw && h >= minh) {
467 pix1 = pixaGetPix(pixam, i, L_CLONE);
468 box = pixaGetBox(pixam, i, L_COPY);
469 pix2 = pixClipRectangle(pixs, box, NULL);
470 pixAnd(pix2, pix2, pix1);
471 pix3 = pixMorphCompSequence(pix2, sequence, 0);
472 pixDestroy(&pix1);
473 pixDestroy(&pix2);
474 if (!pix3) {
475 boxDestroy(&box);
476 pixaDestroy(&pixad);
477 L_ERROR("pix3 not made in iter %d; aborting\n", procName, i);
478 break;
479 }
480 pixaAddPix(pixad, pix3, L_INSERT);
481 pixaAddBox(pixad, box, L_INSERT);
482 }
483 }
484
485 return pixad;
486 }
487
488
489 /*-----------------------------------------------------------------*
490 * Union and intersection of parallel composite operations *
491 *-----------------------------------------------------------------*/
492 /*!
493 * \brief pixUnionOfMorphOps()
494 *
495 * \param[in] pixs binary
496 * \param[in] sela
497 * \param[in] type L_MORPH_DILATE, etc.
498 * \return pixd union of the specified morphological operation
499 * on pixs for each Sel in the Sela, or NULL on error
500 */
501 PIX *
pixUnionOfMorphOps(PIX * pixs,SELA * sela,l_int32 type)502 pixUnionOfMorphOps(PIX *pixs,
503 SELA *sela,
504 l_int32 type)
505 {
506 l_int32 n, i;
507 PIX *pixt, *pixd;
508 SEL *sel;
509
510 PROCNAME("pixUnionOfMorphOps");
511
512 if (!pixs || pixGetDepth(pixs) != 1)
513 return (PIX *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
514 if (!sela)
515 return (PIX *)ERROR_PTR("sela not defined", procName, NULL);
516 n = selaGetCount(sela);
517 if (n == 0)
518 return (PIX *)ERROR_PTR("no sels in sela", procName, NULL);
519 if (type != L_MORPH_DILATE && type != L_MORPH_ERODE &&
520 type != L_MORPH_OPEN && type != L_MORPH_CLOSE &&
521 type != L_MORPH_HMT)
522 return (PIX *)ERROR_PTR("invalid type", procName, NULL);
523
524 pixd = pixCreateTemplate(pixs);
525 for (i = 0; i < n; i++) {
526 sel = selaGetSel(sela, i);
527 if (type == L_MORPH_DILATE)
528 pixt = pixDilate(NULL, pixs, sel);
529 else if (type == L_MORPH_ERODE)
530 pixt = pixErode(NULL, pixs, sel);
531 else if (type == L_MORPH_OPEN)
532 pixt = pixOpen(NULL, pixs, sel);
533 else if (type == L_MORPH_CLOSE)
534 pixt = pixClose(NULL, pixs, sel);
535 else /* type == L_MORPH_HMT */
536 pixt = pixHMT(NULL, pixs, sel);
537 pixOr(pixd, pixd, pixt);
538 pixDestroy(&pixt);
539 }
540
541 return pixd;
542 }
543
544
545 /*!
546 * \brief pixIntersectionOfMorphOps()
547 *
548 * \param[in] pixs binary
549 * \param[in] sela
550 * \param[in] type L_MORPH_DILATE, etc.
551 * \return pixd intersection of the specified morphological operation
552 * on pixs for each Sel in the Sela, or NULL on error
553 */
554 PIX *
pixIntersectionOfMorphOps(PIX * pixs,SELA * sela,l_int32 type)555 pixIntersectionOfMorphOps(PIX *pixs,
556 SELA *sela,
557 l_int32 type)
558 {
559 l_int32 n, i;
560 PIX *pixt, *pixd;
561 SEL *sel;
562
563 PROCNAME("pixIntersectionOfMorphOps");
564
565 if (!pixs || pixGetDepth(pixs) != 1)
566 return (PIX *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
567 if (!sela)
568 return (PIX *)ERROR_PTR("sela not defined", procName, NULL);
569 n = selaGetCount(sela);
570 if (n == 0)
571 return (PIX *)ERROR_PTR("no sels in sela", procName, NULL);
572 if (type != L_MORPH_DILATE && type != L_MORPH_ERODE &&
573 type != L_MORPH_OPEN && type != L_MORPH_CLOSE &&
574 type != L_MORPH_HMT)
575 return (PIX *)ERROR_PTR("invalid type", procName, NULL);
576
577 pixd = pixCreateTemplate(pixs);
578 pixSetAll(pixd);
579 for (i = 0; i < n; i++) {
580 sel = selaGetSel(sela, i);
581 if (type == L_MORPH_DILATE)
582 pixt = pixDilate(NULL, pixs, sel);
583 else if (type == L_MORPH_ERODE)
584 pixt = pixErode(NULL, pixs, sel);
585 else if (type == L_MORPH_OPEN)
586 pixt = pixOpen(NULL, pixs, sel);
587 else if (type == L_MORPH_CLOSE)
588 pixt = pixClose(NULL, pixs, sel);
589 else /* type == L_MORPH_HMT */
590 pixt = pixHMT(NULL, pixs, sel);
591 pixAnd(pixd, pixd, pixt);
592 pixDestroy(&pixt);
593 }
594
595 return pixd;
596 }
597
598
599
600 /*-----------------------------------------------------------------*
601 * Selective connected component filling *
602 *-----------------------------------------------------------------*/
603 /*!
604 * \brief pixSelectiveConnCompFill()
605 *
606 * \param[in] pixs binary
607 * \param[in] connectivity 4 or 8
608 * \param[in] minw minimum width to consider; use 0 or 1 for any width
609 * \param[in] minh minimum height to consider; use 0 or 1 for any height
610 * \return pix with holes filled in selected c.c., or NULL on error
611 */
612 PIX *
pixSelectiveConnCompFill(PIX * pixs,l_int32 connectivity,l_int32 minw,l_int32 minh)613 pixSelectiveConnCompFill(PIX *pixs,
614 l_int32 connectivity,
615 l_int32 minw,
616 l_int32 minh)
617 {
618 l_int32 n, i, x, y, w, h;
619 BOXA *boxa;
620 PIX *pix1, *pix2, *pixd;
621 PIXA *pixa;
622
623 PROCNAME("pixSelectiveConnCompFill");
624
625 if (!pixs)
626 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
627 if (pixGetDepth(pixs) != 1)
628 return (PIX *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
629 if (minw <= 0) minw = 1;
630 if (minh <= 0) minh = 1;
631
632 if ((boxa = pixConnComp(pixs, &pixa, connectivity)) == NULL)
633 return (PIX *)ERROR_PTR("boxa not made", procName, NULL);
634 n = boxaGetCount(boxa);
635 pixd = pixCopy(NULL, pixs);
636 for (i = 0; i < n; i++) {
637 boxaGetBoxGeometry(boxa, i, &x, &y, &w, &h);
638 if (w >= minw && h >= minh) {
639 pix1 = pixaGetPix(pixa, i, L_CLONE);
640 if ((pix2 = pixHolesByFilling(pix1, 12 - connectivity)) == NULL) {
641 L_ERROR("pix2 not made in iter %d\n", procName, i);
642 pixDestroy(&pix1);
643 continue;
644 }
645 pixRasterop(pixd, x, y, w, h, PIX_PAINT, pix2, 0, 0);
646 pixDestroy(&pix1);
647 pixDestroy(&pix2);
648 }
649 }
650 pixaDestroy(&pixa);
651 boxaDestroy(&boxa);
652
653 return pixd;
654 }
655
656
657 /*-----------------------------------------------------------------*
658 * Removal of matched patterns *
659 *-----------------------------------------------------------------*/
660 /*!
661 * \brief pixRemoveMatchedPattern()
662 *
663 * \param[in] pixs input image, 1 bpp
664 * \param[in] pixp pattern to be removed from image, 1 bpp
665 * \param[in] pixe image after erosion by Sel that approximates pixp, 1 bpp
666 * \param[in] x0, y0 center of Sel
667 * \param[in] dsize number of pixels on each side by which pixp is
668 * dilated before being subtracted from pixs;
669 * valid values are {0, 1, 2, 3, 4}
670 * \return 0 if OK, 1 on error
671 *
672 * <pre>
673 * Notes:
674 * (1) This is in-place.
675 * (2) You can use various functions in selgen to create a Sel
676 * that is used to generate pixe from pixs.
677 * (3) This function is applied after pixe has been computed.
678 * It finds the centroid of each c.c., and subtracts
679 * (the appropriately dilated version of) pixp, with the center
680 * of the Sel used to align pixp with pixs.
681 * </pre>
682 */
683 l_int32
pixRemoveMatchedPattern(PIX * pixs,PIX * pixp,PIX * pixe,l_int32 x0,l_int32 y0,l_int32 dsize)684 pixRemoveMatchedPattern(PIX *pixs,
685 PIX *pixp,
686 PIX *pixe,
687 l_int32 x0,
688 l_int32 y0,
689 l_int32 dsize)
690 {
691 l_int32 i, nc, x, y, w, h, xb, yb;
692 BOXA *boxa;
693 PIX *pix1, *pix2;
694 PIXA *pixa;
695 PTA *pta;
696 SEL *sel;
697
698 PROCNAME("pixRemoveMatchedPattern");
699
700 if (!pixs)
701 return ERROR_INT("pixs not defined", procName, 1);
702 if (!pixp)
703 return ERROR_INT("pixp not defined", procName, 1);
704 if (!pixe)
705 return ERROR_INT("pixe not defined", procName, 1);
706 if (pixGetDepth(pixs) != 1 || pixGetDepth(pixp) != 1 ||
707 pixGetDepth(pixe) != 1)
708 return ERROR_INT("all input pix not 1 bpp", procName, 1);
709 if (dsize < 0 || dsize > 4)
710 return ERROR_INT("dsize not in {0,1,2,3,4}", procName, 1);
711
712 /* Find the connected components and their centroids */
713 boxa = pixConnComp(pixe, &pixa, 8);
714 if ((nc = boxaGetCount(boxa)) == 0) {
715 L_WARNING("no matched patterns\n", procName);
716 boxaDestroy(&boxa);
717 pixaDestroy(&pixa);
718 return 0;
719 }
720 pta = pixaCentroids(pixa);
721 pixaDestroy(&pixa);
722
723 /* Optionally dilate the pattern, first adding a border that
724 * is large enough to accommodate the dilated pixels */
725 sel = NULL;
726 if (dsize > 0) {
727 sel = selCreateBrick(2 * dsize + 1, 2 * dsize + 1, dsize, dsize,
728 SEL_HIT);
729 pix1 = pixAddBorder(pixp, dsize, 0);
730 pix2 = pixDilate(NULL, pix1, sel);
731 selDestroy(&sel);
732 pixDestroy(&pix1);
733 } else {
734 pix2 = pixClone(pixp);
735 }
736
737 /* Subtract out each dilated pattern. The centroid of each
738 * component is located at:
739 * (box->x + x, box->y + y)
740 * and the 'center' of the pattern used in making pixe is located at
741 * (x0 + dsize, (y0 + dsize)
742 * relative to the UL corner of the pattern. The center of the
743 * pattern is placed at the center of the component. */
744 pixGetDimensions(pix2, &w, &h, NULL);
745 for (i = 0; i < nc; i++) {
746 ptaGetIPt(pta, i, &x, &y);
747 boxaGetBoxGeometry(boxa, i, &xb, &yb, NULL, NULL);
748 pixRasterop(pixs, xb + x - x0 - dsize, yb + y - y0 - dsize,
749 w, h, PIX_DST & PIX_NOT(PIX_SRC), pix2, 0, 0);
750 }
751
752 boxaDestroy(&boxa);
753 ptaDestroy(&pta);
754 pixDestroy(&pix2);
755 return 0;
756 }
757
758
759 /*-----------------------------------------------------------------*
760 * Display of matched patterns *
761 *-----------------------------------------------------------------*/
762 /*!
763 * \brief pixDisplayMatchedPattern()
764 *
765 * \param[in] pixs input image, 1 bpp
766 * \param[in] pixp pattern to be removed from image, 1 bpp
767 * \param[in] pixe image after erosion by Sel that approximates pixp, 1 bpp
768 * \param[in] x0, y0 center of Sel
769 * \param[in] color to paint the matched patterns; 0xrrggbb00
770 * \param[in] scale reduction factor for output pixd
771 * \param[in] nlevels if scale < 1.0, threshold to this number of levels
772 * \return pixd 8 bpp, colormapped, or NULL on error
773 *
774 * <pre>
775 * Notes:
776 * (1) A 4 bpp colormapped image is generated.
777 * (2) If scale <= 1.0, do scale to gray for the output, and threshold
778 * to nlevels of gray.
779 * (3) You can use various functions in selgen to create a Sel
780 * that will generate pixe from pixs.
781 * (4) This function is applied after pixe has been computed.
782 * It finds the centroid of each c.c., and colors the output
783 * pixels using pixp (appropriately aligned) as a stencil.
784 * Alignment is done using the origin of the Sel and the
785 * centroid of the eroded image to place the stencil pixp.
786 * </pre>
787 */
788 PIX *
pixDisplayMatchedPattern(PIX * pixs,PIX * pixp,PIX * pixe,l_int32 x0,l_int32 y0,l_uint32 color,l_float32 scale,l_int32 nlevels)789 pixDisplayMatchedPattern(PIX *pixs,
790 PIX *pixp,
791 PIX *pixe,
792 l_int32 x0,
793 l_int32 y0,
794 l_uint32 color,
795 l_float32 scale,
796 l_int32 nlevels)
797 {
798 l_int32 i, nc, xb, yb, x, y, xi, yi, rval, gval, bval;
799 BOXA *boxa;
800 PIX *pixd, *pixt, *pixps;
801 PIXA *pixa;
802 PTA *pta;
803 PIXCMAP *cmap;
804
805 PROCNAME("pixDisplayMatchedPattern");
806
807 if (!pixs)
808 return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
809 if (!pixp)
810 return (PIX *)ERROR_PTR("pixp not defined", procName, NULL);
811 if (!pixe)
812 return (PIX *)ERROR_PTR("pixe not defined", procName, NULL);
813 if (pixGetDepth(pixs) != 1 || pixGetDepth(pixp) != 1 ||
814 pixGetDepth(pixe) != 1)
815 return (PIX *)ERROR_PTR("all input pix not 1 bpp", procName, NULL);
816 if (scale > 1.0 || scale <= 0.0) {
817 L_WARNING("scale > 1.0 or < 0.0; setting to 1.0\n", procName);
818 scale = 1.0;
819 }
820
821 /* Find the connected components and their centroids */
822 boxa = pixConnComp(pixe, &pixa, 8);
823 if ((nc = boxaGetCount(boxa)) == 0) {
824 L_WARNING("no matched patterns\n", procName);
825 boxaDestroy(&boxa);
826 pixaDestroy(&pixa);
827 return 0;
828 }
829 pta = pixaCentroids(pixa);
830
831 extractRGBValues(color, &rval, &gval, &bval);
832 if (scale == 1.0) { /* output 4 bpp at full resolution */
833 pixd = pixConvert1To4(NULL, pixs, 0, 1);
834 cmap = pixcmapCreate(4);
835 pixcmapAddColor(cmap, 255, 255, 255);
836 pixcmapAddColor(cmap, 0, 0, 0);
837 pixSetColormap(pixd, cmap);
838
839 /* Paint through pixp for each match location. The centroid of each
840 * component in pixe is located at:
841 * (box->x + x, box->y + y)
842 * and the 'center' of the pattern used in making pixe is located at
843 * (x0, y0)
844 * relative to the UL corner of the pattern. The center of the
845 * pattern is placed at the center of the component. */
846 for (i = 0; i < nc; i++) {
847 ptaGetIPt(pta, i, &x, &y);
848 boxaGetBoxGeometry(boxa, i, &xb, &yb, NULL, NULL);
849 pixSetMaskedCmap(pixd, pixp, xb + x - x0, yb + y - y0,
850 rval, gval, bval);
851 }
852 } else { /* output 4 bpp downscaled */
853 pixt = pixScaleToGray(pixs, scale);
854 pixd = pixThresholdTo4bpp(pixt, nlevels, 1);
855 pixps = pixScaleBySampling(pixp, scale, scale);
856
857 for (i = 0; i < nc; i++) {
858 ptaGetIPt(pta, i, &x, &y);
859 boxaGetBoxGeometry(boxa, i, &xb, &yb, NULL, NULL);
860 xi = (l_int32)(scale * (xb + x - x0));
861 yi = (l_int32)(scale * (yb + y - y0));
862 pixSetMaskedCmap(pixd, pixps, xi, yi, rval, gval, bval);
863 }
864 pixDestroy(&pixt);
865 pixDestroy(&pixps);
866 }
867
868 boxaDestroy(&boxa);
869 pixaDestroy(&pixa);
870 ptaDestroy(&pta);
871 return pixd;
872 }
873
874
875 /*------------------------------------------------------------------------*
876 * Extension of pixa by iterative erosion or dilation (and by scaling) *
877 *------------------------------------------------------------------------*/
878 /*!
879 * \brief pixaExtendByMorph()
880 *
881 * \param[in] pixas
882 * \param[in] type L_MORPH_DILATE, L_MORPH_ERODE
883 * \param[in] niters
884 * \param[in] sel used for dilation, erosion; uses 2x2 if null
885 * \param[in] include 1 to include a copy of the input pixas in pixad;
886 * 0 to omit
887 * \return pixad with derived pix, using all iterations, or NULL on error
888 *
889 * <pre>
890 * Notes:
891 * (1) This dilates or erodes every pix in %pixas, iteratively,
892 * using the input Sel (or, if null, a 2x2 Sel by default),
893 * and puts the results in %pixad.
894 * (2) If %niters <= 0, this is a no-op; it returns a clone of pixas.
895 * (3) If %include == 1, the output %pixad contains all the pix
896 * in %pixas. Otherwise, it doesn't, but pixaJoin() can be
897 * used later to join pixas with pixad.
898 * </pre>
899 */
900 PIXA *
pixaExtendByMorph(PIXA * pixas,l_int32 type,l_int32 niters,SEL * sel,l_int32 include)901 pixaExtendByMorph(PIXA *pixas,
902 l_int32 type,
903 l_int32 niters,
904 SEL *sel,
905 l_int32 include)
906 {
907 l_int32 maxdepth, i, j, n;
908 PIX *pix0, *pix1, *pix2;
909 SEL *selt;
910 PIXA *pixad;
911
912 PROCNAME("pixaExtendByMorph");
913
914 if (!pixas)
915 return (PIXA *)ERROR_PTR("pixas undefined", procName, NULL);
916 if (niters <= 0) {
917 L_INFO("niters = %d; nothing to do\n", procName, niters);
918 return pixaCopy(pixas, L_CLONE);
919 }
920 if (type != L_MORPH_DILATE && type != L_MORPH_ERODE)
921 return (PIXA *)ERROR_PTR("invalid type", procName, NULL);
922 pixaGetDepthInfo(pixas, &maxdepth, NULL);
923 if (maxdepth > 1)
924 return (PIXA *)ERROR_PTR("some pix have bpp > 1", procName, NULL);
925
926 if (!sel)
927 selt = selCreateBrick(2, 2, 0, 0, SEL_HIT); /* default */
928 else
929 selt = selCopy(sel);
930 n = pixaGetCount(pixas);
931 pixad = pixaCreate(n * niters);
932 for (i = 0; i < n; i++) {
933 pix1 = pixaGetPix(pixas, i, L_CLONE);
934 if (include) pixaAddPix(pixad, pix1, L_COPY);
935 pix0 = pix1; /* need to keep the handle to destroy the clone */
936 for (j = 0; j < niters; j++) {
937 if (type == L_MORPH_DILATE) {
938 pix2 = pixDilate(NULL, pix1, selt);
939 } else { /* L_MORPH_ERODE */
940 pix2 = pixErode(NULL, pix1, selt);
941 }
942 pixaAddPix(pixad, pix2, L_INSERT);
943 pix1 = pix2; /* owned by pixad; do not destroy */
944 }
945 pixDestroy(&pix0);
946 }
947
948 selDestroy(&selt);
949 return pixad;
950 }
951
952
953 /*!
954 * \brief pixaExtendByScaling()
955 *
956 * \param[in] pixas
957 * \param[in] nasc numa of scaling factors
958 * \param[in] type L_HORIZ, L_VERT, L_BOTH_DIRECTIONS
959 * \param[in] include 1 to include a copy of the input pixas in pixad;
960 * 0 to omit
961 * \return pixad with derived pix, using all scalings, or NULL on error
962 *
963 * <pre>
964 * Notes:
965 * (1) This scales every pix in %pixas by each factor in %nasc.
966 * and puts the results in %pixad.
967 * (2) If %include == 1, the output %pixad contains all the pix
968 * in %pixas. Otherwise, it doesn't, but pixaJoin() can be
969 * used later to join pixas with pixad.
970 * </pre>
971 */
972 PIXA *
pixaExtendByScaling(PIXA * pixas,NUMA * nasc,l_int32 type,l_int32 include)973 pixaExtendByScaling(PIXA *pixas,
974 NUMA *nasc,
975 l_int32 type,
976 l_int32 include)
977 {
978 l_int32 i, j, n, nsc, w, h, scalew, scaleh;
979 l_float32 scalefact;
980 PIX *pix1, *pix2;
981 PIXA *pixad;
982
983 PROCNAME("pixaExtendByScaling");
984
985 if (!pixas)
986 return (PIXA *)ERROR_PTR("pixas undefined", procName, NULL);
987 if (!nasc || numaGetCount(nasc) == 0)
988 return (PIXA *)ERROR_PTR("nasc undefined or empty", procName, NULL);
989 if (type != L_HORIZ && type != L_VERT && type != L_BOTH_DIRECTIONS)
990 return (PIXA *)ERROR_PTR("invalid type", procName, NULL);
991
992 n = pixaGetCount(pixas);
993 nsc = numaGetCount(nasc);
994 if ((pixad = pixaCreate(n * (nsc + 1))) == NULL) {
995 L_ERROR("pixad not made: n = %d, nsc = %d\n", procName, n, nsc);
996 return NULL;
997 }
998 for (i = 0; i < n; i++) {
999 pix1 = pixaGetPix(pixas, i, L_CLONE);
1000 if (include) pixaAddPix(pixad, pix1, L_COPY);
1001 pixGetDimensions(pix1, &w, &h, NULL);
1002 for (j = 0; j < nsc; j++) {
1003 numaGetFValue(nasc, j, &scalefact);
1004 scalew = w;
1005 scaleh = h;
1006 if (type == L_HORIZ || type == L_BOTH_DIRECTIONS)
1007 scalew = w * scalefact;
1008 if (type == L_VERT || type == L_BOTH_DIRECTIONS)
1009 scaleh = h * scalefact;
1010 pix2 = pixScaleToSize(pix1, scalew, scaleh);
1011 pixaAddPix(pixad, pix2, L_INSERT);
1012 }
1013 pixDestroy(&pix1);
1014 }
1015 return pixad;
1016 }
1017
1018
1019 /*-----------------------------------------------------------------*
1020 * Iterative morphological seed filling *
1021 *-----------------------------------------------------------------*/
1022 /*!
1023 * \brief pixSeedfillMorph()
1024 *
1025 * \param[in] pixs seed
1026 * \param[in] pixm mask
1027 * \param[in] maxiters use 0 to go to completion
1028 * \param[in] connectivity 4 or 8
1029 * \return pixd after filling into the mask or NULL on error
1030 *
1031 * <pre>
1032 * Notes:
1033 * (1) This is in general a very inefficient method for filling
1034 * from a seed into a mask. Use it for a small number of iterations,
1035 * but if you expect more than a few iterations, use
1036 * pixSeedfillBinary().
1037 * (2) We use a 3x3 brick SEL for 8-cc filling and a 3x3 plus SEL for 4-cc.
1038 * </pre>
1039 */
1040 PIX *
pixSeedfillMorph(PIX * pixs,PIX * pixm,l_int32 maxiters,l_int32 connectivity)1041 pixSeedfillMorph(PIX *pixs,
1042 PIX *pixm,
1043 l_int32 maxiters,
1044 l_int32 connectivity)
1045 {
1046 l_int32 same, i;
1047 PIX *pixt, *pixd, *temp;
1048 SEL *sel_3;
1049
1050 PROCNAME("pixSeedfillMorph");
1051
1052 if (!pixs || pixGetDepth(pixs) != 1)
1053 return (PIX *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
1054 if (!pixm)
1055 return (PIX *)ERROR_PTR("mask pix not defined", procName, NULL);
1056 if (connectivity != 4 && connectivity != 8)
1057 return (PIX *)ERROR_PTR("connectivity not in {4,8}", procName, NULL);
1058 if (maxiters <= 0) maxiters = 1000;
1059 if (pixSizesEqual(pixs, pixm) == 0)
1060 return (PIX *)ERROR_PTR("pix sizes unequal", procName, NULL);
1061
1062 if ((sel_3 = selCreateBrick(3, 3, 1, 1, SEL_HIT)) == NULL)
1063 return (PIX *)ERROR_PTR("sel_3 not made", procName, NULL);
1064 if (connectivity == 4) { /* remove corner hits to make a '+' */
1065 selSetElement(sel_3, 0, 0, SEL_DONT_CARE);
1066 selSetElement(sel_3, 2, 2, SEL_DONT_CARE);
1067 selSetElement(sel_3, 2, 0, SEL_DONT_CARE);
1068 selSetElement(sel_3, 0, 2, SEL_DONT_CARE);
1069 }
1070
1071 pixt = pixCopy(NULL, pixs);
1072 pixd = pixCreateTemplate(pixs);
1073 for (i = 1; i <= maxiters; i++) {
1074 pixDilate(pixd, pixt, sel_3);
1075 pixAnd(pixd, pixd, pixm);
1076 pixEqual(pixd, pixt, &same);
1077 if (same || i == maxiters)
1078 break;
1079 else
1080 SWAP(pixt, pixd);
1081 }
1082 fprintf(stderr, " Num iters in binary reconstruction = %d\n", i);
1083
1084 pixDestroy(&pixt);
1085 selDestroy(&sel_3);
1086 return pixd;
1087 }
1088
1089
1090 /*-----------------------------------------------------------------*
1091 * Granulometry on binary images *
1092 *-----------------------------------------------------------------*/
1093 /*!
1094 * \brief pixRunHistogramMorph()
1095 *
1096 * \param[in] pixs
1097 * \param[in] runtype L_RUN_OFF, L_RUN_ON
1098 * \param[in] direction L_HORIZ, L_VERT
1099 * \param[in] maxsize size of largest runlength counted
1100 * \return numa of run-lengths
1101 */
1102 NUMA *
pixRunHistogramMorph(PIX * pixs,l_int32 runtype,l_int32 direction,l_int32 maxsize)1103 pixRunHistogramMorph(PIX *pixs,
1104 l_int32 runtype,
1105 l_int32 direction,
1106 l_int32 maxsize)
1107 {
1108 l_int32 count, i, size;
1109 l_float32 val;
1110 NUMA *na, *nah;
1111 PIX *pix1, *pix2, *pix3;
1112 SEL *sel_2a;
1113
1114 PROCNAME("pixRunHistogramMorph");
1115
1116 if (!pixs)
1117 return (NUMA *)ERROR_PTR("seed pix not defined", procName, NULL);
1118 if (runtype != L_RUN_OFF && runtype != L_RUN_ON)
1119 return (NUMA *)ERROR_PTR("invalid run type", procName, NULL);
1120 if (direction != L_HORIZ && direction != L_VERT)
1121 return (NUMA *)ERROR_PTR("direction not in {L_HORIZ, L_VERT}",
1122 procName, NULL);
1123 if (pixGetDepth(pixs) != 1)
1124 return (NUMA *)ERROR_PTR("pixs must be binary", procName, NULL);
1125
1126 if (direction == L_HORIZ)
1127 sel_2a = selCreateBrick(1, 2, 0, 0, SEL_HIT);
1128 else /* direction == L_VERT */
1129 sel_2a = selCreateBrick(2, 1, 0, 0, SEL_HIT);
1130 if (!sel_2a)
1131 return (NUMA *)ERROR_PTR("sel_2a not made", procName, NULL);
1132
1133 if (runtype == L_RUN_OFF) {
1134 if ((pix1 = pixCopy(NULL, pixs)) == NULL) {
1135 selDestroy(&sel_2a);
1136 return (NUMA *)ERROR_PTR("pix1 not made", procName, NULL);
1137 }
1138 pixInvert(pix1, pix1);
1139 } else { /* runtype == L_RUN_ON */
1140 pix1 = pixClone(pixs);
1141 }
1142
1143 /* Get pixel counts at different stages of erosion */
1144 na = numaCreate(0);
1145 pix2 = pixCreateTemplate(pixs);
1146 pix3 = pixCreateTemplate(pixs);
1147 pixCountPixels(pix1, &count, NULL);
1148 numaAddNumber(na, count);
1149 pixErode(pix2, pix1, sel_2a);
1150 pixCountPixels(pix2, &count, NULL);
1151 numaAddNumber(na, count);
1152 for (i = 0; i < maxsize / 2; i++) {
1153 pixErode(pix3, pix2, sel_2a);
1154 pixCountPixels(pix3, &count, NULL);
1155 numaAddNumber(na, count);
1156 pixErode(pix2, pix3, sel_2a);
1157 pixCountPixels(pix2, &count, NULL);
1158 numaAddNumber(na, count);
1159 }
1160
1161 /* Compute length histogram */
1162 size = numaGetCount(na);
1163 nah = numaCreate(size);
1164 numaAddNumber(nah, 0); /* number at length 0 */
1165 for (i = 1; i < size - 1; i++) {
1166 val = na->array[i+1] - 2 * na->array[i] + na->array[i-1];
1167 numaAddNumber(nah, val);
1168 }
1169
1170 pixDestroy(&pix1);
1171 pixDestroy(&pix2);
1172 pixDestroy(&pix3);
1173 selDestroy(&sel_2a);
1174 numaDestroy(&na);
1175 return nah;
1176 }
1177
1178
1179 /*-----------------------------------------------------------------*
1180 * Composite operations on grayscale images *
1181 *-----------------------------------------------------------------*/
1182 /*!
1183 * \brief pixTophat()
1184 *
1185 * \param[in] pixs
1186 * \param[in] hsize of Sel; must be odd; origin implicitly in center
1187 * \param[in] vsize ditto
1188 * \param[in] type L_TOPHAT_WHITE: image - opening
1189 * L_TOPHAT_BLACK: closing - image
1190 * \return pixd, or NULL on error
1191 *
1192 * <pre>
1193 * Notes:
1194 * (1) Sel is a brick with all elements being hits
1195 * (2) If hsize = vsize = 1, returns an image with all 0 data.
1196 * (3) The L_TOPHAT_WHITE flag emphasizes small bright regions,
1197 * whereas the L_TOPHAT_BLACK flag emphasizes small dark regions.
1198 * The L_TOPHAT_WHITE tophat can be accomplished by doing a
1199 * L_TOPHAT_BLACK tophat on the inverse, or v.v.
1200 * </pre>
1201 */
1202 PIX *
pixTophat(PIX * pixs,l_int32 hsize,l_int32 vsize,l_int32 type)1203 pixTophat(PIX *pixs,
1204 l_int32 hsize,
1205 l_int32 vsize,
1206 l_int32 type)
1207 {
1208 PIX *pixt, *pixd;
1209
1210 PROCNAME("pixTophat");
1211
1212 if (!pixs)
1213 return (PIX *)ERROR_PTR("seed pix not defined", procName, NULL);
1214 if (pixGetDepth(pixs) != 8)
1215 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
1216 if (hsize < 1 || vsize < 1)
1217 return (PIX *)ERROR_PTR("hsize or vsize < 1", procName, NULL);
1218 if ((hsize & 1) == 0 ) {
1219 L_WARNING("horiz sel size must be odd; increasing by 1\n", procName);
1220 hsize++;
1221 }
1222 if ((vsize & 1) == 0 ) {
1223 L_WARNING("vert sel size must be odd; increasing by 1\n", procName);
1224 vsize++;
1225 }
1226 if (type != L_TOPHAT_WHITE && type != L_TOPHAT_BLACK)
1227 return (PIX *)ERROR_PTR("type must be L_TOPHAT_BLACK or L_TOPHAT_WHITE",
1228 procName, NULL);
1229
1230 if (hsize == 1 && vsize == 1)
1231 return pixCreateTemplate(pixs);
1232
1233 switch (type)
1234 {
1235 case L_TOPHAT_WHITE:
1236 if ((pixt = pixOpenGray(pixs, hsize, vsize)) == NULL)
1237 return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
1238 pixd = pixSubtractGray(NULL, pixs, pixt);
1239 pixDestroy(&pixt);
1240 break;
1241 case L_TOPHAT_BLACK:
1242 if ((pixd = pixCloseGray(pixs, hsize, vsize)) == NULL)
1243 return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
1244 pixSubtractGray(pixd, pixd, pixs);
1245 break;
1246 default:
1247 return (PIX *)ERROR_PTR("invalid type", procName, NULL);
1248 }
1249
1250 return pixd;
1251 }
1252
1253
1254 /*!
1255 * \brief pixHDome()
1256 *
1257 * \param[in] pixs 8 bpp, filling mask
1258 * \param[in] height of seed below the filling maskhdome; must be >= 0
1259 * \param[in] connectivity 4 or 8
1260 * \return pixd 8 bpp, or NULL on error
1261 *
1262 * <pre>
1263 * Notes:
1264 * (1) It is more efficient to use a connectivity of 4 for the fill.
1265 * (2) This fills bumps to some level, and extracts the unfilled
1266 * part of the bump. To extract the troughs of basins, first
1267 * invert pixs and then apply pixHDome().
1268 * (3) It is useful to compare the HDome operation with the TopHat.
1269 * The latter extracts peaks or valleys that have a width
1270 * not exceeding the size of the structuring element used
1271 * in the opening or closing, rsp. The height of the peak is
1272 * irrelevant. By contrast, for the HDome, the gray seedfill
1273 * is used to extract all peaks that have a height not exceeding
1274 * a given value, regardless of their width!
1275 * (4) Slightly more precisely, suppose you set 'height' = 40.
1276 * Then all bumps in pixs with a height greater than or equal
1277 * to 40 become, in pixd, bumps with a max value of exactly 40.
1278 * All shorter bumps have a max value in pixd equal to the height
1279 * of the bump.
1280 * (5) The method: the filling mask, pixs, is the image whose peaks
1281 * are to be extracted. The height of a peak is the distance
1282 * between the top of the peak and the highest "leak" to the
1283 * outside -- think of a sombrero, where the leak occurs
1284 * at the highest point on the rim.
1285 * (a) Generate a seed, pixd, by subtracting some value, p, from
1286 * each pixel in the filling mask, pixs. The value p is
1287 * the 'height' input to this function.
1288 * (b) Fill in pixd starting with this seed, clipping by pixs,
1289 * in the way described in seedfillGrayLow(). The filling
1290 * stops before the peaks in pixs are filled.
1291 * For peaks that have a height > p, pixd is filled to
1292 * the level equal to the (top-of-the-peak - p).
1293 * For peaks of height < p, the peak is left unfilled
1294 * from its highest saddle point (the leak to the outside).
1295 * (c) Subtract the filled seed (pixd) from the filling mask (pixs).
1296 * Note that in this procedure, everything is done starting
1297 * with the filling mask, pixs.
1298 * (6) For segmentation, the resulting image, pixd, can be thresholded
1299 * and used as a seed for another filling operation.
1300 * </pre>
1301 */
1302 PIX *
pixHDome(PIX * pixs,l_int32 height,l_int32 connectivity)1303 pixHDome(PIX *pixs,
1304 l_int32 height,
1305 l_int32 connectivity)
1306 {
1307 PIX *pixsd, *pixd;
1308
1309 PROCNAME("pixHDome");
1310
1311 if (!pixs)
1312 return (PIX *)ERROR_PTR("src pix not defined", procName, NULL);
1313 if (pixGetDepth(pixs) != 8)
1314 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
1315 if (height < 0)
1316 return (PIX *)ERROR_PTR("height not >= 0", procName, NULL);
1317 if (height == 0)
1318 return pixCreateTemplate(pixs);
1319
1320 if ((pixsd = pixCopy(NULL, pixs)) == NULL)
1321 return (PIX *)ERROR_PTR("pixsd not made", procName, NULL);
1322 pixAddConstantGray(pixsd, -height);
1323 pixSeedfillGray(pixsd, pixs, connectivity);
1324 pixd = pixSubtractGray(NULL, pixs, pixsd);
1325 pixDestroy(&pixsd);
1326 return pixd;
1327 }
1328
1329
1330 /*!
1331 * \brief pixFastTophat()
1332 *
1333 * \param[in] pixs
1334 * \param[in] xsize width of max/min op, smoothing; any integer >= 1
1335 * \param[in] ysize height of max/min op, smoothing; any integer >= 1
1336 * \param[in] type L_TOPHAT_WHITE: image - min
1337 * L_TOPHAT_BLACK: max - image
1338 * \return pixd, or NULL on error
1339 *
1340 * <pre>
1341 * Notes:
1342 * (1) Don't be fooled. This is NOT a tophat. It is a tophat-like
1343 * operation, where the result is similar to what you'd get
1344 * if you used an erosion instead of an opening, or a dilation
1345 * instead of a closing.
1346 * (2) Instead of opening or closing at full resolution, it does
1347 * a fast downscale/minmax operation, then a quick small smoothing
1348 * at low res, a replicative expansion of the "background"
1349 * to full res, and finally a removal of the background level
1350 * from the input image. The smoothing step may not be important.
1351 * (3) It does not remove noise as well as a tophat, but it is
1352 * 5 to 10 times faster.
1353 * If you need the preciseness of the tophat, don't use this.
1354 * (4) The L_TOPHAT_WHITE flag emphasizes small bright regions,
1355 * whereas the L_TOPHAT_BLACK flag emphasizes small dark regions.
1356 * </pre>
1357 */
1358 PIX *
pixFastTophat(PIX * pixs,l_int32 xsize,l_int32 ysize,l_int32 type)1359 pixFastTophat(PIX *pixs,
1360 l_int32 xsize,
1361 l_int32 ysize,
1362 l_int32 type)
1363 {
1364 PIX *pix1, *pix2, *pix3, *pixd;
1365
1366 PROCNAME("pixFastTophat");
1367
1368 if (!pixs)
1369 return (PIX *)ERROR_PTR("seed pix not defined", procName, NULL);
1370 if (pixGetDepth(pixs) != 8)
1371 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
1372 if (xsize < 1 || ysize < 1)
1373 return (PIX *)ERROR_PTR("size < 1", procName, NULL);
1374 if (type != L_TOPHAT_WHITE && type != L_TOPHAT_BLACK)
1375 return (PIX *)ERROR_PTR("type must be L_TOPHAT_BLACK or L_TOPHAT_WHITE",
1376 procName, NULL);
1377
1378 if (xsize == 1 && ysize == 1)
1379 return pixCreateTemplate(pixs);
1380
1381 switch (type)
1382 {
1383 case L_TOPHAT_WHITE:
1384 if ((pix1 = pixScaleGrayMinMax(pixs, xsize, ysize, L_CHOOSE_MIN))
1385 == NULL)
1386 return (PIX *)ERROR_PTR("pix1 not made", procName, NULL);
1387 pix2 = pixBlockconv(pix1, 1, 1); /* small smoothing */
1388 pix3 = pixScaleBySampling(pix2, xsize, ysize);
1389 pixd = pixSubtractGray(NULL, pixs, pix3);
1390 pixDestroy(&pix3);
1391 break;
1392 case L_TOPHAT_BLACK:
1393 if ((pix1 = pixScaleGrayMinMax(pixs, xsize, ysize, L_CHOOSE_MAX))
1394 == NULL)
1395 return (PIX *)ERROR_PTR("pix1 not made", procName, NULL);
1396 pix2 = pixBlockconv(pix1, 1, 1); /* small smoothing */
1397 pixd = pixScaleBySampling(pix2, xsize, ysize);
1398 pixSubtractGray(pixd, pixd, pixs);
1399 break;
1400 default:
1401 return (PIX *)ERROR_PTR("invalid type", procName, NULL);
1402 }
1403
1404 pixDestroy(&pix1);
1405 pixDestroy(&pix2);
1406 return pixd;
1407 }
1408
1409
1410 /*!
1411 * \brief pixMorphGradient()
1412 *
1413 * \param[in] pixs
1414 * \param[in] hsize of Sel; must be odd; origin implicitly in center
1415 * \param[in] vsize ditto
1416 * \param[in] smoothing half-width of convolution smoothing filter.
1417 * The width is (2 * smoothing + 1, so 0 is no-op.
1418 * \return pixd, or NULL on error
1419 */
1420 PIX *
pixMorphGradient(PIX * pixs,l_int32 hsize,l_int32 vsize,l_int32 smoothing)1421 pixMorphGradient(PIX *pixs,
1422 l_int32 hsize,
1423 l_int32 vsize,
1424 l_int32 smoothing)
1425 {
1426 PIX *pixg, *pixd;
1427
1428 PROCNAME("pixMorphGradient");
1429
1430 if (!pixs)
1431 return (PIX *)ERROR_PTR("seed pix not defined", procName, NULL);
1432 if (pixGetDepth(pixs) != 8)
1433 return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
1434 if (hsize < 1 || vsize < 1)
1435 return (PIX *)ERROR_PTR("hsize or vsize < 1", procName, NULL);
1436 if ((hsize & 1) == 0 ) {
1437 L_WARNING("horiz sel size must be odd; increasing by 1\n", procName);
1438 hsize++;
1439 }
1440 if ((vsize & 1) == 0 ) {
1441 L_WARNING("vert sel size must be odd; increasing by 1\n", procName);
1442 vsize++;
1443 }
1444
1445 /* Optionally smooth first to remove noise.
1446 * If smoothing is 0, just get a copy */
1447 pixg = pixBlockconvGray(pixs, NULL, smoothing, smoothing);
1448
1449 /* This gives approximately the gradient of a transition */
1450 pixd = pixDilateGray(pixg, hsize, vsize);
1451 pixSubtractGray(pixd, pixd, pixg);
1452 pixDestroy(&pixg);
1453 return pixd;
1454 }
1455
1456
1457 /*-----------------------------------------------------------------*
1458 * Centroid of component *
1459 *-----------------------------------------------------------------*/
1460 /*!
1461 * \brief pixaCentroids()
1462 *
1463 * \param[in] pixa of components 1 or 8 bpp
1464 * \return pta of centroids relative to the UL corner of
1465 * each pix, or NULL on error
1466 *
1467 * <pre>
1468 * Notes:
1469 * (1) An error message is returned if any pix has something other
1470 * than 1 bpp or 8 bpp depth, and the centroid from that pix
1471 * is saved as (0, 0).
1472 * </pre>
1473 */
1474 PTA *
pixaCentroids(PIXA * pixa)1475 pixaCentroids(PIXA *pixa)
1476 {
1477 l_int32 i, n;
1478 l_int32 *centtab = NULL;
1479 l_int32 *sumtab = NULL;
1480 l_float32 x, y;
1481 PIX *pix;
1482 PTA *pta;
1483
1484 PROCNAME("pixaCentroids");
1485
1486 if (!pixa)
1487 return (PTA *)ERROR_PTR("pixa not defined", procName, NULL);
1488 if ((n = pixaGetCount(pixa)) == 0)
1489 return (PTA *)ERROR_PTR("no pix in pixa", procName, NULL);
1490
1491 if ((pta = ptaCreate(n)) == NULL)
1492 return (PTA *)ERROR_PTR("pta not defined", procName, NULL);
1493 centtab = makePixelCentroidTab8();
1494 sumtab = makePixelSumTab8();
1495
1496 for (i = 0; i < n; i++) {
1497 pix = pixaGetPix(pixa, i, L_CLONE);
1498 if (pixCentroid(pix, centtab, sumtab, &x, &y) == 1)
1499 L_ERROR("centroid failure for pix %d\n", procName, i);
1500 pixDestroy(&pix);
1501 ptaAddPt(pta, x, y);
1502 }
1503
1504 LEPT_FREE(centtab);
1505 LEPT_FREE(sumtab);
1506 return pta;
1507 }
1508
1509
1510 /*!
1511 * \brief pixCentroid()
1512 *
1513 * \param[in] pix 1 or 8 bpp
1514 * \param[in] centtab [optional] table for finding centroids; can be null
1515 * \param[in] sumtab [optional] table for finding pixel sums; can be null
1516 * \param[out] pxave, pyave coordinates of centroid, relative to
1517 * the UL corner of the pix
1518 * \return 0 if OK, 1 on error
1519 *
1520 * <pre>
1521 * Notes:
1522 * (1) Any table not passed in will be made internally and destroyed
1523 * after use.
1524 * </pre>
1525 */
1526 l_int32
pixCentroid(PIX * pix,l_int32 * centtab,l_int32 * sumtab,l_float32 * pxave,l_float32 * pyave)1527 pixCentroid(PIX *pix,
1528 l_int32 *centtab,
1529 l_int32 *sumtab,
1530 l_float32 *pxave,
1531 l_float32 *pyave)
1532 {
1533 l_int32 w, h, d, i, j, wpl, pixsum, rowsum, val;
1534 l_float32 xsum, ysum;
1535 l_uint32 *data, *line;
1536 l_uint32 word;
1537 l_uint8 byte;
1538 l_int32 *ctab, *stab;
1539
1540 PROCNAME("pixCentroid");
1541
1542 if (!pxave || !pyave)
1543 return ERROR_INT("&pxave and &pyave not defined", procName, 1);
1544 *pxave = *pyave = 0.0;
1545 if (!pix)
1546 return ERROR_INT("pix not defined", procName, 1);
1547 pixGetDimensions(pix, &w, &h, &d);
1548 if (d != 1 && d != 8)
1549 return ERROR_INT("pix not 1 or 8 bpp", procName, 1);
1550
1551 if (!centtab)
1552 ctab = makePixelCentroidTab8();
1553 else
1554 ctab = centtab;
1555 if (!sumtab)
1556 stab = makePixelSumTab8();
1557 else
1558 stab = sumtab;
1559
1560 data = pixGetData(pix);
1561 wpl = pixGetWpl(pix);
1562 xsum = ysum = 0.0;
1563 pixsum = 0;
1564 if (d == 1) {
1565 for (i = 0; i < h; i++) {
1566 /* The body of this loop computes the sum of the set
1567 * (1) bits on this row, weighted by their distance
1568 * from the left edge of pix, and accumulates that into
1569 * xsum; it accumulates their distance from the top
1570 * edge of pix into ysum, and their total count into
1571 * pixsum. It's equivalent to
1572 * for (j = 0; j < w; j++) {
1573 * if (GET_DATA_BIT(line, j)) {
1574 * xsum += j;
1575 * ysum += i;
1576 * pixsum++;
1577 * }
1578 * }
1579 */
1580 line = data + wpl * i;
1581 rowsum = 0;
1582 for (j = 0; j < wpl; j++) {
1583 word = line[j];
1584 if (word) {
1585 byte = word & 0xff;
1586 rowsum += stab[byte];
1587 xsum += ctab[byte] + (j * 32 + 24) * stab[byte];
1588 byte = (word >> 8) & 0xff;
1589 rowsum += stab[byte];
1590 xsum += ctab[byte] + (j * 32 + 16) * stab[byte];
1591 byte = (word >> 16) & 0xff;
1592 rowsum += stab[byte];
1593 xsum += ctab[byte] + (j * 32 + 8) * stab[byte];
1594 byte = (word >> 24) & 0xff;
1595 rowsum += stab[byte];
1596 xsum += ctab[byte] + j * 32 * stab[byte];
1597 }
1598 }
1599 pixsum += rowsum;
1600 ysum += rowsum * i;
1601 }
1602 if (pixsum == 0) {
1603 L_WARNING("no ON pixels in pix\n", procName);
1604 } else {
1605 *pxave = xsum / (l_float32)pixsum;
1606 *pyave = ysum / (l_float32)pixsum;
1607 }
1608 } else { /* d == 8 */
1609 for (i = 0; i < h; i++) {
1610 line = data + wpl * i;
1611 for (j = 0; j < w; j++) {
1612 val = GET_DATA_BYTE(line, j);
1613 xsum += val * j;
1614 ysum += val * i;
1615 pixsum += val;
1616 }
1617 }
1618 if (pixsum == 0) {
1619 L_WARNING("all pixels are 0\n", procName);
1620 } else {
1621 *pxave = xsum / (l_float32)pixsum;
1622 *pyave = ysum / (l_float32)pixsum;
1623 }
1624 }
1625
1626 if (!centtab) LEPT_FREE(ctab);
1627 if (!sumtab) LEPT_FREE(stab);
1628 return 0;
1629 }
1630