/*====================================================================*
- 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.
*====================================================================*/
/*
* pixarith.c
*
* One-image grayscale arithmetic operations (8, 16, 32 bpp)
* l_int32 pixAddConstantGray()
* l_int32 pixMultConstantGray()
*
* Two-image grayscale arithmetic operations (8, 16, 32 bpp)
* PIX *pixAddGray()
* PIX *pixSubtractGray()
*
* Grayscale threshold operation (8, 16, 32 bpp)
* PIX *pixThresholdToValue()
*
* Image accumulator arithmetic operations
* PIX *pixInitAccumulate()
* PIX *pixFinalAccumulate()
* PIX *pixFinalAccumulateThreshold()
* l_int32 pixAccumulate()
* l_int32 pixMultConstAccumulate()
*
* Absolute value of difference
* PIX *pixAbsDifference()
*
* Two-image min and max operations (8 and 16 bpp)
* PIX *pixMinOrMax()
*
* Scale pix for maximum dynamic range in 8 bpp image:
* PIX *pixMaxDynamicRange()
*
* Log base2 lookup
* l_float32 *makeLogBase2Tab()
* l_float32 getLogBase2()
*
* The image accumulator operations are used when you expect
* overflow from 8 bits on intermediate results. For example,
* you might want a tophat contrast operator which is
* 3*I - opening(I,S) - closing(I,S)
* To use these operations, first use the init to generate
* a 16 bpp image, use the accumulate to add or subtract 8 bpp
* images from that, or the multiply constant to multiply
* by a small constant (much less than 256 -- we don't want
* overflow from the 16 bit images!), and when you're finished
* use final to bring the result back to 8 bpp, clipped
* if necessary. There is also a divide function, which
* can be used to divide one image by another, scaling the
* result for maximum dynamic range, and giving back the
* 8 bpp result.
*
* A simpler interface to the arithmetic operations is
* provided in pixacc.c.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "allheaders.h"
/*-------------------------------------------------------------*
* One-image grayscale arithmetic operations *
*-------------------------------------------------------------*/
/*!
* pixAddConstantGray()
*
* Input: pixs (8, 16 or 32 bpp)
* val (amount to add to each pixel)
* Return: 0 if OK, 1 on error
*
* Notes:
* (1) In-place operation.
* (2) No clipping for 32 bpp.
* (3) For 8 and 16 bpp, if val > 0 the result is clipped
* to 0xff and 0xffff, rsp.
* (4) For 8 and 16 bpp, if val < 0 the result is clipped to 0.
*/
l_int32
pixAddConstantGray(PIX *pixs,
l_int32 val)
{
l_int32 w, h, d, wpl;
l_uint32 *data;
PROCNAME("pixAddConstantGray");
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 8 && d != 16 && d != 32)
return ERROR_INT("pixs not 8, 16 or 32 bpp", procName, 1);
data = pixGetData(pixs);
wpl = pixGetWpl(pixs);
addConstantGrayLow(data, w, h, d, wpl, val);
return 0;
}
/*!
* pixMultConstantGray()
*
* Input: pixs (8, 16 or 32 bpp)
* val (>= 0.0; amount to multiply by each pixel)
* Return: 0 if OK, 1 on error
*
* Notes:
* (1) In-place operation; val must be >= 0.
* (2) No clipping for 32 bpp.
* (3) For 8 and 16 bpp, the result is clipped to 0xff and 0xffff, rsp.
*/
l_int32
pixMultConstantGray(PIX *pixs,
l_float32 val)
{
l_int32 w, h, d, wpl;
l_uint32 *data;
PROCNAME("pixMultConstantGray");
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
pixGetDimensions(pixs, &w, &h, &d);
if (d != 8 && d != 16 && d != 32)
return ERROR_INT("pixs not 8, 16 or 32 bpp", procName, 1);
if (val < 0.0)
return ERROR_INT("val < 0.0", procName, 1);
data = pixGetData(pixs);
wpl = pixGetWpl(pixs);
multConstantGrayLow(data, w, h, d, wpl, val);
return 0;
}
/*-------------------------------------------------------------*
* Two-image grayscale arithmetic ops *
*-------------------------------------------------------------*/
/*!
* pixAddGray()
*
* Input: pixd (<optional>; this can be null, equal to pixs1, or
* different from pixs1)
* pixs1 (can be == to pixd)
* pixs2
* Return: pixd always
*
* Notes:
* (1) Arithmetic addition of two 8, 16 or 32 bpp images.
* (2) For 8 and 16 bpp, we do explicit clipping to 0xff and 0xffff,
* respectively.
* (3) Alignment is to UL corner.
* (4) There are 3 cases. The result can go to a new dest,
* in-place to pixs1, or to an existing input dest:
* * pixd == null: (src1 + src2) --> new pixd
* * pixd == pixs1: (src1 + src2) --> src1 (in-place)
* * pixd != pixs1: (src1 + src2) --> input pixd
* (5) pixs2 must be different from both pixd and pixs1.
*/
PIX *
pixAddGray(PIX *pixd,
PIX *pixs1,
PIX *pixs2)
{
l_int32 d, ws, hs, w, h, wpls, wpld;
l_uint32 *datas, *datad;
PROCNAME("pixAddGray");
if (!pixs1)
return (PIX *)ERROR_PTR("pixs1 not defined", procName, pixd);
if (!pixs2)
return (PIX *)ERROR_PTR("pixs2 not defined", procName, pixd);
if (pixs2 == pixs1)
return (PIX *)ERROR_PTR("pixs2 and pixs1 must differ", procName, pixd);
if (pixs2 == pixd)
return (PIX *)ERROR_PTR("pixs2 and pixd must differ", procName, pixd);
d = pixGetDepth(pixs1);
if (d != 8 && d != 16 && d != 32)
return (PIX *)ERROR_PTR("pix are not 8, 16 or 32 bpp", procName, pixd);
if (pixGetDepth(pixs2) != d)
return (PIX *)ERROR_PTR("depths differ (pixs1, pixs2)", procName, pixd);
if (pixd && (pixGetDepth(pixd) != d))
return (PIX *)ERROR_PTR("depths differ (pixs1, pixd)", procName, pixd);
if (!pixSizesEqual(pixs1, pixs2))
L_WARNING("pixs1 and pixs2 not equal in size", procName);
if (pixd && !pixSizesEqual(pixs1, pixd))
L_WARNING("pixs1 and pixd not equal in size", procName);
if (pixs1 != pixd)
pixd = pixCopy(pixd, pixs1);
/* pixd + pixs2 ==> pixd */
datas = pixGetData(pixs2);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs2);
wpld = pixGetWpl(pixd);
pixGetDimensions(pixs2, &ws, &hs, NULL);
pixGetDimensions(pixd, &w, &h, NULL);
w = L_MIN(ws, w);
h = L_MIN(hs, h);
addGrayLow(datad, w, h, d, wpld, datas, wpls);
return pixd;
}
/*!
* pixSubtractGray()
*
* Input: pixd (<optional>; this can be null, equal to pixs1, or
* different from pixs1)
* pixs1 (can be == to pixd)
* pixs2
* Return: pixd always
*
* Notes:
* (1) Arithmetic subtraction of two 8, 16 or 32 bpp images.
* (2) Source pixs2 is always subtracted from source pixs1.
* (3) Do explicit clipping to 0.
* (4) Alignment is to UL corner.
* (5) There are 3 cases. The result can go to a new dest,
* in-place to pixs1, or to an existing input dest:
* (a) pixd == null (src1 - src2) --> new pixd
* (b) pixd == pixs1 (src1 - src2) --> src1 (in-place)
* (d) pixd != pixs1 (src1 - src2) --> input pixd
* (6) pixs2 must be different from both pixd and pixs1.
*/
PIX *
pixSubtractGray(PIX *pixd,
PIX *pixs1,
PIX *pixs2)
{
l_int32 w, h, ws, hs, d, wpls, wpld;
l_uint32 *datas, *datad;
PROCNAME("pixSubtractGray");
if (!pixs1)
return (PIX *)ERROR_PTR("pixs1 not defined", procName, pixd);
if (!pixs2)
return (PIX *)ERROR_PTR("pixs2 not defined", procName, pixd);
if (pixs2 == pixs1)
return (PIX *)ERROR_PTR("pixs2 and pixs1 must differ", procName, pixd);
if (pixs2 == pixd)
return (PIX *)ERROR_PTR("pixs2 and pixd must differ", procName, pixd);
d = pixGetDepth(pixs1);
if (d != 8 && d != 16 && d != 32)
return (PIX *)ERROR_PTR("pix are not 8, 16 or 32 bpp", procName, pixd);
if (pixGetDepth(pixs2) != d)
return (PIX *)ERROR_PTR("depths differ (pixs1, pixs2)", procName, pixd);
if (pixd && (pixGetDepth(pixd) != d))
return (PIX *)ERROR_PTR("depths differ (pixs1, pixd)", procName, pixd);
if (!pixSizesEqual(pixs1, pixs2))
L_WARNING("pixs1 and pixs2 not equal in size", procName);
if (pixd && !pixSizesEqual(pixs1, pixd))
L_WARNING("pixs1 and pixd not equal in size", procName);
if (pixs1 != pixd)
pixd = pixCopy(pixd, pixs1);
/* pixd - pixs2 ==> pixd */
datas = pixGetData(pixs2);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs2);
wpld = pixGetWpl(pixd);
pixGetDimensions(pixs2, &ws, &hs, NULL);
pixGetDimensions(pixd, &w, &h, NULL);
w = L_MIN(ws, w);
h = L_MIN(hs, h);
subtractGrayLow(datad, w, h, d, wpld, datas, wpls);
return pixd;
}
/*-------------------------------------------------------------*
* Grayscale threshold operation *
*-------------------------------------------------------------*/
/*!
* pixThresholdToValue()
*
* Input: pixd (<optional>; if not null, must be equal to pixs)
* pixs (8, 16, 32 bpp)
* threshval
* setval
* Return: pixd always
*
* Notes:
* - operation can be in-place (pixs == pixd) or to a new pixd
* - if setval > threshval, sets pixels with a value >= threshval to setval
* - if setval < threshval, sets pixels with a value <= threshval to setval
* - if setval == threshval, no-op
*/
PIX *
pixThresholdToValue(PIX *pixd,
PIX *pixs,
l_int32 threshval,
l_int32 setval)
{
l_int32 w, h, d, wpld;
l_uint32 *datad;
PROCNAME("pixThresholdToValue");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
d = pixGetDepth(pixs);
if (d != 8 && d != 16 && d != 32)
return (PIX *)ERROR_PTR("pixs not 8, 16 or 32 bpp", procName, pixd);
if (pixd && (pixs != pixd))
return (PIX *)ERROR_PTR("pixd exists and is not pixs", procName, pixd);
if (threshval < 0 || setval < 0)
return (PIX *)ERROR_PTR("threshval & setval not < 0", procName, pixd);
if (d == 8 && setval > 255)
return (PIX *)ERROR_PTR("setval > 255 for 8 bpp", procName, pixd);
if (d == 16 && setval > 0xffff)
return (PIX *)ERROR_PTR("setval > 0xffff for 16 bpp", procName, pixd);
if (!pixd)
pixd = pixCopy(NULL, pixs);
if (setval == threshval) {
L_WARNING("setval == threshval; no operation", procName);
return pixd;
}
datad = pixGetData(pixd);
pixGetDimensions(pixd, &w, &h, NULL);
wpld = pixGetWpl(pixd);
thresholdToValueLow(datad, w, h, d, wpld, threshval, setval);
return pixd;
}
/*-------------------------------------------------------------*
* Image accumulator arithmetic operations *
*-------------------------------------------------------------*/
/*!
* pixInitAccumulate()
*
* Input: w, h (of accumulate array)
* offset (initialize the 32 bpp to have this
* value; not more than 0x40000000)
* Return: pixd (32 bpp), or null on error
*
* Notes:
* (1) The offset must be >= 0.
* (2) The offset is used so that we can do arithmetic
* with negative number results on l_uint32 data; it
* prevents the l_uint32 data from going negative.
* (3) Because we use l_int32 intermediate data results,
* these should never exceed the max of l_int32 (0x7fffffff).
* We do not permit the offset to be above 0x40000000,
* which is half way between 0 and the max of l_int32.
* (4) The same offset should be used for initialization,
* multiplication by a constant, and final extraction!
* (5) If you're only adding positive values, offset can be 0.
*/
PIX *
pixInitAccumulate(l_int32 w,
l_int32 h,
l_uint32 offset)
{
PIX *pixd;
PROCNAME("pixInitAccumulate");
if ((pixd = pixCreate(w, h, 32)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
if (offset > 0x40000000)
offset = 0x40000000;
pixSetAllArbitrary(pixd, offset);
return pixd;
}
/*!
* pixFinalAccumulate()
*
* Input: pixs (32 bpp)
* offset (same as used for initialization)
* depth (8, 16 or 32 bpp, of destination)
* Return: pixd (8, 16 or 32 bpp), or null on error
*
* Notes:
* (1) The offset must be >= 0 and should not exceed 0x40000000.
* (2) The offset is subtracted from the src 32 bpp image
* (3) For 8 bpp dest, the result is clipped to [0, 0xff]
* (4) For 16 bpp dest, the result is clipped to [0, 0xffff]
*/
PIX *
pixFinalAccumulate(PIX *pixs,
l_uint32 offset,
l_int32 depth)
{
l_int32 w, h, wpls, wpld;
l_uint32 *datas, *datad;
PIX *pixd;
PROCNAME("pixFinalAccumulate");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (pixGetDepth(pixs) != 32)
return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
if (depth != 8 && depth != 16 && depth != 32)
return (PIX *)ERROR_PTR("dest depth not 8, 16, 32 bpp", procName, NULL);
if (offset > 0x40000000)
offset = 0x40000000;
pixGetDimensions(pixs, &w, &h, NULL);
if ((pixd = pixCreate(w, h, depth)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
pixCopyResolution(pixd, pixs); /* but how did pixs get it initially? */
datas = pixGetData(pixs);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs);
wpld = pixGetWpl(pixd);
finalAccumulateLow(datad, w, h, depth, wpld, datas, wpls, offset);
return pixd;
}
/*!
* pixFinalAccumulateThreshold()
*
* Input: pixs (32 bpp)
* offset (same as used for initialization)
* threshold (values less than this are set in the destination)
* Return: pixd (1 bpp), or null on error
*
* Notes:
* (1) The offset must be >= 0 and should not exceed 0x40000000.
* (2) The offset is subtracted from the src 32 bpp image
*/
PIX *
pixFinalAccumulateThreshold(PIX *pixs,
l_uint32 offset,
l_uint32 threshold)
{
l_int32 w, h, wpls, wpld;
l_uint32 *datas, *datad;
PIX *pixd;
PROCNAME("pixFinalAccumulateThreshold");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
if (pixGetDepth(pixs) != 32)
return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
if (offset > 0x40000000)
offset = 0x40000000;
pixGetDimensions(pixs, &w, &h, NULL);
if ((pixd = pixCreate(w, h, 1)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
pixCopyResolution(pixd, pixs); /* but how did pixs get it initially? */
datas = pixGetData(pixs);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs);
wpld = pixGetWpl(pixd);
finalAccumulateThreshLow(datad, w, h, wpld, datas, wpls, offset, threshold);
return pixd;
}
/*!
* pixAccumulate()
*
* Input: pixd (32 bpp)
* pixs (1, 8, 16 or 32 bpp)
* op (L_ARITH_ADD or L_ARITH_SUBTRACT)
* Return: 0 if OK; 1 on error
*
* Notes:
* (1) This adds or subtracts each pixs value from pixd.
* (2) This clips to the minimum of pixs and pixd, so they
* do not need to be the same size.
* (3) The alignment is to the origin (UL corner) of pixs & pixd.
*/
l_int32
pixAccumulate(PIX *pixd,
PIX *pixs,
l_int32 op)
{
l_int32 w, h, d, wd, hd, wpls, wpld;
l_uint32 *datas, *datad;
PROCNAME("pixAccumulate");
if (!pixd || (pixGetDepth(pixd) != 32))
return ERROR_INT("pixd not defined or not 32 bpp", procName, 1);
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
d = pixGetDepth(pixs);
if (d != 1 && d != 8 && d != 16 && d != 32)
return ERROR_INT("pixs not 1, 8, 16 or 32 bpp", procName, 1);
if (op != L_ARITH_ADD && op != L_ARITH_SUBTRACT)
return ERROR_INT("op must be in {L_ARITH_ADD, L_ARITH_SUBTRACT}",
procName, 1);
datas = pixGetData(pixs);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs);
wpld = pixGetWpl(pixd);
pixGetDimensions(pixs, &w, &h, NULL);
pixGetDimensions(pixd, &wd, &hd, NULL);
w = L_MIN(w, wd);
h = L_MIN(h, hd);
accumulateLow(datad, w, h, wpld, datas, d, wpls, op);
return 0;
}
/*!
* pixMultConstAccumulate()
*
* Input: pixs (32 bpp)
* factor
* offset (same as used for initialization)
* Return: 0 if OK; 1 on error
*
* Notes:
* (1) The offset must be >= 0 and should not exceed 0x40000000.
* (2) This multiplies each pixel, relative to offset, by the input factor
* (3) The result is returned with the offset back in place.
*/
l_int32
pixMultConstAccumulate(PIX *pixs,
l_float32 factor,
l_uint32 offset)
{
l_int32 w, h, wpl;
l_uint32 *data;
PROCNAME("pixMultConstAccumulate");
if (!pixs)
return ERROR_INT("pixs not defined", procName, 1);
if (pixGetDepth(pixs) != 32)
return ERROR_INT("pixs not 32 bpp", procName, 1);
if (offset > 0x40000000)
offset = 0x40000000;
pixGetDimensions(pixs, &w, &h, NULL);
data = pixGetData(pixs);
wpl = pixGetWpl(pixs);
multConstAccumulateLow(data, w, h, wpl, factor, offset);
return 0;
}
/*-----------------------------------------------------------------------*
* Absolute value of difference *
*-----------------------------------------------------------------------*/
/*!
* pixAbsDifference()
*
* Input: pixs1, pixs2 (both either 8 or 16 bpp gray, or 32 bpp RGB)
* Return: pixd, or null on error
*
* Notes:
* (1) The depth of pixs1 and pixs2 must be equal.
* (2) Clips computation to the min size, aligning the UL corners
* (3) For 8 and 16 bpp, assumes one gray component.
* (4) For 32 bpp, assumes 3 color components, and ignores the
* LSB of each word (the alpha channel)
* (5) Computes the absolute value of the difference between
* each component value.
*/
PIX *
pixAbsDifference(PIX *pixs1,
PIX *pixs2)
{
l_int32 w, h, w2, h2, d, wpls, wpld;
l_uint32 *datas1, *datas2, *datad;
PIX *pixd;
PROCNAME("pixAbsDifference");
if (!pixs1)
return (PIX *)ERROR_PTR("pixs1 not defined", procName, NULL);
if (!pixs2)
return (PIX *)ERROR_PTR("pixs2 not defined", procName, NULL);
d = pixGetDepth(pixs1);
if (d != pixGetDepth(pixs2))
return (PIX *)ERROR_PTR("src1 and src2 depths unequal", procName, NULL);
if (d != 8 && d != 16 && d != 32)
return (PIX *)ERROR_PTR("depths not in {8, 16, 32}", procName, NULL);
pixGetDimensions(pixs1, &w, &h, NULL);
pixGetDimensions(pixs2, &w2, &h2, NULL);
w = L_MIN(w, w2);
h = L_MIN(h, h2);
if ((pixd = pixCreate(w, h, d)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
pixCopyResolution(pixd, pixs1);
datas1 = pixGetData(pixs1);
datas2 = pixGetData(pixs2);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs1);
wpld = pixGetWpl(pixd);
absDifferenceLow(datad, w, h, wpld, datas1, datas2, d, wpls);
return pixd;
}
/*-----------------------------------------------------------------------*
* Two-image min and max operations (8 and 16 bpp) *
*-----------------------------------------------------------------------*/
/*!
* pixMinOrMax()
*
* Input: pixd (<optional> destination: this can be null,
* equal to pixs1, or different from pixs1)
* pixs1 (can be == to pixd)
* pixs2
* type (L_CHOOSE_MIN, L_CHOOSE_MAX)
* Return: pixd always
*
* Notes:
* (1) This gives the min or max of two images.
* (2) The depth can be 8 or 16 bpp.
* (3) There are 3 cases:
* - if pixd == null, Min(src1, src2) --> new pixd
* - if pixd == pixs1, Min(src1, src2) --> src1 (in-place)
* - if pixd != pixs1, Min(src1, src2) --> input pixd
*/
PIX *
pixMinOrMax(PIX *pixd,
PIX *pixs1,
PIX *pixs2,
l_int32 type)
{
l_int32 d, ws, hs, w, h, wpls, wpld, i, j;
l_int32 vals, vald, val;
l_uint32 *datas, *datad, *lines, *lined;
PROCNAME("pixMinOrMax");
if (!pixs1)
return (PIX *)ERROR_PTR("pixs1 not defined", procName, pixd);
if (!pixs2)
return (PIX *)ERROR_PTR("pixs2 not defined", procName, pixd);
if (pixs1 == pixs2)
return (PIX *)ERROR_PTR("pixs1 and pixs2 must differ", procName, pixd);
if (type != L_CHOOSE_MIN && type != L_CHOOSE_MAX)
return (PIX *)ERROR_PTR("invalid type", procName, pixd);
d = pixGetDepth(pixs1);
if (pixGetDepth(pixs2) != d)
return (PIX *)ERROR_PTR("depths unequal", procName, pixd);
if (d != 8 && d != 16)
return (PIX *)ERROR_PTR("depth not 8 or 16 bpp", procName, pixd);
if (pixs1 != pixd)
pixd = pixCopy(pixd, pixs1);
pixGetDimensions(pixs2, &ws, &hs, NULL);
pixGetDimensions(pixd, &w, &h, NULL);
w = L_MIN(w, ws);
h = L_MIN(h, hs);
datas = pixGetData(pixs2);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs2);
wpld = pixGetWpl(pixd);
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
if (d == 8) {
if (type == L_CHOOSE_MIN) {
for (j = 0; j < w; j++) {
vals = GET_DATA_BYTE(lines, j);
vald = GET_DATA_BYTE(lined, j);
val = L_MIN(vals, vald);
SET_DATA_BYTE(lined, j, val);
}
} else { /* type == L_CHOOSE_MAX */
for (j = 0; j < w; j++) {
vals = GET_DATA_BYTE(lines, j);
vald = GET_DATA_BYTE(lined, j);
val = L_MAX(vals, vald);
SET_DATA_BYTE(lined, j, val);
}
}
} else { /* d == 16 */
if (type == L_CHOOSE_MIN) {
for (j = 0; j < w; j++) {
vals = GET_DATA_TWO_BYTES(lines, j);
vald = GET_DATA_TWO_BYTES(lined, j);
val = L_MIN(vals, vald);
SET_DATA_TWO_BYTES(lined, j, val);
}
} else { /* type == L_CHOOSE_MAX */
for (j = 0; j < w; j++) {
vals = GET_DATA_TWO_BYTES(lines, j);
vald = GET_DATA_TWO_BYTES(lined, j);
val = L_MAX(vals, vald);
SET_DATA_TWO_BYTES(lined, j, val);
}
}
}
}
return pixd;
}
/*-----------------------------------------------------------------------*
* Scale for maximum dynamic range in 8 bpp image *
*-----------------------------------------------------------------------*/
/*!
* pixMaxDynamicRange()
*
* Input: pixs (4, 8, 16 or 32 bpp source)
* type (L_LINEAR_SCALE or L_LOG_SCALE)
* Return: pixd (8 bpp), or null on error
*
* Notes:
* (1) Scales pixel values to fit maximally within the dest 8 bpp pixd
* (2) Uses a LUT for log scaling
*/
PIX *
pixMaxDynamicRange(PIX *pixs,
l_int32 type)
{
l_uint8 dval;
l_int32 i, j, w, h, d, wpls, wpld, max, sval;
l_uint32 *datas, *datad;
l_uint32 word;
l_uint32 *lines, *lined;
l_float32 factor;
l_float32 *tab;
PIX *pixd;
PROCNAME("pixMaxDynamicRange");
if (!pixs)
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
d = pixGetDepth(pixs);
if (d != 4 && d != 8 && d != 16 && d != 32)
return (PIX *)ERROR_PTR("pixs not in {4,8,16,32} bpp", procName, NULL);
if (type != L_LINEAR_SCALE && type != L_LOG_SCALE)
return (PIX *)ERROR_PTR("invalid type", procName, NULL);
pixGetDimensions(pixs, &w, &h, NULL);
if ((pixd = pixCreate(w, h, 8)) == NULL)
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
pixCopyResolution(pixd, pixs);
datas = pixGetData(pixs);
datad = pixGetData(pixd);
wpls = pixGetWpl(pixs);
wpld = pixGetWpl(pixd);
/* Get max */
max = 0;
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
for (j = 0; j < wpls; j++) {
word = *(lines + j);
if (d == 4) {
max = L_MAX(max, word >> 28);
max = L_MAX(max, (word >> 24) & 0xf);
max = L_MAX(max, (word >> 20) & 0xf);
max = L_MAX(max, (word >> 16) & 0xf);
max = L_MAX(max, (word >> 12) & 0xf);
max = L_MAX(max, (word >> 8) & 0xf);
max = L_MAX(max, (word >> 4) & 0xf);
max = L_MAX(max, word & 0xf);
} else if (d == 8) {
max = L_MAX(max, word >> 24);
max = L_MAX(max, (word >> 16) & 0xff);
max = L_MAX(max, (word >> 8) & 0xff);
max = L_MAX(max, word & 0xff);
} else if (d == 16) {
max = L_MAX(max, word >> 16);
max = L_MAX(max, word & 0xffff);
} else { /* d == 32 */
max = L_MAX(max, word);
}
}
}
/* Map to the full dynamic range of 8 bpp output */
if (d == 4) {
if (type == L_LINEAR_SCALE) {
factor = 255. / (l_float32)max;
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = GET_DATA_QBIT(lines, j);
dval = (l_uint8)(factor * (l_float32)sval + 0.5);
SET_DATA_QBIT(lined, j, dval);
}
}
} else { /* type == L_LOG_SCALE) */
tab = makeLogBase2Tab();
factor = 255. / getLogBase2(max, tab);
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = GET_DATA_QBIT(lines, j);
dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
SET_DATA_BYTE(lined, j, dval);
}
}
FREE(tab);
}
} else if (d == 8) {
if (type == L_LINEAR_SCALE) {
factor = 255. / (l_float32)max;
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = GET_DATA_BYTE(lines, j);
dval = (l_uint8)(factor * (l_float32)sval + 0.5);
SET_DATA_BYTE(lined, j, dval);
}
}
} else { /* type == L_LOG_SCALE) */
tab = makeLogBase2Tab();
factor = 255. / getLogBase2(max, tab);
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = GET_DATA_BYTE(lines, j);
dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
SET_DATA_BYTE(lined, j, dval);
}
}
FREE(tab);
}
} else if (d == 16) {
if (type == L_LINEAR_SCALE) {
factor = 255. / (l_float32)max;
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = GET_DATA_TWO_BYTES(lines, j);
dval = (l_uint8)(factor * (l_float32)sval + 0.5);
SET_DATA_BYTE(lined, j, dval);
}
}
} else { /* type == L_LOG_SCALE) */
tab = makeLogBase2Tab();
factor = 255. / getLogBase2(max, tab);
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = GET_DATA_TWO_BYTES(lines, j);
dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
SET_DATA_BYTE(lined, j, dval);
}
}
FREE(tab);
}
} else { /* d == 32 */
if (type == L_LINEAR_SCALE) {
factor = 255. / (l_float32)max;
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = lines[j];
dval = (l_uint8)(factor * (l_float32)sval + 0.5);
SET_DATA_BYTE(lined, j, dval);
}
}
} else { /* type == L_LOG_SCALE) */
tab = makeLogBase2Tab();
factor = 255. / getLogBase2(max, tab);
for (i = 0; i < h; i++) {
lines = datas + i * wpls;
lined = datad + i * wpld;
for (j = 0; j < w; j++) {
sval = lines[j];
dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
SET_DATA_BYTE(lined, j, dval);
}
}
FREE(tab);
}
}
return pixd;
}
/*-----------------------------------------------------------------------*
* Log base2 lookup *
*-----------------------------------------------------------------------*/
/*
* makeLogBase2Tab()
*
* Input: void
* Return: table (giving the log[base 2] of val)
*/
l_float32 *
makeLogBase2Tab(void)
{
l_int32 i;
l_float32 log2;
l_float32 *tab;
PROCNAME("makeLogBase2Tab");
if ((tab = (l_float32 *)CALLOC(256, sizeof(l_float32))) == NULL)
return (l_float32 *)ERROR_PTR("tab not made", procName, NULL);
log2 = (l_float32)log((l_float32)2);
for (i = 0; i < 256; i++)
tab[i] = (l_float32)log((l_float32)i) / log2;
return tab;
}
/*
* getLogBase2()
*
* Input: val
* logtab (256-entry table of logs)
* Return: logdist, or 0 on error
*/
l_float32
getLogBase2(l_int32 val,
l_float32 *logtab)
{
PROCNAME("getLogBase2");
if (!logtab)
return ERROR_INT("logtab not defined", procName, 0);
if (val < 0x100)
return logtab[val];
else if (val < 0x10000)
return 8.0 + logtab[val >> 8];
else if (val < 0x1000000)
return 16.0 + logtab[val >> 16];
else
return 24.0 + logtab[val >> 24];
}