#include "allheaders.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif /* M_PI */
/*---------------------------------------------------------------------*
* Simple rearrangements *
*---------------------------------------------------------------------*/
/*!
* \brief ptaSubsample()
*
* \param[in] ptas
* \param[in] subfactor subsample factor, >= 1
* \return ptad evenly sampled pt values from ptas, or NULL on error
*/
PTA *
ptaSubsample(PTA *ptas,
l_int32 subfactor)
{
l_int32 n, i;
l_float32 x, y;
PTA *ptad;
PROCNAME("pixSubsample");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
if (subfactor < 1)
return (PTA *)ERROR_PTR("subfactor < 1", procName, NULL);
ptad = ptaCreate(0);
n = ptaGetCount(ptas);
for (i = 0; i < n; i++) {
if (i % subfactor != 0) continue;
ptaGetPt(ptas, i, &x, &y);
ptaAddPt(ptad, x, y);
}
return ptad;
}
/*!
* \brief ptaJoin()
*
* \param[in] ptad dest pta; add to this one
* \param[in] ptas source pta; add from this one
* \param[in] istart starting index in ptas
* \param[in] iend ending index in ptas; use -1 to cat all
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) istart < 0 is taken to mean 'read from the start' (istart = 0)
* (2) iend < 0 means 'read to the end'
* (3) if ptas == NULL, this is a no-op
*
*/
l_ok
ptaJoin(PTA *ptad,
PTA *ptas,
l_int32 istart,
l_int32 iend)
{
l_int32 n, i, x, y;
PROCNAME("ptaJoin");
if (!ptad)
return ERROR_INT("ptad not defined", procName, 1);
if (!ptas)
return 0;
if (istart < 0)
istart = 0;
n = ptaGetCount(ptas);
if (iend < 0 || iend >= n)
iend = n - 1;
if (istart > iend)
return ERROR_INT("istart > iend; no pts", procName, 1);
for (i = istart; i <= iend; i++) {
ptaGetIPt(ptas, i, &x, &y);
ptaAddPt(ptad, x, y);
}
return 0;
}
/*!
* \brief ptaaJoin()
*
* \param[in] ptaad dest ptaa; add to this one
* \param[in] ptaas source ptaa; add from this one
* \param[in] istart starting index in ptaas
* \param[in] iend ending index in ptaas; use -1 to cat all
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) istart < 0 is taken to mean 'read from the start' (istart = 0)
* (2) iend < 0 means 'read to the end'
* (3) if ptas == NULL, this is a no-op
*
*/
l_ok
ptaaJoin(PTAA *ptaad,
PTAA *ptaas,
l_int32 istart,
l_int32 iend)
{
l_int32 n, i;
PTA *pta;
PROCNAME("ptaaJoin");
if (!ptaad)
return ERROR_INT("ptaad not defined", procName, 1);
if (!ptaas)
return 0;
if (istart < 0)
istart = 0;
n = ptaaGetCount(ptaas);
if (iend < 0 || iend >= n)
iend = n - 1;
if (istart > iend)
return ERROR_INT("istart > iend; no pts", procName, 1);
for (i = istart; i <= iend; i++) {
pta = ptaaGetPta(ptaas, i, L_CLONE);
ptaaAddPta(ptaad, pta, L_INSERT);
}
return 0;
}
/*!
* \brief ptaReverse()
*
* \param[in] ptas
* \param[in] type 0 for float values; 1 for integer values
* \return ptad reversed pta, or NULL on error
*/
PTA *
ptaReverse(PTA *ptas,
l_int32 type)
{
l_int32 n, i, ix, iy;
l_float32 x, y;
PTA *ptad;
PROCNAME("ptaReverse");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
n = ptaGetCount(ptas);
if ((ptad = ptaCreate(n)) == NULL)
return (PTA *)ERROR_PTR("ptad not made", procName, NULL);
for (i = n - 1; i >= 0; i--) {
if (type == 0) {
ptaGetPt(ptas, i, &x, &y);
ptaAddPt(ptad, x, y);
} else { /* type == 1 */
ptaGetIPt(ptas, i, &ix, &iy);
ptaAddPt(ptad, ix, iy);
}
}
return ptad;
}
/*!
* \brief ptaTranspose()
*
* \param[in] ptas
* \return ptad with x and y values swapped, or NULL on error
*/
PTA *
ptaTranspose(PTA *ptas)
{
l_int32 n, i;
l_float32 x, y;
PTA *ptad;
PROCNAME("ptaTranspose");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
n = ptaGetCount(ptas);
if ((ptad = ptaCreate(n)) == NULL)
return (PTA *)ERROR_PTR("ptad not made", procName, NULL);
for (i = 0; i < n; i++) {
ptaGetPt(ptas, i, &x, &y);
ptaAddPt(ptad, y, x);
}
return ptad;
}
/*!
* \brief ptaCyclicPerm()
*
* \param[in] ptas
* \param[in] xs, ys start point; must be in ptas
* \return ptad cyclic permutation, starting and ending at (xs, ys,
* or NULL on error
*
*
* Notes:
* (1) Check to insure that (a) ptas is a closed path where
* the first and last points are identical, and (b) the
* resulting pta also starts and ends on the same point
* (which in this case is (xs, ys).
*
*/
PTA *
ptaCyclicPerm(PTA *ptas,
l_int32 xs,
l_int32 ys)
{
l_int32 n, i, x, y, j, index, state;
l_int32 x1, y1, x2, y2;
PTA *ptad;
PROCNAME("ptaCyclicPerm");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
n = ptaGetCount(ptas);
/* Verify input data */
ptaGetIPt(ptas, 0, &x1, &y1);
ptaGetIPt(ptas, n - 1, &x2, &y2);
if (x1 != x2 || y1 != y2)
return (PTA *)ERROR_PTR("start and end pts not same", procName, NULL);
state = L_NOT_FOUND;
for (i = 0; i < n; i++) {
ptaGetIPt(ptas, i, &x, &y);
if (x == xs && y == ys) {
state = L_FOUND;
break;
}
}
if (state == L_NOT_FOUND)
return (PTA *)ERROR_PTR("start pt not in ptas", procName, NULL);
if ((ptad = ptaCreate(n)) == NULL)
return (PTA *)ERROR_PTR("ptad not made", procName, NULL);
for (j = 0; j < n - 1; j++) {
if (i + j < n - 1)
index = i + j;
else
index = (i + j + 1) % n;
ptaGetIPt(ptas, index, &x, &y);
ptaAddPt(ptad, x, y);
}
ptaAddPt(ptad, xs, ys);
return ptad;
}
/*!
* \brief ptaSelectRange()
*
* \param[in] ptas
* \param[in] first use 0 to select from the beginning
* \param[in] last use -1 to select to the end
* \return ptad, or NULL on error
*/
PTA *
ptaSelectRange(PTA *ptas,
l_int32 first,
l_int32 last)
{
l_int32 n, npt, i;
l_float32 x, y;
PTA *ptad;
PROCNAME("ptaSelectRange");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
if ((n = ptaGetCount(ptas)) == 0) {
L_WARNING("ptas is empty\n", procName);
return ptaCopy(ptas);
}
first = L_MAX(0, first);
if (last < 0) last = n - 1;
if (first >= n)
return (PTA *)ERROR_PTR("invalid first", procName, NULL);
if (last >= n) {
L_WARNING("last = %d is beyond max index = %d; adjusting\n",
procName, last, n - 1);
last = n - 1;
}
if (first > last)
return (PTA *)ERROR_PTR("first > last", procName, NULL);
npt = last - first + 1;
ptad = ptaCreate(npt);
for (i = first; i <= last; i++) {
ptaGetPt(ptas, i, &x, &y);
ptaAddPt(ptad, x, y);
}
return ptad;
}
/*---------------------------------------------------------------------*
* Geometric *
*---------------------------------------------------------------------*/
/*!
* \brief ptaGetBoundingRegion()
*
* \param[in] pta
* \return box, or NULL on error
*
*
* Notes:
* (1) This is used when the pta represents a set of points in
* a two-dimensional image. It returns the box of minimum
* size containing the pts in the pta.
*
*/
BOX *
ptaGetBoundingRegion(PTA *pta)
{
l_int32 n, i, x, y, minx, maxx, miny, maxy;
PROCNAME("ptaGetBoundingRegion");
if (!pta)
return (BOX *)ERROR_PTR("pta not defined", procName, NULL);
minx = 10000000;
miny = 10000000;
maxx = -10000000;
maxy = -10000000;
n = ptaGetCount(pta);
for (i = 0; i < n; i++) {
ptaGetIPt(pta, i, &x, &y);
if (x < minx) minx = x;
if (x > maxx) maxx = x;
if (y < miny) miny = y;
if (y > maxy) maxy = y;
}
return boxCreate(minx, miny, maxx - minx + 1, maxy - miny + 1);
}
/*!
* \brief ptaGetRange()
*
* \param[in] pta
* \param[out] pminx [optional] min value of x
* \param[out] pmaxx [optional] max value of x
* \param[out] pminy [optional] min value of y
* \param[out] pmaxy [optional] max value of y
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) We can use pts to represent pairs of floating values, that
* are not necessarily tied to a two-dimension region. For
* example, the pts can represent a general function y(x).
*
*/
l_ok
ptaGetRange(PTA *pta,
l_float32 *pminx,
l_float32 *pmaxx,
l_float32 *pminy,
l_float32 *pmaxy)
{
l_int32 n, i;
l_float32 x, y, minx, maxx, miny, maxy;
PROCNAME("ptaGetRange");
if (!pminx && !pmaxx && !pminy && !pmaxy)
return ERROR_INT("no output requested", procName, 1);
if (pminx) *pminx = 0;
if (pmaxx) *pmaxx = 0;
if (pminy) *pminy = 0;
if (pmaxy) *pmaxy = 0;
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if ((n = ptaGetCount(pta)) == 0)
return ERROR_INT("no points in pta", procName, 1);
ptaGetPt(pta, 0, &x, &y);
minx = x;
maxx = x;
miny = y;
maxy = y;
for (i = 1; i < n; i++) {
ptaGetPt(pta, i, &x, &y);
if (x < minx) minx = x;
if (x > maxx) maxx = x;
if (y < miny) miny = y;
if (y > maxy) maxy = y;
}
if (pminx) *pminx = minx;
if (pmaxx) *pmaxx = maxx;
if (pminy) *pminy = miny;
if (pmaxy) *pmaxy = maxy;
return 0;
}
/*!
* \brief ptaGetInsideBox()
*
* \param[in] ptas input pts
* \param[in] box
* \return ptad of pts in ptas that are inside the box, or NULL on error
*/
PTA *
ptaGetInsideBox(PTA *ptas,
BOX *box)
{
PTA *ptad;
l_int32 n, i, contains;
l_float32 x, y;
PROCNAME("ptaGetInsideBox");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
if (!box)
return (PTA *)ERROR_PTR("box not defined", procName, NULL);
n = ptaGetCount(ptas);
ptad = ptaCreate(0);
for (i = 0; i < n; i++) {
ptaGetPt(ptas, i, &x, &y);
boxContainsPt(box, x, y, &contains);
if (contains)
ptaAddPt(ptad, x, y);
}
return ptad;
}
/*!
* \brief pixFindCornerPixels()
*
* \param[in] pixs 1 bpp
* \return pta, or NULL on error
*
*
* Notes:
* (1) Finds the 4 corner-most pixels, as defined by a search
* inward from each corner, using a 45 degree line.
*
*/
PTA *
pixFindCornerPixels(PIX *pixs)
{
l_int32 i, j, x, y, w, h, wpl, mindim, found;
l_uint32 *data, *line;
PTA *pta;
PROCNAME("pixFindCornerPixels");
if (!pixs)
return (PTA *)ERROR_PTR("pixs not defined", procName, NULL);
if (pixGetDepth(pixs) != 1)
return (PTA *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
w = pixGetWidth(pixs);
h = pixGetHeight(pixs);
mindim = L_MIN(w, h);
data = pixGetData(pixs);
wpl = pixGetWpl(pixs);
if ((pta = ptaCreate(4)) == NULL)
return (PTA *)ERROR_PTR("pta not made", procName, NULL);
for (found = FALSE, i = 0; i < mindim; i++) {
for (j = 0; j <= i; j++) {
y = i - j;
line = data + y * wpl;
if (GET_DATA_BIT(line, j)) {
ptaAddPt(pta, j, y);
found = TRUE;
break;
}
}
if (found == TRUE)
break;
}
for (found = FALSE, i = 0; i < mindim; i++) {
for (j = 0; j <= i; j++) {
y = i - j;
line = data + y * wpl;
x = w - 1 - j;
if (GET_DATA_BIT(line, x)) {
ptaAddPt(pta, x, y);
found = TRUE;
break;
}
}
if (found == TRUE)
break;
}
for (found = FALSE, i = 0; i < mindim; i++) {
for (j = 0; j <= i; j++) {
y = h - 1 - i + j;
line = data + y * wpl;
if (GET_DATA_BIT(line, j)) {
ptaAddPt(pta, j, y);
found = TRUE;
break;
}
}
if (found == TRUE)
break;
}
for (found = FALSE, i = 0; i < mindim; i++) {
for (j = 0; j <= i; j++) {
y = h - 1 - i + j;
line = data + y * wpl;
x = w - 1 - j;
if (GET_DATA_BIT(line, x)) {
ptaAddPt(pta, x, y);
found = TRUE;
break;
}
}
if (found == TRUE)
break;
}
return pta;
}
/*!
* \brief ptaContainsPt()
*
* \param[in] pta
* \param[in] x, y point
* \return 1 if contained, 0 otherwise or on error
*/
l_int32
ptaContainsPt(PTA *pta,
l_int32 x,
l_int32 y)
{
l_int32 i, n, ix, iy;
PROCNAME("ptaContainsPt");
if (!pta)
return ERROR_INT("pta not defined", procName, 0);
n = ptaGetCount(pta);
for (i = 0; i < n; i++) {
ptaGetIPt(pta, i, &ix, &iy);
if (x == ix && y == iy)
return 1;
}
return 0;
}
/*!
* \brief ptaTestIntersection()
*
* \param[in] pta1, pta2
* \return bval which is 1 if they have any elements in common;
* 0 otherwise or on error.
*/
l_int32
ptaTestIntersection(PTA *pta1,
PTA *pta2)
{
l_int32 i, j, n1, n2, x1, y1, x2, y2;
PROCNAME("ptaTestIntersection");
if (!pta1)
return ERROR_INT("pta1 not defined", procName, 0);
if (!pta2)
return ERROR_INT("pta2 not defined", procName, 0);
n1 = ptaGetCount(pta1);
n2 = ptaGetCount(pta2);
for (i = 0; i < n1; i++) {
ptaGetIPt(pta1, i, &x1, &y1);
for (j = 0; j < n2; j++) {
ptaGetIPt(pta2, i, &x2, &y2);
if (x1 == x2 && y1 == y2)
return 1;
}
}
return 0;
}
/*!
* \brief ptaTransform()
*
* \param[in] ptas
* \param[in] shiftx, shifty
* \param[in] scalex, scaley
* \return pta, or NULL on error
*
*
* Notes:
* (1) Shift first, then scale.
*
*/
PTA *
ptaTransform(PTA *ptas,
l_int32 shiftx,
l_int32 shifty,
l_float32 scalex,
l_float32 scaley)
{
l_int32 n, i, x, y;
PTA *ptad;
PROCNAME("ptaTransform");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
n = ptaGetCount(ptas);
ptad = ptaCreate(n);
for (i = 0; i < n; i++) {
ptaGetIPt(ptas, i, &x, &y);
x = (l_int32)(scalex * (x + shiftx) + 0.5);
y = (l_int32)(scaley * (y + shifty) + 0.5);
ptaAddPt(ptad, x, y);
}
return ptad;
}
/*!
* \brief ptaPtInsidePolygon()
*
* \param[in] pta vertices of a polygon
* \param[in] x, y point to be tested
* \param[out] pinside 1 if inside; 0 if outside or on boundary
* \return 1 if OK, 0 on error
*
* The abs value of the sum of the angles subtended from a point by
* the sides of a polygon, when taken in order traversing the polygon,
* is 0 if the point is outside the polygon and 2*pi if inside.
* The sign will be positive if traversed cw and negative if ccw.
*/
l_int32
ptaPtInsidePolygon(PTA *pta,
l_float32 x,
l_float32 y,
l_int32 *pinside)
{
l_int32 i, n;
l_float32 sum, x1, y1, x2, y2, xp1, yp1, xp2, yp2;
PROCNAME("ptaPtInsidePolygon");
if (!pinside)
return ERROR_INT("&inside not defined", procName, 1);
*pinside = 0;
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
/* Think of (x1,y1) as the end point of a vector that starts
* from the origin (0,0), and ditto for (x2,y2). */
n = ptaGetCount(pta);
sum = 0.0;
for (i = 0; i < n; i++) {
ptaGetPt(pta, i, &xp1, &yp1);
ptaGetPt(pta, (i + 1) % n, &xp2, &yp2);
x1 = xp1 - x;
y1 = yp1 - y;
x2 = xp2 - x;
y2 = yp2 - y;
sum += l_angleBetweenVectors(x1, y1, x2, y2);
}
if (L_ABS(sum) > M_PI)
*pinside = 1;
return 0;
}
/*!
* \brief l_angleBetweenVectors()
*
* \param[in] x1, y1 end point of first vector
* \param[in] x2, y2 end point of second vector
* \return angle radians, or 0.0 on error
*
*
* Notes:
* (1) This gives the angle between two vectors, going between
* vector1 (x1,y1) and vector2 (x2,y2). The angle is swept
* out from 1 --> 2. If this is clockwise, the angle is
* positive, but the result is folded into the interval [-pi, pi].
*
*/
l_float32
l_angleBetweenVectors(l_float32 x1,
l_float32 y1,
l_float32 x2,
l_float32 y2)
{
l_float64 ang;
ang = atan2(y2, x2) - atan2(y1, x1);
if (ang > M_PI) ang -= 2.0 * M_PI;
if (ang < -M_PI) ang += 2.0 * M_PI;
return ang;
}
/*---------------------------------------------------------------------*
* Min/max and filtering *
*---------------------------------------------------------------------*/
/*!
* \brief ptaGetMinMax()
*
* \param[in] pta
* \param[out] pxmin [optional] min of x
* \param[out] pymin [optional] min of y
* \param[out] pxmax [optional] max of x
* \param[out] pymax [optional] max of y
* \return 0 if OK, 1 on error. If pta is empty, requested
* values are returned as -1.0.
*/
l_ok
ptaGetMinMax(PTA *pta,
l_float32 *pxmin,
l_float32 *pymin,
l_float32 *pxmax,
l_float32 *pymax)
{
l_int32 i, n;
l_float32 x, y, xmin, ymin, xmax, ymax;
PROCNAME("ptaGetMinMax");
if (pxmin) *pxmin = -1.0;
if (pymin) *pymin = -1.0;
if (pxmax) *pxmax = -1.0;
if (pymax) *pymax = -1.0;
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if (!pxmin && !pxmax && !pymin && !pymax)
return ERROR_INT("no output requested", procName, 1);
if ((n = ptaGetCount(pta)) == 0) {
L_WARNING("pta is empty\n", procName);
return 0;
}
xmin = ymin = 1.0e20f;
xmax = ymax = -1.0e20f;
for (i = 0; i < n; i++) {
ptaGetPt(pta, i, &x, &y);
if (x < xmin) xmin = x;
if (y < ymin) ymin = y;
if (x > xmax) xmax = x;
if (y > ymax) ymax = y;
}
if (pxmin) *pxmin = xmin;
if (pymin) *pymin = ymin;
if (pxmax) *pxmax = xmax;
if (pymax) *pymax = ymax;
return 0;
}
/*!
* \brief ptaSelectByValue()
*
* \param[in] ptas
* \param[in] xth, yth threshold values
* \param[in] type L_SELECT_XVAL, L_SELECT_YVAL,
* L_SELECT_IF_EITHER, L_SELECT_IF_BOTH
* \param[in] relation L_SELECT_IF_LT, L_SELECT_IF_GT,
* L_SELECT_IF_LTE, L_SELECT_IF_GTE
* \return ptad filtered set, or NULL on error
*/
PTA *
ptaSelectByValue(PTA *ptas,
l_float32 xth,
l_float32 yth,
l_int32 type,
l_int32 relation)
{
l_int32 i, n;
l_float32 x, y;
PTA *ptad;
PROCNAME("ptaSelectByValue");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
if (ptaGetCount(ptas) == 0) {
L_WARNING("ptas is empty\n", procName);
return ptaCopy(ptas);
}
if (type != L_SELECT_XVAL && type != L_SELECT_YVAL &&
type != L_SELECT_IF_EITHER && type != L_SELECT_IF_BOTH)
return (PTA *)ERROR_PTR("invalid type", procName, NULL);
if (relation != L_SELECT_IF_LT && relation != L_SELECT_IF_GT &&
relation != L_SELECT_IF_LTE && relation != L_SELECT_IF_GTE)
return (PTA *)ERROR_PTR("invalid relation", procName, NULL);
n = ptaGetCount(ptas);
ptad = ptaCreate(n);
for (i = 0; i < n; i++) {
ptaGetPt(ptas, i, &x, &y);
if (type == L_SELECT_XVAL) {
if ((relation == L_SELECT_IF_LT && x < xth) ||
(relation == L_SELECT_IF_GT && x > xth) ||
(relation == L_SELECT_IF_LTE && x <= xth) ||
(relation == L_SELECT_IF_GTE && x >= xth))
ptaAddPt(ptad, x, y);
} else if (type == L_SELECT_YVAL) {
if ((relation == L_SELECT_IF_LT && y < yth) ||
(relation == L_SELECT_IF_GT && y > yth) ||
(relation == L_SELECT_IF_LTE && y <= yth) ||
(relation == L_SELECT_IF_GTE && y >= yth))
ptaAddPt(ptad, x, y);
} else if (type == L_SELECT_IF_EITHER) {
if (((relation == L_SELECT_IF_LT) && (x < xth || y < yth)) ||
((relation == L_SELECT_IF_GT) && (x > xth || y > yth)) ||
((relation == L_SELECT_IF_LTE) && (x <= xth || y <= yth)) ||
((relation == L_SELECT_IF_GTE) && (x >= xth || y >= yth)))
ptaAddPt(ptad, x, y);
} else { /* L_SELECT_IF_BOTH */
if (((relation == L_SELECT_IF_LT) && (x < xth && y < yth)) ||
((relation == L_SELECT_IF_GT) && (x > xth && y > yth)) ||
((relation == L_SELECT_IF_LTE) && (x <= xth && y <= yth)) ||
((relation == L_SELECT_IF_GTE) && (x >= xth && y >= yth)))
ptaAddPt(ptad, x, y);
}
}
return ptad;
}
/*!
* \brief ptaCropToMask()
*
* \param[in] ptas input pta
* \param[in] pixm 1 bpp mask
* \return ptad with only pts under the mask fg, or NULL on error
*/
PTA *
ptaCropToMask(PTA *ptas,
PIX *pixm)
{
l_int32 i, n, x, y;
l_uint32 val;
PTA *ptad;
PROCNAME("ptaCropToMask");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
if (!pixm || pixGetDepth(pixm) != 1)
return (PTA *)ERROR_PTR("pixm undefined or not 1 bpp", procName, NULL);
if (ptaGetCount(ptas) == 0) {
L_INFO("ptas is empty\n", procName);
return ptaCopy(ptas);
}
n = ptaGetCount(ptas);
ptad = ptaCreate(n);
for (i = 0; i < n; i++) {
ptaGetIPt(ptas, i, &x, &y);
pixGetPixel(pixm, x, y, &val);
if (val == 1)
ptaAddPt(ptad, x, y);
}
return ptad;
}
/*---------------------------------------------------------------------*
* Least Squares Fit *
*---------------------------------------------------------------------*/
/*!
* \brief ptaGetLinearLSF()
*
* \param[in] pta
* \param[out] pa [optional] slope a of least square fit: y = ax + b
* \param[out] pb [optional] intercept b of least square fit
* \param[out] pnafit [optional] numa of least square fit
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) Either or both &a and &b must be input. They determine the
* type of line that is fit.
* (2) If both &a and &b are defined, this returns a and b that minimize:
*
* sum (yi - axi -b)^2
* i
*
* The method is simple: differentiate this expression w/rt a and b,
* and solve the resulting two equations for a and b in terms of
* various sums over the input data (xi, yi).
* (3) We also allow two special cases, where either a = 0 or b = 0:
* (a) If &a is given and &b = null, find the linear LSF that
* goes through the origin (b = 0).
* (b) If &b is given and &a = null, find the linear LSF with
* zero slope (a = 0).
* (4) If &nafit is defined, this returns an array of fitted values,
* corresponding to the two implicit Numa arrays (nax and nay) in pta.
* Thus, just as you can plot the data in pta as nay vs. nax,
* you can plot the linear least square fit as nafit vs. nax.
* Get the nax array using ptaGetArrays(pta, &nax, NULL);
*
*/
l_ok
ptaGetLinearLSF(PTA *pta,
l_float32 *pa,
l_float32 *pb,
NUMA **pnafit)
{
l_int32 n, i;
l_float32 a, b, factor, sx, sy, sxx, sxy, val;
l_float32 *xa, *ya;
PROCNAME("ptaGetLinearLSF");
if (pa) *pa = 0.0;
if (pb) *pb = 0.0;
if (pnafit) *pnafit = NULL;
if (!pa && !pb && !pnafit)
return ERROR_INT("no output requested", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if ((n = ptaGetCount(pta)) < 2)
return ERROR_INT("less than 2 pts found", procName, 1);
xa = pta->x; /* not a copy */
ya = pta->y; /* not a copy */
sx = sy = sxx = sxy = 0.;
if (pa && pb) { /* general line */
for (i = 0; i < n; i++) {
sx += xa[i];
sy += ya[i];
sxx += xa[i] * xa[i];
sxy += xa[i] * ya[i];
}
factor = n * sxx - sx * sx;
if (factor == 0.0)
return ERROR_INT("no solution found", procName, 1);
factor = 1. / factor;
a = factor * ((l_float32)n * sxy - sx * sy);
b = factor * (sxx * sy - sx * sxy);
} else if (pa) { /* b = 0; line through origin */
for (i = 0; i < n; i++) {
sxx += xa[i] * xa[i];
sxy += xa[i] * ya[i];
}
if (sxx == 0.0)
return ERROR_INT("no solution found", procName, 1);
a = sxy / sxx;
b = 0.0;
} else { /* a = 0; horizontal line */
for (i = 0; i < n; i++)
sy += ya[i];
a = 0.0;
b = sy / (l_float32)n;
}
if (pnafit) {
*pnafit = numaCreate(n);
for (i = 0; i < n; i++) {
val = a * xa[i] + b;
numaAddNumber(*pnafit, val);
}
}
if (pa) *pa = a;
if (pb) *pb = b;
return 0;
}
/*!
* \brief ptaGetQuadraticLSF()
*
* \param[in] pta
* \param[out] pa [optional] coeff a of LSF: y = ax^2 + bx + c
* \param[out] pb [optional] coeff b of LSF: y = ax^2 + bx + c
* \param[out] pc [optional] coeff c of LSF: y = ax^2 + bx + c
* \param[out] pnafit [optional] numa of least square fit
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) This does a quadratic least square fit to the set of points
* in %pta. That is, it finds coefficients a, b and c that minimize:
*
* sum (yi - a*xi*xi -b*xi -c)^2
* i
*
* The method is simple: differentiate this expression w/rt
* a, b and c, and solve the resulting three equations for these
* coefficients in terms of various sums over the input data (xi, yi).
* The three equations are in the form:
* f[0][0]a + f[0][1]b + f[0][2]c = g[0]
* f[1][0]a + f[1][1]b + f[1][2]c = g[1]
* f[2][0]a + f[2][1]b + f[2][2]c = g[2]
* (2) If &nafit is defined, this returns an array of fitted values,
* corresponding to the two implicit Numa arrays (nax and nay) in pta.
* Thus, just as you can plot the data in pta as nay vs. nax,
* you can plot the linear least square fit as nafit vs. nax.
* Get the nax array using ptaGetArrays(pta, &nax, NULL);
*
*/
l_ok
ptaGetQuadraticLSF(PTA *pta,
l_float32 *pa,
l_float32 *pb,
l_float32 *pc,
NUMA **pnafit)
{
l_int32 n, i, ret;
l_float32 x, y, sx, sy, sx2, sx3, sx4, sxy, sx2y;
l_float32 *xa, *ya;
l_float32 *f[3];
l_float32 g[3];
PROCNAME("ptaGetQuadraticLSF");
if (pa) *pa = 0.0;
if (pb) *pb = 0.0;
if (pc) *pc = 0.0;
if (pnafit) *pnafit = NULL;
if (!pa && !pb && !pc && !pnafit)
return ERROR_INT("no output requested", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if ((n = ptaGetCount(pta)) < 3)
return ERROR_INT("less than 3 pts found", procName, 1);
xa = pta->x; /* not a copy */
ya = pta->y; /* not a copy */
sx = sy = sx2 = sx3 = sx4 = sxy = sx2y = 0.;
for (i = 0; i < n; i++) {
x = xa[i];
y = ya[i];
sx += x;
sy += y;
sx2 += x * x;
sx3 += x * x * x;
sx4 += x * x * x * x;
sxy += x * y;
sx2y += x * x * y;
}
for (i = 0; i < 3; i++)
f[i] = (l_float32 *)LEPT_CALLOC(3, sizeof(l_float32));
f[0][0] = sx4;
f[0][1] = sx3;
f[0][2] = sx2;
f[1][0] = sx3;
f[1][1] = sx2;
f[1][2] = sx;
f[2][0] = sx2;
f[2][1] = sx;
f[2][2] = n;
g[0] = sx2y;
g[1] = sxy;
g[2] = sy;
/* Solve for the unknowns, also putting f-inverse into f */
ret = gaussjordan(f, g, 3);
for (i = 0; i < 3; i++)
LEPT_FREE(f[i]);
if (ret)
return ERROR_INT("quadratic solution failed", procName, 1);
if (pa) *pa = g[0];
if (pb) *pb = g[1];
if (pc) *pc = g[2];
if (pnafit) {
*pnafit = numaCreate(n);
for (i = 0; i < n; i++) {
x = xa[i];
y = g[0] * x * x + g[1] * x + g[2];
numaAddNumber(*pnafit, y);
}
}
return 0;
}
/*!
* \brief ptaGetCubicLSF()
*
* \param[in] pta
* \param[out] pa [optional] coeff a of LSF: y = ax^3 + bx^2 + cx + d
* \param[out] pb [optional] coeff b of LSF
* \param[out] pc [optional] coeff c of LSF
* \param[out] pd [optional] coeff d of LSF
* \param[out] pnafit [optional] numa of least square fit
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) This does a cubic least square fit to the set of points
* in %pta. That is, it finds coefficients a, b, c and d
* that minimize:
*
* sum (yi - a*xi*xi*xi -b*xi*xi -c*xi - d)^2
* i
*
* Differentiate this expression w/rt a, b, c and d, and solve
* the resulting four equations for these coefficients in
* terms of various sums over the input data (xi, yi).
* The four equations are in the form:
* f[0][0]a + f[0][1]b + f[0][2]c + f[0][3] = g[0]
* f[1][0]a + f[1][1]b + f[1][2]c + f[1][3] = g[1]
* f[2][0]a + f[2][1]b + f[2][2]c + f[2][3] = g[2]
* f[3][0]a + f[3][1]b + f[3][2]c + f[3][3] = g[3]
* (2) If &nafit is defined, this returns an array of fitted values,
* corresponding to the two implicit Numa arrays (nax and nay) in pta.
* Thus, just as you can plot the data in pta as nay vs. nax,
* you can plot the linear least square fit as nafit vs. nax.
* Get the nax array using ptaGetArrays(pta, &nax, NULL);
*
*/
l_ok
ptaGetCubicLSF(PTA *pta,
l_float32 *pa,
l_float32 *pb,
l_float32 *pc,
l_float32 *pd,
NUMA **pnafit)
{
l_int32 n, i, ret;
l_float32 x, y, sx, sy, sx2, sx3, sx4, sx5, sx6, sxy, sx2y, sx3y;
l_float32 *xa, *ya;
l_float32 *f[4];
l_float32 g[4];
PROCNAME("ptaGetCubicLSF");
if (pa) *pa = 0.0;
if (pb) *pb = 0.0;
if (pc) *pc = 0.0;
if (pd) *pd = 0.0;
if (pnafit) *pnafit = NULL;
if (!pa && !pb && !pc && !pd && !pnafit)
return ERROR_INT("no output requested", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if ((n = ptaGetCount(pta)) < 4)
return ERROR_INT("less than 4 pts found", procName, 1);
xa = pta->x; /* not a copy */
ya = pta->y; /* not a copy */
sx = sy = sx2 = sx3 = sx4 = sx5 = sx6 = sxy = sx2y = sx3y = 0.;
for (i = 0; i < n; i++) {
x = xa[i];
y = ya[i];
sx += x;
sy += y;
sx2 += x * x;
sx3 += x * x * x;
sx4 += x * x * x * x;
sx5 += x * x * x * x * x;
sx6 += x * x * x * x * x * x;
sxy += x * y;
sx2y += x * x * y;
sx3y += x * x * x * y;
}
for (i = 0; i < 4; i++)
f[i] = (l_float32 *)LEPT_CALLOC(4, sizeof(l_float32));
f[0][0] = sx6;
f[0][1] = sx5;
f[0][2] = sx4;
f[0][3] = sx3;
f[1][0] = sx5;
f[1][1] = sx4;
f[1][2] = sx3;
f[1][3] = sx2;
f[2][0] = sx4;
f[2][1] = sx3;
f[2][2] = sx2;
f[2][3] = sx;
f[3][0] = sx3;
f[3][1] = sx2;
f[3][2] = sx;
f[3][3] = n;
g[0] = sx3y;
g[1] = sx2y;
g[2] = sxy;
g[3] = sy;
/* Solve for the unknowns, also putting f-inverse into f */
ret = gaussjordan(f, g, 4);
for (i = 0; i < 4; i++)
LEPT_FREE(f[i]);
if (ret)
return ERROR_INT("cubic solution failed", procName, 1);
if (pa) *pa = g[0];
if (pb) *pb = g[1];
if (pc) *pc = g[2];
if (pd) *pd = g[3];
if (pnafit) {
*pnafit = numaCreate(n);
for (i = 0; i < n; i++) {
x = xa[i];
y = g[0] * x * x * x + g[1] * x * x + g[2] * x + g[3];
numaAddNumber(*pnafit, y);
}
}
return 0;
}
/*!
* \brief ptaGetQuarticLSF()
*
* \param[in] pta
* \param[out] pa [optional] coeff a of LSF:
* y = ax^4 + bx^3 + cx^2 + dx + e
* \param[out] pb [optional] coeff b of LSF
* \param[out] pc [optional] coeff c of LSF
* \param[out] pd [optional] coeff d of LSF
* \param[out] pe [optional] coeff e of LSF
* \param[out] pnafit [optional] numa of least square fit
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) This does a quartic least square fit to the set of points
* in %pta. That is, it finds coefficients a, b, c, d and 3
* that minimize:
*
* sum (yi - a*xi*xi*xi*xi -b*xi*xi*xi -c*xi*xi - d*xi - e)^2
* i
*
* Differentiate this expression w/rt a, b, c, d and e, and solve
* the resulting five equations for these coefficients in
* terms of various sums over the input data (xi, yi).
* The five equations are in the form:
* f[0][0]a + f[0][1]b + f[0][2]c + f[0][3] + f[0][4] = g[0]
* f[1][0]a + f[1][1]b + f[1][2]c + f[1][3] + f[1][4] = g[1]
* f[2][0]a + f[2][1]b + f[2][2]c + f[2][3] + f[2][4] = g[2]
* f[3][0]a + f[3][1]b + f[3][2]c + f[3][3] + f[3][4] = g[3]
* f[4][0]a + f[4][1]b + f[4][2]c + f[4][3] + f[4][4] = g[4]
* (2) If &nafit is defined, this returns an array of fitted values,
* corresponding to the two implicit Numa arrays (nax and nay) in pta.
* Thus, just as you can plot the data in pta as nay vs. nax,
* you can plot the linear least square fit as nafit vs. nax.
* Get the nax array using ptaGetArrays(pta, &nax, NULL);
*
*/
l_ok
ptaGetQuarticLSF(PTA *pta,
l_float32 *pa,
l_float32 *pb,
l_float32 *pc,
l_float32 *pd,
l_float32 *pe,
NUMA **pnafit)
{
l_int32 n, i, ret;
l_float32 x, y, sx, sy, sx2, sx3, sx4, sx5, sx6, sx7, sx8;
l_float32 sxy, sx2y, sx3y, sx4y;
l_float32 *xa, *ya;
l_float32 *f[5];
l_float32 g[5];
PROCNAME("ptaGetQuarticLSF");
if (pa) *pa = 0.0;
if (pb) *pb = 0.0;
if (pc) *pc = 0.0;
if (pd) *pd = 0.0;
if (pe) *pe = 0.0;
if (pnafit) *pnafit = NULL;
if (!pa && !pb && !pc && !pd && !pe && !pnafit)
return ERROR_INT("no output requested", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if ((n = ptaGetCount(pta)) < 5)
return ERROR_INT("less than 5 pts found", procName, 1);
xa = pta->x; /* not a copy */
ya = pta->y; /* not a copy */
sx = sy = sx2 = sx3 = sx4 = sx5 = sx6 = sx7 = sx8 = 0;
sxy = sx2y = sx3y = sx4y = 0.;
for (i = 0; i < n; i++) {
x = xa[i];
y = ya[i];
sx += x;
sy += y;
sx2 += x * x;
sx3 += x * x * x;
sx4 += x * x * x * x;
sx5 += x * x * x * x * x;
sx6 += x * x * x * x * x * x;
sx7 += x * x * x * x * x * x * x;
sx8 += x * x * x * x * x * x * x * x;
sxy += x * y;
sx2y += x * x * y;
sx3y += x * x * x * y;
sx4y += x * x * x * x * y;
}
for (i = 0; i < 5; i++)
f[i] = (l_float32 *)LEPT_CALLOC(5, sizeof(l_float32));
f[0][0] = sx8;
f[0][1] = sx7;
f[0][2] = sx6;
f[0][3] = sx5;
f[0][4] = sx4;
f[1][0] = sx7;
f[1][1] = sx6;
f[1][2] = sx5;
f[1][3] = sx4;
f[1][4] = sx3;
f[2][0] = sx6;
f[2][1] = sx5;
f[2][2] = sx4;
f[2][3] = sx3;
f[2][4] = sx2;
f[3][0] = sx5;
f[3][1] = sx4;
f[3][2] = sx3;
f[3][3] = sx2;
f[3][4] = sx;
f[4][0] = sx4;
f[4][1] = sx3;
f[4][2] = sx2;
f[4][3] = sx;
f[4][4] = n;
g[0] = sx4y;
g[1] = sx3y;
g[2] = sx2y;
g[3] = sxy;
g[4] = sy;
/* Solve for the unknowns, also putting f-inverse into f */
ret = gaussjordan(f, g, 5);
for (i = 0; i < 5; i++)
LEPT_FREE(f[i]);
if (ret)
return ERROR_INT("quartic solution failed", procName, 1);
if (pa) *pa = g[0];
if (pb) *pb = g[1];
if (pc) *pc = g[2];
if (pd) *pd = g[3];
if (pe) *pe = g[4];
if (pnafit) {
*pnafit = numaCreate(n);
for (i = 0; i < n; i++) {
x = xa[i];
y = g[0] * x * x * x * x + g[1] * x * x * x + g[2] * x * x
+ g[3] * x + g[4];
numaAddNumber(*pnafit, y);
}
}
return 0;
}
/*!
* \brief ptaNoisyLinearLSF()
*
* \param[in] pta
* \param[in] factor reject outliers with error greater than this
* number of medians; typically ~ 3
* \param[out] pptad [optional] with outliers removed
* \param[out] pa [optional] slope a of least square fit: y = ax + b
* \param[out] pb [optional] intercept b of least square fit
* \param[out] pmederr [optional] median error
* \param[out] pnafit [optional] numa of least square fit to ptad
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) This does a linear least square fit to the set of points
* in %pta. It then evaluates the errors and removes points
* whose error is >= factor * median_error. It then re-runs
* the linear LSF on the resulting points.
* (2) Either or both &a and &b must be input. They determine the
* type of line that is fit.
* (3) The median error can give an indication of how good the fit
* is likely to be.
*
*/
l_ok
ptaNoisyLinearLSF(PTA *pta,
l_float32 factor,
PTA **pptad,
l_float32 *pa,
l_float32 *pb,
l_float32 *pmederr,
NUMA **pnafit)
{
l_int32 n, i, ret;
l_float32 x, y, yf, val, mederr;
NUMA *nafit, *naerror;
PTA *ptad;
PROCNAME("ptaNoisyLinearLSF");
if (pptad) *pptad = NULL;
if (pa) *pa = 0.0;
if (pb) *pb = 0.0;
if (pmederr) *pmederr = 0.0;
if (pnafit) *pnafit = NULL;
if (!pptad && !pa && !pb && !pnafit)
return ERROR_INT("no output requested", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if (factor <= 0.0)
return ERROR_INT("factor must be > 0.0", procName, 1);
if ((n = ptaGetCount(pta)) < 3)
return ERROR_INT("less than 2 pts found", procName, 1);
if (ptaGetLinearLSF(pta, pa, pb, &nafit) != 0)
return ERROR_INT("error in linear LSF", procName, 1);
/* Get the median error */
naerror = numaCreate(n);
for (i = 0; i < n; i++) {
ptaGetPt(pta, i, &x, &y);
numaGetFValue(nafit, i, &yf);
numaAddNumber(naerror, L_ABS(y - yf));
}
numaGetMedian(naerror, &mederr);
if (pmederr) *pmederr = mederr;
numaDestroy(&nafit);
/* Remove outliers */
ptad = ptaCreate(n);
for (i = 0; i < n; i++) {
ptaGetPt(pta, i, &x, &y);
numaGetFValue(naerror, i, &val);
if (val <= factor * mederr) /* <= in case mederr = 0 */
ptaAddPt(ptad, x, y);
}
numaDestroy(&naerror);
/* Do LSF again */
ret = ptaGetLinearLSF(ptad, pa, pb, pnafit);
if (pptad)
*pptad = ptad;
else
ptaDestroy(&ptad);
return ret;
}
/*!
* \brief ptaNoisyQuadraticLSF()
*
* \param[in] pta
* \param[in] factor reject outliers with error greater than this
* number of medians; typically ~ 3
* \param[out] pptad [optional] with outliers removed
* \param[out] pa [optional] coeff a of LSF: y = ax^2 + bx + c
* \param[out] pb [optional] coeff b of LSF: y = ax^2 + bx + c
* \param[out] pc [optional] coeff c of LSF: y = ax^2 + bx + c
* \param[out] pmederr [optional] median error
* \param[out] pnafit [optional] numa of least square fit to ptad
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) This does a quadratic least square fit to the set of points
* in %pta. It then evaluates the errors and removes points
* whose error is >= factor * median_error. It then re-runs
* a quadratic LSF on the resulting points.
*
*/
l_ok
ptaNoisyQuadraticLSF(PTA *pta,
l_float32 factor,
PTA **pptad,
l_float32 *pa,
l_float32 *pb,
l_float32 *pc,
l_float32 *pmederr,
NUMA **pnafit)
{
l_int32 n, i, ret;
l_float32 x, y, yf, val, mederr;
NUMA *nafit, *naerror;
PTA *ptad;
PROCNAME("ptaNoisyQuadraticLSF");
if (pptad) *pptad = NULL;
if (pa) *pa = 0.0;
if (pb) *pb = 0.0;
if (pc) *pc = 0.0;
if (pmederr) *pmederr = 0.0;
if (pnafit) *pnafit = NULL;
if (!pptad && !pa && !pb && !pc && !pnafit)
return ERROR_INT("no output requested", procName, 1);
if (factor <= 0.0)
return ERROR_INT("factor must be > 0.0", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if ((n = ptaGetCount(pta)) < 3)
return ERROR_INT("less than 3 pts found", procName, 1);
if (ptaGetQuadraticLSF(pta, NULL, NULL, NULL, &nafit) != 0)
return ERROR_INT("error in quadratic LSF", procName, 1);
/* Get the median error */
naerror = numaCreate(n);
for (i = 0; i < n; i++) {
ptaGetPt(pta, i, &x, &y);
numaGetFValue(nafit, i, &yf);
numaAddNumber(naerror, L_ABS(y - yf));
}
numaGetMedian(naerror, &mederr);
if (pmederr) *pmederr = mederr;
numaDestroy(&nafit);
/* Remove outliers */
ptad = ptaCreate(n);
for (i = 0; i < n; i++) {
ptaGetPt(pta, i, &x, &y);
numaGetFValue(naerror, i, &val);
if (val <= factor * mederr) /* <= in case mederr = 0 */
ptaAddPt(ptad, x, y);
}
numaDestroy(&naerror);
n = ptaGetCount(ptad);
if ((n = ptaGetCount(ptad)) < 3) {
ptaDestroy(&ptad);
return ERROR_INT("less than 3 pts found", procName, 1);
}
/* Do LSF again */
ret = ptaGetQuadraticLSF(ptad, pa, pb, pc, pnafit);
if (pptad)
*pptad = ptad;
else
ptaDestroy(&ptad);
return ret;
}
/*!
* \brief applyLinearFit()
*
* \param[in] a, b linear fit coefficients
* \param[in] x
* \param[out] py y = a * x + b
* \return 0 if OK, 1 on error
*/
l_ok
applyLinearFit(l_float32 a,
l_float32 b,
l_float32 x,
l_float32 *py)
{
PROCNAME("applyLinearFit");
if (!py)
return ERROR_INT("&y not defined", procName, 1);
*py = a * x + b;
return 0;
}
/*!
* \brief applyQuadraticFit()
*
* \param[in] a, b, c quadratic fit coefficients
* \param[in] x
* \param[out] py y = a * x^2 + b * x + c
* \return 0 if OK, 1 on error
*/
l_ok
applyQuadraticFit(l_float32 a,
l_float32 b,
l_float32 c,
l_float32 x,
l_float32 *py)
{
PROCNAME("applyQuadraticFit");
if (!py)
return ERROR_INT("&y not defined", procName, 1);
*py = a * x * x + b * x + c;
return 0;
}
/*!
* \brief applyCubicFit()
*
* \param[in] a, b, c, d cubic fit coefficients
* \param[in] x
* \param[out] py y = a * x^3 + b * x^2 + c * x + d
* \return 0 if OK, 1 on error
*/
l_ok
applyCubicFit(l_float32 a,
l_float32 b,
l_float32 c,
l_float32 d,
l_float32 x,
l_float32 *py)
{
PROCNAME("applyCubicFit");
if (!py)
return ERROR_INT("&y not defined", procName, 1);
*py = a * x * x * x + b * x * x + c * x + d;
return 0;
}
/*!
* \brief applyQuarticFit()
*
* \param[in] a, b, c, d, e quartic fit coefficients
* \param[in] x
* \param[out] py y = a * x^4 + b * x^3 + c * x^2 + d * x + e
* \return 0 if OK, 1 on error
*/
l_ok
applyQuarticFit(l_float32 a,
l_float32 b,
l_float32 c,
l_float32 d,
l_float32 e,
l_float32 x,
l_float32 *py)
{
l_float32 x2;
PROCNAME("applyQuarticFit");
if (!py)
return ERROR_INT("&y not defined", procName, 1);
x2 = x * x;
*py = a * x2 * x2 + b * x2 * x + c * x2 + d * x + e;
return 0;
}
/*---------------------------------------------------------------------*
* Interconversions with Pix *
*---------------------------------------------------------------------*/
/*!
* \brief pixPlotAlongPta()
*
* \param[in] pixs any depth
* \param[in] pta set of points on which to plot
* \param[in] outformat GPLOT_PNG, GPLOT_PS, GPLOT_EPS, GPLOT_LATEX
* \param[in] title [optional] for plot; can be null
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) This is a debugging function.
* (2) Removes existing colormaps and clips the pta to the input %pixs.
* (3) If the image is RGB, three separate plots are generated.
*
*/
l_ok
pixPlotAlongPta(PIX *pixs,
PTA *pta,
l_int32 outformat,
const char *title)
{
char buffer[128];
char *rtitle, *gtitle, *btitle;
static l_int32 count = 0; /* require separate temp files for each call */
l_int32 i, x, y, d, w, h, npts, rval, gval, bval;
l_uint32 val;
NUMA *na, *nar, *nag, *nab;
PIX *pixt;
PROCNAME("pixPlotAlongPta");
lept_mkdir("lept/plot");
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
if (outformat != GPLOT_PNG && outformat != GPLOT_PS &&
outformat != GPLOT_EPS && outformat != GPLOT_LATEX) {
L_WARNING("outformat invalid; using GPLOT_PNG\n", procName);
outformat = GPLOT_PNG;
}
pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
d = pixGetDepth(pixt);
w = pixGetWidth(pixt);
h = pixGetHeight(pixt);
npts = ptaGetCount(pta);
if (d == 32) {
nar = numaCreate(npts);
nag = numaCreate(npts);
nab = numaCreate(npts);
for (i = 0; i < npts; i++) {
ptaGetIPt(pta, i, &x, &y);
if (x < 0 || x >= w)
continue;
if (y < 0 || y >= h)
continue;
pixGetPixel(pixt, x, y, &val);
rval = GET_DATA_BYTE(&val, COLOR_RED);
gval = GET_DATA_BYTE(&val, COLOR_GREEN);
bval = GET_DATA_BYTE(&val, COLOR_BLUE);
numaAddNumber(nar, rval);
numaAddNumber(nag, gval);
numaAddNumber(nab, bval);
}
snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
rtitle = stringJoin("Red: ", title);
gplotSimple1(nar, outformat, buffer, rtitle);
snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
gtitle = stringJoin("Green: ", title);
gplotSimple1(nag, outformat, buffer, gtitle);
snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
btitle = stringJoin("Blue: ", title);
gplotSimple1(nab, outformat, buffer, btitle);
numaDestroy(&nar);
numaDestroy(&nag);
numaDestroy(&nab);
LEPT_FREE(rtitle);
LEPT_FREE(gtitle);
LEPT_FREE(btitle);
} else {
na = numaCreate(npts);
for (i = 0; i < npts; i++) {
ptaGetIPt(pta, i, &x, &y);
if (x < 0 || x >= w)
continue;
if (y < 0 || y >= h)
continue;
pixGetPixel(pixt, x, y, &val);
numaAddNumber(na, (l_float32)val);
}
snprintf(buffer, sizeof(buffer), "/tmp/lept/plot/%03d", count++);
gplotSimple1(na, outformat, buffer, title);
numaDestroy(&na);
}
pixDestroy(&pixt);
return 0;
}
/*!
* \brief ptaGetPixelsFromPix()
*
* \param[in] pixs 1 bpp
* \param[in] box [optional] can be null
* \return pta, or NULL on error
*
*
* Notes:
* (1) Generates a pta of fg pixels in the pix, within the box.
* If box == NULL, it uses the entire pix.
*
*/
PTA *
ptaGetPixelsFromPix(PIX *pixs,
BOX *box)
{
l_int32 i, j, w, h, wpl, xstart, xend, ystart, yend, bw, bh;
l_uint32 *data, *line;
PTA *pta;
PROCNAME("ptaGetPixelsFromPix");
if (!pixs || (pixGetDepth(pixs) != 1))
return (PTA *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
pixGetDimensions(pixs, &w, &h, NULL);
data = pixGetData(pixs);
wpl = pixGetWpl(pixs);
xstart = ystart = 0;
xend = w - 1;
yend = h - 1;
if (box) {
boxGetGeometry(box, &xstart, &ystart, &bw, &bh);
xend = xstart + bw - 1;
yend = ystart + bh - 1;
}
if ((pta = ptaCreate(0)) == NULL)
return (PTA *)ERROR_PTR("pta not made", procName, NULL);
for (i = ystart; i <= yend; i++) {
line = data + i * wpl;
for (j = xstart; j <= xend; j++) {
if (GET_DATA_BIT(line, j))
ptaAddPt(pta, j, i);
}
}
return pta;
}
/*!
* \brief pixGenerateFromPta()
*
* \param[in] pta
* \param[in] w, h of pix
* \return pix 1 bpp, or NULL on error
*
*
* Notes:
* (1) Points are rounded to nearest ints.
* (2) Any points outside (w,h) are silently discarded.
* (3) Output 1 bpp pix has values 1 for each point in the pta.
*
*/
PIX *
pixGenerateFromPta(PTA *pta,
l_int32 w,
l_int32 h)
{
l_int32 n, i, x, y;
PIX *pix;
PROCNAME("pixGenerateFromPta");
if (!pta)
return (PIX *)ERROR_PTR("pta not defined", procName, NULL);
if ((pix = pixCreate(w, h, 1)) == NULL)
return (PIX *)ERROR_PTR("pix not made", procName, NULL);
n = ptaGetCount(pta);
for (i = 0; i < n; i++) {
ptaGetIPt(pta, i, &x, &y);
if (x < 0 || x >= w || y < 0 || y >= h)
continue;
pixSetPixel(pix, x, y, 1);
}
return pix;
}
/*!
* \brief ptaGetBoundaryPixels()
*
* \param[in] pixs 1 bpp
* \param[in] type L_BOUNDARY_FG, L_BOUNDARY_BG
* \return pta, or NULL on error
*
*
* Notes:
* (1) This generates a pta of either fg or bg boundary pixels.
* (2) See also pixGeneratePtaBoundary() for rendering of
* fg boundary pixels.
*
*/
PTA *
ptaGetBoundaryPixels(PIX *pixs,
l_int32 type)
{
PIX *pixt;
PTA *pta;
PROCNAME("ptaGetBoundaryPixels");
if (!pixs || (pixGetDepth(pixs) != 1))
return (PTA *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
if (type != L_BOUNDARY_FG && type != L_BOUNDARY_BG)
return (PTA *)ERROR_PTR("invalid type", procName, NULL);
if (type == L_BOUNDARY_FG)
pixt = pixMorphSequence(pixs, "e3.3", 0);
else
pixt = pixMorphSequence(pixs, "d3.3", 0);
pixXor(pixt, pixt, pixs);
pta = ptaGetPixelsFromPix(pixt, NULL);
pixDestroy(&pixt);
return pta;
}
/*!
* \brief ptaaGetBoundaryPixels()
*
* \param[in] pixs 1 bpp
* \param[in] type L_BOUNDARY_FG, L_BOUNDARY_BG
* \param[in] connectivity 4 or 8
* \param[out] pboxa [optional] bounding boxes of the c.c.
* \param[out] ppixa [optional] pixa of the c.c.
* \return ptaa, or NULL on error
*
*
* Notes:
* (1) This generates a ptaa of either fg or bg boundary pixels,
* where each pta has the boundary pixels for a connected
* component.
* (2) We can't simply find all the boundary pixels and then select
* those within the bounding box of each component, because
* bounding boxes can overlap. It is necessary to extract and
* dilate or erode each component separately. Note also that
* special handling is required for bg pixels when the
* component touches the pix boundary.
*
*/
PTAA *
ptaaGetBoundaryPixels(PIX *pixs,
l_int32 type,
l_int32 connectivity,
BOXA **pboxa,
PIXA **ppixa)
{
l_int32 i, n, w, h, x, y, bw, bh, left, right, top, bot;
BOXA *boxa;
PIX *pixt1, *pixt2;
PIXA *pixa;
PTA *pta1, *pta2;
PTAA *ptaa;
PROCNAME("ptaaGetBoundaryPixels");
if (pboxa) *pboxa = NULL;
if (ppixa) *ppixa = NULL;
if (!pixs || (pixGetDepth(pixs) != 1))
return (PTAA *)ERROR_PTR("pixs undefined or not 1 bpp", procName, NULL);
if (type != L_BOUNDARY_FG && type != L_BOUNDARY_BG)
return (PTAA *)ERROR_PTR("invalid type", procName, NULL);
if (connectivity != 4 && connectivity != 8)
return (PTAA *)ERROR_PTR("connectivity not 4 or 8", procName, NULL);
pixGetDimensions(pixs, &w, &h, NULL);
boxa = pixConnComp(pixs, &pixa, connectivity);
n = boxaGetCount(boxa);
ptaa = ptaaCreate(0);
for (i = 0; i < n; i++) {
pixt1 = pixaGetPix(pixa, i, L_CLONE);
boxaGetBoxGeometry(boxa, i, &x, &y, &bw, &bh);
left = right = top = bot = 0;
if (type == L_BOUNDARY_BG) {
if (x > 0) left = 1;
if (y > 0) top = 1;
if (x + bw < w) right = 1;
if (y + bh < h) bot = 1;
pixt2 = pixAddBorderGeneral(pixt1, left, right, top, bot, 0);
} else {
pixt2 = pixClone(pixt1);
}
pta1 = ptaGetBoundaryPixels(pixt2, type);
pta2 = ptaTransform(pta1, x - left, y - top, 1.0, 1.0);
ptaaAddPta(ptaa, pta2, L_INSERT);
ptaDestroy(&pta1);
pixDestroy(&pixt1);
pixDestroy(&pixt2);
}
if (pboxa)
*pboxa = boxa;
else
boxaDestroy(&boxa);
if (ppixa)
*ppixa = pixa;
else
pixaDestroy(&pixa);
return ptaa;
}
/*!
* \brief ptaaIndexLabeledPixels()
*
* \param[in] pixs 32 bpp, of indices of c.c.
* \param[out] pncc [optional] number of connected components
* \return ptaa, or NULL on error
*
*
* Notes:
* (1) The pixel values in %pixs are the index of the connected component
* to which the pixel belongs; %pixs is typically generated from
* a 1 bpp pix by pixConnCompTransform(). Background pixels in
* the generating 1 bpp pix are represented in %pixs by 0.
* We do not check that the pixel values are correctly labelled.
* (2) Each pta in the returned ptaa gives the pixel locations
* correspnding to a connected component, with the label of each
* given by the index of the pta into the ptaa.
* (3) Initialize with the first pta in ptaa being empty and
* representing the background value (index 0) in the pix.
*
*/
PTAA *
ptaaIndexLabeledPixels(PIX *pixs,
l_int32 *pncc)
{
l_int32 wpl, index, i, j, w, h;
l_uint32 maxval;
l_uint32 *data, *line;
PTA *pta;
PTAA *ptaa;
PROCNAME("ptaaIndexLabeledPixels");
if (pncc) *pncc = 0;
if (!pixs || (pixGetDepth(pixs) != 32))
return (PTAA *)ERROR_PTR("pixs undef or not 32 bpp", procName, NULL);
/* The number of c.c. is the maximum pixel value. Use this to
* initialize ptaa with sufficient pta arrays */
pixGetMaxValueInRect(pixs, NULL, &maxval, NULL, NULL);
if (pncc) *pncc = maxval;
pta = ptaCreate(1);
ptaa = ptaaCreate(maxval + 1);
ptaaInitFull(ptaa, pta);
ptaDestroy(&pta);
/* Sweep over %pixs, saving the pixel coordinates of each pixel
* with nonzero value in the appropriate pta, indexed by that value. */
pixGetDimensions(pixs, &w, &h, NULL);
data = pixGetData(pixs);
wpl = pixGetWpl(pixs);
for (i = 0; i < h; i++) {
line = data + wpl * i;
for (j = 0; j < w; j++) {
index = line[j];
if (index > 0)
ptaaAddPt(ptaa, index, j, i);
}
}
return ptaa;
}
/*!
* \brief ptaGetNeighborPixLocs()
*
* \param[in] pixs any depth
* \param[in] x, y pixel from which we search for nearest neighbors
* \param[in] conn 4 or 8 connectivity
* \return pta, or NULL on error
*
*
* Notes:
* (1) Generates a pta of all valid neighbor pixel locations,
* or NULL on error.
*
*/
PTA *
ptaGetNeighborPixLocs(PIX *pixs,
l_int32 x,
l_int32 y,
l_int32 conn)
{
l_int32 w, h;
PTA *pta;
PROCNAME("ptaGetNeighborPixLocs");
if (!pixs)
return (PTA *)ERROR_PTR("pixs not defined", procName, NULL);
pixGetDimensions(pixs, &w, &h, NULL);
if (x < 0 || x >= w || y < 0 || y >= h)
return (PTA *)ERROR_PTR("(x,y) not in pixs", procName, NULL);
if (conn != 4 && conn != 8)
return (PTA *)ERROR_PTR("conn not 4 or 8", procName, NULL);
pta = ptaCreate(conn);
if (x > 0)
ptaAddPt(pta, x - 1, y);
if (x < w - 1)
ptaAddPt(pta, x + 1, y);
if (y > 0)
ptaAddPt(pta, x, y - 1);
if (y < h - 1)
ptaAddPt(pta, x, y + 1);
if (conn == 8) {
if (x > 0) {
if (y > 0)
ptaAddPt(pta, x - 1, y - 1);
if (y < h - 1)
ptaAddPt(pta, x - 1, y + 1);
}
if (x < w - 1) {
if (y > 0)
ptaAddPt(pta, x + 1, y - 1);
if (y < h - 1)
ptaAddPt(pta, x + 1, y + 1);
}
}
return pta;
}
/*---------------------------------------------------------------------*
* Interconversion with Numa *
*---------------------------------------------------------------------*/
/*!
* \brief numaConvertToPta1()
*
* \param[in] na numa with implicit y(x)
* \return pta if OK; null on error
*/
PTA *
numaConvertToPta1(NUMA *na)
{
l_int32 i, n;
l_float32 startx, delx, val;
PTA *pta;
PROCNAME("numaConvertToPta1");
if (!na)
return (PTA *)ERROR_PTR("na not defined", procName, NULL);
n = numaGetCount(na);
pta = ptaCreate(n);
numaGetParameters(na, &startx, &delx);
for (i = 0; i < n; i++) {
numaGetFValue(na, i, &val);
ptaAddPt(pta, startx + i * delx, val);
}
return pta;
}
/*!
* \brief numaConvertToPta2()
*
* \param[in] nax
* \param[in] nay
* \return pta if OK; null on error
*/
PTA *
numaConvertToPta2(NUMA *nax,
NUMA *nay)
{
l_int32 i, n, nx, ny;
l_float32 valx, valy;
PTA *pta;
PROCNAME("numaConvertToPta2");
if (!nax || !nay)
return (PTA *)ERROR_PTR("nax and nay not both defined", procName, NULL);
nx = numaGetCount(nax);
ny = numaGetCount(nay);
n = L_MIN(nx, ny);
if (nx != ny)
L_WARNING("nx = %d does not equal ny = %d\n", procName, nx, ny);
pta = ptaCreate(n);
for (i = 0; i < n; i++) {
numaGetFValue(nax, i, &valx);
numaGetFValue(nay, i, &valy);
ptaAddPt(pta, valx, valy);
}
return pta;
}
/*!
* \brief ptaConvertToNuma()
*
* \param[in] pta
* \param[out] pnax addr of nax
* \param[out] pnay addr of nay
* \return 0 if OK, 1 on error
*/
l_ok
ptaConvertToNuma(PTA *pta,
NUMA **pnax,
NUMA **pnay)
{
l_int32 i, n;
l_float32 valx, valy;
PROCNAME("ptaConvertToNuma");
if (pnax) *pnax = NULL;
if (pnay) *pnay = NULL;
if (!pnax || !pnay)
return ERROR_INT("&nax and &nay not both defined", procName, 1);
if (!pta)
return ERROR_INT("pta not defined", procName, 1);
n = ptaGetCount(pta);
*pnax = numaCreate(n);
*pnay = numaCreate(n);
for (i = 0; i < n; i++) {
ptaGetPt(pta, i, &valx, &valy);
numaAddNumber(*pnax, valx);
numaAddNumber(*pnay, valy);
}
return 0;
}
/*---------------------------------------------------------------------*
* Display Pta and Ptaa *
*---------------------------------------------------------------------*/
/*!
* \brief pixDisplayPta()
*
* \param[in] pixd can be same as pixs or NULL; 32 bpp if in-place
* \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp
* \param[in] pta of path to be plotted
* \return pixd 32 bpp RGB version of pixs, with path in green.
*
*
* Notes:
* (1) To write on an existing pixs, pixs must be 32 bpp and
* call with pixd == pixs:
* pixDisplayPta(pixs, pixs, pta);
* To write to a new pix, use pixd == NULL and call:
* pixd = pixDisplayPta(NULL, pixs, pta);
* (2) On error, returns pixd to avoid losing pixs if called as
* pixs = pixDisplayPta(pixs, pixs, pta);
*
*/
PIX *
pixDisplayPta(PIX *pixd,
PIX *pixs,
PTA *pta)
{
l_int32 i, n, w, h, x, y;
l_uint32 rpixel, gpixel, bpixel;
PROCNAME("pixDisplayPta");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
if (!pta)
return (PIX *)ERROR_PTR("pta not defined", procName, pixd);
if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
return (PIX *)ERROR_PTR("invalid pixd", procName, pixd);
if (!pixd)
pixd = pixConvertTo32(pixs);
pixGetDimensions(pixd, &w, &h, NULL);
composeRGBPixel(255, 0, 0, &rpixel); /* start point */
composeRGBPixel(0, 255, 0, &gpixel);
composeRGBPixel(0, 0, 255, &bpixel); /* end point */
n = ptaGetCount(pta);
for (i = 0; i < n; i++) {
ptaGetIPt(pta, i, &x, &y);
if (x < 0 || x >= w || y < 0 || y >= h)
continue;
if (i == 0)
pixSetPixel(pixd, x, y, rpixel);
else if (i < n - 1)
pixSetPixel(pixd, x, y, gpixel);
else
pixSetPixel(pixd, x, y, bpixel);
}
return pixd;
}
/*!
* \brief pixDisplayPtaaPattern()
*
* \param[in] pixd 32 bpp
* \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp; 32 bpp if in place
* \param[in] ptaa giving locations at which the pattern is displayed
* \param[in] pixp 1 bpp pattern to be placed such that its reference
* point co-locates with each point in pta
* \param[in] cx, cy reference point in pattern
* \return pixd 32 bpp RGB version of pixs.
*
*
* Notes:
* (1) To write on an existing pixs, pixs must be 32 bpp and
* call with pixd == pixs:
* pixDisplayPtaPattern(pixs, pixs, pta, ...);
* To write to a new pix, use pixd == NULL and call:
* pixd = pixDisplayPtaPattern(NULL, pixs, pta, ...);
* (2) Puts a random color on each pattern associated with a pta.
* (3) On error, returns pixd to avoid losing pixs if called as
* pixs = pixDisplayPtaPattern(pixs, pixs, pta, ...);
* (4) A typical pattern to be used is a circle, generated with
* generatePtaFilledCircle()
*
*/
PIX *
pixDisplayPtaaPattern(PIX *pixd,
PIX *pixs,
PTAA *ptaa,
PIX *pixp,
l_int32 cx,
l_int32 cy)
{
l_int32 i, n;
l_uint32 color;
PIXCMAP *cmap;
PTA *pta;
PROCNAME("pixDisplayPtaaPattern");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
if (!ptaa)
return (PIX *)ERROR_PTR("ptaa not defined", procName, pixd);
if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
return (PIX *)ERROR_PTR("invalid pixd", procName, pixd);
if (!pixp)
return (PIX *)ERROR_PTR("pixp not defined", procName, pixd);
if (!pixd)
pixd = pixConvertTo32(pixs);
/* Use 256 random colors */
cmap = pixcmapCreateRandom(8, 0, 0);
n = ptaaGetCount(ptaa);
for (i = 0; i < n; i++) {
pixcmapGetColor32(cmap, i % 256, &color);
pta = ptaaGetPta(ptaa, i, L_CLONE);
pixDisplayPtaPattern(pixd, pixd, pta, pixp, cx, cy, color);
ptaDestroy(&pta);
}
pixcmapDestroy(&cmap);
return pixd;
}
/*!
* \brief pixDisplayPtaPattern()
*
* \param[in] pixd can be same as pixs or NULL; 32 bpp if in-place
* \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp
* \param[in] pta giving locations at which the pattern is displayed
* \param[in] pixp 1 bpp pattern to be placed such that its reference
* point co-locates with each point in pta
* \param[in] cx, cy reference point in pattern
* \param[in] color in 0xrrggbb00 format
* \return pixd 32 bpp RGB version of pixs.
*
*
* Notes:
* (1) To write on an existing pixs, pixs must be 32 bpp and
* call with pixd == pixs:
* pixDisplayPtaPattern(pixs, pixs, pta, ...);
* To write to a new pix, use pixd == NULL and call:
* pixd = pixDisplayPtaPattern(NULL, pixs, pta, ...);
* (2) On error, returns pixd to avoid losing pixs if called as
* pixs = pixDisplayPtaPattern(pixs, pixs, pta, ...);
* (3) A typical pattern to be used is a circle, generated with
* generatePtaFilledCircle()
*
*/
PIX *
pixDisplayPtaPattern(PIX *pixd,
PIX *pixs,
PTA *pta,
PIX *pixp,
l_int32 cx,
l_int32 cy,
l_uint32 color)
{
l_int32 i, n, w, h, x, y;
PTA *ptat;
PROCNAME("pixDisplayPtaPattern");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
if (!pta)
return (PIX *)ERROR_PTR("pta not defined", procName, pixd);
if (pixd && (pixd != pixs || pixGetDepth(pixd) != 32))
return (PIX *)ERROR_PTR("invalid pixd", procName, pixd);
if (!pixp)
return (PIX *)ERROR_PTR("pixp not defined", procName, pixd);
if (!pixd)
pixd = pixConvertTo32(pixs);
pixGetDimensions(pixs, &w, &h, NULL);
ptat = ptaReplicatePattern(pta, pixp, NULL, cx, cy, w, h);
n = ptaGetCount(ptat);
for (i = 0; i < n; i++) {
ptaGetIPt(ptat, i, &x, &y);
if (x < 0 || x >= w || y < 0 || y >= h)
continue;
pixSetPixel(pixd, x, y, color);
}
ptaDestroy(&ptat);
return pixd;
}
/*!
* \brief ptaReplicatePattern()
*
* \param[in] ptas "sparse" input pta
* \param[in] pixp [optional] 1 bpp pattern, to be replicated
* in output pta
* \param[in] ptap [optional] set of pts, to be replicated in output pta
* \param[in] cx, cy reference point in pattern
* \param[in] w, h clipping sizes for output pta
* \return ptad with all points of replicated pattern, or NULL on error
*
*
* Notes:
* (1) You can use either the image %pixp or the set of pts %ptap.
* (2) The pattern is placed with its reference point at each point
* in ptas, and all the fg pixels are colleced into ptad.
* For %pixp, this is equivalent to blitting pixp at each point
* in ptas, and then converting the resulting pix to a pta.
*
*/
PTA *
ptaReplicatePattern(PTA *ptas,
PIX *pixp,
PTA *ptap,
l_int32 cx,
l_int32 cy,
l_int32 w,
l_int32 h)
{
l_int32 i, j, n, np, x, y, xp, yp, xf, yf;
PTA *ptat, *ptad;
PROCNAME("ptaReplicatePattern");
if (!ptas)
return (PTA *)ERROR_PTR("ptas not defined", procName, NULL);
if (!pixp && !ptap)
return (PTA *)ERROR_PTR("no pattern is defined", procName, NULL);
if (pixp && ptap)
L_WARNING("pixp and ptap defined; using ptap\n", procName);
n = ptaGetCount(ptas);
ptad = ptaCreate(n);
if (ptap)
ptat = ptaClone(ptap);
else
ptat = ptaGetPixelsFromPix(pixp, NULL);
np = ptaGetCount(ptat);
for (i = 0; i < n; i++) {
ptaGetIPt(ptas, i, &x, &y);
for (j = 0; j < np; j++) {
ptaGetIPt(ptat, j, &xp, &yp);
xf = x - cx + xp;
yf = y - cy + yp;
if (xf >= 0 && xf < w && yf >= 0 && yf < h)
ptaAddPt(ptad, xf, yf);
}
}
ptaDestroy(&ptat);
return ptad;
}
/*!
* \brief pixDisplayPtaa()
*
* \param[in] pixs 1, 2, 4, 8, 16 or 32 bpp
* \param[in] ptaa array of paths to be plotted
* \return pixd 32 bpp RGB version of pixs, with paths plotted
* in different colors, or NULL on error
*/
PIX *
pixDisplayPtaa(PIX *pixs,
PTAA *ptaa)
{
l_int32 i, j, w, h, npta, npt, x, y, rv, gv, bv;
l_uint32 *pixela;
NUMA *na1, *na2, *na3;
PIX *pixd;
PTA *pta;
PROCNAME("pixDisplayPtaa");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (!ptaa)
return (PIX *)ERROR_PTR("ptaa not defined", procName, NULL);
npta = ptaaGetCount(ptaa);
if (npta == 0)
return (PIX *)ERROR_PTR("no pta", procName, NULL);
if ((pixd = pixConvertTo32(pixs)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
pixGetDimensions(pixd, &w, &h, NULL);
/* Make a colormap for the paths */
if ((pixela = (l_uint32 *)LEPT_CALLOC(npta, sizeof(l_uint32))) == NULL) {
pixDestroy(&pixd);
return (PIX *)ERROR_PTR("calloc fail for pixela", procName, NULL);
}
na1 = numaPseudorandomSequence(256, 14657);
na2 = numaPseudorandomSequence(256, 34631);
na3 = numaPseudorandomSequence(256, 54617);
for (i = 0; i < npta; i++) {
numaGetIValue(na1, i % 256, &rv);
numaGetIValue(na2, i % 256, &gv);
numaGetIValue(na3, i % 256, &bv);
composeRGBPixel(rv, gv, bv, &pixela[i]);
}
numaDestroy(&na1);
numaDestroy(&na2);
numaDestroy(&na3);
for (i = 0; i < npta; i++) {
pta = ptaaGetPta(ptaa, i, L_CLONE);
npt = ptaGetCount(pta);
for (j = 0; j < npt; j++) {
ptaGetIPt(pta, j, &x, &y);
if (x < 0 || x >= w || y < 0 || y >= h)
continue;
pixSetPixel(pixd, x, y, pixela[i]);
}
ptaDestroy(&pta);
}
LEPT_FREE(pixela);
return pixd;
}