/*====================================================================* - 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]; }