#include "allheaders.h"
/* Array size must be > 0 and not larger than this */
static const l_uint32 MaxArraySize = 100000;
/*------------------------------------------------------------------------*
* Create / Destroy *
*------------------------------------------------------------------------*/
/*!
* \brief kernelCreate()
*
* \param[in] height, width
* \return kernel, or NULL on error
*
*
* Notes:
* (1) kernelCreate() initializes all values to 0.
* (2) After this call, (cy,cx) and nonzero data values must be
* assigned.
* (2) The number of kernel elements must be less than 2^29.
*
*/
L_KERNEL *
kernelCreate(l_int32 height,
l_int32 width)
{
l_uint64 size64;
L_KERNEL *kel;
PROCNAME("kernelCreate");
if (width <= 0)
return (L_KERNEL *)ERROR_PTR("width must be > 0", procName, NULL);
if (height <= 0)
return (L_KERNEL *)ERROR_PTR("height must be > 0", procName, NULL);
/* Avoid overflow in malloc arg */
size64 = (l_uint64)width * (l_uint64)height;
if (size64 >= (1LL << 29)) {
L_ERROR("requested width = %d, height = %d\n", procName, width, height);
return (L_KERNEL *)ERROR_PTR("size >= 2^29", procName, NULL);
}
kel = (L_KERNEL *)LEPT_CALLOC(1, sizeof(L_KERNEL));
kel->sy = height;
kel->sx = width;
if ((kel->data = create2dFloatArray(height, width)) == NULL) {
LEPT_FREE(kel);
return (L_KERNEL *)ERROR_PTR("data not allocated", procName, NULL);
}
return kel;
}
/*!
* \brief kernelDestroy()
*
* \param[in,out] pkel will be set to null before returning
* \return void
*/
void
kernelDestroy(L_KERNEL **pkel)
{
l_int32 i;
L_KERNEL *kel;
PROCNAME("kernelDestroy");
if (pkel == NULL) {
L_WARNING("ptr address is NULL!\n", procName);
return;
}
if ((kel = *pkel) == NULL)
return;
for (i = 0; i < kel->sy; i++)
LEPT_FREE(kel->data[i]);
LEPT_FREE(kel->data);
LEPT_FREE(kel);
*pkel = NULL;
}
/*!
* \brief kernelCopy()
*
* \param[in] kels source kernel
* \return keld copy of kels, or NULL on error
*/
L_KERNEL *
kernelCopy(L_KERNEL *kels)
{
l_int32 i, j, sx, sy, cx, cy;
L_KERNEL *keld;
PROCNAME("kernelCopy");
if (!kels)
return (L_KERNEL *)ERROR_PTR("kels not defined", procName, NULL);
kernelGetParameters(kels, &sy, &sx, &cy, &cx);
if ((keld = kernelCreate(sy, sx)) == NULL)
return (L_KERNEL *)ERROR_PTR("keld not made", procName, NULL);
keld->cy = cy;
keld->cx = cx;
for (i = 0; i < sy; i++)
for (j = 0; j < sx; j++)
keld->data[i][j] = kels->data[i][j];
return keld;
}
/*----------------------------------------------------------------------*
* Accessors *
*----------------------------------------------------------------------*/
/*!
* \brief kernelGetElement()
*
* \param[in] kel
* \param[in] row
* \param[in] col
* \param[out] pval
* \return 0 if OK; 1 on error
*/
l_ok
kernelGetElement(L_KERNEL *kel,
l_int32 row,
l_int32 col,
l_float32 *pval)
{
PROCNAME("kernelGetElement");
if (!pval)
return ERROR_INT("&val not defined", procName, 1);
*pval = 0;
if (!kel)
return ERROR_INT("kernel not defined", procName, 1);
if (row < 0 || row >= kel->sy)
return ERROR_INT("kernel row out of bounds", procName, 1);
if (col < 0 || col >= kel->sx)
return ERROR_INT("kernel col out of bounds", procName, 1);
*pval = kel->data[row][col];
return 0;
}
/*!
* \brief kernelSetElement()
*
* \param[in] kel kernel
* \param[in] row
* \param[in] col
* \param[in] val
* \return 0 if OK; 1 on error
*/
l_ok
kernelSetElement(L_KERNEL *kel,
l_int32 row,
l_int32 col,
l_float32 val)
{
PROCNAME("kernelSetElement");
if (!kel)
return ERROR_INT("kel not defined", procName, 1);
if (row < 0 || row >= kel->sy)
return ERROR_INT("kernel row out of bounds", procName, 1);
if (col < 0 || col >= kel->sx)
return ERROR_INT("kernel col out of bounds", procName, 1);
kel->data[row][col] = val;
return 0;
}
/*!
* \brief kernelGetParameters()
*
* \param[in] kel kernel
* \param[out] psy, psx, pcy, pcx [optional] each can be null
* \return 0 if OK, 1 on error
*/
l_ok
kernelGetParameters(L_KERNEL *kel,
l_int32 *psy,
l_int32 *psx,
l_int32 *pcy,
l_int32 *pcx)
{
PROCNAME("kernelGetParameters");
if (psy) *psy = 0;
if (psx) *psx = 0;
if (pcy) *pcy = 0;
if (pcx) *pcx = 0;
if (!kel)
return ERROR_INT("kernel not defined", procName, 1);
if (psy) *psy = kel->sy;
if (psx) *psx = kel->sx;
if (pcy) *pcy = kel->cy;
if (pcx) *pcx = kel->cx;
return 0;
}
/*!
* \brief kernelSetOrigin()
*
* \param[in] kel kernel
* \param[in] cy, cx
* \return 0 if OK; 1 on error
*/
l_ok
kernelSetOrigin(L_KERNEL *kel,
l_int32 cy,
l_int32 cx)
{
PROCNAME("kernelSetOrigin");
if (!kel)
return ERROR_INT("kel not defined", procName, 1);
kel->cy = cy;
kel->cx = cx;
return 0;
}
/*!
* \brief kernelGetSum()
*
* \param[in] kel kernel
* \param[out] psum sum of all kernel values
* \return 0 if OK, 1 on error
*/
l_ok
kernelGetSum(L_KERNEL *kel,
l_float32 *psum)
{
l_int32 sx, sy, i, j;
PROCNAME("kernelGetSum");
if (!psum)
return ERROR_INT("&sum not defined", procName, 1);
*psum = 0.0;
if (!kel)
return ERROR_INT("kernel not defined", procName, 1);
kernelGetParameters(kel, &sy, &sx, NULL, NULL);
for (i = 0; i < sy; i++) {
for (j = 0; j < sx; j++) {
*psum += kel->data[i][j];
}
}
return 0;
}
/*!
* \brief kernelGetMinMax()
*
* \param[in] kel kernel
* \param[out] pmin [optional] minimum value
* \param[out] pmax [optional] maximum value
* \return 0 if OK, 1 on error
*/
l_ok
kernelGetMinMax(L_KERNEL *kel,
l_float32 *pmin,
l_float32 *pmax)
{
l_int32 sx, sy, i, j;
l_float32 val, minval, maxval;
PROCNAME("kernelGetMinmax");
if (!pmin && !pmax)
return ERROR_INT("neither &min nor &max defined", procName, 1);
if (pmin) *pmin = 0.0;
if (pmax) *pmax = 0.0;
if (!kel)
return ERROR_INT("kernel not defined", procName, 1);
kernelGetParameters(kel, &sy, &sx, NULL, NULL);
minval = 10000000.0;
maxval = -10000000.0;
for (i = 0; i < sy; i++) {
for (j = 0; j < sx; j++) {
val = kel->data[i][j];
if (val < minval)
minval = val;
if (val > maxval)
maxval = val;
}
}
if (pmin)
*pmin = minval;
if (pmax)
*pmax = maxval;
return 0;
}
/*----------------------------------------------------------------------*
* Normalize/Invert *
*----------------------------------------------------------------------*/
/*!
* \brief kernelNormalize()
*
* \param[in] kels source kel, to be normalized
* \param[in] normsum desired sum of elements in keld
* \return keld normalized version of kels, or NULL on error
* or if sum of elements is very close to 0)
*
*
* Notes:
* (1) If the sum of kernel elements is close to 0, do not
* try to calculate the normalized kernel. Instead,
* return a copy of the input kernel, with a warning.
*
*/
L_KERNEL *
kernelNormalize(L_KERNEL *kels,
l_float32 normsum)
{
l_int32 i, j, sx, sy, cx, cy;
l_float32 sum, factor;
L_KERNEL *keld;
PROCNAME("kernelNormalize");
if (!kels)
return (L_KERNEL *)ERROR_PTR("kels not defined", procName, NULL);
kernelGetSum(kels, &sum);
if (L_ABS(sum) < 0.00001) {
L_WARNING("null sum; not normalizing; returning a copy\n", procName);
return kernelCopy(kels);
}
kernelGetParameters(kels, &sy, &sx, &cy, &cx);
if ((keld = kernelCreate(sy, sx)) == NULL)
return (L_KERNEL *)ERROR_PTR("keld not made", procName, NULL);
keld->cy = cy;
keld->cx = cx;
factor = normsum / sum;
for (i = 0; i < sy; i++)
for (j = 0; j < sx; j++)
keld->data[i][j] = factor * kels->data[i][j];
return keld;
}
/*!
* \brief kernelInvert()
*
* \param[in] kels source kel, to be inverted
* \return keld spatially inverted, about the origin, or NULL on error
*
*
* Notes:
* (1) For convolution, the kernel is spatially inverted before
* a "correlation" operation is done between the kernel and the image.
*
*/
L_KERNEL *
kernelInvert(L_KERNEL *kels)
{
l_int32 i, j, sx, sy, cx, cy;
L_KERNEL *keld;
PROCNAME("kernelInvert");
if (!kels)
return (L_KERNEL *)ERROR_PTR("kels not defined", procName, NULL);
kernelGetParameters(kels, &sy, &sx, &cy, &cx);
if ((keld = kernelCreate(sy, sx)) == NULL)
return (L_KERNEL *)ERROR_PTR("keld not made", procName, NULL);
keld->cy = sy - 1 - cy;
keld->cx = sx - 1 - cx;
for (i = 0; i < sy; i++)
for (j = 0; j < sx; j++)
keld->data[i][j] = kels->data[sy - 1 - i][sx - 1 - j];
return keld;
}
/*----------------------------------------------------------------------*
* Helper function *
*----------------------------------------------------------------------*/
/*!
* \brief create2dFloatArray()
*
* \param[in] sy rows == height
* \param[in] sx columns == width
* \return doubly indexed array i.e., an array of sy row pointers,
* each of which points to an array of sx floats
*
*
* Notes:
* (1) The array[%sy][%sx] is indexed in standard "matrix notation",
* with the row index first.
* (2) The caller kernelCreate() limits the size to < 2^29 pixels.
*
*/
l_float32 **
create2dFloatArray(l_int32 sy,
l_int32 sx)
{
l_int32 i;
l_float32 **array;
PROCNAME("create2dFloatArray");
if (sx <= 0 || sx > MaxArraySize)
return (l_float32 **)ERROR_PTR("sx out of bounds", procName, NULL);
if (sy <= 0 || sy > MaxArraySize)
return (l_float32 **)ERROR_PTR("sy out of bounds", procName, NULL);
if ((array = (l_float32 **)LEPT_CALLOC(sy, sizeof(l_float32 *))) == NULL)
return (l_float32 **)ERROR_PTR("ptr array not made", procName, NULL);
for (i = 0; i < sy; i++)
array[i] = (l_float32 *)LEPT_CALLOC(sx, sizeof(l_float32));
return array;
}
/*----------------------------------------------------------------------*
* Kernel serialized I/O *
*----------------------------------------------------------------------*/
/*!
* \brief kernelRead()
*
* \param[in] fname filename
* \return kernel, or NULL on error
*/
L_KERNEL *
kernelRead(const char *fname)
{
FILE *fp;
L_KERNEL *kel;
PROCNAME("kernelRead");
if (!fname)
return (L_KERNEL *)ERROR_PTR("fname not defined", procName, NULL);
if ((fp = fopenReadStream(fname)) == NULL)
return (L_KERNEL *)ERROR_PTR("stream not opened", procName, NULL);
if ((kel = kernelReadStream(fp)) == NULL) {
fclose(fp);
return (L_KERNEL *)ERROR_PTR("kel not returned", procName, NULL);
}
fclose(fp);
return kel;
}
/*!
* \brief kernelReadStream()
*
* \param[in] fp file stream
* \return kernel, or NULL on error
*/
L_KERNEL *
kernelReadStream(FILE *fp)
{
l_int32 sy, sx, cy, cx, i, j, ret, version, ignore;
L_KERNEL *kel;
PROCNAME("kernelReadStream");
if (!fp)
return (L_KERNEL *)ERROR_PTR("stream not defined", procName, NULL);
ret = fscanf(fp, " Kernel Version %d\n", &version);
if (ret != 1)
return (L_KERNEL *)ERROR_PTR("not a kernel file", procName, NULL);
if (version != KERNEL_VERSION_NUMBER)
return (L_KERNEL *)ERROR_PTR("invalid kernel version", procName, NULL);
if (fscanf(fp, " sy = %d, sx = %d, cy = %d, cx = %d\n",
&sy, &sx, &cy, &cx) != 4)
return (L_KERNEL *)ERROR_PTR("dimensions not read", procName, NULL);
if (sx > MaxArraySize || sy > MaxArraySize) {
L_ERROR("sx = %d or sy = %d > %d\n", procName, sx, sy, MaxArraySize);
return NULL;
}
if ((kel = kernelCreate(sy, sx)) == NULL)
return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
kernelSetOrigin(kel, cy, cx);
for (i = 0; i < sy; i++) {
for (j = 0; j < sx; j++)
ignore = fscanf(fp, "%15f", &kel->data[i][j]);
ignore = fscanf(fp, "\n");
}
ignore = fscanf(fp, "\n");
return kel;
}
/*!
* \brief kernelWrite()
*
* \param[in] fname output file
* \param[in] kel kernel
* \return 0 if OK, 1 on error
*/
l_ok
kernelWrite(const char *fname,
L_KERNEL *kel)
{
FILE *fp;
PROCNAME("kernelWrite");
if (!fname)
return ERROR_INT("fname not defined", procName, 1);
if (!kel)
return ERROR_INT("kel not defined", procName, 1);
if ((fp = fopenWriteStream(fname, "wb")) == NULL)
return ERROR_INT("stream not opened", procName, 1);
kernelWriteStream(fp, kel);
fclose(fp);
return 0;
}
/*!
* \brief kernelWriteStream()
*
* \param[in] fp file stream
* \param[in] kel
* \return 0 if OK, 1 on error
*/
l_ok
kernelWriteStream(FILE *fp,
L_KERNEL *kel)
{
l_int32 sx, sy, cx, cy, i, j;
PROCNAME("kernelWriteStream");
if (!fp)
return ERROR_INT("stream not defined", procName, 1);
if (!kel)
return ERROR_INT("kel not defined", procName, 1);
kernelGetParameters(kel, &sy, &sx, &cy, &cx);
fprintf(fp, " Kernel Version %d\n", KERNEL_VERSION_NUMBER);
fprintf(fp, " sy = %d, sx = %d, cy = %d, cx = %d\n", sy, sx, cy, cx);
for (i = 0; i < sy; i++) {
for (j = 0; j < sx; j++)
fprintf(fp, "%15.4f", kel->data[i][j]);
fprintf(fp, "\n");
}
fprintf(fp, "\n");
return 0;
}
/*----------------------------------------------------------------------*
* Making a kernel from a compiled string *
*----------------------------------------------------------------------*/
/*!
* \brief kernelCreateFromString()
*
* \param[in] h, w height, width
* \param[in] cy, cx origin
* \param[in] kdata
* \return kernel of the given size, or NULL on error
*
*
* Notes:
* (1) The data is an array of chars, in row-major order, giving
* space separated integers in the range [-255 ... 255].
* (2) The only other formatting limitation is that you must
* leave space between the last number in each row and
* the double-quote. If possible, it's also nice to have each
* line in the string represent a line in the kernel; e.g.,
* static const char *kdata =
* " 20 50 20 "
* " 70 140 70 "
* " 20 50 20 ";
*
*/
L_KERNEL *
kernelCreateFromString(l_int32 h,
l_int32 w,
l_int32 cy,
l_int32 cx,
const char *kdata)
{
l_int32 n, i, j, index;
l_float32 val;
L_KERNEL *kel;
NUMA *na;
PROCNAME("kernelCreateFromString");
if (h < 1)
return (L_KERNEL *)ERROR_PTR("height must be > 0", procName, NULL);
if (w < 1)
return (L_KERNEL *)ERROR_PTR("width must be > 0", procName, NULL);
if (cy < 0 || cy >= h)
return (L_KERNEL *)ERROR_PTR("cy invalid", procName, NULL);
if (cx < 0 || cx >= w)
return (L_KERNEL *)ERROR_PTR("cx invalid", procName, NULL);
kel = kernelCreate(h, w);
kernelSetOrigin(kel, cy, cx);
na = parseStringForNumbers(kdata, " \t\n");
n = numaGetCount(na);
if (n != w * h) {
kernelDestroy(&kel);
numaDestroy(&na);
lept_stderr("w = %d, h = %d, num ints = %d\n", w, h, n);
return (L_KERNEL *)ERROR_PTR("invalid integer data", procName, NULL);
}
index = 0;
for (i = 0; i < h; i++) {
for (j = 0; j < w; j++) {
numaGetFValue(na, index, &val);
kernelSetElement(kel, i, j, val);
index++;
}
}
numaDestroy(&na);
return kel;
}
/*----------------------------------------------------------------------*
* Making a kernel from a simple file format *
*----------------------------------------------------------------------*/
/*!
* \brief kernelCreateFromFile()
*
* \param[in] filename
* \return kernel, or NULL on error
*
*
* Notes:
* (1) The file contains, in the following order:
* ~ Any number of comment lines starting with '#' are ignored
* ~ The height and width of the kernel
* ~ The y and x values of the kernel origin
* ~ The kernel data, formatted as lines of numbers (integers
* or floats) for the kernel values in row-major order,
* and with no other punctuation.
* (Note: this differs from kernelCreateFromString(),
* where each line must begin and end with a double-quote
* to tell the compiler it's part of a string.)
* ~ The kernel specification ends when a blank line,
* a comment line, or the end of file is reached.
* (2) All lines must be left-justified.
* (3) See kernelCreateFromString() for a description of the string
* format for the kernel data. As an example, here are the lines
* of a valid kernel description file In the file, all lines
* are left-justified:
* \code
* # small 3x3 kernel
* 3 3
* 1 1
* 25.5 51 24.3
* 70.2 146.3 73.4
* 20 50.9 18.4
* \endcode
*
*/
L_KERNEL *
kernelCreateFromFile(const char *filename)
{
char *filestr, *line;
l_int32 nlines, i, j, first, index, w, h, cx, cy, n;
l_float32 val;
size_t size;
NUMA *na, *nat;
SARRAY *sa;
L_KERNEL *kel;
PROCNAME("kernelCreateFromFile");
if (!filename)
return (L_KERNEL *)ERROR_PTR("filename not defined", procName, NULL);
if ((filestr = (char *)l_binaryRead(filename, &size)) == NULL)
return (L_KERNEL *)ERROR_PTR("file not found", procName, NULL);
if (size == 0) {
LEPT_FREE(filestr);
return (L_KERNEL *)ERROR_PTR("file is empty", procName, NULL);
}
sa = sarrayCreateLinesFromString(filestr, 1);
LEPT_FREE(filestr);
nlines = sarrayGetCount(sa);
/* Find the first data line. */
for (i = 0, first = 0; i < nlines; i++) {
line = sarrayGetString(sa, i, L_NOCOPY);
if (line[0] != '#') {
first = i;
break;
}
}
/* Find the kernel dimensions and origin location. */
line = sarrayGetString(sa, first, L_NOCOPY);
if (sscanf(line, "%d %d", &h, &w) != 2) {
sarrayDestroy(&sa);
return (L_KERNEL *)ERROR_PTR("error reading h,w", procName, NULL);
}
if (h > MaxArraySize || w > MaxArraySize) {
L_ERROR("h = %d or w = %d > %d\n", procName, h, w, MaxArraySize);
sarrayDestroy(&sa);
return NULL;
}
line = sarrayGetString(sa, first + 1, L_NOCOPY);
if (sscanf(line, "%d %d", &cy, &cx) != 2) {
sarrayDestroy(&sa);
return (L_KERNEL *)ERROR_PTR("error reading cy,cx", procName, NULL);
}
/* Extract the data. This ends when we reach eof, or when we
* encounter a line of data that is either a null string or
* contains just a newline. */
na = numaCreate(0);
for (i = first + 2; i < nlines; i++) {
line = sarrayGetString(sa, i, L_NOCOPY);
if (line[0] == '\0' || line[0] == '\n' || line[0] == '#')
break;
nat = parseStringForNumbers(line, " \t\n");
numaJoin(na, nat, 0, -1);
numaDestroy(&nat);
}
sarrayDestroy(&sa);
n = numaGetCount(na);
if (n != w * h) {
numaDestroy(&na);
lept_stderr("w = %d, h = %d, num ints = %d\n", w, h, n);
return (L_KERNEL *)ERROR_PTR("invalid integer data", procName, NULL);
}
kel = kernelCreate(h, w);
kernelSetOrigin(kel, cy, cx);
index = 0;
for (i = 0; i < h; i++) {
for (j = 0; j < w; j++) {
numaGetFValue(na, index, &val);
kernelSetElement(kel, i, j, val);
index++;
}
}
numaDestroy(&na);
return kel;
}
/*----------------------------------------------------------------------*
* Making a kernel from a Pix *
*----------------------------------------------------------------------*/
/*!
* \brief kernelCreateFromPix()
*
* \param[in] pix
* \param[in] cy, cx origin of kernel
* \return kernel, or NULL on error
*
*
* Notes:
* (1) The origin must be positive and within the dimensions of the pix.
*
*/
L_KERNEL *
kernelCreateFromPix(PIX *pix,
l_int32 cy,
l_int32 cx)
{
l_int32 i, j, w, h, d;
l_uint32 val;
L_KERNEL *kel;
PROCNAME("kernelCreateFromPix");
if (!pix)
return (L_KERNEL *)ERROR_PTR("pix not defined", procName, NULL);
pixGetDimensions(pix, &w, &h, &d);
if (d != 8)
return (L_KERNEL *)ERROR_PTR("pix not 8 bpp", procName, NULL);
if (cy < 0 || cx < 0 || cy >= h || cx >= w)
return (L_KERNEL *)ERROR_PTR("(cy, cx) invalid", procName, NULL);
kel = kernelCreate(h, w);
kernelSetOrigin(kel, cy, cx);
for (i = 0; i < h; i++) {
for (j = 0; j < w; j++) {
pixGetPixel(pix, j, i, &val);
kernelSetElement(kel, i, j, (l_float32)val);
}
}
return kel;
}
/*----------------------------------------------------------------------*
* Display a kernel in a pix *
*----------------------------------------------------------------------*/
/*!
* \brief kernelDisplayInPix()
*
* \param[in] kel kernel
* \param[in] size of grid interiors; odd; either 1 or a minimum size
* of 17 is enforced
* \param[in] gthick grid thickness; either 0 or a minimum size of 2
* is enforced
* \return pix display of kernel, or NULL on error
*
*
* Notes:
* (1) This gives a visual representation of a kernel.
* (2) There are two modes of display:
* (a) Grid lines of minimum width 2, surrounding regions
* representing kernel elements of minimum size 17,
* with a "plus" mark at the kernel origin, or
* (b) A pix without grid lines and using 1 pixel per kernel element.
* (3) For both cases, the kernel absolute value is displayed,
* normalized such that the maximum absolute value is 255.
* (4) Large 2D separable kernels should be used for convolution
* with two 1D kernels. However, for the bilateral filter,
* the computation time is independent of the size of the
* 2D content kernel.
*
*/
PIX *
kernelDisplayInPix(L_KERNEL *kel,
l_int32 size,
l_int32 gthick)
{
l_int32 i, j, w, h, sx, sy, cx, cy, width, x0, y0;
l_int32 normval;
l_float32 minval, maxval, max, val, norm;
PIX *pixd, *pixt0, *pixt1;
PROCNAME("kernelDisplayInPix");
if (!kel)
return (PIX *)ERROR_PTR("kernel not defined", procName, NULL);
/* Normalize the max value to be 255 for display */
kernelGetParameters(kel, &sy, &sx, &cy, &cx);
kernelGetMinMax(kel, &minval, &maxval);
max = L_MAX(maxval, -minval);
if (max == 0.0)
return (PIX *)ERROR_PTR("kernel elements all 0.0", procName, NULL);
norm = 255. / (l_float32)max;
/* Handle the 1 element/pixel case; typically with large kernels */
if (size == 1 && gthick == 0) {
pixd = pixCreate(sx, sy, 8);
for (i = 0; i < sy; i++) {
for (j = 0; j < sx; j++) {
kernelGetElement(kel, i, j, &val);
normval = (l_int32)(norm * L_ABS(val));
pixSetPixel(pixd, j, i, normval);
}
}
return pixd;
}
/* Enforce the constraints for the grid line version */
if (size < 17) {
L_WARNING("size < 17; setting to 17\n", procName);
size = 17;
}
if (size % 2 == 0)
size++;
if (gthick < 2) {
L_WARNING("grid thickness < 2; setting to 2\n", procName);
gthick = 2;
}
w = size * sx + gthick * (sx + 1);
h = size * sy + gthick * (sy + 1);
pixd = pixCreate(w, h, 8);
/* Generate grid lines */
for (i = 0; i <= sy; i++)
pixRenderLine(pixd, 0, gthick / 2 + i * (size + gthick),
w - 1, gthick / 2 + i * (size + gthick),
gthick, L_SET_PIXELS);
for (j = 0; j <= sx; j++)
pixRenderLine(pixd, gthick / 2 + j * (size + gthick), 0,
gthick / 2 + j * (size + gthick), h - 1,
gthick, L_SET_PIXELS);
/* Generate mask for each element */
pixt0 = pixCreate(size, size, 1);
pixSetAll(pixt0);
/* Generate crossed lines for origin pattern */
pixt1 = pixCreate(size, size, 1);
width = size / 8;
pixRenderLine(pixt1, size / 2, (l_int32)(0.12 * size),
size / 2, (l_int32)(0.88 * size),
width, L_SET_PIXELS);
pixRenderLine(pixt1, (l_int32)(0.15 * size), size / 2,
(l_int32)(0.85 * size), size / 2,
width, L_FLIP_PIXELS);
pixRasterop(pixt1, size / 2 - width, size / 2 - width,
2 * width, 2 * width, PIX_NOT(PIX_DST), NULL, 0, 0);
/* Paste the patterns in */
y0 = gthick;
for (i = 0; i < sy; i++) {
x0 = gthick;
for (j = 0; j < sx; j++) {
kernelGetElement(kel, i, j, &val);
normval = (l_int32)(norm * L_ABS(val));
pixSetMaskedGeneral(pixd, pixt0, normval, x0, y0);
if (i == cy && j == cx)
pixPaintThroughMask(pixd, pixt1, x0, y0, 255 - normval);
x0 += size + gthick;
}
y0 += size + gthick;
}
pixDestroy(&pixt0);
pixDestroy(&pixt1);
return pixd;
}
/*------------------------------------------------------------------------*
* Parse string to extract numbers *
*------------------------------------------------------------------------*/
/*!
* \brief parseStringForNumbers()
*
* \param[in] str string containing numbers; not changed
* \param[in] seps string of characters that can be used between ints
* \return numa of numbers found, or NULL on error
*
*
* Notes:
* (1) The numbers can be ints or floats.
*
*/
NUMA *
parseStringForNumbers(const char *str,
const char *seps)
{
char *newstr, *head;
char *tail = NULL;
l_float32 val;
NUMA *na;
PROCNAME("parseStringForNumbers");
if (!str)
return (NUMA *)ERROR_PTR("str not defined", procName, NULL);
newstr = stringNew(str); /* to enforce const-ness of str */
na = numaCreate(0);
head = strtokSafe(newstr, seps, &tail);
val = atof(head);
numaAddNumber(na, val);
LEPT_FREE(head);
while ((head = strtokSafe(NULL, seps, &tail)) != NULL) {
val = atof(head);
numaAddNumber(na, val);
LEPT_FREE(head);
}
LEPT_FREE(newstr);
return na;
}
/*------------------------------------------------------------------------*
* Simple parametric kernels *
*------------------------------------------------------------------------*/
/*!
* \brief makeFlatKernel()
*
* \param[in] height, width
* \param[in] cy, cx origin of kernel
* \return kernel, or NULL on error
*
*
* Notes:
* (1) This is the same low-pass filtering kernel that is used
* in the block convolution functions.
* (2) The kernel origin (%cy, %cx) is typically placed as near
* the center of the kernel as possible. If height and
* width are odd, then using %cy = height / 2 and
* %cx = width / 2 places the origin at the exact center.
* (3) This returns a normalized kernel.
*
*/
L_KERNEL *
makeFlatKernel(l_int32 height,
l_int32 width,
l_int32 cy,
l_int32 cx)
{
l_int32 i, j;
l_float32 normval;
L_KERNEL *kel;
PROCNAME("makeFlatKernel");
if ((kel = kernelCreate(height, width)) == NULL)
return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
kernelSetOrigin(kel, cy, cx);
normval = 1.0 / (l_float32)(height * width);
for (i = 0; i < height; i++) {
for (j = 0; j < width; j++) {
kernelSetElement(kel, i, j, normval);
}
}
return kel;
}
/*!
* \brief makeGaussianKernel()
*
* \param[in] halfh sy = 2 * halfh + 1
* \param[in] halfw sx = 2 * halfw + 1
* \param[in] stdev standard deviation
* \param[in] max value at (cx,cy)
* \return kernel, or NULL on error
*
*
* Notes:
* (1) The kernel size (sx, sy) = (2 * %halfw + 1, 2 * %halfh + 1)
* (2) The kernel center (cx, cy) = (%halfw, %halfh).
* (3) %halfw and %halfh are typically equal, and
* are typically several times larger than the standard deviation.
* (4) If pixConvolve() is invoked with normalization (the sum of
* kernel elements = 1.0), use 1.0 for max (or any number that's
* not too small or too large).
*
*/
L_KERNEL *
makeGaussianKernel(l_int32 halfh,
l_int32 halfw,
l_float32 stdev,
l_float32 max)
{
l_int32 sx, sy, i, j;
l_float32 val;
L_KERNEL *kel;
PROCNAME("makeGaussianKernel");
sx = 2 * halfw + 1;
sy = 2 * halfh + 1;
if ((kel = kernelCreate(sy, sx)) == NULL)
return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
kernelSetOrigin(kel, halfh, halfw);
for (i = 0; i < sy; i++) {
for (j = 0; j < sx; j++) {
val = expf(-(l_float32)((i - halfh) * (i - halfh) +
(j - halfw) * (j - halfw)) /
(2. * stdev * stdev));
kernelSetElement(kel, i, j, max * val);
}
}
return kel;
}
/*!
* \brief makeGaussianKernelSep()
*
* \param[in] halfh sy = 2 * halfh + 1
* \param[in] halfw sx = 2 * halfw + 1
* \param[in] stdev standard deviation
* \param[in] max value at (cx,cy)
* \param[out] pkelx x part of kernel
* \param[out] pkely y part of kernel
* \return 0 if OK, 1 on error
*
*
* Notes:
* (1) See makeGaussianKernel() for description of input parameters.
* (2) These kernels are constructed so that the result of both
* normalized and un-normalized convolution will be the same
* as when convolving with pixConvolve() using the full kernel.
* (3) The trick for the un-normalized convolution is to have the
* product of the two kernel elemets at (cx,cy) be equal to %max,
* not max**2. That's why %max for kely is 1.0. If instead
* we use sqrt(%max) for both, the results are slightly less
* accurate, when compared to using the full kernel in
* makeGaussianKernel().
*
*/
l_ok
makeGaussianKernelSep(l_int32 halfh,
l_int32 halfw,
l_float32 stdev,
l_float32 max,
L_KERNEL **pkelx,
L_KERNEL **pkely)
{
PROCNAME("makeGaussianKernelSep");
if (!pkelx || !pkely)
return ERROR_INT("&kelx and &kely not defined", procName, 1);
*pkelx = makeGaussianKernel(0, halfw, stdev, max);
*pkely = makeGaussianKernel(halfh, 0, stdev, 1.0);
return 0;
}
/*!
* \brief makeDoGKernel()
*
* \param[in] halfh sy = 2 * halfh + 1
* \param[in] halfw sx = 2 * halfw + 1
* \param[in] stdev standard deviation of narrower gaussian
* \param[in] ratio of stdev for wide filter to stdev for narrow one
* \return kernel, or NULL on error
*
*
* Notes:
* (1) The DoG (difference of gaussians) is a wavelet mother
* function with null total sum. By subtracting two blurred
* versions of the image, it acts as a bandpass filter for
* frequencies passed by the narrow gaussian but stopped
* by the wide one.See:
* http://en.wikipedia.org/wiki/Difference_of_Gaussians
* (2) The kernel size (sx, sy) = (2 * halfw + 1, 2 * halfh + 1).
* (3) The kernel center (cx, cy) = (halfw, halfh).
* (4) %halfw and %halfh are typically equal, and are typically
* several times larger than the standard deviation.
* (5) %ratio is the ratio of standard deviations of the wide
* to narrow gaussian. It must be >= 1.0; 1.0 is a no-op.
* (6) Because the kernel is a null sum, it must be invoked without
* normalization in pixConvolve().
*
*/
L_KERNEL *
makeDoGKernel(l_int32 halfh,
l_int32 halfw,
l_float32 stdev,
l_float32 ratio)
{
l_int32 sx, sy, i, j;
l_float32 pi, squaredist, highnorm, lownorm, val;
L_KERNEL *kel;
PROCNAME("makeDoGKernel");
sx = 2 * halfw + 1;
sy = 2 * halfh + 1;
if ((kel = kernelCreate(sy, sx)) == NULL)
return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
kernelSetOrigin(kel, halfh, halfw);
pi = 3.1415926535f;
for (i = 0; i < sy; i++) {
for (j = 0; j < sx; j++) {
squaredist = (l_float32)((i - halfh) * (i - halfh) +
(j - halfw) * (j - halfw));
highnorm = 1. / (2 * stdev * stdev);
lownorm = highnorm / (ratio * ratio);
val = (highnorm / pi) * expf(-(highnorm * squaredist))
- (lownorm / pi) * expf(-(lownorm * squaredist));
kernelSetElement(kel, i, j, val);
}
}
return kel;
}