Logo Search packages:      
Sourcecode: leptonlib version File versions  Download package

pixarith.c

/*====================================================================*
 -  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()
 *
 *      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()
 *           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.
 */


#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);
    d = pixGetDepth(pixs);
    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);
    w = pixGetWidth(pixs);
    h = pixGetHeight(pixs);
    addConstantGrayLow(data, w, h, d, wpl, val);

    return 0;
}

          
/*-------------------------------------------------------------*
 *             Two-image grayscale arithmetic ops              *
 *-------------------------------------------------------------*/
/*!
 *  pixAddGray()
 *
 *      Input:  pixd (<optional>)
 *              pixs1 (can be == to pixd)
 *              pixs2 (can be == to pixd)
 *      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) 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 && pixd != pixs2: (src1 + src2) --> input pixd
 */
PIX *
pixAddGray(PIX  *pixd,
           PIX  *pixs1,
           PIX  *pixs2)
{
l_int32    d, 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 (pixs1 == pixs2)
      return (PIX *)ERROR_PTR("pixs1 and pixs2 must differ", procName, pixd);
    if (!pixSizesEqual(pixs1, pixs2))
      return (PIX *)ERROR_PTR("pixs1 and pixs2 not equal size",
                               procName, pixd);
    if (pixd) {
      if (!pixSizesEqual(pixs1, pixd))
          return (PIX *)ERROR_PTR("pixs1 and pixd not equal size",
                                   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 (pixs1 != pixd) {
      if ((pixd = pixCopy(pixd, pixs1)) == NULL)
          return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    }

    w = pixGetWidth(pixs2);
    h = pixGetHeight(pixs2);
    datas = pixGetData(pixs2);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs2);
    wpld = pixGetWpl(pixd);

    addGrayLow(datad, w, h, d, wpld, datas, wpls);

    return pixd;
}


/*!
 *  pixSubtractGray()
 *
 *      Input:  pixd (<optional>)
 *              pixs1 (can be == to pixd)
 *              pixs2 (can be == to pixd)
 *      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) There are 4 cases.  The result can go to a new dest,
 *          in-place to either pixs1 or pixs2, or to an
 *          existing input dest:
 *          (a) pixd == null   (src1 - src2) --> new pixd
 *          (b) pixd == pixs1  (src1 - src2) --> src1  (in-place)
 *          (c) pixd == pixs2  (src1 - src2) --> src2  (in-place)
 *          (d) pixd != pixs1 && pixd != pixs2)   (src1 - src2) --> input pixd
 *      (5) Case (c) requires swap = 1. For the other three, we can
 *          arrange for the data in pixs1 to be copied to pixd before
 *          calling the low-level code.
 */
PIX *
pixSubtractGray(PIX  *pixd,
                PIX  *pixs1,
                PIX  *pixs2)
{
l_int32    w, h, d, wpls, wpld, swap;
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 (pixs1 == pixs2)
      return (PIX *)ERROR_PTR("pixs1 and pixs2 must differ", procName, pixd);
    if (!pixSizesEqual(pixs1, pixs2))
      return (PIX *)ERROR_PTR("pixs1 and pixs2 not equal size",
                               procName, pixd);
    if (pixd) {
      if (!pixSizesEqual(pixs1, pixd))
          return (PIX *)ERROR_PTR("pixs1 and pixd not equal size",
                                   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 (!pixd || (pixd != pixs1 && pixd != pixs2)) {  /* not in-place */
      if ((pixd = pixCopy(pixd, pixs1)) == NULL)
          return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    }

    w = pixGetWidth(pixs2);
    h = pixGetHeight(pixs2);
    if (pixd != pixs2) {   /* straightforward:  pixd - pixs2 ==> pixd  */
      datas = pixGetData(pixs2);
      datad = pixGetData(pixd);
      wpls = pixGetWpl(pixs2);
      wpld = pixGetWpl(pixd);
      swap = 0;
    }
    else {  /* swap around:  pixs1 - pixd ==> pixd */
      datas = pixGetData(pixs1);
      datad = pixGetData(pixd);
      wpls = pixGetWpl(pixs1);
      wpld = pixGetWpl(pixd);
      swap = 1;
    }
    subtractGrayLow(datad, w, h, d, wpld, datas, wpls, swap);

    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);
    w = pixGetWidth(pixd);
    h = pixGetHeight(pixd);
    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;

    w = pixGetWidth(pixs);
    h = pixGetHeight(pixs);
    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;
}


/*!
 *  pixAccumulate() 
 *
 *      Input:  pixd (32 bpp)
 *              pixs (1, 8, 16 or 32 bpp)
 *              op  (ARITH_ADD or 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
 *          so 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 != ARITH_ADD && op != ARITH_SUBTRACT)
      return ERROR_INT("op must be in {ARITH_ADD, ARITH_SUBTRACT}",
                         procName, 1);

    w = pixGetWidth(pixs);
    h = pixGetHeight(pixs);
    wd = pixGetWidth(pixd);
    hd = pixGetHeight(pixd);
    w = L_MIN(w, wd);
    h = L_MIN(h, hd);

    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    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;

    w = pixGetWidth(pixs);
    h = pixGetHeight(pixs);
    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, 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);

    w = pixGetWidth(pixs1);
    w = L_MIN(w, pixGetWidth(pixs2));
    h = pixGetHeight(pixs1);
    h = L_MIN(h, pixGetHeight(pixs1));
    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, wd, hd, 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);

    if (pixs1 != pixd) {
        if ((pixd = pixCopy(pixd, pixs1)) == NULL)
          return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    }
    d = pixGetDepth(pixd);
    if (pixGetDepth(pixs2) != d)
      return (PIX *)ERROR_PTR("depths unequal", procName, pixd);

    ws = pixGetWidth(pixs2);
    hs = pixGetHeight(pixs2);
    wd = pixGetWidth(pixd);
    hd = pixGetHeight(pixd);
    w = L_MIN(ws, wd);
    h = L_MIN(hs, hd);
    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 or 16 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)
        return (PIX *)ERROR_PTR("pixs must be 4, 8 or 16 bpp", procName, NULL);
    if (type != L_LINEAR_SCALE && type != L_LOG_SCALE)
        return (PIX *)ERROR_PTR("invalid type", procName, NULL);

    w = pixGetWidth(pixs);
    h = pixGetHeight(pixs);
    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 {   /* d == 16 */
            max = L_MAX(max, word >> 16);
            max = L_MAX(max, word & 0xffff);
          }
      }
    }

        /* 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_BYTE(lines, j);
                dval = (l_uint8)(factor * getLogBase2(sval, tab) + 0.5);
                SET_DATA_BYTE(lined, j, dval);
            }
          }
            FREE((void *)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((void *)tab);
      }
    } else {   /* 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((void *)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(2);
    for (i = 0; i < 256; i++)
        tab[i] = (l_float32)log(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];
    else if (val < 0x1000000)
        return 16.0 + logtab[val];
    else
        return 24.0 + logtab[val];
}


Generated by  Doxygen 1.6.0   Back to index