/*====================================================================*
- Copyright (C) 2001 Leptonica. All rights reserved.
- This software is distributed in the hope that it will be
- useful, but with NO WARRANTY OF ANY KIND.
- No author or distributor accepts responsibility to anyone for the
- consequences of using this software, or for whether it serves any
- particular purpose or works at all, unless he or she says so in
- writing. Everyone is granted permission to copy, modify and
- redistribute this source code, for commercial or non-commercial
- purposes, with the following restrictions: (1) the origin of this
- source code must not be misrepresented; (2) modified versions must
- be plainly marked as such; and (3) this notice may not be removed
- or altered from any source or modified source distribution.
*====================================================================*/
/*
* morphapp.c
*
* These are some useful and/or interesting composite
* image processing operations, of the type that are often
* useful in applications. Most are morphological in
* nature.
*
* Selective morph sequence operation on each component
* PIX *pixMorphSequenceByComponent()
* PIXA *pixaMorphSequenceByComponent()
*
* Selective morph sequence operation on each region
* PIX *pixMorphSequenceByRegion()
* PIXA *pixaMorphSequenceByRegion()
*
* Union and intersection of parallel composite operations
* PIX *pixUnionOfMorphOps()
* PIX *pixIntersectionOfMorphOps()
*
* Selective connected component filling
* PIX *pixSelectiveConnCompFill()
*
* Removal of matched patterns
* PIX *pixRemoveMatchedPattern()
*
* Display of matched patterns
* PIX *pixDisplayMatchedPattern()
*
* Iterative morphological seed filling (don't use for real work)
* PIX *pixSeedfillMorph()
*
* Granulometry on binary images
* NUMA *pixRunHistogramMorph()
*
* Composite operations on grayscale images
* PIX *pixTophat()
* PIX *pixHDome()
* PIX *pixFastTophat()
* PIX *pixMorphGradient()
*
* Centroids of PIXA
* PTA *pixaCentroids()
*/
#include <stdio.h>
#include "allheaders.h"
#define SWAP(x, y) {temp = (x); (x) = (y); (y) = temp;}
/*-----------------------------------------------------------------*
* Morph sequence operation on each component *
*-----------------------------------------------------------------*/
/*!
* pixMorphSequenceByComponent()
*
* Input: pixs (1 bpp)
* sequence (string specifying sequence)
* connectivity (4 or 8)
* minw (minimum width to consider; use 0 or 1 for any width)
* minh (minimum height to consider; use 0 or 1 for any height)
* &boxa (<optional> return boxa of c.c. in pixs)
* Return: pixd, or null on error
*
* Notes:
* (1) See pixMorphSequence() for composing operation sequences.
* (2) This operates separately on each c.c. in the input pix.
* (3) The dilation does NOT increase the c.c. size; it is clipped
* to the size of the original c.c. This is necessary to
* keep the c.c. independent after the operation.
* (4) You can specify that the width and/or height must equal
* or exceed a minimum size for the operation to take place.
* (5) Use NULL for boxa to avoid returning the boxa.
*/
PIX *
pixMorphSequenceByComponent(PIX *pixs,
const char *sequence,
l_int32 connectivity,
l_int32 minw,
l_int32 minh,
BOXA **pboxa)
{
l_int32 n, i, x, y, w, h;
BOXA *boxa;
PIX *pix, *pixd;
PIXA *pixas, *pixad;
PROCNAME("pixMorphSequenceByComponent");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (!sequence)
return (PIX *)ERROR_PTR("sequence not defined", procName, NULL);
if (minw <= 0) minw = 1;
if (minh <= 0) minh = 1;
/* Get the c.c. */
if ((boxa = pixConnComp(pixs, &pixas, connectivity)) == NULL)
return (PIX *)ERROR_PTR("boxa not made", procName, NULL);
/* Operate on each c.c. independently */
pixad = pixaMorphSequenceByComponent(pixas, sequence, minw, minh);
pixaDestroy(&pixas);
boxaDestroy(&boxa);
if (!pixad)
return (PIX *)ERROR_PTR("pixad not made", procName, NULL);
/* Display the result out into pixd */
pixd = pixCreateTemplate(pixs);
n = pixaGetCount(pixad);
for (i = 0; i < n; i++) {
pixaGetBoxGeometry(pixad, i, &x, &y, &w, &h);
pix = pixaGetPix(pixad, i, L_CLONE);
pixRasterop(pixd, x, y, w, h, PIX_PAINT, pix, 0, 0);
pixDestroy(&pix);
}
if (pboxa)
*pboxa = pixaGetBoxa(pixad, L_CLONE);
pixaDestroy(&pixad);
return pixd;
}
/*!
* pixaMorphSequenceByComponent()
*
* Input: pixas (of 1 bpp pix)
* sequence (string specifying sequence)
* minw (minimum width to consider; use 0 or 1 for any width)
* minh (minimum height to consider; use 0 or 1 for any height)
* Return: pixad, or null on error
*
* Notes:
* (1) See pixMorphSequence() for composing operation sequences.
* (2) This operates separately on each c.c. in the input pixa.
* (3) You can specify that the width and/or height must equal
* or exceed a minimum size for the operation to take place.
* (4) The input pixa should have a boxa giving the locations
* of the pix components.
*/
PIXA *
pixaMorphSequenceByComponent(PIXA *pixas,
const char *sequence,
l_int32 minw,
l_int32 minh)
{
l_int32 n, i, w, h, d;
BOX *box;
PIX *pixt1, *pixt2;
PIXA *pixad;
PROCNAME("pixaMorphSequenceByComponent");
if (!pixas)
return (PIXA *)ERROR_PTR("pixas not defined", procName, NULL);
if ((n = pixaGetCount(pixas)) == 0)
return (PIXA *)ERROR_PTR("no pix in pixas", procName, NULL);
if (n != pixaGetBoxaCount(pixas))
L_WARNING("boxa size != n", procName);
pixaGetPixDimensions(pixas, 0, NULL, NULL, &d);
if (d != 1)
return (PIXA *)ERROR_PTR("depth not 1 bpp", procName, NULL);
if (!sequence)
return (PIXA *)ERROR_PTR("sequence not defined", procName, NULL);
if (minw <= 0) minw = 1;
if (minh <= 0) minh = 1;
if ((pixad = pixaCreate(n)) == NULL)
return (PIXA *)ERROR_PTR("pixad not made", procName, NULL);
for (i = 0; i < n; i++) {
pixaGetPixDimensions(pixas, i, &w, &h, NULL);
if (w >= minw && h >= minh) {
if ((pixt1 = pixaGetPix(pixas, i, L_CLONE)) == NULL)
return (PIXA *)ERROR_PTR("pixt1 not found", procName, NULL);
if ((pixt2 = pixMorphCompSequence(pixt1, sequence, 0)) == NULL)
return (PIXA *)ERROR_PTR("pixt2 not made", procName, NULL);
pixaAddPix(pixad, pixt2, L_INSERT);
box = pixaGetBox(pixas, i, L_COPY);
pixaAddBox(pixad, box, L_INSERT);
pixDestroy(&pixt1);
}
}
return pixad;
}
/*-----------------------------------------------------------------*
* Morph sequence operation on each region *
*-----------------------------------------------------------------*/
/*!
* pixMorphSequenceByRegion()
*
* Input: pixs (1 bpp)
* pixm (mask specifying regions)
* sequence (string specifying sequence)
* connectivity (4 or 8, used on mask)
* minw (minimum width to consider; use 0 or 1 for any width)
* minh (minimum height to consider; use 0 or 1 for any height)
* &boxa (<optional> return boxa of c.c. in pixm)
* Return: pixd, or null on error
*
* Notes:
* (1) See pixMorphCompSequence() for composing operation sequences.
* (2) This operates separately on the region in pixs corresponding
* to each c.c. in the mask pixm. It differs from
* pixMorphSequenceByComponent() in that the latter does not have
* a pixm (mask), but instead operates independently on each
* component in pixs.
* (3) Dilation will NOT increase the region size; the result
* is clipped to the size of the mask region. This is necessary
* to make regions independent after the operation.
* (4) You can specify that the width and/or height of a region must
* equal or exceed a minimum size for the operation to take place.
* (5) Use NULL for @pboxa to avoid returning the boxa.
*/
PIX *
pixMorphSequenceByRegion(PIX *pixs,
PIX *pixm,
const char *sequence,
l_int32 connectivity,
l_int32 minw,
l_int32 minh,
BOXA **pboxa)
{
l_int32 n, i, x, y, w, h;
BOXA *boxa;
PIX *pix, *pixd;
PIXA *pixam, *pixad;
PROCNAME("pixMorphSequenceByRegion");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (!pixm)
return (PIX *)ERROR_PTR("pixm not defined", procName, NULL);
if (pixGetDepth(pixs) != 1 || pixGetDepth(pixm) != 1)
return (PIX *)ERROR_PTR("pixs and pixm not both 1 bpp", procName, NULL);
if (!sequence)
return (PIX *)ERROR_PTR("sequence not defined", procName, NULL);
if (minw <= 0) minw = 1;
if (minh <= 0) minh = 1;
/* Get the c.c. of the mask */
if ((boxa = pixConnComp(pixm, &pixam, connectivity)) == NULL)
return (PIX *)ERROR_PTR("boxa not made", procName, NULL);
/* Operate on each region in pixs independently */
pixad = pixaMorphSequenceByRegion(pixs, pixam, sequence, minw, minh);
pixaDestroy(&pixam);
boxaDestroy(&boxa);
if (!pixad)
return (PIX *)ERROR_PTR("pixad not made", procName, NULL);
/* Display the result out into pixd */
pixd = pixCreateTemplate(pixs);
n = pixaGetCount(pixad);
for (i = 0; i < n; i++) {
pixaGetBoxGeometry(pixad, i, &x, &y, &w, &h);
pix = pixaGetPix(pixad, i, L_CLONE);
pixRasterop(pixd, x, y, w, h, PIX_PAINT, pix, 0, 0);
pixDestroy(&pix);
}
if (pboxa)
*pboxa = pixaGetBoxa(pixad, L_CLONE);
pixaDestroy(&pixad);
return pixd;
}
/*!
* pixaMorphSequenceByRegion()
*
* Input: pixs (1 bpp)
* pixam (of 1 bpp mask elements)
* sequence (string specifying sequence)
* minw (minimum width to consider; use 0 or 1 for any width)
* minh (minimum height to consider; use 0 or 1 for any height)
* Return: pixad, or null on error
*
* Notes:
* (1) See pixMorphSequence() for composing operation sequences.
* (2) This operates separately on each region in the input pixs
* defined by the components in pixam.
* (3) You can specify that the width and/or height of a mask
* component must equal or exceed a minimum size for the
* operation to take place.
* (4) The input pixam should have a boxa giving the locations
* of the regions in pixs.
*/
PIXA *
pixaMorphSequenceByRegion(PIX *pixs,
PIXA *pixam,
const char *sequence,
l_int32 minw,
l_int32 minh)
{
l_int32 n, i, w, h, d;
BOX *box;
PIX *pixt1, *pixt2, *pixt3;
PIXA *pixad;
PROCNAME("pixaMorphSequenceByRegion");
if (!pixs)
return (PIXA *)ERROR_PTR("pixs not defined", procName, NULL);
if (pixGetDepth(pixs) != 1)
return (PIXA *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
if (!pixam)
return (PIXA *)ERROR_PTR("pixam not defined", procName, NULL);
pixaGetPixDimensions(pixam, 0, NULL, NULL, &d);
if (d != 1)
return (PIXA *)ERROR_PTR("mask depth not 1 bpp", procName, NULL);
if ((n = pixaGetCount(pixam)) == 0)
return (PIXA *)ERROR_PTR("no regions specified", procName, NULL);
if (n != pixaGetBoxaCount(pixam))
L_WARNING("boxa size != n", procName);
if (!sequence)
return (PIXA *)ERROR_PTR("sequence not defined", procName, NULL);
if (minw <= 0) minw = 1;
if (minh <= 0) minh = 1;
if ((pixad = pixaCreate(n)) == NULL)
return (PIXA *)ERROR_PTR("pixad not made", procName, NULL);
/* Use the rectangle to remove the appropriate part of pixs;
* then AND with the mask component to get the actual fg
* of pixs that is under the mask component. */
for (i = 0; i < n; i++) {
pixaGetPixDimensions(pixam, i, &w, &h, NULL);
if (w >= minw && h >= minh) {
if ((pixt1 = pixaGetPix(pixam, i, L_CLONE)) == NULL)
return (PIXA *)ERROR_PTR("pixt1 not found", procName, NULL);
box = pixaGetBox(pixam, i, L_COPY);
pixt2 = pixClipRectangle(pixs, box, NULL);
pixAnd(pixt2, pixt2, pixt1);
if ((pixt3 = pixMorphCompSequence(pixt2, sequence, 0)) == NULL)
return (PIXA *)ERROR_PTR("pixt3 not made", procName, NULL);
pixaAddPix(pixad, pixt3, L_INSERT);
pixaAddBox(pixad, box, L_INSERT);
pixDestroy(&pixt1);
pixDestroy(&pixt2);
}
}
return pixad;
}
/*-----------------------------------------------------------------*
* Union and intersection of parallel composite operations *
*-----------------------------------------------------------------*/
/*!
* pixUnionOfMorphOps()
*
* Input: pixs (binary)
* sela
* type (L_MORPH_DILATE, etc.)
* Return: pixd (union of the specified morphological operation
* on pixs for each Sel in the Sela), or null on error
*/
PIX *
pixUnionOfMorphOps(PIX *pixs,
SELA *sela,
l_int32 type)
{
l_int32 n, i;
PIX *pixt, *pixd;
SEL *sel;
PROCNAME("pixUnionOfMorphOps");
if (!pixs || pixGetDepth(pixs) != 1)
return (PIX *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
if (!sela)
return (PIX *)ERROR_PTR("sela not defined", procName, NULL);
n = selaGetCount(sela);
if (n == 0)
return (PIX *)ERROR_PTR("no sels in sela", procName, NULL);
if (type != L_MORPH_DILATE && type != L_MORPH_ERODE &&
type != L_MORPH_OPEN && type != L_MORPH_CLOSE &&
type != L_MORPH_HMT)
return (PIX *)ERROR_PTR("invalid type", procName, NULL);
pixd = pixCreateTemplate(pixs);
for (i = 0; i < n; i++) {
sel = selaGetSel(sela, i);
if (type == L_MORPH_DILATE)
pixt = pixDilate(NULL, pixs, sel);
else if (type == L_MORPH_ERODE)
pixt = pixErode(NULL, pixs, sel);
else if (type == L_MORPH_OPEN)
pixt = pixOpen(NULL, pixs, sel);
else if (type == L_MORPH_CLOSE)
pixt = pixClose(NULL, pixs, sel);
else /* type == L_MORPH_HMT */
pixt = pixHMT(NULL, pixs, sel);
pixOr(pixd, pixd, pixt);
pixDestroy(&pixt);
}
return pixd;
}
/*!
* pixIntersectionOfMorphOps()
*
* Input: pixs (binary)
* sela
* type (L_MORPH_DILATE, etc.)
* Return: pixd (intersection of the specified morphological operation
* on pixs for each Sel in the Sela), or null on error
*/
PIX *
pixIntersectionOfMorphOps(PIX *pixs,
SELA *sela,
l_int32 type)
{
l_int32 n, i;
PIX *pixt, *pixd;
SEL *sel;
PROCNAME("pixIntersectionOfMorphOps");
if (!pixs || pixGetDepth(pixs) != 1)
return (PIX *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
if (!sela)
return (PIX *)ERROR_PTR("sela not defined", procName, NULL);
n = selaGetCount(sela);
if (n == 0)
return (PIX *)ERROR_PTR("no sels in sela", procName, NULL);
if (type != L_MORPH_DILATE && type != L_MORPH_ERODE &&
type != L_MORPH_OPEN && type != L_MORPH_CLOSE &&
type != L_MORPH_HMT)
return (PIX *)ERROR_PTR("invalid type", procName, NULL);
pixd = pixCreateTemplate(pixs);
pixSetAll(pixd);
for (i = 0; i < n; i++) {
sel = selaGetSel(sela, i);
if (type == L_MORPH_DILATE)
pixt = pixDilate(NULL, pixs, sel);
else if (type == L_MORPH_ERODE)
pixt = pixErode(NULL, pixs, sel);
else if (type == L_MORPH_OPEN)
pixt = pixOpen(NULL, pixs, sel);
else if (type == L_MORPH_CLOSE)
pixt = pixClose(NULL, pixs, sel);
else /* type == L_MORPH_HMT */
pixt = pixHMT(NULL, pixs, sel);
pixAnd(pixd, pixd, pixt);
pixDestroy(&pixt);
}
return pixd;
}
/*-----------------------------------------------------------------*
* Selective connected component filling *
*-----------------------------------------------------------------*/
/*!
* pixSelectiveConnCompFill()
*
* Input: pixs (binary)
* connectivity (4 or 8)
* minw (minimum width to consider; use 0 or 1 for any width)
* minh (minimum height to consider; use 0 or 1 for any height)
* Return: pix (with holes filled in selected c.c.), or null on error
*/
PIX *
pixSelectiveConnCompFill(PIX *pixs,
l_int32 connectivity,
l_int32 minw,
l_int32 minh)
{
l_int32 n, i, x, y, w, h;
BOXA *boxa;
PIX *pixt1, *pixt2, *pixd;
PIXA *pixa;
PROCNAME("pixSelectiveConnCompFill");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (pixGetDepth(pixs) != 1)
return (PIX *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
if (minw <= 0) minw = 1;
if (minh <= 0) minh = 1;
if ((pixd = pixCopy(NULL, pixs)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
if ((boxa = pixConnComp(pixs, &pixa, connectivity)) == NULL)
return (PIX *)ERROR_PTR("boxa not made", procName, NULL);
n = boxaGetCount(boxa);
for (i = 0; i < n; i++) {
boxaGetBoxGeometry(boxa, i, &x, &y, &w, &h);
if (w >= minw && h >= minh) {
if ((pixt1 = pixaGetPix(pixa, i, L_CLONE)) == NULL)
return (PIX *)ERROR_PTR("pixt1 not found", procName, NULL);
if ((pixt2 = pixHolesByFilling(pixt1, 12 - connectivity)) == NULL)
return (PIX *)ERROR_PTR("pixt2 not made", procName, NULL);
pixRasterop(pixd, x, y, w, h, PIX_PAINT, pixt2, 0, 0);
pixDestroy(&pixt1);
pixDestroy(&pixt2);
}
}
pixaDestroy(&pixa);
boxaDestroy(&boxa);
return pixd;
}
/*-----------------------------------------------------------------*
* Removal of matched patterns *
*-----------------------------------------------------------------*/
/*!
* pixRemoveMatchedPattern()
*
* Input: pixs (input image, 1 bpp)
* pixp (pattern to be removed from image, 1 bpp)
* pixe (image after erosion by Sel that approximates pixp, 1 bpp)
* x0, y0 (center of Sel)
* dsize (number of pixels on each side by which pixp is
* dilated before being subtracted from pixs;
* valid values are {0, 1, 2, 3, 4})
* Return: 0 if OK, 1 on error
*
* Notes:
* (1) This is in-place.
* (2) You can use various functions in selgen to create a Sel
* that is used to generate pixe from pixs.
* (3) This function is applied after pixe has been computed.
* It finds the centroid of each c.c., and subtracts
* (the appropriately dilated version of) pixp, with the center
* of the Sel used to align pixp with pixs.
*/
l_int32
pixRemoveMatchedPattern(PIX *pixs,
PIX *pixp,
PIX *pixe,
l_int32 x0,
l_int32 y0,
l_int32 dsize)
{
l_int32 i, nc, x, y, w, h, xb, yb;
BOXA *boxa;
PIX *pixt1, *pixt2;
PIXA *pixa;
PTA *pta;
SEL *sel;
PROCNAME("pixRemoveMatchedPattern");
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
if (!pixp)
return ERROR_INT("pixp not defined", procName, 1);
if (!pixe)
return ERROR_INT("pixe not defined", procName, 1);
if (pixGetDepth(pixs) != 1 || pixGetDepth(pixp) != 1 ||
pixGetDepth(pixe) != 1)
return ERROR_INT("all input pix not 1 bpp", procName, 1);
if (dsize < 0 || dsize > 4)
return ERROR_INT("dsize not in {0,1,2,3,4}", procName, 1);
/* Find the connected components and their centroids */
boxa = pixConnComp(pixe, &pixa, 8);
if ((nc = boxaGetCount(boxa)) == 0) {
L_WARNING("no matched patterns", procName);
boxaDestroy(&boxa);
pixaDestroy(&pixa);
return 0;
}
pta = pixaCentroids(pixa);
/* Optionally dilate the pattern, first adding a border that
* is large enough to accommodate the dilated pixels */
sel = NULL;
if (dsize > 0) {
sel = selCreateBrick(2 * dsize + 1, 2 * dsize + 1, dsize, dsize,
SEL_HIT);
pixt1 = pixAddBorder(pixp, dsize, 0);
pixt2 = pixDilate(NULL, pixt1, sel);
selDestroy(&sel);
pixDestroy(&pixt1);
}
else
pixt2 = pixClone(pixp);
/* Subtract out each dilated pattern. The centroid of each
* component is located at:
* (box->x + x, box->y + y)
* and the 'center' of the pattern used in making pixe is located at
* (x0 + dsize, (y0 + dsize)
* relative to the UL corner of the pattern. The center of the
* pattern is placed at the center of the component. */
w = pixGetWidth(pixt2);
h = pixGetHeight(pixt2);
for (i = 0; i < nc; i++) {
ptaGetIPt(pta, i, &x, &y);
boxaGetBoxGeometry(boxa, i, &xb, &yb, NULL, NULL);
pixRasterop(pixs, xb + x - x0 - dsize, yb + y - y0 - dsize,
w, h, PIX_DST & PIX_NOT(PIX_SRC), pixt2, 0, 0);
}
boxaDestroy(&boxa);
pixaDestroy(&pixa);
ptaDestroy(&pta);
pixDestroy(&pixt2);
return 0;
}
/*-----------------------------------------------------------------*
* Display of matched patterns *
*-----------------------------------------------------------------*/
/*!
* pixDisplayMatchedPattern()
*
* Input: pixs (input image, 1 bpp)
* pixp (pattern to be removed from image, 1 bpp)
* pixe (image after erosion by Sel that approximates pixp, 1 bpp)
* x0, y0 (center of Sel)
* color (to paint the matched patterns; 0xrrggbb00)
* scale (reduction factor for output pixd)
* nlevels (if scale < 1.0, threshold to this number of levels)
* Return: pixd (8 bpp, colormapped), or null on error
*
* Notes:
* (1) A 4 bpp colormapped image is generated.
* (2) If scale <= 1.0, do scale to gray for the output, and threshold
* to nlevels of gray.
* (3) You can use various functions in selgen to create a Sel
* that will generate pixe from pixs.
* (4) This function is applied after pixe has been computed.
* It finds the centroid of each c.c., and colors the output
* pixels using pixp (appropriately aligned) as a stencil.
* Alignment is done using the origin of the Sel and the
* centroid of the eroded image to place the stencil pixp.
*/
PIX *
pixDisplayMatchedPattern(PIX *pixs,
PIX *pixp,
PIX *pixe,
l_int32 x0,
l_int32 y0,
l_uint32 color,
l_float32 scale,
l_int32 nlevels)
{
l_int32 i, nc, xb, yb, x, y, xi, yi, rval, gval, bval;
BOXA *boxa;
PIX *pixd, *pixt, *pixps;
PIXA *pixa;
PTA *pta;
PIXCMAP *cmap;
PROCNAME("pixDisplayMatchedPattern");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (!pixp)
return (PIX *)ERROR_PTR("pixp not defined", procName, NULL);
if (!pixe)
return (PIX *)ERROR_PTR("pixe not defined", procName, NULL);
if (pixGetDepth(pixs) != 1 || pixGetDepth(pixp) != 1 ||
pixGetDepth(pixe) != 1)
return (PIX *)ERROR_PTR("all input pix not 1 bpp", procName, NULL);
if (scale > 1.0 || scale <= 0.0) {
L_WARNING("scale > 1.0 or < 0.0; setting to 1.0", procName);
scale = 1.0;
}
/* Find the connected components and their centroids */
boxa = pixConnComp(pixe, &pixa, 8);
if ((nc = boxaGetCount(boxa)) == 0) {
L_WARNING("no matched patterns", procName);
boxaDestroy(&boxa);
pixaDestroy(&pixa);
return 0;
}
pta = pixaCentroids(pixa);
rval = GET_DATA_BYTE(&color, COLOR_RED);
gval = GET_DATA_BYTE(&color, COLOR_GREEN);
bval = GET_DATA_BYTE(&color, COLOR_BLUE);
if (scale == 1.0) { /* output 4 bpp at full resolution */
pixd = pixConvert1To4(NULL, pixs, 0, 1);
cmap = pixcmapCreate(4);
pixcmapAddColor(cmap, 255, 255, 255);
pixcmapAddColor(cmap, 0, 0, 0);
pixSetColormap(pixd, cmap);
/* Paint through pixp for each match location. The centroid of each
* component in pixe is located at:
* (box->x + x, box->y + y)
* and the 'center' of the pattern used in making pixe is located at
* (x0, y0)
* relative to the UL corner of the pattern. The center of the
* pattern is placed at the center of the component. */
for (i = 0; i < nc; i++) {
ptaGetIPt(pta, i, &x, &y);
boxaGetBoxGeometry(boxa, i, &xb, &yb, NULL, NULL);
pixSetMaskedCmap(pixd, pixp, xb + x - x0, yb + y - y0,
rval, gval, bval);
}
}
else { /* output 4 bpp downscaled */
pixt = pixScaleToGray(pixs, scale);
pixd = pixThresholdTo4bpp(pixt, nlevels, 1);
pixps = pixScaleBySampling(pixp, scale, scale);
for (i = 0; i < nc; i++) {
ptaGetIPt(pta, i, &x, &y);
boxaGetBoxGeometry(boxa, i, &xb, &yb, NULL, NULL);
xi = (l_int32)(scale * (xb + x - x0));
yi = (l_int32)(scale * (yb + y - y0));
pixSetMaskedCmap(pixd, pixps, xi, yi, rval, gval, bval);
}
pixDestroy(&pixt);
pixDestroy(&pixps);
}
boxaDestroy(&boxa);
pixaDestroy(&pixa);
ptaDestroy(&pta);
return pixd;
}
/*-----------------------------------------------------------------*
* Iterative morphological seed filling *
*-----------------------------------------------------------------*/
/*!
* pixSeedfillMorph()
*
* Input: pixs (seed)
* pixm (mask)
* connectivity (4 or 8)
* Return: pix where seed has been grown to completion
* into the mask, or null on error
*
* Notes:
* (1) This is in general a very inefficient method for filling
* from a seed into a mask. I've included it here for
* pedagogical reasons, but it should NEVER be used if
* efficiency is any consideration -- use pixSeedfillBinary()!
* (2) We use a 3x3 brick SEL for 8-cc filling and a 3x3 plus SEL for 4-cc.
*/
PIX *
pixSeedfillMorph(PIX *pixs,
PIX *pixm,
l_int32 connectivity)
{
l_int32 same, iter;
PIX *pixt1, *pixd, *temp;
SEL *sel_3;
PROCNAME("pixSeedfillMorph");
if (!pixs)
return (PIX *)ERROR_PTR("seed pix not defined", procName, NULL);
if (!pixm)
return (PIX *)ERROR_PTR("mask pix not defined", procName, NULL);
if (connectivity != 4 && connectivity != 8)
return (PIX *)ERROR_PTR("connectivity not in {4,8}", procName, NULL);
if (pixSizesEqual(pixs, pixm) == 0)
return (PIX *)ERROR_PTR("pix sizes unequal", procName, NULL);
if (pixGetDepth(pixs) != 1)
return (PIX *)ERROR_PTR("pix not binary", procName, NULL);
if ((sel_3 = selCreateBrick(3, 3, 1, 1, 1)) == NULL)
return (PIX *)ERROR_PTR("sel_3 not made", procName, NULL);
if (connectivity == 4) { /* remove corner hits to make a '+' */
selSetElement(sel_3, 0, 0, SEL_DONT_CARE);
selSetElement(sel_3, 2, 2, SEL_DONT_CARE);
selSetElement(sel_3, 2, 0, SEL_DONT_CARE);
selSetElement(sel_3, 0, 2, SEL_DONT_CARE);
}
if ((pixt1 = pixCopy(NULL, pixs)) == NULL)
return (PIX *)ERROR_PTR("pixt1 not made", procName, NULL);
if ((pixd = pixCreateTemplate(pixs)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
iter = 0;
while (1) {
iter++;
pixDilate(pixd, pixt1, sel_3);
pixAnd(pixd, pixd, pixm);
pixEqual(pixd, pixt1, &same);
if (same)
break;
else
SWAP(pixt1, pixd);
}
fprintf(stderr, " Num iters in binary reconstruction = %d\n", iter);
pixDestroy(&pixt1);
selDestroy(&sel_3);
return pixd;
}
/*-----------------------------------------------------------------*
* Granulometry on binary images *
*-----------------------------------------------------------------*/
/*!
* pixRunHistogramMorph()
*
* Input: pixs
* runtype (L_RUN_OFF, L_RUN_ON)
* direction (L_HORIZ, L_VERT)
* maxsize (size of largest runlength counted)
* Return: numa of run-lengths
*/
NUMA *
pixRunHistogramMorph(PIX *pixs,
l_int32 runtype,
l_int32 direction,
l_int32 maxsize)
{
l_int32 count, i;
l_float32 val;
NUMA *na, *nah;
PIX *pixt1, *pixt2, *pixt3;
SEL *sel_2a;
PROCNAME("pixRunHistogramMorph");
if (!pixs)
return (NUMA *)ERROR_PTR("seed pix not defined", procName, NULL);
if (runtype != L_RUN_OFF && runtype != L_RUN_ON)
return (NUMA *)ERROR_PTR("invalid run type", procName, NULL);
if (direction != L_HORIZ && direction != L_VERT)
return (NUMA *)ERROR_PTR("direction not in {L_HORIZ, L_VERT}",
procName, NULL);
if (pixGetDepth(pixs) != 1)
return (NUMA *)ERROR_PTR("pixs must be binary", procName, NULL);
if ((na = numaCreate(0)) == NULL)
return (NUMA *)ERROR_PTR("na not made", procName, NULL);
if (direction == L_HORIZ)
sel_2a = selCreateBrick(1, 2, 0, 0, 1);
else /* direction == L_VERT */
sel_2a = selCreateBrick(2, 1, 0, 0, 1);
if (!sel_2a)
return (NUMA *)ERROR_PTR("sel_2a not made", procName, NULL);
if (runtype == L_RUN_OFF) {
if ((pixt1 = pixCopy(NULL, pixs)) == NULL)
return (NUMA *)ERROR_PTR("pix1 not made", procName, NULL);
pixInvert(pixt1, pixt1);
}
else /* runtype == L_RUN_ON */
pixt1 = pixClone(pixs);
if ((pixt2 = pixCreateTemplate(pixs)) == NULL)
return (NUMA *)ERROR_PTR("pix2 not made", procName, NULL);
if ((pixt3 = pixCreateTemplate(pixs)) == NULL)
return (NUMA *)ERROR_PTR("pix3 not made", procName, NULL);
/* Get pixel counts at different stages of erosion */
pixCountPixels(pixt1, &count, NULL);
numaAddNumber(na, count);
pixErode(pixt2, pixt1, sel_2a);
pixCountPixels(pixt2, &count, NULL);
numaAddNumber(na, count);
for (i = 0; i < maxsize / 2; i++) {
pixErode(pixt3, pixt2, sel_2a);
pixCountPixels(pixt3, &count, NULL);
numaAddNumber(na, count);
pixErode(pixt2, pixt3, sel_2a);
pixCountPixels(pixt2, &count, NULL);
numaAddNumber(na, count);
}
/* Compute length histogram */
if ((nah = numaCreate(na->n)) == NULL)
return (NUMA *)ERROR_PTR("nah not made", procName, NULL);
numaAddNumber(nah, 0); /* number at length 0 */
for (i = 1; i < na->n - 1; i++) {
val = na->array[i+1] - 2 * na->array[i] + na->array[i-1];
numaAddNumber(nah, val);
}
pixDestroy(&pixt1);
pixDestroy(&pixt2);
pixDestroy(&pixt3);
selDestroy(&sel_2a);
numaDestroy(&na);
return nah;
}
/*-----------------------------------------------------------------*
* Composite operations on grayscale images *
*-----------------------------------------------------------------*/
/*!
* pixTophat()
*
* Input: pixs
* hsize (of Sel; must be odd; origin implicitly in center)
* vsize (ditto)
* type (L_TOPHAT_WHITE: image - opening
* L_TOPHAT_BLACK: closing - image)
* Return: pixd, or null on error
*
* Notes:
* (1) Sel is a brick with all elements being hits
* (2) If hsize = vsize = 1, returns an image with all 0 data.
* (3) The L_TOPHAT_WHITE flag emphasizes small bright regions,
* whereas the L_TOPHAT_BLACK flag emphasizes small dark regions.
* The L_TOPHAT_WHITE tophat can be accomplished by doing a
* L_TOPHAT_BLACK tophat on the inverse, or v.v.
*/
PIX *
pixTophat(PIX *pixs,
l_int32 hsize,
l_int32 vsize,
l_int32 type)
{
PIX *pixt, *pixd;
PROCNAME("pixTophat");
if (!pixs)
return (PIX *)ERROR_PTR("seed pix not defined", procName, NULL);
if (pixGetDepth(pixs) != 8)
return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
if (hsize < 1 || vsize < 1)
return (PIX *)ERROR_PTR("hsize or vsize < 1", procName, NULL);
if ((hsize & 1) == 0 ) {
L_WARNING("horiz sel size must be odd; increasing by 1", procName);
hsize++;
}
if ((vsize & 1) == 0 ) {
L_WARNING("vert sel size must be odd; increasing by 1", procName);
vsize++;
}
if (type != L_TOPHAT_WHITE && type != L_TOPHAT_BLACK)
return (PIX *)ERROR_PTR("type must be L_TOPHAT_BLACK or L_TOPHAT_WHITE",
procName, NULL);
if (hsize == 1 && vsize == 1)
return pixCreateTemplate(pixs);
switch (type)
{
case L_TOPHAT_WHITE:
if ((pixt = pixOpenGray(pixs, hsize, vsize)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
pixd = pixSubtractGray(NULL, pixs, pixt);
pixDestroy(&pixt);
break;
case L_TOPHAT_BLACK:
if ((pixd = pixCloseGray(pixs, hsize, vsize)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
pixSubtractGray(pixd, pixd, pixs);
break;
default:
return (PIX *)ERROR_PTR("invalid type", procName, NULL);
}
return pixd;
}
/*!
* pixHDome()
*
* Input: pixs (8 bpp, filling mask)
* height (of seed below the filling maskhdome; must be >= 0)
* connectivity (4 or 8)
* Return: pixd (8 bpp), or null on error
*
* Notes:
* (1) It is more efficient to use a connectivity of 4 for the fill.
* (2) This fills bumps to some level, and extracts the unfilled
* part of the bump. To extract the troughs of basins, first
* invert pixs and then apply pixHDome().
* (3) It is useful to compare the HDome operation with the TopHat.
* The latter extracts peaks or valleys that have a width
* not exceeding the size of the structuring element used
* in the opening or closing, rsp. The height of the peak is
* irrelevant. By contrast, for the HDome, the gray seedfill
* is used to extract all peaks that have a height not exceeding
* a given value, regardless of their width!
* (4) Slightly more precisely, suppose you set 'height' = 40.
* Then all bumps in pixs with a height greater than or equal
* to 40 become, in pixd, bumps with a max value of exactly 40.
* All shorter bumps have a max value in pixd equal to the height
* of the bump.
* (5) The method: the filling mask, pixs, is the image whose peaks
* are to be extracted. The height of a peak is the distance
* between the top of the peak and the highest "leak" to the
* outside -- think of a sombrero, where the leak occurs
* at the highest point on the rim.
* (a) Generate a seed, pixd, by subtracting some value, p, from
* each pixel in the filling mask, pixs. The value p is
* the 'height' input to this function.
* (b) Fill in pixd starting with this seed, clipping by pixs,
* in the way described in seedfillGrayLow(). The filling
* stops before the peaks in pixs are filled.
* For peaks that have a height > p, pixd is filled to
* the level equal to the (top-of-the-peak - p).
* For peaks of height < p, the peak is left unfilled
* from its highest saddle point (the leak to the outside).
* (c) Subtract the filled seed (pixd) from the filling mask (pixs).
* Note that in this procedure, everything is done starting
* with the filling mask, pixs.
* (6) For segmentation, the resulting image, pixd, can be thresholded
* and used as a seed for another filling operation.
*/
PIX *
pixHDome(PIX *pixs,
l_int32 height,
l_int32 connectivity)
{
PIX *pixsd, *pixd;
PROCNAME("pixHDome");
if (!pixs)
return (PIX *)ERROR_PTR("src pix not defined", procName, NULL);
if (pixGetDepth(pixs) != 8)
return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
if (height < 0)
return (PIX *)ERROR_PTR("height not >= 0", procName, NULL);
if (height == 0)
return pixCreateTemplate(pixs);
if ((pixsd = pixCopy(NULL, pixs)) == NULL)
return (PIX *)ERROR_PTR("pixsd not made", procName, NULL);
pixAddConstantGray(pixsd, -height);
pixSeedfillGray(pixsd, pixs, connectivity);
pixd = pixSubtractGray(NULL, pixs, pixsd);
pixDestroy(&pixsd);
return pixd;
}
/*!
* pixFastTophat()
*
* Input: pixs
* xsize (width of max/min op, smoothing; any integer >= 1)
* ysize (height of max/min op, smoothing; any integer >= 1)
* type (L_TOPHAT_WHITE: image - min
* L_TOPHAT_BLACK: max - image)
* Return: pixd, or null on error
*
* Notes:
* (1) Don't be fooled. This is NOT a tophat. It is a tophat-like
* operation, where the result is similar to what you'd get
* if you used an erosion instead of an opening, or a dilation
* instead of a closing.
* (2) Instead of opening or closing at full resolution, it does
* a fast downscale/minmax operation, then a quick small smoothing
* at low res, a replicative expansion of the "background"
* to full res, and finally a removal of the background level
* from the input image. The smoothing step may not be important.
* (3) It does not remove noise as well as a tophat, but it is
* 5 to 10 times faster.
* If you need the preciseness of the tophat, don't use this.
* (4) The L_TOPHAT_WHITE flag emphasizes small bright regions,
* whereas the L_TOPHAT_BLACK flag emphasizes small dark regions.
*/
PIX *
pixFastTophat(PIX *pixs,
l_int32 xsize,
l_int32 ysize,
l_int32 type)
{
PIX *pixt1, *pixt2, *pixt3, *pixd;
PROCNAME("pixFastTophat");
if (!pixs)
return (PIX *)ERROR_PTR("seed pix not defined", procName, NULL);
if (pixGetDepth(pixs) != 8)
return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
if (xsize < 1 || ysize < 1)
return (PIX *)ERROR_PTR("size < 1", procName, NULL);
if (type != L_TOPHAT_WHITE && type != L_TOPHAT_BLACK)
return (PIX *)ERROR_PTR("type must be L_TOPHAT_BLACK or L_TOPHAT_WHITE",
procName, NULL);
if (xsize == 1 && ysize == 1)
return pixCreateTemplate(pixs);
switch (type)
{
case L_TOPHAT_WHITE:
if ((pixt1 = pixScaleGrayMinMax(pixs, xsize, ysize, L_CHOOSE_MIN))
== NULL)
return (PIX *)ERROR_PTR("pixt1 not made", procName, NULL);
pixt2 = pixBlockconv(pixt1, 1, 1); /* small smoothing */
pixt3 = pixScaleBySampling(pixt2, xsize, ysize);
pixd = pixSubtractGray(NULL, pixs, pixt3);
pixDestroy(&pixt3);
break;
case L_TOPHAT_BLACK:
if ((pixt1 = pixScaleGrayMinMax(pixs, xsize, ysize, L_CHOOSE_MAX))
== NULL)
return (PIX *)ERROR_PTR("pixt1 not made", procName, NULL);
pixt2 = pixBlockconv(pixt1, 1, 1); /* small smoothing */
pixd = pixScaleBySampling(pixt2, xsize, ysize);
pixSubtractGray(pixd, pixd, pixs);
break;
default:
return (PIX *)ERROR_PTR("invalid type", procName, NULL);
}
pixDestroy(&pixt1);
pixDestroy(&pixt2);
return pixd;
}
/*!
* pixMorphGradient()
*
* Input: pixs
* hsize (of Sel; must be odd; origin implicitly in center)
* vsize (ditto)
* smoothing (half-width of convolution smoothing filter.
* The width is (2 * smoothing + 1), so 0 is no-op.
* Return: pixd, or null on error
*/
PIX *
pixMorphGradient(PIX *pixs,
l_int32 hsize,
l_int32 vsize,
l_int32 smoothing)
{
PIX *pixg, *pixd;
PROCNAME("pixMorphGradient");
if (!pixs)
return (PIX *)ERROR_PTR("seed pix not defined", procName, NULL);
if (pixGetDepth(pixs) != 8)
return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
if (hsize < 1 || vsize < 1)
return (PIX *)ERROR_PTR("hsize or vsize < 1", procName, NULL);
if ((hsize & 1) == 0 ) {
L_WARNING("horiz sel size must be odd; increasing by 1", procName);
hsize++;
}
if ((vsize & 1) == 0 ) {
L_WARNING("vert sel size must be odd; increasing by 1", procName);
vsize++;
}
/* Optionally smooth first to remove noise.
* If smoothing is 0, just get a copy */
pixg = pixBlockconvGray(pixs, NULL, smoothing, smoothing);
/* This gives approximately the gradient of a transition */
pixd = pixDilateGray(pixg, hsize, vsize);
pixSubtractGray(pixd, pixd, pixg);
pixDestroy(&pixg);
return pixd;
}
/*-----------------------------------------------------------------*
* Center of mass *
*-----------------------------------------------------------------*/
/*!
* pixaCentroids()
*
* Input: pixa of components
* Return: pta of centroids relative to the UL corner of
* each pix, or null on error
*
* Notes:
* (1) It is assumed that all pix are the same depth.
* (2) Only depths of 1 and 8 bpp are allowed
*/
PTA *
pixaCentroids(PIXA *pixa)
{
l_int32 d, i, j, k, n, w, h, wpl, pixsum, rowsum, val;
l_float32 xsum, ysum, xave, yave;
l_uint32 *data, *line;
l_uint32 word;
l_uint8 byte;
PIX *pix;
PTA *pta;
static l_int32 *centtab = NULL;
static l_int32 *sumtab = NULL;
PROCNAME("pixaCentroids");
if (!pixa)
return (PTA *)ERROR_PTR("pixa not defined", procName, NULL);
if ((n = pixaGetCount(pixa)) == 0)
return (PTA *)ERROR_PTR("no pix in pixa", procName, NULL);
pix = pixaGetPix(pixa, 0, L_CLONE);
d = pixGetDepth(pix);
pixDestroy(&pix);
if (d != 1 && d != 8)
return (PTA *)ERROR_PTR("depth not 1 or 8 bpp", procName, NULL);
if ((pta = ptaCreate(n)) == NULL)
return (PTA *)ERROR_PTR("pta not defined", procName, NULL);
if ((centtab == NULL) &&
((centtab = makePixelCentroidTab8()) == NULL))
return (PTA *)ERROR_PTR("couldn't make centtab", procName, NULL);
if ((sumtab == NULL) &&
((sumtab = makePixelSumTab8()) == NULL))
return (PTA *)ERROR_PTR("couldn't make sumtab", procName, NULL);
for (k = 0; k < n; k++) {
pix = pixaGetPix(pixa, k, L_CLONE);
w = pixGetWidth(pix);
h = pixGetHeight(pix);
data = pixGetData(pix);
wpl = pixGetWpl(pix);
xsum = ysum = 0.0;
pixsum = 0;
if (d == 1) {
for (i = 0; i < h; i++) {
/* The body of this loop computes the sum of the set
* (1) bits on this row, weighted by their distance
* from the left edge of pix, and accumulates that into
* xsum; it accumulates their distance from the top
* edge of pix into ysum, and their total count into
* pixsum. It's equivalent to
* for (j = 0; j < w; j++) {
* if (GET_DATA_BIT(line, j)) {
* xsum += j;
* ysum += i;
* pixsum++;
* }
* }
*/
line = data + wpl * i;
rowsum = 0;
for (j = 0; j < wpl; j++) {
word = line[j];
if (word) {
byte = word & 0xff;
rowsum += sumtab[byte];
xsum += centtab[byte] + (j * 32 + 24) * sumtab[byte];
byte = (word >> 8) & 0xff;
rowsum += sumtab[byte];
xsum += centtab[byte] + (j * 32 + 16) * sumtab[byte];
byte = (word >> 16) & 0xff;
rowsum += sumtab[byte];
xsum += centtab[byte] + (j * 32 + 8) * sumtab[byte];
byte = (word >> 24) & 0xff;
rowsum += sumtab[byte];
xsum += centtab[byte] + j * 32 * sumtab[byte];
}
}
pixsum += rowsum;
ysum += rowsum * i;
}
if (pixsum == 0) {
L_WARNING("no ON pixels in pix", procName);
ptaAddPt(pta, 0.0, 0.0); /* this shouldn't happen */
}
else {
xave = xsum / (l_float32)pixsum;
yave = ysum / (l_float32)pixsum;
ptaAddPt(pta, xave, yave);
}
}
else { /* d == 8 */
for (i = 0; i < h; i++) {
line = data + wpl * i;
for (j = 0; j < w; j++) {
val = GET_DATA_BYTE(line, j);
xsum += val * j;
ysum += val * i;
pixsum += val;
}
}
if (pixsum == 0) {
L_WARNING("all pixels are 0", procName);
ptaAddPt(pta, 0.0, 0.0); /* this shouldn't happen */
}
else {
xave = xsum / (l_float32)pixsum;
yave = ysum / (l_float32)pixsum;
ptaAddPt(pta, xave, yave);
}
}
pixDestroy(&pix);
}
return pta;
}