/*====================================================================*
- 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.
*====================================================================*/
/*
* enhance.c
*
* Gamma TRC (tone reproduction curve) mapping
* PIX *pixGammaTRC()
* PIX *pixGammaTRCMasked()
* NUMA *numaGammaTRC()
*
* Contrast enhancement
* PIX *pixContrastTRC()
* PIX *pixContrastTRCMasked()
* NUMA *numaContrastTRC()
*
* Histogram equalization
* PIX *pixEqualizeTRC()
* NUMA *numaEqualizeTRC()
*
* Generic TRC mapper
* PIX *pixTRCMap()
*
* Unsharp-masking
* PIX *pixUnsharpMasking()
* PIX *pixUnsharpMaskingGray()
* PIX *pixUnsharpMaskingFast()
* PIX *pixUnsharpMaskingGrayFast()
* PIX *pixUnsharpMaskingGray1D()
* PIX *pixUnsharpMaskingGray2D()
*
* Hue and saturation modification
* PIX *pixModifyHue()
* PIX *pixModifySaturation()
* l_int32 pixMeasureSaturation()
*
* General multiplicative constant color transform
* PIX *pixMultConstantColor()
* PIX *pixMultMatrixColor()
*
* Edge by bandpass
* PIX *pixHalfEdgeByBandpass()
*
* Gamma correction, contrast enhancement and histogram equalization
* apply a simple mapping function to each pixel (or, for color
* images, to each sample (i.e., r,g,b) of the pixel).
*
* - Gamma correction either lightens the image or darkens
* it, depending on whether the gamma factor is greater
* or less than 1.0, respectively.
*
* - Contrast enhancement darkens the pixels that are already
* darker than the middle of the dynamic range (128)
* and lightens pixels that are lighter than 128.
*
* - Histogram equalization remaps to have the same number
* of image pixels at each of 256 intensity values. This is
* a quick and dirty method of adjusting contrast and brightness
* to bring out details in both light and dark regions.
*
* Unsharp masking is a more complicated enhancement.
* A "high frequency" image, generated by subtracting
* the smoothed ("low frequency") part of the image from
* itself, has all the energy at the edges. This "edge image"
* has 0 average value. A fraction of the edge image is
* then added to the original, enhancing the differences
* between pixel values at edges. Because we represent
* images as l_uint8 arrays, we preserve dynamic range and
* handle negative values by doing all the arithmetic on
* shifted l_uint16 arrays; the l_uint8 values are recovered
* at the end.
*
* Hue and saturation modification work in HSV space. Because
* this is too large for efficient table lookup, each pixel value
* is transformed to HSV, modified, and transformed back.
* It's not the fastest way to do this, but the method is
* easily understood.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "allheaders.h"
/* Scales contrast enhancement factor to have a useful range
* between 0.0 and 1.0 */
static const l_float32 ENHANCE_SCALE_FACTOR = 5.;
/* Default number of pixels sampled to determine histogram */
static const l_int32 DEFAULT_HISTO_SAMPLES = 100000;
/*-------------------------------------------------------------*
* Gamma TRC (tone reproduction curve) mapping *
*-------------------------------------------------------------*/
/*!
* pixGammaTRC()
*
* Input: pixd (<optional> null or equal to pixs)
* pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
* gamma (gamma correction; must be > 0.0)
* minval (input value that gives 0 for output; can be < 0)
* maxval (input value that gives 255 for output; can be > 255)
* Return: pixd always
*
* Notes:
* (1) pixd must either be null or equal to pixs.
* Set pixd == pixs to get in-place operation;
* set pixd == null to get new image.
* (2) If pixs is colormapped, the colormap is transformed,
* either in-place or in a copy of pixs.
* (3) We use a gamma mapping between minval and maxval.
* (4) If gamma < 1.0, the image will appear darker;
* if gamma > 1.0, the image will appear lighter;
* if gamma == 1.0 and minval == 0 and maxval == 255, return a clone
* (5) For color images that are not colormapped, the mapping
* is applied to each component.
* (6) minval and maxval are not restricted to the interval [0, 255].
* If minval < 0, an input value of 0 is mapped to a
* nonzero output. This will turn black to gray.
* If maxval > 255, an input value of 255 is mapped to
* an output value less than 255. This will turn
* white (e.g., in the background) to gray.
* (7) Increasing minval darkens the image.
* (8) Decreasing maxval bleaches the image.
* (9) Simultaneously increasing minval and decreasing maxval
* will darken the image and make the colors more intense;
* e.g., minval = 50, maxval = 200.
* (10) See numaGammaTRC() for further examples of use.
*/
PIX *
pixGammaTRC(PIX *pixd,
PIX *pixs,
l_float32 gamma,
l_int32 minval,
l_int32 maxval)
{
l_int32 d;
NUMA *nag;
PIXCMAP *cmap;
PROCNAME("pixGammaTRC");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
if (pixd && (pixd != pixs))
return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
if (gamma <= 0.0) {
L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
gamma = 1.0;
}
if (minval >= maxval)
return (PIX *)ERROR_PTR("minval not < maxval", procName, pixd);
if (gamma == 1.0 && minval == 0 && maxval == 255)
return pixClone(pixs);
cmap = pixGetColormap(pixs);
d = pixGetDepth(pixs);
if (!cmap && d != 8 && d != 32)
return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
if (!pixd) /* start with a copy if not in-place */
pixd = pixCopy(NULL, pixs);
if (cmap) {
pixcmapGammaTRC(pixGetColormap(pixd), gamma, minval, maxval);
return pixd;
}
/* pixd is 8 or 32 bpp */
if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
return (PIX *)ERROR_PTR("nag not made", procName, pixd);
pixTRCMap(pixd, NULL, nag);
numaDestroy(&nag);
return pixd;
}
/*!
* pixGammaTRCMasked()
*
* Input: pixd (<optional> null or equal to pixs)
* pixs (8 or 32 bpp; not colormapped)
* pixm (<optional> null or 1 bpp)
* gamma (gamma correction; must be > 0.0)
* minval (input value that gives 0 for output; can be < 0)
* maxval (input value that gives 255 for output; can be > 255)
* Return: pixd always
*
* Notes:
* (1) Same as pixGammaTRC() except mapping is optionally over
* a subset of pixels described by pixm.
* (2) Masking does not work for colormapped images.
* (3) See pixGammaTRC() for details on how to use the parameters.
*/
PIX *
pixGammaTRCMasked(PIX *pixd,
PIX *pixs,
PIX *pixm,
l_float32 gamma,
l_int32 minval,
l_int32 maxval)
{
l_int32 d;
NUMA *nag;
PROCNAME("pixGammaTRCMasked");
if (!pixm)
return pixGammaTRC(pixd, pixs, gamma, minval, maxval);
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
if (pixGetColormap(pixs))
return (PIX *)ERROR_PTR("invalid: pixs has a colormap", procName, pixd);
if (pixd && (pixd != pixs))
return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
d = pixGetDepth(pixs);
if (d != 8 && d != 32)
return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
if (minval >= maxval)
return (PIX *)ERROR_PTR("minval not < maxval", procName, pixd);
if (gamma <= 0.0) {
L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
gamma = 1.0;
}
if (!pixd) /* start with a copy if not in-place */
pixd = pixCopy(NULL, pixs);
if ((nag = numaGammaTRC(gamma, minval, maxval)) == NULL)
return (PIX *)ERROR_PTR("nag not made", procName, pixd);
pixTRCMap(pixd, pixm, nag);
numaDestroy(&nag);
return pixd;
}
/*!
* numaGammaTRC()
*
* Input: gamma (gamma factor; must be > 0.0)
* minval (input value that gives 0 for output)
* maxval (input value that gives 255 for output)
* Return: na, or null on error
*
* Notes:
* (1) The map is returned as a numa; values are clipped to [0, 255].
* (2) To force all intensities into a range within fraction delta
* of white, use: minval = -256 * (1 - delta) / delta
* maxval = 255
* (3) To force all intensities into a range within fraction delta
* of black, use: minval = 0
* maxval = 256 * (1 - delta) / delta
*/
NUMA *
numaGammaTRC(l_float32 gamma,
l_int32 minval,
l_int32 maxval)
{
l_int32 i, val;
l_float32 x, invgamma;
NUMA *na;
PROCNAME("numaGammaTRC");
if (minval >= maxval)
return (NUMA *)ERROR_PTR("minval not < maxval", procName, NULL);
if (gamma <= 0.0) {
L_WARNING("gamma must be > 0.0; setting to 1.0", procName);
gamma = 1.0;
}
invgamma = 1. / gamma;
na = numaCreate(256);
for (i = 0; i < minval; i++)
numaAddNumber(na, 0);
for (i = minval; i <= maxval; i++) {
if (i < 0) continue;
if (i > 255) continue;
x = (l_float32)(i - minval) / (l_float32)(maxval - minval);
val = (l_int32)(255. * powf(x, invgamma) + 0.5);
val = L_MAX(val, 0);
val = L_MIN(val, 255);
numaAddNumber(na, val);
}
for (i = maxval + 1; i < 256; i++)
numaAddNumber(na, 255);
return na;
}
/*-------------------------------------------------------------*
* Contrast enhancement *
*-------------------------------------------------------------*/
/*!
* pixContrastTRC()
*
* Input: pixd (<optional> null or equal to pixs)
* pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
* factor (0.0 is no enhancement)
* Return: pixd always
*
* Notes:
* (1) pixd must either be null or equal to pixs.
* Set pixd == pixs to get in-place operation;
* set pixd == null to get new image.
* (2) If pixs is colormapped, the colormap is transformed,
* either in-place or in a copy of pixs.
* (3) Contrast is enhanced by mapping each color component
* using an atan function with maximum slope at 127.
* Pixels below 127 are lowered in intensity and pixels
* above 127 are increased.
* (4) The useful range for the contrast factor is scaled to
* be in (0.0 to 1.0), but larger values can also be used.
* 0.0 corresponds to no enhancement; return a clone.
* (5) For color images that are not colormapped, the mapping
* is applied to each component.
*/
PIX *
pixContrastTRC(PIX *pixd,
PIX *pixs,
l_float32 factor)
{
l_int32 d;
NUMA *nac;
PIXCMAP *cmap;
PROCNAME("pixContrastTRC");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
if (pixd && (pixd != pixs))
return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
if (factor < 0.0) {
L_WARNING("factor must be >= 0.0; using 0.0", procName);
factor = 0.0;
}
if (factor == 0.0)
return pixClone(pixs);
cmap = pixGetColormap(pixs);
d = pixGetDepth(pixs);
if (!cmap && d != 8 && d != 32)
return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
if (!pixd) /* start with a copy if not in-place */
pixd = pixCopy(NULL, pixs);
if (cmap) {
pixcmapContrastTRC(pixGetColormap(pixd), factor);
return pixd;
}
/* pixd is 8 or 32 bpp */
if ((nac = numaContrastTRC(factor)) == NULL)
return (PIX *)ERROR_PTR("nac not made", procName, pixd);
pixTRCMap(pixd, NULL, nac);
numaDestroy(&nac);
return pixd;
}
/*!
* pixContrastTRCMasked()
*
* Input: pixd (<optional> null or equal to pixs)
* pixs (8 or 32 bpp; or 2, 4 or 8 bpp with colormap)
* pixm (<optional> null or 1 bpp)
* factor (0.0 is no enhancement)
* Return: pixd always
*
* Notes:
* (1) Same as pixContrastTRC() except mapping is optionally over
* a subset of pixels described by pixm.
* (2) Masking does not work for colormapped images.
* (3) See pixContrastTRC() for details on how to use the parameters.
*/
PIX *
pixContrastTRCMasked(PIX *pixd,
PIX *pixs,
PIX *pixm,
l_float32 factor)
{
l_int32 d;
NUMA *nac;
PROCNAME("pixContrastTRCMasked");
if (!pixm)
return pixContrastTRC(pixd, pixs, factor);
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
if (pixGetColormap(pixs))
return (PIX *)ERROR_PTR("invalid: pixs has a colormap", procName, pixd);
if (pixd && (pixd != pixs))
return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
d = pixGetDepth(pixs);
if (d != 8 && d != 32)
return (PIX *)ERROR_PTR("depth not 8 or 32 bpp", procName, pixd);
if (factor < 0.0) {
L_WARNING("factor must be >= 0.0; using 0.0", procName);
factor = 0.0;
}
if (factor == 0.0)
return pixClone(pixs);
if (!pixd) /* start with a copy if not in-place */
pixd = pixCopy(NULL, pixs);
if ((nac = numaContrastTRC(factor)) == NULL)
return (PIX *)ERROR_PTR("nac not made", procName, pixd);
pixTRCMap(pixd, pixm, nac);
numaDestroy(&nac);
return pixd;
}
/*!
* numaContrastTRC()
*
* Input: factor (generally between 0.0 (no enhancement)
* and 1.0, but can be larger than 1.0)
* Return: na, or null on error
*
* Notes:
* (1) The mapping is monotonic increasing, where 0 is mapped
* to 0 and 255 is mapped to 255.
* (2) As 'factor' is increased from 0.0 (where the mapping is linear),
* the map gets closer to its limit as a step function that
* jumps from 0 to 255 at the center (input value = 127).
*/
NUMA *
numaContrastTRC(l_float32 factor)
{
l_int32 i, val;
l_float64 x, ymax, ymin, dely, scale;
NUMA *na;
PROCNAME("numaContrastTRC");
if (factor < 0.0) {
L_WARNING("factor must be >= 0.0; using 0.0; no enhancement", procName);
factor = 0.0;
}
if (factor == 0.0)
return numaMakeSequence(0, 1, 256); /* linear map */
scale = ENHANCE_SCALE_FACTOR;
ymax = atan((l_float64)(1.0 * factor * scale));
ymin = atan((l_float64)(-127. * factor * scale / 128.));
dely = ymax - ymin;
na = numaCreate(256);
for (i = 0; i < 256; i++) {
x = (l_float64)i;
val = (l_int32)((255. / dely) *
(-ymin + atan((l_float64)(factor * scale * (x - 127.) / 128.))) +
0.5);
numaAddNumber(na, val);
}
return na;
}
/*-------------------------------------------------------------*
* Histogram equalization *
*-------------------------------------------------------------*/
/*!
* pixEqualizeTRC()
*
* Input: pixd (<optional> null or equal to pixs)
* pixs (8 bpp, or colormapped)
* fract (fraction of equalization movement of pixel values)
* factor (subsampling factor; integer >= 1)
* Return: pixd, or null on error
*
* Notes:
* (1) pixd must either be null or equal to pixs.
* Set pixd == pixs to get in-place operation;
* set pixd == null to get new image.
* (2) In histogram equalization, a tone reproduction curve
* mapping is used to make the number of pixels at each
* intensity equal.
* (3) If fract == 0.0, no equalization is performed; return a copy.
* If fract == 1.0, equalization is complete.
* (4) Set the subsampling factor > 1 to reduce the amount of computation.
* (5) If the pix is colormapped, the colormap is removed and
* it is converted to grayscale.
* (6) Note that even if there is a colormap, we can get an
* in-place operation because the intermediate image pixt
* is copied back to pixs (which for in-place is the same
* as pixd).
*/
PIX *
pixEqualizeTRC(PIX *pixd,
PIX *pixs,
l_float32 fract,
l_int32 factor)
{
NUMA *na;
PIX *pixt;
PIXCMAP *cmap;
PROCNAME("pixEqualizeTRC");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (pixd && (pixd != pixs))
return (PIX *)ERROR_PTR("pixd not null or pixs", procName, pixd);
cmap = pixGetColormap(pixs);
if ((pixGetDepth(pixs) != 8) && !cmap)
return (PIX *)ERROR_PTR("pixs not 8 bpp or cmapped", procName, NULL);
if (fract < 0.0 || fract > 1.0)
return (PIX *)ERROR_PTR("fract not in [0.0 ... 1.0]", procName, NULL);
if (factor < 1)
return (PIX *)ERROR_PTR("sampling factor < 1", procName, NULL);
if (fract == 0.0)
return pixCopy(pixd, pixs);
/* If there is a colormap, remove it. */
if (cmap)
pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_GRAYSCALE);
else
pixt = pixClone(pixs);
/* Make a copy if necessary */
pixd = pixCopy(pixd, pixt);
pixDestroy(&pixt);
if ((na = numaEqualizeTRC(pixd, fract, factor)) == NULL)
return (PIX *)ERROR_PTR("na not made", procName, pixd);
pixTRCMap(pixd, NULL, na);
numaDestroy(&na);
return pixd;
}
/*!
* numaEqualizeTRC()
*
* Input: pix (8 bpp, no colormap)
* fract (fraction of equalization movement of pixel values)
* factor (subsampling factor; integer >= 1)
* Return: nad, or null on error
*
* Notes:
* (1) If fract == 0.0, no equalization will be performed.
* If fract == 1.0, equalization is complete.
* (2) Set the subsampling factor > 1 to reduce the amount of computation.
* (3) The map is returned as a numa with 256 values, specifying
* the equalized value (array value) for every input value
* (the array index).
*/
NUMA *
numaEqualizeTRC(PIX *pix,
l_float32 fract,
l_int32 factor)
{
l_int32 iin, iout, itarg;
l_float32 val, sum;
NUMA *nah, *nasum, *nad;
PROCNAME("numaEqualizeTRC");
if (!pix)
return (NUMA *)ERROR_PTR("pix not defined", procName, NULL);
if (pixGetDepth(pix) != 8)
return (NUMA *)ERROR_PTR("pix not 8 bpp", procName, NULL);
if (fract < 0.0 || fract > 1.0)
return (NUMA *)ERROR_PTR("fract not in [0.0 ... 1.0]", procName, NULL);
if (factor < 1)
return (NUMA *)ERROR_PTR("sampling factor < 1", procName, NULL);
if (fract == 0.0)
L_WARNING("fract = 0.0; no equalization requested", procName);
if ((nah = pixGetGrayHistogram(pix, factor)) == NULL)
return (NUMA *)ERROR_PTR("histogram not made", procName, NULL);
numaGetSum(nah, &sum);
nasum = numaGetPartialSums(nah);
nad = numaCreate(256);
for (iin = 0; iin < 256; iin++) {
numaGetFValue(nasum, iin, &val);
itarg = (l_int32)(255. * val / sum + 0.5);
iout = iin + (l_int32)(fract * (itarg - iin));
iout = L_MIN(iout, 255); /* to be safe */
numaAddNumber(nad, iout);
}
numaDestroy(&nah);
numaDestroy(&nasum);
return nad;
}
/*-------------------------------------------------------------*
* Generic TRC mapping *
*-------------------------------------------------------------*/
/*!
* pixTRCMap()
*
* Input: pixs (8 grayscale or 32 bpp rgb; not colormapped)
* pixm (<optional> 1 bpp mask)
* na (mapping array)
* Return: pixd, or null on error
*
* Notes:
* (1) This operation is in-place on pixs.
* (2) For 32 bpp, this applies the same map to each of the r,g,b
* components.
* (3) The mapping array is of size 256, and it maps the input
* index into values in the range [0, 255].
* (4) If defined, the optional 1 bpp mask pixm has its origin
* aligned with pixs, and the map function is applied only
* to pixels in pixs under the fg of pixm.
*/
l_int32
pixTRCMap(PIX *pixs,
PIX *pixm,
NUMA *na)
{
l_int32 w, h, d, wm, hm, wpl, wplm, i, j, sval8, dval8;
l_int32 *tab;
l_uint32 sval32, dval32;
l_uint32 *data, *datam, *line, *linem;
PROCNAME("pixTRCMap");
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
if (pixGetColormap(pixs))
return ERROR_INT("pixs is colormapped", procName, 1);
if (!na)
return ERROR_INT("na not defined", procName, 1);
if (numaGetCount(na) != 256)
return ERROR_INT("na not of size 256", procName, 1);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 8 && d != 32)
return ERROR_INT("pixs not 8 or 32 bpp", procName, 1);
if (pixm) {
if (pixGetDepth(pixm) != 1)
return ERROR_INT("pixm not 1 bpp", procName, 1);
}
tab = numaGetIArray(na); /* get the array for efficiency */
wpl = pixGetWpl(pixs);
data = pixGetData(pixs);
if (!pixm) {
if (d == 8) {
for (i = 0; i < h; i++) {
line = data + i * wpl;
for (j = 0; j < w; j++) {
sval8 = GET_DATA_BYTE(line, j);
dval8 = tab[sval8];
SET_DATA_BYTE(line, j, dval8);
}
}
}
else { /* d == 32 */
for (i = 0; i < h; i++) {
line = data + i * wpl;
for (j = 0; j < w; j++) {
sval32 = *(line + j);
dval32 =
tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
*(line + j) = dval32;
}
}
}
}
else {
datam = pixGetData(pixm);
wplm = pixGetWpl(pixm);
pixGetDimensions(pixm, &wm, &hm, NULL);
if (d == 8) {
for (i = 0; i < h; i++) {
if (i >= hm)
break;
line = data + i * wpl;
linem = datam + i * wplm;
for (j = 0; j < w; j++) {
if (j >= wm)
break;
if (GET_DATA_BIT(linem, j) == 0)
continue;
sval8 = GET_DATA_BYTE(line, j);
dval8 = tab[sval8];
SET_DATA_BYTE(line, j, dval8);
}
}
}
else { /* d == 32 */
for (i = 0; i < h; i++) {
if (i >= hm)
break;
line = data + i * wpl;
linem = datam + i * wplm;
for (j = 0; j < w; j++) {
if (j >= wm)
break;
if (GET_DATA_BIT(linem, j) == 0)
continue;
sval32 = *(line + j);
dval32 =
tab[(sval32 >> L_RED_SHIFT) & 0xff] << L_RED_SHIFT |
tab[(sval32 >> L_GREEN_SHIFT) & 0xff] << L_GREEN_SHIFT |
tab[(sval32 >> L_BLUE_SHIFT) & 0xff] << L_BLUE_SHIFT;
*(line + j) = dval32;
}
}
}
}
FREE(tab);
return 0;
}
/*-----------------------------------------------------------------------*
* Unsharp masking *
*-----------------------------------------------------------------------*/
/*!
* pixUnsharpMasking()
*
* Input: pixs (all depths except 1 bpp; with or without colormaps)
* halfwidth ("half-width" of smoothing filter)
* fract (fraction of edge added back into image)
* Return: pixd, or null on error
*
* Notes:
* (1) We use symmetric smoothing filters of odd dimension,
* typically use sizes of 3, 5, 7, etc. The @halfwidth parameter
* for these is (size - 1)/2; i.e., 1, 2, 3, etc.
* (2) The fract parameter is typically taken in the
* range: 0.2 < fract < 0.7
*/
PIX *
pixUnsharpMasking(PIX *pixs,
l_int32 halfwidth,
l_float32 fract)
{
l_int32 d;
PIX *pixt, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;
PROCNAME("pixUnsharpMasking");
if (!pixs || (pixGetDepth(pixs) == 1))
return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", procName, NULL);
if (fract <= 0.0 || halfwidth <= 0) {
L_WARNING("no sharpening requested; clone returned", procName);
return pixClone(pixs);
}
if (halfwidth == 1 || halfwidth == 2)
return pixUnsharpMaskingFast(pixs, halfwidth, fract, L_BOTH_DIRECTIONS);
/* Remove colormap; clone if possible; result is either 8 or 32 bpp */
if ((pixt = pixConvertTo8Or32(pixs, 0, 1)) == NULL)
return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
/* Sharpen */
d = pixGetDepth(pixt);
if (d == 8)
pixd = pixUnsharpMaskingGray(pixt, halfwidth, fract);
else { /* d == 32 */
pixr = pixGetRGBComponent(pixs, COLOR_RED);
pixrs = pixUnsharpMaskingGray(pixr, halfwidth, fract);
pixDestroy(&pixr);
pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
pixgs = pixUnsharpMaskingGray(pixg, halfwidth, fract);
pixDestroy(&pixg);
pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
pixbs = pixUnsharpMaskingGray(pixb, halfwidth, fract);
pixDestroy(&pixb);
pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
pixDestroy(&pixrs);
pixDestroy(&pixgs);
pixDestroy(&pixbs);
}
pixDestroy(&pixt);
return pixd;
}
/*!
* pixUnsharpMaskingGray()
*
* Input: pixs (8 bpp; no colormap)
* halfwidth ("half-width" of smoothing filter)
* fract (fraction of edge added back into image)
* Return: pixd, or null on error
*
* Notes:
* (1) We use symmetric smoothing filters of odd dimension,
* typically use sizes of 3, 5, 7, etc. The @halfwidth parameter
* for these is (size - 1)/2; i.e., 1, 2, 3, etc.
* (2) The fract parameter is typically taken in the range:
* 0.2 < fract < 0.7
*/
PIX *
pixUnsharpMaskingGray(PIX *pixs,
l_int32 halfwidth,
l_float32 fract)
{
l_int32 w, h, d;
PIX *pixc, *pixd;
PIXACC *pixacc;
PROCNAME("pixUnsharpMaskingGray");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 8 || pixGetColormap(pixs) != NULL)
return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
if (fract <= 0.0 || halfwidth <= 0) {
L_WARNING("no sharpening requested; clone returned", procName);
return pixClone(pixs);
}
if (halfwidth == 1 || halfwidth == 2)
return pixUnsharpMaskingGrayFast(pixs, halfwidth, fract,
L_BOTH_DIRECTIONS);
if ((pixc = pixBlockconvGray(pixs, NULL, halfwidth, halfwidth)) == NULL)
return (PIX *)ERROR_PTR("pixc not made", procName, NULL);
/* Steps:
* (1) edge image is pixs - pixc (this is highpass part)
* (2) multiply edge image by fract
* (3) add fraction of edge to pixs
*
* To show how this is done with both interfaces to arithmetic
* on integer Pix, here is the implementation in the lower-level
* function calls:
* pixt = pixInitAccumulate(w, h, 0x10000000)) == NULL)
* pixAccumulate(pixt, pixs, L_ARITH_ADD);
* pixAccumulate(pixt, pixc, L_ARITH_SUBTRACT);
* pixMultConstAccumulate(pixt, fract, 0x10000000);
* pixAccumulate(pixt, pixs, L_ARITH_ADD);
* pixd = pixFinalAccumulate(pixt, 0x10000000, 8)) == NULL)
* pixDestroy(&pixt);
*
* The code below does the same thing using the Pixacc accumulator,
* hiding the details of the offset that is needed for subtraction.
*/
pixacc = pixaccCreate(w, h, 1);
pixaccAdd(pixacc, pixs);
pixaccSubtract(pixacc, pixc);
pixaccMultConst(pixacc, fract);
pixaccAdd(pixacc, pixs);
pixd = pixaccFinal(pixacc, 8);
pixaccDestroy(&pixacc);
pixDestroy(&pixc);
return pixd;
}
/*!
* pixUnsharpMaskingFast()
*
* Input: pixs (all depths except 1 bpp; with or without colormaps)
* halfwidth ("half-width" of smoothing filter; 1 and 2 only)
* fract (fraction of high frequency added to image)
* direction (L_HORIZ, L_VERT, L_BOTH_DIRECTIONS)
* Return: pixd, or null on error
*
* Notes:
* (1) The fast version uses separable 1-D filters directly on
* the input image. The halfwidth is either 1 (full width = 3)
* or 2 (full width = 5).
* (2) The fract parameter is typically taken in the
* range: 0.2 < fract < 0.7
* (3) To skip horizontal sharpening, use @fracth = 0.0; ditto for @fractv
* (4) For one dimensional filtering (as an example):
* For @halfwidth = 1, the low-pass filter is
* L: 1/3 1/3 1/3
* and the high-pass filter is
* H = I - L: -1/3 2/3 -1/3
* For @halfwidth = 2, the low-pass filter is
* L: 1/5 1/5 1/5 1/5 1/5
* and the high-pass filter is
* H = I - L: -1/5 -1/5 4/5 -1/5 -1/5
* The new sharpened pixel value is found by adding some fraction
* of the high-pass filter value (which sums to 0) to the
* initial pixel value:
* N = I + fract * H
* (5) For 2D, the sharpening filter is not separable, because the
* vertical filter depends on the horizontal location relative
* to the filter origin, and v.v. So we either do the full
* 2D filter (for @halfwidth == 1) or do the low-pass
* convolution separably and then compose with the original pix.
*/
PIX *
pixUnsharpMaskingFast(PIX *pixs,
l_int32 halfwidth,
l_float32 fract,
l_int32 direction)
{
l_int32 d;
PIX *pixt, *pixd, *pixr, *pixrs, *pixg, *pixgs, *pixb, *pixbs;
PROCNAME("pixUnsharpMaskingFast");
if (!pixs || (pixGetDepth(pixs) == 1))
return (PIX *)ERROR_PTR("pixs not defined or 1 bpp", procName, NULL);
if (fract <= 0.0 || halfwidth <= 0) {
L_WARNING("no sharpening requested; clone returned", procName);
return pixClone(pixs);
}
if (halfwidth != 1 && halfwidth != 2)
return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
if (direction != L_HORIZ && direction != L_VERT &&
direction != L_BOTH_DIRECTIONS)
return (PIX *)ERROR_PTR("invalid direction", procName, NULL);
/* Remove colormap; clone if possible; result is either 8 or 32 bpp */
if ((pixt = pixConvertTo8Or32(pixs, 0, 1)) == NULL)
return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
/* Sharpen */
d = pixGetDepth(pixt);
if (d == 8)
pixd = pixUnsharpMaskingGrayFast(pixt, halfwidth, fract, direction);
else { /* d == 32 */
pixr = pixGetRGBComponent(pixs, COLOR_RED);
pixrs = pixUnsharpMaskingGrayFast(pixr, halfwidth, fract, direction);
pixDestroy(&pixr);
pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
pixgs = pixUnsharpMaskingGrayFast(pixg, halfwidth, fract, direction);
pixDestroy(&pixg);
pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
pixbs = pixUnsharpMaskingGrayFast(pixb, halfwidth, fract, direction);
pixDestroy(&pixb);
pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
pixDestroy(&pixrs);
pixDestroy(&pixgs);
pixDestroy(&pixbs);
}
pixDestroy(&pixt);
return pixd;
}
/*!
* pixUnsharpMaskingGrayFast()
*
* Input: pixs (8 bpp; no colormap)
* halfwidth ("half-width" of smoothing filter: 1 or 2)
* fract (fraction of high frequency added to image)
* direction (L_HORIZ, L_VERT, L_BOTH_DIRECTIONS)
* Return: pixd, or null on error
*
* Notes:
* (1) For usage and explanation of the algorithm, see notes
* in pixUnsharpMaskingFast().
*/
PIX *
pixUnsharpMaskingGrayFast(PIX *pixs,
l_int32 halfwidth,
l_float32 fract,
l_int32 direction)
{
PIX *pixd;
PROCNAME("pixUnsharpMaskingGrayFast");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (pixGetDepth(pixs) != 8 || pixGetColormap(pixs) != NULL)
return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
if (fract <= 0.0 || halfwidth <= 0) {
L_WARNING("no sharpening requested; clone returned", procName);
return pixClone(pixs);
}
if (halfwidth != 1 && halfwidth != 2)
return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
if (direction != L_HORIZ && direction != L_VERT &&
direction != L_BOTH_DIRECTIONS)
return (PIX *)ERROR_PTR("invalid direction", procName, NULL);
if (direction != L_BOTH_DIRECTIONS)
pixd = pixUnsharpMaskingGray1D(pixs, halfwidth, fract, direction);
else /* 2D sharpening */
pixd = pixUnsharpMaskingGray2D(pixs, halfwidth, fract);
return pixd;
}
/*!
* pixUnsharpMaskingGray1D()
*
* Input: pixs (8 bpp; no colormap)
* halfwidth ("half-width" of smoothing filter: 1 or 2)
* fract (fraction of high frequency added to image)
* direction (of filtering; use L_HORIZ or L_VERT)
* Return: pixd, or null on error
*
* Notes:
* (1) For usage and explanation of the algorithm, see notes
* in pixUnsharpMaskingFast().
*/
PIX *
pixUnsharpMaskingGray1D(PIX *pixs,
l_int32 halfwidth,
l_float32 fract,
l_int32 direction)
{
l_int32 w, h, d, wpls, wpld, i, j, ival;
l_uint32 *datas, *datad;
l_uint32 *lines, *lines0, *lines1, *lines2, *lines3, *lines4, *lined;
l_float32 val, a[5];
PIX *pixd;
PROCNAME("pixUnsharpMaskingGray1D");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 8 || pixGetColormap(pixs) != NULL)
return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
if (fract <= 0.0 || halfwidth <= 0) {
L_WARNING("no sharpening requested; clone returned", procName);
return pixClone(pixs);
}
if (halfwidth != 1 && halfwidth != 2)
return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
/* Initialize pixd with pixels from pixs that will not be
* set when computing the sharpened values. */
pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
halfwidth, halfwidth);
datas = pixGetData(pixs);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs);
wpld = pixGetWpl(pixd);
if (halfwidth == 1) {
a[0] = -fract / 3.0;
a[1] = 1.0 + fract * 2.0 / 3.0;
a[2] = a[0];
}
else { /* halfwidth == 2 */
a[0] = -fract / 5.0;
a[1] = a[0];
a[2] = 1.0 + fract * 4.0 / 5.0;
a[3] = a[0];
a[4] = a[0];
}
if (direction == L_HORIZ) {
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
if (halfwidth == 1) {
for (j = 1; j < w - 1; j++) {
val = a[0] * GET_DATA_BYTE(lines, j - 1) +
a[1] * GET_DATA_BYTE(lines, j) +
a[2] * GET_DATA_BYTE(lines, j + 1);
ival = (l_int32)val;
ival = L_MAX(0, ival);
ival = L_MIN(255, ival);
SET_DATA_BYTE(lined, j, ival);
}
}
else { /* halfwidth == 2 */
for (j = 2; j < w - 2; j++) {
val = a[0] * GET_DATA_BYTE(lines, j - 2) +
a[1] * GET_DATA_BYTE(lines, j - 1) +
a[2] * GET_DATA_BYTE(lines, j) +
a[3] * GET_DATA_BYTE(lines, j + 1) +
a[4] * GET_DATA_BYTE(lines, j + 2);
ival = (l_int32)val;
ival = L_MAX(0, ival);
ival = L_MIN(255, ival);
SET_DATA_BYTE(lined, j, ival);
}
}
}
}
else { /* direction == L_VERT */
if (halfwidth == 1) {
for (i = 1; i < h - 1; i++) {
lines0 = datas + (i - 1) * wpls;
lines1 = datas + i * wpls;
lines2 = datas + (i + 1) * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
val = a[0] * GET_DATA_BYTE(lines0, j) +
a[1] * GET_DATA_BYTE(lines1, j) +
a[2] * GET_DATA_BYTE(lines2, j);
ival = (l_int32)val;
ival = L_MAX(0, ival);
ival = L_MIN(255, ival);
SET_DATA_BYTE(lined, j, ival);
}
}
}
else { /* halfwidth == 2 */
for (i = 2; i < h - 2; i++) {
lines0 = datas + (i - 2) * wpls;
lines1 = datas + (i - 1) * wpls;
lines2 = datas + i * wpls;
lines3 = datas + (i + 1) * wpls;
lines4 = datas + (i + 2) * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
val = a[0] * GET_DATA_BYTE(lines0, j) +
a[1] * GET_DATA_BYTE(lines1, j) +
a[2] * GET_DATA_BYTE(lines2, j) +
a[3] * GET_DATA_BYTE(lines3, j) +
a[4] * GET_DATA_BYTE(lines4, j);
ival = (l_int32)val;
ival = L_MAX(0, ival);
ival = L_MIN(255, ival);
SET_DATA_BYTE(lined, j, ival);
}
}
}
}
return pixd;
}
/*!
* pixUnsharpMaskingGray2D()
*
* Input: pixs (8 bpp; no colormap)
* halfwidth ("half-width" of smoothing filter: 1 or 2)
* fract (fraction of high frequency added to image)
* Return: pixd, or null on error
*
* Notes:
* (1) For halfwidth == 1, we implement the full sharpening filter
* directly. For halfwidth == 2, we implement the the lowpass
* filter separably and then compute the sharpening result locally.
*/
PIX *
pixUnsharpMaskingGray2D(PIX *pixs,
l_int32 halfwidth,
l_float32 fract)
{
l_int32 w, h, d, wpls, wpld, wplf, i, j, ival, sval;
l_uint32 *datas, *datad, *lines, *lines0, *lines1, *lines2, *lined;
l_float32 val, a[9];
l_float32 *dataf, *linef, *linef0, *linef1, *linef2, *linef3, *linef4;
PIX *pixd;
FPIX *fpix;
PROCNAME("pixUnsharpMaskingGray2D");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 8 || pixGetColormap(pixs) != NULL)
return (PIX *)ERROR_PTR("pixs not 8 bpp or has cmap", procName, NULL);
if (fract <= 0.0 || halfwidth <= 0) {
L_WARNING("no sharpening requested; clone returned", procName);
return pixClone(pixs);
}
if (halfwidth != 1 && halfwidth != 2)
return (PIX *)ERROR_PTR("halfwidth must be 1 or 2", procName, NULL);
pixd = pixCopyBorder(NULL, pixs, halfwidth, halfwidth,
halfwidth, halfwidth);
datad = pixGetData(pixd);
wpld = pixGetWpl(pixd);
datas = pixGetData(pixs);
wpls = pixGetWpl(pixs);
if (halfwidth == 1) {
for (i = 0; i < 9; i++)
a[i] = -fract / 9.0;
a[4] = 1.0 + fract * 8.0 / 9.0;
for (i = 1; i < h - 1; i++) {
lines0 = datas + (i - 1) * wpls;
lines1 = datas + i * wpls;
lines2 = datas + (i + 1) * wpls;
lined = datad + i * wpld;
for (j = 1; j < w - 1; j++) {
val = a[0] * GET_DATA_BYTE(lines0, j - 1) +
a[1] * GET_DATA_BYTE(lines0, j) +
a[2] * GET_DATA_BYTE(lines0, j + 1) +
a[3] * GET_DATA_BYTE(lines1, j - 1) +
a[4] * GET_DATA_BYTE(lines1, j) +
a[5] * GET_DATA_BYTE(lines1, j + 1) +
a[6] * GET_DATA_BYTE(lines2, j - 1) +
a[7] * GET_DATA_BYTE(lines2, j) +
a[8] * GET_DATA_BYTE(lines2, j + 1);
ival = (l_int32)(val + 0.5);
ival = L_MAX(0, ival);
ival = L_MIN(255, ival);
SET_DATA_BYTE(lined, j, ival);
}
}
return pixd;
}
/* For halfwidth == 2, do the low pass separably. Store
* the result of horizontal smoothing in an intermediate fpix. */
fpix = fpixCreate(w, h);
dataf = fpixGetData(fpix);
wplf = fpixGetWpl(fpix);
for (i = 2; i < h - 2; i++) {
lines = datas + i * wpls;
linef = dataf + i * wplf;
for (j = 2; j < w - 2; j++) {
val = GET_DATA_BYTE(lines, j - 2) +
GET_DATA_BYTE(lines, j - 1) +
GET_DATA_BYTE(lines, j) +
GET_DATA_BYTE(lines, j + 1) +
GET_DATA_BYTE(lines, j + 2);
linef[j] = val;
}
}
/* Do vertical smoothing to finish the low-pass filter.
* At each pixel, if L is the lowpass value, I is the
* src pixel value and f is the fraction of highpass to
* be added to I, then the highpass filter value is
* H = I - L
* and the new sharpened value is
* N = I + f * H.
*/
for (i = 2; i < h - 2; i++) {
linef0 = dataf + (i - 2) * wplf;
linef1 = dataf + (i - 1) * wplf;
linef2 = dataf + i * wplf;
linef3 = dataf + (i + 1) * wplf;
linef4 = dataf + (i + 2) * wplf;
lined = datad + i * wpld;
lines = datas + i * wpls;
for (j = 2; j < w - 2; j++) {
val = 0.04 * (linef0[j] + linef1[j] + linef2[j] +
linef3[j] + linef4[j]); /* L: lowpass filter value */
sval = GET_DATA_BYTE(lines, j); /* I: source pixel */
ival = (l_int32)(sval + fract * (sval - val) + 0.5);
ival = L_MAX(0, ival);
ival = L_MIN(255, ival);
SET_DATA_BYTE(lined, j, ival);
}
}
fpixDestroy(&fpix);
return pixd;
}
/*-----------------------------------------------------------------------*
* Hue and saturation modification *
*-----------------------------------------------------------------------*/
/*!
* pixModifyHue()
*
* Input: pixd (<optional> can be null, existing or equal to pixs)
* pixs (32 bpp rgb)
* fract (between -1.0 and 1.0)
* Return: pixd, or null on error
*
* Notes:
* (1) Use fract > 0.0 to increase hue value; < 0.0 to decrease it.
* 1.0 (or -1.0) represents a 360 degree rotation; i.e., no change.
*/
PIX *
pixModifyHue(PIX *pixd,
PIX *pixs,
l_float32 fract)
{
l_int32 w, h, d, i, j, wpl, delhue;
l_int32 rval, gval, bval, hval, sval, vval;
l_uint32 *data, *line;
PROCNAME("pixModifyHue");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 32)
return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
if (L_ABS(fract) > 1.0)
return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", procName, NULL);
pixd = pixCopy(pixd, pixs);
delhue = (l_int32)(240 * fract);
if (delhue == 0 || delhue == 240 || delhue == -240) {
L_WARNING("no change requested in hue", procName);
return pixd;
}
if (delhue < 0)
delhue += 240;
data = pixGetData(pixd);
wpl = pixGetWpl(pixd);
for (i = 0; i < h; i++) {
line = data + i * wpl;
for (j = 0; j < w; j++) {
extractRGBValues(line[j], &rval, &gval, &bval);
convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
hval = (hval + delhue) % 240;
convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
composeRGBPixel(rval, gval, bval, line + j);
}
}
return pixd;
}
/*!
* pixModifySaturation()
*
* Input: pixd (<optional> can be null, existing or equal to pixs)
* pixs (32 bpp rgb)
* fract (between -1.0 and 1.0)
* Return: pixd, or null on error
*
* Notes:
* (1) If fract > 0.0, it gives the fraction that the pixel
* saturation is moved from its initial value toward 255.
* If fract < 0.0, it gives the fraction that the pixel
* saturation is moved from its initial value toward 0.
* The limiting values for fract = -1.0 (1.0) thus set the
* saturation to 0 (255).
*/
PIX *
pixModifySaturation(PIX *pixd,
PIX *pixs,
l_float32 fract)
{
l_int32 w, h, d, i, j, wpl;
l_int32 rval, gval, bval, hval, sval, vval;
l_uint32 *data, *line;
PROCNAME("pixModifySaturation");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 32)
return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
if (L_ABS(fract) > 1.0)
return (PIX *)ERROR_PTR("fract not in [-1.0 ... 1.0]", procName, NULL);
pixd = pixCopy(pixd, pixs);
if (fract == 0.0) {
L_WARNING("no change requested in saturation", procName);
return pixd;
}
data = pixGetData(pixd);
wpl = pixGetWpl(pixd);
for (i = 0; i < h; i++) {
line = data + i * wpl;
for (j = 0; j < w; j++) {
extractRGBValues(line[j], &rval, &gval, &bval);
convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
if (fract < 0.0)
sval = (l_int32)(sval * (1.0 + fract));
else
sval = (l_int32)(sval + fract * (255 - sval));
convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
composeRGBPixel(rval, gval, bval, line + j);
}
}
return pixd;
}
/*!
* pixMeasureSaturation()
*
* Input: pixs (32 bpp rgb)
* factor (subsampling factor; integer >= 1)
* &sat (<return> average saturation)
* Return: pixd, or null on error
*/
l_int32
pixMeasureSaturation(PIX *pixs,
l_int32 factor,
l_float32 *psat)
{
l_int32 w, h, d, i, j, wpl, sum, count;
l_int32 rval, gval, bval, hval, sval, vval;
l_uint32 *data, *line;
PROCNAME("pixMeasureSaturation");
if (!psat)
return ERROR_INT("pixs not defined", procName, 1);
*psat = 0.0;
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 32)
return ERROR_INT("pixs not 32 bpp", procName, 1);
if (factor < 1)
return ERROR_INT("subsampling factor < 1", procName, 1);
data = pixGetData(pixs);
wpl = pixGetWpl(pixs);
for (i = 0, sum = 0, count = 0; i < h; i += factor) {
line = data + i * wpl;
for (j = 0; j < w; j += factor) {
extractRGBValues(line[j], &rval, &gval, &bval);
convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
sum += sval;
count++;
}
}
if (count > 0)
*psat = (l_float32)sum / (l_float32)count;
return 0;
}
/*-----------------------------------------------------------------------*
* General multiplicative constant color transform *
*-----------------------------------------------------------------------*/
/*
* pixMultConstantColor()
*
* Input: pixs (colormapped or rgb)
* rfact, gfact, bfact (multiplicative factors on each component)
* Return: pixd (colormapped or rgb, with colors scaled), or null on error
*
* Notes:
* (1) rfact, gfact and bfact can only have non-negative values.
* They can be greater than 1.0. All transformed component
* values are clipped to the interval [0, 255].
* (2) For multiplication with a general 3x3 matrix of constants,
* use pixMultMatrixColor().
*/
PIX *
pixMultConstantColor(PIX *pixs,
l_float32 rfact,
l_float32 gfact,
l_float32 bfact)
{
l_int32 i, j, w, h, d, wpls, wpld;
l_int32 ncolors, rval, gval, bval, nrval, ngval, nbval;
l_uint32 nval;
l_uint32 *datas, *datad, *lines, *lined;
PIX *pixd;
PIXCMAP *cmap;
PROCNAME("pixMultConstantColor");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
pixGetDimensions(pixs, &w, &h, &d);
cmap = pixGetColormap(pixs);
if (!cmap && d != 32)
return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
rfact = L_MAX(0.0, rfact);
gfact = L_MAX(0.0, gfact);
bfact = L_MAX(0.0, bfact);
if (cmap) {
if ((pixd = pixCopy(NULL, pixs)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
cmap = pixGetColormap(pixd);
ncolors = pixcmapGetCount(cmap);
for (i = 0; i < ncolors; i++) {
pixcmapGetColor(cmap, i, &rval, &gval, &bval);
nrval = (l_int32)(rfact * rval);
ngval = (l_int32)(gfact * gval);
nbval = (l_int32)(bfact * bval);
nrval = L_MIN(255, nrval);
ngval = L_MIN(255, ngval);
nbval = L_MIN(255, nbval);
pixcmapResetColor(cmap, i, nrval, ngval, nbval);
}
return pixd;
}
if ((pixd = pixCreateTemplateNoInit(pixs)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
datas = pixGetData(pixs);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs);
wpld = pixGetWpl(pixd);
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
extractRGBValues(lines[j], &rval, &gval, &bval);
nrval = (l_int32)(rfact * rval);
ngval = (l_int32)(gfact * gval);
nbval = (l_int32)(bfact * bval);
nrval = L_MIN(255, nrval);
ngval = L_MIN(255, ngval);
nbval = L_MIN(255, nbval);
composeRGBPixel(nrval, ngval, nbval, &nval);
*(lined + j) = nval;
}
}
return pixd;
}
/*
* pixMultMatrixColor()
*
* Input: pixs (colormapped or rgb)
* kernel (3x3 matrix of floats)
* Return: pixd (colormapped or rgb), or null on error
*
* Notes:
* (1) The kernel is a data structure used mostly for floating point
* convolution. Here it is a 3x3 matrix of floats that are used
* to transform the pixel values by matrix multiplication:
* nrval = a[0,0] * rval + a[0,1] * gval + a[0,2] * bval
* ngval = a[1,0] * rval + a[1,1] * gval + a[1,2] * bval
* nbval = a[2,0] * rval + a[2,1] * gval + a[2,2] * bval
* (2) The matrix can be generated in several ways.
* See kernel.c for details. Here are two of them:
* (a) kel = kernelCreate(3, 3);
* kernelSetElement(kel, 0, 0, val00);
* kernelSetElement(kel, 0, 1, val01);
* ...
* (b) from a static string; e.g.,:
* const char *kdata = " 0.6 0.3 -0.2 "
* " 0.1 1.2 0.4 "
* " -0.4 0.2 0.9 ";
* kel = kernelCreateFromString(3, 3, 0, 0, kdata);
* (3) For the special case where the matrix is diagonal, it is easier
* to use pixMultConstantColor().
* (4) Matrix entries can have positive and negative values, and can
* be larger than 1.0. All transformed component values
* are clipped to [0, 255].
*/
PIX *
pixMultMatrixColor(PIX *pixs,
L_KERNEL *kel)
{
l_int32 i, j, index, kw, kh, w, h, d, wpls, wpld;
l_int32 ncolors, rval, gval, bval, nrval, ngval, nbval;
l_uint32 nval;
l_uint32 *datas, *datad, *lines, *lined;
l_float32 v[9]; /* use linear array for convenience */
PIX *pixd;
PIXCMAP *cmap;
PROCNAME("pixMultMatrixColor");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (!kel)
return (PIX *)ERROR_PTR("kel not defined", procName, NULL);
kernelGetParameters(kel, &kw, &kh, NULL, NULL);
if (kw != 3 || kh != 3)
return (PIX *)ERROR_PTR("matrix not 3x3", procName, NULL);
pixGetDimensions(pixs, &w, &h, &d);
cmap = pixGetColormap(pixs);
if (!cmap && d != 32)
return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
for (i = 0, index = 0; i < 3; i++)
for (j = 0; j < 3; j++, index++)
kernelGetElement(kel, i, j, v + index);
if (cmap) {
if ((pixd = pixCopy(NULL, pixs)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
cmap = pixGetColormap(pixd);
ncolors = pixcmapGetCount(cmap);
for (i = 0; i < ncolors; i++) {
pixcmapGetColor(cmap, i, &rval, &gval, &bval);
nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
nrval = L_MAX(0, L_MIN(255, nrval));
ngval = L_MAX(0, L_MIN(255, ngval));
nbval = L_MAX(0, L_MIN(255, nbval));
pixcmapResetColor(cmap, i, nrval, ngval, nbval);
}
return pixd;
}
if ((pixd = pixCreateTemplateNoInit(pixs)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
datas = pixGetData(pixs);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs);
wpld = pixGetWpl(pixd);
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
extractRGBValues(lines[j], &rval, &gval, &bval);
nrval = (l_int32)(v[0] * rval + v[1] * gval + v[2] * bval);
ngval = (l_int32)(v[3] * rval + v[4] * gval + v[5] * bval);
nbval = (l_int32)(v[6] * rval + v[7] * gval + v[8] * bval);
nrval = L_MAX(0, L_MIN(255, nrval));
ngval = L_MAX(0, L_MIN(255, ngval));
nbval = L_MAX(0, L_MIN(255, nbval));
composeRGBPixel(nrval, ngval, nbval, &nval);
*(lined + j) = nval;
}
}
return pixd;
}
/*-------------------------------------------------------------*
* Half-edge by bandpass *
*-------------------------------------------------------------*/
/*!
* pixHalfEdgeByBandpass()
*
* Input: pixs (8 bpp gray or 32 bpp rgb)
* sm1h, sm1v ("half-widths" of smoothing filter sm1)
* sm2h, sm2v ("half-widths" of smoothing filter sm2)
* (require sm2 != sm1)
* Return: pixd, or null on error
*
* Notes:
* (1) We use symmetric smoothing filters of odd dimension,
* typically use 3, 5, 7, etc. The smoothing parameters
* for these are 1, 2, 3, etc. The filter size is related
* to the smoothing parameter by
* size = 2 * smoothing + 1
* (2) Because we take the difference of two lowpass filters,
* this is actually a bandpass filter.
* (3) We allow both filters to be anisotropic.
* (4) Consider either the h or v component of the 2 filters.
* Depending on whether sm1 > sm2 or sm2 > sm1, we get
* different halves of the smoothed gradients (or "edges").
* This difference of smoothed signals looks more like
* a second derivative of a transition, which we rectify
* by not allowing the signal to go below zero. If sm1 < sm2,
* the sm2 transition is broader, so the difference between
* sm1 and sm2 signals is positive on the upper half of
* the transition. Likewise, if sm1 > sm2, the sm1 - sm2
* signal difference is positive on the lower half of
* the transition.
*/
PIX *
pixHalfEdgeByBandpass(PIX *pixs,
l_int32 sm1h,
l_int32 sm1v,
l_int32 sm2h,
l_int32 sm2v)
{
l_int32 d;
PIX *pixg, *pixacc, *pixc1, *pixc2;
PROCNAME("pixHalfEdgeByBandpass");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (sm1h == sm2h && sm1v == sm2v)
return (PIX *)ERROR_PTR("sm2 = sm1", procName, NULL);
d = pixGetDepth(pixs);
if (d != 8 && d != 32)
return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
if (d == 32)
pixg = pixConvertRGBToLuminance(pixs);
else /* d == 8 */
pixg = pixClone(pixs);
/* Make a convolution accumulator and use it twice */
if ((pixacc = pixBlockconvAccum(pixg)) == NULL)
return (PIX *)ERROR_PTR("pixacc not made", procName, NULL);
if ((pixc1 = pixBlockconvGray(pixg, pixacc, sm1h, sm1v)) == NULL)
return (PIX *)ERROR_PTR("pixc1 not made", procName, NULL);
if ((pixc2 = pixBlockconvGray(pixg, pixacc, sm2h, sm2v)) == NULL)
return (PIX *)ERROR_PTR("pixc2 not made", procName, NULL);
pixDestroy(&pixacc);
/* Compute the half-edge using pixc1 - pixc2. */
pixSubtractGray(pixc1, pixc1, pixc2);
pixDestroy(&pixg);
pixDestroy(&pixc2);
return pixc1;
}