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


/*
 *  scale.c
 *
 *         Top-level scaling
 *               PIX    *pixScale()     ***
 *
 *         24-bit Color (linearly interpolated) scaling
 *               PIX    *pixScaleColorLI()     ***
 *               PIX    *pixScaleColor2xLI()   ***
 *               PIX    *pixScaleColor4xLI()   ***
 *
 *         Grayscale (linearly interpolated) scaling
 *               PIX    *pixScaleGrayLI()
 *               PIX    *pixScaleGray2xLI()
 *               PIX    *pixScaleGray4xLI()
 *
 *         General scaling by closest pixel sampling
 *               PIX    *pixScaleBySampling()
 *
 *         Fast integer factor subsampling RGB to gray and to binary
 *               PIX    *pixScaleRGBToGrayFast()       ***
 *               PIX    *pixScaleRGBToBinaryFast()       ***
 *
 *         Downscaling with (antialias) smoothing
 *               PIX    *pixScaleSmooth()      ***
 *               PIX    *pixScaleRGBToGray2()  ***  special 2x reduction to gray
 *
 *         Downscaling with (antialias) area mapping
 *               PIX    *pixScaleAreaMap()     ***
 *
 *         Binary scaling by closest pixel sampling
 *               PIX    *pixScaleBinary()
 *
 *         Scale-to-gray (1 bpp --> 8 bpp, arbitrary reduction)
 *               PIX    *pixScaleToGray()
 *
 *         Scale-to-gray (1 bpp --> 8 bpp, 2x reduction)
 *               PIX    *pixScaleToGray2()
 *
 *         Scale-to-gray (1 bpp --> 8 bpp, 3x reduction)
 *               PIX    *pixScaleToGray3()
 *
 *         Scale-to-gray (1 bpp --> 8 bpp, 4x reduction)
 *               PIX    *pixScaleToGray4()
 *
 *         Scale-to-gray (1 bpp --> 8 bpp, 8x reduction)
 *               PIX    *pixScaleToGray8()
 *
 *         Scale-to-gray (1 bpp --> 8 bpp, 16x reduction)
 *               PIX    *pixScaleToGray16()
 *
 *         Scale-to-gray by mipmap(1 bpp --> 8 bpp, arbitrary reduction)
 *               PIX    *pixScaleToGrayMipmap()
 *
 *         Grayscale scaling using mipmap
 *               PIX    *pixScaleMipmap()
 *
 *         Upscale 2x followed by binarization
 *               PIX    *pixScaleGray2xLIThresh()
 *               PIX    *pixScaleGray2xLIDither()
 *
 *         Upscale 4x followed by binarization
 *               PIX    *pixScaleGray4xLIThresh()
 *               PIX    *pixScaleGray4xLIDither()
 *
 *
 *  *** Note: these functions make an implicit assumption about RGB
 *            component ordering.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "allheaders.h"



/*------------------------------------------------------------------*
 *                    Top level scaling dispatcher                  *
 *------------------------------------------------------------------*/
/*!
 *  pixScale()
 *
 *      Input:  pixs (1, 2, 4, 8, 16 and 32 bpp)
 *              scalex, scaley
 *      Return: pixd, or null on error
 * 
 *  This function scales 24 bpp RGB; 2, 4 or 8 bpp palette color;
 *  2, 4, 8 or 16 bpp gray; and binary images.
 *
 *  When the input has palette color, the colormap is removed and
 *  a 24 bpp full color pix is made, which is then scaled with
 *  either pixScaleAreaMap() or pixScaleColorLI(), depending on
 *  the scale factor.
 *
 *  Images with 2, 4 or 16 bpp are converted to 8 bpp.
 *
 *  Grayscale and color images are scaled using either antialiased
 *  subsampling (lowpass filtering followed by subsampling), or
 *  by linear interpolation.  For scale factors less than 0.7,
 *  antialiased area mapping is used.
 *
 *  Binary image are scaled by sampling the closest pixel, without
 *  any low-pass filtering (averaging of neighboring pixels).
 *  This will introduce aliasing for reductions, which can be
 *  prevented by using pixScaleToGray() instead.
 *
 *  See notes below regarding use of filtering when subsampling.
 *  
 *  *** Warning: implicit assumption about RGB component order
 *               for LI color scaling
 */
PIX *
pixScale(PIX       *pixs,
	 l_float32  scalex,
	 l_float32  scaley)
{
l_int32    d;
l_float32  maxscale;
PIX       *pixt, *pixd;

    PROCNAME("pixScale");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (scalex == 1.0 && scaley == 1.0)
        return pixCopy(NULL, pixs);

    maxscale = L_MAX(scalex, scaley);
    d = pixGetDepth(pixs);
    if (d == 1)
	return pixScaleBinary(pixs, scalex, scaley);
    else if (d == 16) {
        L_WARNING("pix has 16 bpp; converting to 8 bpp with MSB", procName);
        pixt = pixConvert16To8(pixs, 1);
        if (maxscale >= 0.7) {
            pixd = pixScaleGrayLI(pixt, scalex, scaley);
        } 
        else {  /* maxscale < 0.7 */
            pixd = pixScaleAreaMap(pixt, scalex, scaley);
        } 
	pixDestroy(&pixt);
	return pixd;
    }
    else if (pixGetColormap(pixs)) {   /* 2, 4 or 8 bpp */
        L_WARNING("pixs has colormap; removing", procName);
	if ((pixt = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC)) == NULL)
            return (PIX *)ERROR_PTR("colormap not removed", procName, NULL);
        d = pixGetDepth(pixt);
        if (maxscale >= 0.7) {
            if (d == 8)
                pixd = pixScaleGrayLI(pixt, scalex, scaley);
            else  /* d == 32 */
                pixd = pixScaleColorLI(pixt, scalex, scaley);
	}
        else {  /* maxscale < 0.7 */
            pixd = pixScaleAreaMap(pixt, scalex, scaley);
        } 
	pixDestroy(&pixt);
	return pixd;
    }
    else if (d == 2 || d == 4) {
        L_WARNING("pix has 2 or 4 bpp without colormap; converting to 8 bpp",
                procName);
        pixt = pixConvertTo8(pixs);
        if (maxscale >= 0.7) {
            pixd = pixScaleGrayLI(pixt, scalex, scaley);
        } 
        else {  /* maxscale < 0.7 */
            pixd = pixScaleAreaMap(pixt, scalex, scaley);
        } 
	pixDestroy(&pixt);
	return pixd;
    }
    else if (d == 8) {
        if (maxscale >= 0.7) {
            return pixScaleGrayLI(pixs, scalex, scaley);
        } 
        else {  /* maxscale < 0.7 */
            return pixScaleAreaMap(pixs, scalex, scaley);
        } 
    }
    else if (d == 32) {
        if (maxscale >= 0.7) {
	    return pixScaleColorLI(pixs, scalex, scaley);
        } 
        else {  /* maxscale < 0.7 */
            return pixScaleAreaMap(pixs, scalex, scaley);
        } 
    }
    else
        return (PIX *)ERROR_PTR("pixs not {1,2,4,8,16,32} bpp", procName, NULL);
}


/*------------------------------------------------------------------*
 *            24-bit color scaling by linear interpolation          *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleColorLI()
 *
 *      Input:  pixs  (32 bpp, representing 24 bpp 3-color)
 *              scalex, scaley
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) If this is used for scale factors less than about 0.7,
 *          it will suffer from antialiasing.  Particularly for
 *          document images with sharp edges, use pixScaleSmooth()
 *          or pixScaleAreaMap() instead.
 *      (2) For the general case, it's 3-4x faster to manipulate
 *          the color pixels directly, rather than to make component
 *          images, do each 8 bpp component separately, and combine
 *          the result.  The speed on intel hardware for the
 *          general case (not 2x or 4x LI) is about
 *          4.8 * 10^6 dest-pixels/sec/GHz.
 *      (3) The slow method does about 2 * 10^6 dest-pixels/sec/GHz.
 *          It is implemented using:
 *              pixr = pixGetRGBComponent(pixs, COLOR_RED);
 *              pixrs = pixScaleGrayLI(pixr, scalex, scaley);
 *              pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
 *              pixgs = pixScaleGrayLI(pixg, scalex, scaley);
 *              pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
 *              pixbs = pixScaleGrayLI(pixb, scalex, scaley);
 *              pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleColorLI(PIX      *pixs,
	       l_float32  scalex,
	       l_float32  scaley)
{
l_int32    ws, hs, wpls, wd, hd, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleColorLI");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs must be 32 bpp", procName, NULL);

        /* do fast special cases if possible */
    if (scalex == 1.0 && scaley == 1.0)
	return pixCopy(NULL, pixs);
    if (scalex == 2.0 && scaley == 2.0)
	return pixScaleColor2xLI(pixs);
    if (scalex == 4.0 && scaley == 4.0)
	return pixScaleColor4xLI(pixs);

        /* general case */
    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    wd = (l_int32)(scalex * (l_float32)ws + 0.5);
    hd = (l_int32)(scaley * (l_float32)hs + 0.5);
    if ((pixd = pixCreate(wd, hd, 32)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, scalex, scaley);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleColorLILow(datad, wd, hd, wpld, datas, ws, hs, wpls);
    return pixd;
}


/*!
 *  pixScaleColor2xLI()
 * 
 *      Input:  pixs  (32 bpp, representing 24 bpp 3-color)
 *      Return: pixd, or null on error
 *
 *  This is a special case of linear interpolated scaling
 *  for 2x upscaling.  It is about 4x faster than using
 *  pixScaleGray2xLI() on each component separately.
 *  The speed on intel hardware for this special case
 *  is about 64 * 10^6 dest-pixels/sec/GHz (!!)
 *
 *  The slow method, operating on each component separately, has 
 *  a speed of about 16 * 10^6 dest-pixels/sec/GHz:
 *       pixr = pixGetRGBComponent(pixs, COLOR_RED);
 *       pixrs = pixScaleGray2xLI(pixr);
 *       pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
 *       pixgs = pixScaleGray2xLI(pixg);
 *       pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
 *       pixbs = pixScaleGray2xLI(pixb);
 *       pixd = pixCreateRGBImage(pixrs, pixgs, pixbs);
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleColor2xLI(PIX  *pixs)
{
l_int32    d, ws, hs, wpls, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleColor2xLI");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if ((d = pixGetDepth(pixs)) != 32)
        return (PIX *)ERROR_PTR("pixs must be 32 bpp", procName, NULL);
    
    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    if ((pixd = pixCreate(2 * ws, 2 * hs, 32)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 2.0, 2.0);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleColor2xLILow(datad, wpld, datas, ws, hs, wpls);

    return pixd;
}


/*!
 *  pixScaleColor4xLI()
 *
 *      Input:  pixs  (32 bpp, representing 24 bpp 3-color)
 *      Return: pixd, or null on error
 *
 *  Special version using separate component 4x upscaling
 *  It would be faster to inline results analogously to
 *  scaleColor4xLILow(), but is this used enough to make
 *  it worth the effort?
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleColor4xLI(PIX  *pixs)
{
PIX  *pixr, *pixg, *pixb;
PIX  *pixrs, *pixgs, *pixbs;
PIX  *pixd;

    PROCNAME("pixScaleColor4xLI");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs must be 32 bpp", procName, NULL);

    pixr = pixGetRGBComponent(pixs, COLOR_RED);
    pixrs = pixScaleGray4xLI(pixr);
    pixDestroy(&pixr);
    pixg = pixGetRGBComponent(pixs, COLOR_GREEN);
    pixgs = pixScaleGray4xLI(pixg);
    pixDestroy(&pixg);
    pixb = pixGetRGBComponent(pixs, COLOR_BLUE);
    pixbs = pixScaleGray4xLI(pixb);
    pixDestroy(&pixb);

    if ((pixd = pixCreateRGBImage(pixrs, pixgs, pixbs)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);

    pixDestroy(&pixrs);
    pixDestroy(&pixgs);
    pixDestroy(&pixbs);

    return pixd;
}


/*------------------------------------------------------------------*
 *              General linear interpolated scaling                 *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleGrayLI()
 *
 *      Input:  pixs (8 bpp grayscale)
 *              scalex 
 *              scaley 
 *      Return: pixd, or null on error
 *
 *  This function is appropriate for upscaling
 *  (magnification: scale factors > 1), and for a
 *  small amount of downscaling (reduction: scale
 *  factors > 0.5).   For scale factors less than 0.5,
 *  the best result is obtained by area mapping,
 *  but this is very expensive.  So for such large
 *  reductions, it is more appropriate to do low pass
 *  filtering followed by subsampling, a combination
 *  which is effectively a cheap form of area mapping.
 *
 *  Some details follow.
 *
 *  For each pixel in the dest, this does a linear
 *  interpolation of 4 neighboring pixels in the src.
 *  Specifically, consider the UL corner of src and
 *  dest pixels.  The UL corner of the dest falls within
 *  a src pixel, whose four corners are the UL corners
 *  of 4 adjacent src pixels.  The value of the dest
 *  is taken by linear interpolation using the values of
 *  the four src pixels and the distance of the UL corner
 *  of the dest from each corner.
 *
 *  If the image is expanded so that the dest pixel is
 *  smaller than the src pixel, such interpolation
 *  is a reasonable approach.  This interpolation is
 *  also good for a small image reduction factor that
 *  is not more than a 2x reduction.
 *
 *  Note that the linear interpolation algorithm for scaling
 *  is identical in form to the area-mapping algorithm
 *  for grayscale rotation.  The latter corresponds to a
 *  translation of each pixel without scaling.
 * 
 *  This function is NOT optimal if the scaling involves
 *  a large reduction.    If the image is significantly
 *  reduced, so that the dest pixel is much larger than
 *  the src pixels, this interpolation, which is over src
 *  pixels only near the UL corner of the dest pixel,
 *  is not going to give a good area-mapping average.
 *  Because area mapping for image scaling is considerably
 *  more computationally intensive than linear interpolation,
 *  we choose not to use it.   For large image reduction,
 *  linear interpolation over adjacent src pixels
 *  degenerates asymptotically to subsampling.  But
 *  subsampling without a low-pass pre-filter causes
 *  aliasing by the nyquist theorem.  To avoid aliasing,
 *  a low-pass filter (e.g., an averaging filter) of
 *  size roughly equal to the dest pixel (i.e., the
 *  reduction factor) should be applied to the src before
 *  subsampling.
 *
 *  As an alternative to low-pass filtering and subsampling
 *  for large reduction factors, linear interpolation can
 *  also be done between the (widely separated) src pixels in
 *  which the corners of the dest pixel lie.  This also is
 *  not optimal, as it samples src pixels only near the
 *  corners of the dest pixel, and it is not implemented.
 */
PIX *
pixScaleGrayLI(PIX       *pixs,
	       l_float32  scalex,
	       l_float32  scaley)
{
l_int32    d, ws, hs, wpls, wd, hd, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleGrayLI");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if ((d = pixGetDepth(pixs)) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (pixGetColormap(pixs))
        L_WARNING("pix has colormap", procName);

        /* do fast special cases if possible */
    if (scalex == 1.0 && scaley == 1.0)
	return pixCopy(NULL, pixs);
    if (scalex == 2.0 && scaley == 2.0)
	return pixScaleGray2xLI(pixs);
    if (scalex == 4.0 && scaley == 4.0)
	return pixScaleGray4xLI(pixs);
    
        /* general case */
    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    wd = (l_int32)(scalex * (l_float32)ws + 0.5);
    hd = (l_int32)(scaley * (l_float32)hs + 0.5);
    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, scalex, scaley);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleGrayLILow(datad, wd, hd, wpld, datas, ws, hs, wpls);

    return pixd;
}


/*!
 *  pixScaleGray2xLI()
 * 
 *      Input:  pixs (8 bpp grayscale)
 *      Return: pixd, or null on error
 *
 *  This is a special case of linear interpolated scaling for 2x upscaling.
 */
PIX *
pixScaleGray2xLI(PIX  *pixs)
{
l_int32    d, ws, hs, wpls, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleGray2xLI");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if ((d = pixGetDepth(pixs)) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (pixGetColormap(pixs))
        L_WARNING("pix has colormap", procName);
    
    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    if ((pixd = pixCreate(2 * ws, 2 * hs, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 2.0, 2.0);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleGray2xLILow(datad, wpld, datas, ws, hs, wpls);

    return pixd;
}


/*!
 *  pixScaleGray4xLI()
 * 
 *      Input:  pixs (8 bpp grayscale)
 *      Return: pixd, or null on error
 *
 *  This is a special case of linear interpolated scaling for 4x upscaling.
 */
PIX *
pixScaleGray4xLI(PIX  *pixs)
{
l_int32    d, ws, hs, wpls, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleGray4xLI");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if ((d = pixGetDepth(pixs)) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (pixGetColormap(pixs))
        L_WARNING("pix has colormap", procName);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    if ((pixd = pixCreate(4 * ws, 4 * hs, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 4.0, 4.0);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleGray4xLILow(datad, wpld, datas, ws, hs, wpls);

    return pixd;
}



/*------------------------------------------------------------------*
 *              General scaling by closest pixel sampling           *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleBySampling()
 *
 *      Input:  pixs (1, 2, 4, 8, 16, 32 bpp)
 *              scalex, scaley
 *      Return: pixd, or null on error
 * 
 *  Note: this function samples from the source without
 *        filtering, so that aliasing will result for
 *        subsampling (scalex and scaley < 1.0).
 */
PIX *
pixScaleBySampling(PIX       *pixs,
	           l_float32  scalex,
	           l_float32  scaley)
{
l_int32    ws, hs, d, wpls, wd, hd, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleBySampling");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (scalex == 1.0 && scaley == 1.0)
        return pixCopy(NULL, pixs);
    if ((d = pixGetDepth(pixs)) == 1)
        return pixScaleBinary(pixs, scalex, scaley);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    wd = (l_int32)(scalex * (l_float32)ws + 0.5);
    hd = (l_int32)(scaley * (l_float32)hs + 0.5);
    if ((pixd = pixCreate(wd, hd, d)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, scalex, scaley);
    pixCopyColormap(pixd, pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleBySamplingLow(datad, wd, hd, wpld, datas, ws, hs, d, wpls);

    return pixd;
}


/*------------------------------------------------------------------*
 *            Fast integer factor subsampling RGB to gray           *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleRGBToGrayFast()
 *
 *      Input:  pixs (32 bpp RGB)
 *              factor (integer reduction factor)
 *              color (one of COLOR_RED, COLOR_GREEN, COLOR_BLUE)
 *      Return: pixd (8 bpp), or null on error
 * 
 *  Notes:
 *    (1) This does simultaneous subsampling by an integer factor and
 *        extraction of the color from the RGB pix.
 *    (2) It is designed for maximum speed, and is used for quickly
 *        generating a downsized grayscale image from a high resolution
 *        RGB image.  This would typically be used for image analysis.
 *    (3) The standard color byte order (RGBA) is assumed.
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleRGBToGrayFast(PIX     *pixs,
                      l_int32  factor,
                      l_int32  color)
{
l_int32    byteval, shift;
l_int32    i, j, ws, hs, wd, hd, wpls, wpld;
l_uint32  *datas, *words, *datad, *lined;
l_float32  scale;
PIX       *pixd;

    PROCNAME("pixScaleRGBToGrayFast");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("depth not 32 bpp", procName, NULL);
    if (factor < 1)
        return (PIX *)ERROR_PTR("factor must be >= 1", procName, NULL);

    if (color == COLOR_RED)
        shift = 24;
    else if (color == COLOR_GREEN)
        shift = 16;
    else if (color == COLOR_BLUE)
        shift = 8;
    else
        return (PIX *)ERROR_PTR("invalid color", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);

    wd = ws / factor;
    hd = hs / factor;
    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
	return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    scale = 1. / (l_float32) factor;
    pixScaleResolution(pixd, scale, scale);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    
    for (i = 0; i < hd; i++) {
	words = datas + i * factor * wpls;
	lined = datad + i * wpld;
        for (j = 0; j < wd; j++, words += factor) {
            byteval = ((*words) >> shift) & 0xff;
            SET_DATA_BYTE(lined, j, byteval);
        }
    }

    return pixd;
}


/*!
 *  pixScaleRGBToBinaryFast()
 *
 *      Input:  pixs (32 bpp RGB)
 *              factor (integer reduction factor)
 *              thresh (binarization threshold)
 *      Return: pixd (1 bpp), or null on error
 * 
 *  Notes:
 *    (1) This does simultaneous subsampling by an integer factor and
 *        conversion from RGB to gray to binary.
 *    (2) It is designed for maximum speed, and is used for quickly
 *        generating a downsized binary image from a high resolution
 *        RGB image.  This would typically be used for image analysis.
 *    (3) It uses the green channel to represent the RGB pixel intensity.
 *    (4) The standard color byte order (RGBA) is assumed.
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleRGBToBinaryFast(PIX     *pixs,
                        l_int32  factor,
                        l_int32  thresh)
{
l_int32    byteval;
l_int32    i, j, ws, hs, wd, hd, wpls, wpld;
l_uint32  *datas, *words, *datad, *lined;
l_float32  scale;
PIX       *pixd;

    PROCNAME("pixScaleRGBToBinaryFast");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (factor < 1)
        return (PIX *)ERROR_PTR("factor must be >= 1", procName, NULL);
    if (pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("depth not 32 bpp", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);

    wd = ws / factor;
    hd = hs / factor;
    if ((pixd = pixCreate(wd, hd, 1)) == NULL)
	return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    scale = 1. / (l_float32) factor;
    pixScaleResolution(pixd, scale, scale);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    for (i = 0; i < hd; i++) {
	words = datas + i * factor * wpls;
	lined = datad + i * wpld;
	for (j = 0; j < wd; j++, words += factor) {
	    byteval = ((*words) >> 16) & 0xff;
            if (byteval < thresh)
                SET_DATA_BIT(lined, j);
	}
    }

    return pixd;
}


/*------------------------------------------------------------------*
 *               Downscaling with (antialias) smoothing             *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleSmooth()
 *
 *      Input:  pixs (2, 4, 8 or 32 bpp; and 2, 4, 8 bpp with colormap)
 *              scalex, scaley (must both be <= 0.7)
 *      Return: pixd, or null on error
 * 
 *  Notes:
 *    (1) This function should only be used when the scale factors are less
 *        than or equal to 0.7 (i.e., more than about 1.42x reduction).
 *        If either scale factor is larger than 0.7, we issue a warning
 *        and invoke pixScale().
 *    (2) This works only on 2, 4, 8 and 32 bpp images, and if there is
 *        a colormap, it is removed by converting to RGB.  In other
 *        cases, we issue a warning and invoke pixScale().
 *    (3) It does simple (flat filter) convolution, with a filter size
 *        commensurate with the amount of reduction, to avoid antialiasing.
 *    (4) It does simple subsampling after smoothing, which is appropriate
 *        for this range of scaling.  Linear interpolation gives essentially
 *        the same result with more computation for these scale factors,
 *        so we don't use it.
 *    (5) The result is the same as doing a full block convolution followed by
 *        subsampling, but this is faster because the results of the block
 *        convolution are only computed at the subsampling locations.
 *        In fact, the computation time is approximately independent of
 *        the scale factor, because the convolution kernel is adjusted
 *        so that each source pixel is summed approximately once.
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleSmooth(PIX       *pix,
               l_float32  scalex,
               l_float32  scaley)
{
l_int32    ws, hs, d, wd, hd, wpls, wpld, isize;
l_uint32  *datas, *datad;
l_float32  minscale, size;
PIX    	  *pixs, *pixd;

    PROCNAME("pixScaleSmooth");

    if (!pix)
        return (PIX *)ERROR_PTR("pix not defined", procName, NULL);
    if (scalex > 0.7 || scaley > 0.7) {
        L_WARNING("scaling factor not <= 0.7; doing regular scaling", procName);
        return pixScale(pix, scalex, scaley);
    }

        /* remove colormap if necessary;
	 * if 2 bpp or 4 bpp gray, convert to 8 bpp */
    d = pixGetDepth(pix);
    if ((d == 2 || d == 4 || d == 8) && pixGetColormap(pix)) {
        L_WARNING("pix has colormap; removing", procName);
	pixs = pixRemoveColormap(pix, REMOVE_CMAP_BASED_ON_SRC);
	d = pixGetDepth(pixs);
    }
    else if (d == 2 || d == 4) {
        pixs = pixConvertTo8(pix);
	d = 8;
    }
    else
        pixs = pixClone(pix);

    if (d != 8 && d != 32) {   /* d == 1 or d == 16 */
        L_WARNING("depth not 8 or 32 bpp; doing regular scaling", procName);
	pixDestroy(&pixs);
        return pixScale(pix, scalex, scaley);
    }

        /* If 1.42 < 1/minscale < 2.5, use isize = 2
         * If 2.5 =< 1/minscale < 3.5, use isize = 3, etc.
         * Under no conditions use isize < 2  */
    minscale = L_MIN(scalex, scaley);
    size = 1.0 / minscale;   /* ideal filter full width */
    isize = L_MAX(2, (l_int32)(size + 0.5));

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    if ((ws < isize) || (hs < isize)) {
	pixDestroy(&pixs);
        return (PIX *)ERROR_PTR("pixs too small", procName, NULL);
    }
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    wd = (l_int32)(scalex * (l_float32)ws + 0.5);
    hd = (l_int32)(scaley * (l_float32)hs + 0.5);
    if (wd < 1 || hd < 1) {
	pixDestroy(&pixs);
        return (PIX *)ERROR_PTR("pixd too small", procName, NULL);
    }
    if ((pixd = pixCreate(wd, hd, d)) == NULL) {
	pixDestroy(&pixs);
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    }
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, scalex, scaley);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleSmoothLow(datad, wd, hd, wpld, datas, ws, hs, d, wpls, isize);

    pixDestroy(&pixs);
    return pixd;
}


/*!
 *  pixScaleRGBToGray2()
 *
 *      Input:  pixs (32 bpp)
 *              rwt, gwt, bwt (must sum to 1.0)
 *      Return: pixd, (8 bpp, 2x reduced), or null on error
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleRGBToGray2(PIX       *pixs,
                   l_float32  rwt,
                   l_float32  gwt,
                   l_float32  bwt)
{
l_int32    ws, hs, wd, hd, wpls, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleRGBToGray2");

    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 (rwt + gwt + bwt < 0.98 || rwt + gwt + bwt > 1.02)
        return (PIX *)ERROR_PTR("sum of wts should be 1.0", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wpls = pixGetWpl(pixs);
    datas = pixGetData(pixs);
    wd = ws / 2;
    hd = hs / 2;
    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 0.5, 0.5);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);

    scaleRGBToGray2Low(datad, wd, hd, wpld, datas, ws, hs, wpls, rwt, gwt, bwt);
    return pixd;
}


/*------------------------------------------------------------------*
 *             Downscaling with (antialias) area mapping            *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleAreaMap()
 *
 *      Input:  pixs (2, 4, 8 or 32 bpp; and 2, 4, 8 bpp with colormap)
 *              scalex, scaley (must both be <= 0.7)
 *      Return: pixd, or null on error
 * 
 *  Notes:
 *    (1) This function should only be used when the scale factors are less
 *        than or equal to 0.7 (i.e., more than about 1.42x reduction).
 *        If either scale factor is larger than 0.7, we issue a warning
 *        and invoke pixScale().
 *    (2) This works only on 2, 4, 8 and 32 bpp images.  If there is
 *        a colormap, it is removed by converting to RGB.  In other
 *        cases, we issue a warning and invoke pixScale().
 *    (3) It does a relatively expensive area mapping computation, to
 *        avoid antialiasing.  It is about 2x slower than pixScaleSmooth(),
 *        but the results are much better on fine text.
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixScaleAreaMap(PIX       *pix,
                l_float32  scalex,
                l_float32  scaley)
{
l_int32    ws, hs, d, wd, hd, wpls, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixs, *pixd;

    PROCNAME("pixScaleAreaMap");

    if (!pix)
        return (PIX *)ERROR_PTR("pix not defined", procName, NULL);
    d = pixGetDepth(pix);
    if (d == 1 || d == 16)
        return (PIX *)ERROR_PTR("pix not 1 or 16 bpp", procName, NULL);
    if (scalex > 0.7 || scaley > 0.7) {
        L_WARNING("scaling factor not <= 0.7; doing regular scaling", procName);
        return pixScale(pix, scalex, scaley);
    }

        /* remove colormap if necessary;
	 * if 2 bpp or 4 bpp gray, convert to 8 bpp */
    if ((d == 2 || d == 4 || d == 8) && pixGetColormap(pix)) {
        L_WARNING("pix has colormap; removing", procName);
	pixs = pixRemoveColormap(pix, REMOVE_CMAP_BASED_ON_SRC);
	d = pixGetDepth(pixs);
    }
    else if (d == 2 || d == 4) {
        pixs = pixConvertTo8(pix);
	d = 8;
    }
    else
        pixs = pixClone(pix);

    if (d != 8 && d != 32) {   /* d == 1 or d == 16 */
        L_WARNING("depth not 8 or 32 bpp; doing regular scaling", procName);
	pixDestroy(&pixs);
        return pixScale(pix, scalex, scaley);
    }

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    wd = (l_int32)(scalex * (l_float32)ws + 0.5);
    hd = (l_int32)(scaley * (l_float32)hs + 0.5);
    if (wd < 1 || hd < 1) {
	pixDestroy(&pixs);
        return (PIX *)ERROR_PTR("pixd too small", procName, NULL);
    }
    if ((pixd = pixCreate(wd, hd, d)) == NULL) {
	pixDestroy(&pixs);
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    }
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, scalex, scaley);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    if (d == 8)
        scaleGrayAreaMapLow(datad, wd, hd, wpld, datas, ws, hs, wpls);
    else  /* RGB, d == 32 */
        scaleColorAreaMapLow(datad, wd, hd, wpld, datas, ws, hs, wpls);

    pixDestroy(&pixs);
    return pixd;
}


/*------------------------------------------------------------------*
 *               Binary scaling by closest pixel sampling           *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleBinary()
 *
 *      Input:  pixs (1 bpp)
 *              scalex, scaley
 *      Return: pixd, or null on error
 * 
 *  Note: this function samples from the source without
 *        filtering, so that aliasing will result for
 *        subsampling (scalex and scaley < 1.0).
 */
PIX *
pixScaleBinary(PIX       *pixs,
	       l_float32  scalex,
	       l_float32  scaley)
{
l_int32    ws, hs, wpls, wd, hd, wpld;
l_uint32  *datas, *datad;
PIX    	  *pixd;

    PROCNAME("pixScaleBinary");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs must be 1 bpp", procName, NULL);
    if (scalex == 1.0 && scaley == 1.0)
        return pixCopy(NULL, pixs);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    wd = (l_int32)(scalex * (l_float32)ws + 0.5);
    hd = (l_int32)(scaley * (l_float32)hs + 0.5);
    if ((pixd = pixCreate(wd, hd, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyColormap(pixd, pixs);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, scalex, scaley);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleBinaryLow(datad, wd, hd, wpld, datas, ws, hs, wpls);

    return pixd;
}



/*------------------------------------------------------------------*
 *       Scale-to-gray (1 bpp --> 8 bpp, arbitrary reduction)       *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleToGray()
 *
 *      Input:  1 bpp pix
 *              scalefactor (reduction, < 1.0)
 *      Return: 8 bpp pix, scaled down by scalefactor in each direction,
 *              or NULL on error.
 *
 *  Notes:
 *
 *  Binary images have sharp edges, so they intrinsically have very
 *  high frequency content.  To avoid aliasing, they must be low-pass
 *  filtered, which tends to blur the edges.  How can we keep relatively
 *  crisp edges without aliasing?  The trick is to do binary upscaling
 *  followed by a power-of-2 scaleToGray.  For large reductions, where
 *  you don't end up with much detail, some corners can be cut.
 *
 *  The intent here is to get high quality reduced grayscale
 *  images with relatively little computation.  We do binary
 *  pre-scaling followed by scaleToGrayN() for best results,
 *  esp. to avoid excess blur when the scale factor is near
 *  an inverse power of 2.  Where a low-pass filter is required,
 *  we use simple convolution kernels: either the hat filter for
 *  linear interpolation or a flat filter for larger downscaling.
 *  Other choices, such as a perfect bandpass filter with infinite extent
 *  (the sinc) or various approximations to it (e.g., lanczos), are
 *  unnecessarily expensive.
 *
 *  The choices made are as follows:
 *      (1) Do binary upscaling before scaleToGrayN() for scalefactors > 1/8
 *      (2) Do binary downscaling before scaleToGray8() for scalefactors
 *          between 1/16 and 1/8.
 *      (3) Use scaleToGray16() before grayscale downscaling for
 *          scalefactors less than 1/16
 *  Another reasonable choice would be to start binary downscaling
 *  for scalefactors below 1/4, rather than below 1/8 as we do here.
 *
 *  The general scaling rules, not all of which are used here, go as follows:
 *      (1) For grayscale upscaling, use pixScaleGrayLI().  However,
 *          note that edges will be visibly blurred for scalefactors
 *          near (but above) 1.0.  Replication will avoid edge blur,
 *          and should be considered for factors very near 1.0.
 *      (2) For grayscale downscaling with a scale factor larger than
 *          about 0.7, use pixScaleGrayLI().  For scalefactors near
 *          (but below) 1.0, you tread between Scylla and Charybdis.
 *          pixScaleGrayLI() again gives edge blurring, but
 *          pixScaleBySampling() gives visible aliasing.
 *      (3) For grayscale downscaling with a scale factor smaller than
 *          about 0.7, use pixScaleSmooth()
 *      (4) For binary input images, do as much scale to gray as possible
 *          using the special integer functions (2, 3, 4, 8 and 16).
 *      (5) It is better to upscale in binary, followed by scaleToGrayN()
 *          than to do scaleToGrayN() followed by an upscale using either
 *          LI or oversampling.
 *      (6) It may be better to downscale in binary, followed by
 *          scaleToGrayN() than to first use scaleToGrayN() followed by
 *          downscaling.  For downscaling between 8x and 16x, this is
 *          a reasonable option.
 *      (7) For reductions greater than 16x, it's reasonable to use
 *          scaleToGray16() followed by further grayscale downscaling.
 */
PIX *
pixScaleToGray(PIX       *pixs,
               l_float32  scalefactor)
{
l_int32    w, h, minsrc, mindest;
l_float32  mag, red;
PIX       *pixt, *pixd;

    PROCNAME("pixScaleToGray");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
    if (scalefactor >= 1.0)
        return (PIX *)ERROR_PTR("scalefactor not < 1.0", procName, NULL);
    w = pixGetWidth(pixs);
    h = pixGetHeight(pixs);
    minsrc = L_MIN(w, h);
    mindest = (l_int32)((l_float32)minsrc * scalefactor);
    if (mindest < 2)
        return (PIX *)ERROR_PTR("scalefactor too small", procName, NULL);

    if (scalefactor > 0.5) {   /* see note (5) */
        mag = 2.0 * scalefactor;  /* will be < 2.0 */
/*	fprintf(stderr, "2x with mag %7.3f\n", mag);  */
	if ((pixt = pixScaleBinary(pixs, mag, mag)) == NULL)
	    return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
	pixd = pixScaleToGray2(pixt);
    }
    else if (scalefactor == 0.5)
	return pixd = pixScaleToGray2(pixs);
    else if (scalefactor > 0.334) {   /* see note (5) */
        mag = 3.0 * scalefactor;   /* will be < 1.5 */
/*	fprintf(stderr, "3x with mag %7.3f\n", mag);  */
	if ((pixt = pixScaleBinary(pixs, mag, mag)) == NULL)
	    return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
	pixd = pixScaleToGray3(pixt);
    }
    else if (scalefactor > 0.25) {  /* see note (5) */
	mag = 4.0 * scalefactor;   /* will be < 1.333 */
/*	fprintf(stderr, "4x with mag %7.3f\n", mag);  */
	if ((pixt = pixScaleBinary(pixs, mag, mag)) == NULL)
	    return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
	pixd = pixScaleToGray4(pixt);
    }
    else if (scalefactor == 0.25)
	return pixd = pixScaleToGray4(pixs);
    else if (scalefactor > 0.125) {  /* see note (5) */
	mag = 8.0 * scalefactor;   /*  will be < 2.0  */
/*	fprintf(stderr, "8x with mag %7.3f\n", mag);  */
	if ((pixt = pixScaleBinary(pixs, mag, mag)) == NULL)
	    return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
	pixd = pixScaleToGray8(pixt);
    }
    else if (scalefactor == 0.125)
	return pixd = pixScaleToGray8(pixs);
    else if (scalefactor > 0.0625) {  /* see note (6) */
	red = 8.0 * scalefactor;   /* will be > 0.5 */
/*	fprintf(stderr, "8x with red %7.3f\n", red);  */
	if ((pixt = pixScaleBinary(pixs, red, red)) == NULL)
	    return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
	pixd = pixScaleToGray8(pixt);
    }
    else if (scalefactor == 0.0625)
	return pixd = pixScaleToGray16(pixs);
    else {  /* see note (7) */
	red = 16.0 * scalefactor;  /* will be <= 1.0 */
/*	fprintf(stderr, "16x with red %7.3f\n", red);  */
	if ((pixt = pixScaleToGray16(pixs)) == NULL)
	    return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
        if (red < 0.7)
	    pixd = pixScaleSmooth(pixt, red, red);  /* see note (3) */
        else	
	    pixd = pixScaleGrayLI(pixt, red, red);  /* see note (2) */
    }

    pixDestroy(&pixt);
    if (!pixd)
	return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    else
        return pixd;
}



/*------------------------------------------------------------------*
 *           Scale-to-gray (1 bpp --> 8 bpp, 2x reduction)          *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleToGray2()
 *
 *      Input:  1 bpp pix
 *      Return: 8 bpp pix, scaled down by 2x in each direction,
 *              or NULL on error.
 */
PIX *
pixScaleToGray2(PIX  *pixs)
{
l_uint8   *valtab;
l_int32    ws, hs, wd, hd;
l_int32    wpld, wpls;
l_uint32  *sumtab;
l_uint32  *datas, *datad;
PIX       *pixd;

    PROCNAME("pixScaleToGray2");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs must be 1 bpp", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = (ws / 2) & 0xfffffffc;    /* truncate to factor of 4 */
    hd = hs / 2;
    if (wd == 0 || hd == 0)
        return (PIX *)ERROR_PTR("pixs too small", procName, NULL);

    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 0.5, 0.5);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    if ((sumtab = makeSumTabSG2()) == NULL)
	return (PIX *)ERROR_PTR("sumtab not made", procName, NULL);

    if ((valtab = makeValTabSG2()) == NULL)
	return (PIX *)ERROR_PTR("valtab not made", procName, NULL);

    scaleToGray2Low(datad, wd, hd, wpld, datas, wpls, sumtab, valtab);

    FREE((void *)sumtab);
    FREE((void *)valtab);

    return pixd;
}


/*------------------------------------------------------------------*
 *           Scale-to-gray (1 bpp --> 8 bpp, 3x reduction)          *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleToGray3()
 *
 *      Input:  1 bpp pix
 *      Return: 8 bpp pix, scaled down by 3x in each direction,
 *              or NULL on error.
 *
 *  Time: about 0.1 sec on 300 ppi page image (866 MHz P3)
 *  Speed: about 1 src pix per 10 cycles
 */
PIX *
pixScaleToGray3(PIX  *pixs)
{
l_uint8   *valtab;
l_int32    ws, hs, wd, hd;
l_int32    wpld, wpls;
l_uint32  *sumtab;
l_uint32  *datas, *datad;
PIX       *pixd;

    PROCNAME("pixScaleToGray3");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs not 1 bpp", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = (ws / 3) & 0xfffffff8;    /* truncate to factor of 8 */
    hd = hs / 3;
    if (wd == 0 || hd == 0)
        return (PIX *)ERROR_PTR("pixs too small", procName, NULL);

    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 0.3333, 0.3333);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    if ((sumtab = makeSumTabSG3()) == NULL)
	return (PIX *)ERROR_PTR("sumtab not made", procName, NULL);

    if ((valtab = makeValTabSG3()) == NULL)
	return (PIX *)ERROR_PTR("valtab not made", procName, NULL);

    scaleToGray3Low(datad, wd, hd, wpld, datas, wpls, sumtab, valtab);

    FREE((void *)sumtab);
    FREE((void *)valtab);

    return pixd;
}


/*------------------------------------------------------------------*
 *                         Scale-to-gray 4x                         *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleToGray4()
 *
 *      Input:  1 bpp pix
 *      Return: 8 bpp pix, scaled down by 4x in each direction,
 *              or NULL on error.
 */
PIX *
pixScaleToGray4(PIX  *pixs)
{
l_uint8   *valtab;
l_int32    ws, hs, wd, hd;
l_int32    wpld, wpls;
l_uint32  *sumtab;
l_uint32  *datas, *datad;
PIX       *pixd;

    PROCNAME("pixScaleToGray4");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs must be 1 bpp", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = (ws / 4) & 0xfffffffe;    /* truncate to factor of 2 */
    hd = hs / 4;
    if (wd == 0 || hd == 0)
        return (PIX *)ERROR_PTR("pixs too small", procName, NULL);

    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 0.25, 0.25);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    if ((sumtab = makeSumTabSG4()) == NULL)
	return (PIX *)ERROR_PTR("sumtab not made", procName, NULL);

    if ((valtab = makeValTabSG4()) == NULL)
	return (PIX *)ERROR_PTR("valtab not made", procName, NULL);

    scaleToGray4Low(datad, wd, hd, wpld, datas, wpls, sumtab, valtab);

    FREE((void *)sumtab);
    FREE((void *)valtab);

    return pixd;
}


/*------------------------------------------------------------------*
 *                         Scale-to-gray 8x                         *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleToGray8()
 *
 *      Input:  1 bpp pix
 *      Return: 8 bpp pix, scaled down by 8x in each direction,
 *              or NULL on error.
 */
PIX *
pixScaleToGray8(PIX  *pixs)
{
l_uint8   *valtab;
l_int32    ws, hs, wd, hd;
l_int32    wpld, wpls;
l_int32   *tab8;
l_uint32  *datas, *datad;
PIX       *pixd;

    PROCNAME("pixScaleToGray8");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs must be 1 bpp", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = ws / 8;  /* truncate to nearest dest byte */
    hd = hs / 8;
    if (wd == 0 || hd == 0)
        return (PIX *)ERROR_PTR("pixs too small", procName, NULL);

    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 0.125, 0.125);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    if ((tab8 = makePixelSumTab8()) == NULL)
	return (PIX *)ERROR_PTR("tab8 not made", procName, NULL);

    if ((valtab = makeValTabSG8()) == NULL)
	return (PIX *)ERROR_PTR("valtab not made", procName, NULL);

    scaleToGray8Low(datad, wd, hd, wpld, datas, wpls, tab8, valtab);

    FREE((void *)tab8);
    FREE((void *)valtab);

    return pixd;
}


/*------------------------------------------------------------------*
 *                         Scale-to-gray 16x                         *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleToGray16()
 *
 *      Input:  1 bpp pix
 *      Return: 8 bpp pix, scaled down by 16x in each direction,
 *              or null on error.
 */
PIX *
pixScaleToGray16(PIX  *pixs)
{
l_int32    ws, hs, wd, hd;
l_int32    wpld, wpls;
l_int32   *tab8;
l_uint32  *datas, *datad;
PIX       *pixd;

    PROCNAME("pixScaleToGray16");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs must be 1 bpp", procName, NULL);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = ws / 16;
    hd = hs / 16;
    if (wd == 0 || hd == 0)
        return (PIX *)ERROR_PTR("pixs too small", procName, NULL);

    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 0.0625, 0.0625);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wpld = pixGetWpl(pixd);

    if ((tab8 = makePixelSumTab8()) == NULL)
	return (PIX *)ERROR_PTR("tab8 not made", procName, NULL);

    scaleToGray16Low(datad, wd, hd, wpld, datas, wpls, tab8);

    FREE((void *)tab8);

    return pixd;
}


/*------------------------------------------------------------------*
 *    Scale-to-gray mipmap(1 bpp --> 8 bpp, arbitrary reduction)    *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleToGrayMipmap()
 *
 *      Input:  1 bpp pix
 *              scalefactor (reduction, < 1.0)
 *      Return: 8 bpp pix, scaled down by scalefactor in each direction,
 *              or NULL on error.
 *
 *  Notes:
 *
 *  This function is here mainly for pedagogical reasons.
 *  Mip-mapping is widely used in graphics for texture mapping, because
 *  the texture changes smoothly with scale.  This is accomplished by
 *  constructing a multiresolution pyramid and, for each pixel,
 *  doing a linear interpolation between corresponding pixels in
 *  the two planes of the pyramid that bracket the desired resolution.
 *  The computation is very efficient, and is implemented in hardware
 *  in high-end graphics cards.
 *
 *  We can use mip-mapping for scale-to-gray by using two scale-to-gray
 *  reduced images (we don't need the entire pyramid) selected from
 *  the set {2x, 4x, ... 16x}, and interpolating.  However, we get
 *  severe aliasing, probably because we are subsampling from the
 *  higher resolution image.  The method is very fast, but the result
 *  is very poor.  In fact, the results don't look any better than
 *  either subsampling off the higher-res grayscale image or oversampling
 *  on the lower-res image.  Consequently, this method should NOT be used
 *  for generating reduced images, scale-to-gray or otherwise.
 */
PIX *
pixScaleToGrayMipmap(PIX       *pixs,
                     l_float32  scalefactor)
{
l_int32    w, h, minsrc, mindest;
l_float32  red;
PIX       *pixs1, *pixs2, *pixt, *pixd;

    PROCNAME("pixScaleToGrayMipmap");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 1)
        return (PIX *)ERROR_PTR("pixs not 1 bpp", procName, NULL);
    if (scalefactor >= 1.0)
        return (PIX *)ERROR_PTR("scalefactor not < 1.0", procName, NULL);
    w = pixGetWidth(pixs);
    h = pixGetHeight(pixs);
    minsrc = L_MIN(w, h);
    mindest = (l_int32)((l_float32)minsrc * scalefactor);
    if (mindest < 2)
        return (PIX *)ERROR_PTR("scalefactor too small", procName, NULL);

    if (scalefactor > 0.5) {
        pixs1 = pixConvert1To8(NULL, pixs, 255, 0);
	pixs2 = pixScaleToGray2(pixs);
	red = scalefactor;
    }
    else if (scalefactor == 0.5) {
        return pixScaleToGray2(pixs);
    }
    else if (scalefactor > 0.25) {
	pixs1 = pixScaleToGray2(pixs);
	pixs2 = pixScaleToGray4(pixs);
	red = 2. * scalefactor;
    }
    else if (scalefactor == 0.25) {
        return pixScaleToGray4(pixs);
    }
    else if (scalefactor > 0.125) {
	pixs1 = pixScaleToGray4(pixs);
	pixs2 = pixScaleToGray8(pixs);
	red = 4. * scalefactor;
    }
    else if (scalefactor == 0.125) {
        return pixScaleToGray8(pixs);
    }
    else if (scalefactor > 0.0625) {
	pixs1 = pixScaleToGray8(pixs);
	pixs2 = pixScaleToGray16(pixs);
	red = 8. * scalefactor;
    }
    else if (scalefactor == 0.0625) {
        return pixScaleToGray16(pixs);
    }
    else {  /* end of the pyramid; just do it */
	red = 16.0 * scalefactor;  /* will be <= 1.0 */
	if ((pixt = pixScaleToGray16(pixs)) == NULL)
	    return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
        if (red < 0.7)
	    pixd = pixScaleSmooth(pixt, red, red);
        else	
	    pixd = pixScaleGrayLI(pixt, red, red);
        pixDestroy(&pixt);
	return pixd;
    }

    pixd = pixScaleMipmap(pixs1, pixs2, red);

    pixDestroy(&pixs1);
    pixDestroy(&pixs2);
    return pixd;
}


/*------------------------------------------------------------------*
 *                  Grayscale scaling using mipmap                  *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleMipmap()
 *
 *      Input:  pixs1 (high res 8 bpp)
 *              pixs2 (low res -- 2x reduced -- 8 bpp)
 *              scale (reduction with respect to high res image, > 0.5)
 *      Return: 8 bpp pix, scaled down by reduction in each direction,
 *              or NULL on error.
 *
 *  Notes:
 *      (1) See notes in pixScaleToGrayMipmap().
 *      (2) This function suffers from aliasing effects that are
 *          easily seen in document images.
 */
PIX *
pixScaleMipmap(PIX       *pixs1,
               PIX       *pixs2,
               l_float32  scale)
{
l_int32    ws1, hs1, ws2, hs2, wd, hd, wpls1, wpls2, wpld;
l_uint32  *datas1, *datas2, *datad;
PIX       *pixd;

    PROCNAME("pixScaleMipmap");

    if (!pixs1)
	return (PIX *)ERROR_PTR("pixs1 not defined", procName, NULL);
    if (!pixs2)
	return (PIX *)ERROR_PTR("pixs2 not defined", procName, NULL);
    if (pixGetDepth(pixs1) != 8 || pixGetDepth(pixs2) != 8)
        return (PIX *)ERROR_PTR("pixs1, pixs2 not both 8 bpp", procName, NULL);
    if (scale > 1.0 || scale < 0.5)
        return (PIX *)ERROR_PTR("scale not in [0.5, 1.0]", procName, NULL);
    if (pixGetColormap(pixs1) || pixGetColormap(pixs2))
        L_WARNING("pixs1 or pixs2 has colormap", procName);
    ws1 = pixGetWidth(pixs1);
    ws2 = pixGetWidth(pixs2);
    if (ws1 < 2 * ws2)
        return (PIX *)ERROR_PTR("invalid width ratio", procName, NULL);
    hs1 = pixGetHeight(pixs1);
    hs2 = pixGetHeight(pixs2);
    if (hs1 < 2 * hs2)
        return (PIX *)ERROR_PTR("invalid height ratio", procName, NULL);

        /* generate wd and hd from the lower resolution dimensions,
         * to guarantee staying within both src images */
    datas1 = pixGetData(pixs1);
    wpls1 = pixGetWpl(pixs1);
    datas2 = pixGetData(pixs2);
    wpls2 = pixGetWpl(pixs2);
    wd = (l_int32)(2. * scale * pixGetWidth(pixs2));
    hd = (l_int32)(2. * scale * pixGetHeight(pixs2));
    if ((pixd = pixCreate(wd, hd, 8)) == NULL)
	return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs1);
    pixScaleResolution(pixd, scale, scale);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    scaleMipmapLow(datad, wd, hd, wpld, datas1, wpls1, datas2, wpls2, scale);
    return pixd;
}



/*------------------------------------------------------------------*
 *                Scale 2x followed by binarization                 *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleGray2xLIThresh()
 *
 *      Input:  pixs (8 bpp)
 *              thresh  (between 0 and 256)
 *      Return: pixd (1 bpp), or null on error
 *
 *  Notes:
 *      (1) This does 2x upscale on pixs, using linear interpolation,
 *          followed by thresholding to binary.
 *      (2) Buffers are used to avoid making a large grayscale image.
 */
PIX *
pixScaleGray2xLIThresh(PIX     *pixs,
                       l_int32  thresh)
{
l_int32    i, ws, hs, hsm, wd, hd, wpls, wplb, wpld;
l_uint32  *datas, *datad, *lines, *lined, *lineb;
PIX       *pixd;

    PROCNAME("pixScaleGray2xLIThresh");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (thresh < 0 || thresh > 256)
        return (PIX *)ERROR_PTR("thresh must be in [0, ... 256]",
	    procName, NULL);
    if (pixGetColormap(pixs))
        L_WARNING("pixs has colormap", procName);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = 2 * ws;
    hd = 2 * hs;
    hsm = hs - 1;
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);

        /* make line buffer for 2 lines of virtual intermediate image */
    wplb = (wd + 3) / 4;
    if ((lineb = (l_uint32 *)CALLOC(2 * wplb, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("lineb not made", procName, NULL);

	/* make dest binary image */
    if ((pixd = pixCreate(wd, hd, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 2.0, 2.0);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);

        /* do all but last src line */
    for (i = 0; i < hsm; i++) {
	lines = datas + i * wpls;
	lined = datad + 2 * i * wpld;  /* do 2 dest lines at a time */
	scaleGray2xLILineLow(lineb, wplb, lines, ws, wpls, 0);
        thresholdToBinaryLineLow(lined, wd, lineb, 8, thresh);
        thresholdToBinaryLineLow(lined + wpld, wd, lineb + wplb, 8, thresh);
    }

        /* do last src line */
    lines = datas + hsm * wpls;
    lined = datad + 2 * hsm * wpld;
    scaleGray2xLILineLow(lineb, wplb, lines, ws, wpls, 1);
    thresholdToBinaryLineLow(lined, wd, lineb, 8, thresh);
    thresholdToBinaryLineLow(lined + wpld, wd, lineb + wplb, 8, thresh);

    FREE((void *)lineb);

    return pixd;
}


/*!
 *  pixScaleGray2xLIDither()
 *
 *      Input:  pixs (8 bpp)
 *      Return: pixd (1 bpp), or null on error
 *
 *  Notes:
 *      (1) This does 2x upscale on pixs, using linear interpolation,
 *          followed by Floyd-Steinberg dithering to binary.
 *      (2) Buffers are used to avoid making a large grayscale image.
 *          - Two line buffers are used for the src, required for the 2x
 *            LI upscale.
 *          - Three line buffers are used for the intermediate image.
 *            Two are filled with each 2xLI row operation; the third is
 *            needed because the upscale and dithering ops are out of sync.
 */
PIX *
pixScaleGray2xLIDither(PIX  *pixs)
{
l_int32    i, ws, hs, hsm, wd, hd, wpls, wplb, wpld;
l_uint32  *datas, *datad;
l_uint32  *lined;
l_uint32  *lineb;   /* 2 intermediate buffer lines */
l_uint32  *linebp;  /* 1 intermediate buffer line */
l_uint32  *bufs;    /* 2 source buffer lines */
PIX       *pixd;

    PROCNAME("pixScaleGray2xLIDither");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (pixGetColormap(pixs))
        L_WARNING("pixs has colormap", procName);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = 2 * ws;
    hd = 2 * hs;
    hsm = hs - 1;
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);

        /* make line buffers for 2 lines of src image */
    if ((bufs = (l_uint32 *)CALLOC(2 * wpls, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("bufs not made", procName, NULL);

        /* make line buffer for 2 lines of virtual intermediate image */
    wplb = (wd + 3) / 4;
    if ((lineb = (l_uint32 *)CALLOC(2 * wplb, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("lineb not made", procName, NULL);

        /* make line buffer for 1 line of virtual intermediate image */
    if ((linebp = (l_uint32 *)CALLOC(wplb, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("linebp not made", procName, NULL);

	/* make dest binary image */
    if ((pixd = pixCreate(wd, hd, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 2.0, 2.0);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);

        /* start with the first src and the first dest line */
    memcpy(bufs, datas, 4 * wpls);   /* first src line */
    memcpy(bufs + wpls, datas + wpls, 4 * wpls);  /* 2nd src line */
    scaleGray2xLILineLow(lineb, wplb, bufs, ws, wpls, 0);  /* 2 i lines */
    lined = datad;
    ditherToBinaryLineLow(lined, wd, lineb, lineb + wplb,
                          DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
			                            /* 1st d line */

        /* do all but last src line */
    for (i = 1; i < hsm; i++) {
	memcpy(bufs, datas + i * wpls, 4 * wpls);  /* i-th src line */
	memcpy(bufs + wpls, datas + (i + 1) * wpls, 4 * wpls);
	memcpy(linebp, lineb + wplb, 4 * wplb);
	scaleGray2xLILineLow(lineb, wplb, bufs, ws, wpls, 0);  /* 2 i lines */
	lined = datad + 2 * i * wpld;
        ditherToBinaryLineLow(lined - wpld, wd, linebp, lineb,
                              DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
	                                           /* odd dest line */
        ditherToBinaryLineLow(lined, wd, lineb, lineb + wplb,
                              DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
	                                           /* even dest line */
    }

       /* do the last src line and the last 3 dest lines */
    memcpy(bufs, datas + hsm * wpls, 4 * wpls);  /* hsm-th src line */
    memcpy(linebp, lineb + wplb, 4 * wplb);   /* 1 i line */
    scaleGray2xLILineLow(lineb, wplb, bufs, ws, wpls, 1);  /* 2 i lines */
    ditherToBinaryLineLow(lined + wpld, wd, linebp, lineb,
                          DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
                                                   /* odd dest line */
    ditherToBinaryLineLow(lined + 2 * wpld, wd, lineb, lineb + wplb,
                          DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
                                                   /* even dest line */
    ditherToBinaryLineLow(lined + 3 * wpld, wd, lineb + wplb, NULL,
                          DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 1);
                                                   /* last dest line */

    FREE((void *)bufs);
    FREE((void *)lineb);
    FREE((void *)linebp);

    return pixd;
}


/*------------------------------------------------------------------*
 *                Scale 4x followed by binarization                 *
 *------------------------------------------------------------------*/
/*!
 *  pixScaleGray4xLIThresh()
 *
 *      Input:  pixs (8 bpp)
 *              thresh  (between 0 and 256)
 *      Return: pixd (1 bpp), or null on error
 *
 *  Notes:
 *      (1) This does 4x upscale on pixs, using linear interpolation,
 *          followed by thresholding to binary.
 *      (2) Buffers are used to avoid making a large grayscale image.
 *      (3) If a full 4x expanded grayscale image can be kept in memory,
 *          this function is only about 10% faster than separately doing
 *          a linear interpolation to a large grayscale image, followed
 *          by thresholding to binary.
 */
PIX *
pixScaleGray4xLIThresh(PIX     *pixs,
                       l_int32  thresh)
{
l_int32    i, j, ws, hs, hsm, wd, hd, wpls, wplb, wpld;
l_uint32  *datas, *datad, *lines, *lined, *lineb;
PIX       *pixd;

    PROCNAME("pixScaleGray4xLIThresh");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (thresh < 0 || thresh > 256)
        return (PIX *)ERROR_PTR("thresh must be in [0, ... 256]",
	    procName, NULL);
    if (pixGetColormap(pixs))
        L_WARNING("pixs has colormap", procName);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = 4 * ws;
    hd = 4 * hs;
    hsm = hs - 1;
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);

        /* make line buffer for 4 lines of virtual intermediate image */
    wplb = (wd + 3) / 4;
    if ((lineb = (l_uint32 *)CALLOC(4 * wplb, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("lineb not made", procName, NULL);

	/* make dest binary image */
    if ((pixd = pixCreate(wd, hd, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 4.0, 4.0);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);

        /* do all but last src line */
    for (i = 0; i < hsm; i++) {
	lines = datas + i * wpls;
	lined = datad + 4 * i * wpld;  /* do 4 dest lines at a time */
	scaleGray4xLILineLow(lineb, wplb, lines, ws, wpls, 0);
	for (j = 0; j < 4; j++) {
	    thresholdToBinaryLineLow(lined + j * wpld, wd,
	                             lineb + j * wplb, 8, thresh);
	}
    }

        /* do last src line */
    lines = datas + hsm * wpls;
    lined = datad + 4 * hsm * wpld;
    scaleGray4xLILineLow(lineb, wplb, lines, ws, wpls, 1);
    for (j = 0; j < 4; j++) {
	thresholdToBinaryLineLow(lined + j * wpld, wd,
				 lineb + j * wplb, 8, thresh);
    }

    FREE((void *)lineb);

    return pixd;
}


/*!
 *  pixScaleGray4xLIDither()
 *
 *      Input:  pixs (8 bpp)
 *      Return: pixd (1 bpp), or null on error
 *
 *  Notes:
 *      (1) This does 4x upscale on pixs, using linear interpolation,
 *          followed by Floyd-Steinberg dithering to binary.
 *      (2) Buffers are used to avoid making a large grayscale image.
 *          - Two line buffers are used for the src, required for the
 *            4xLI upscale.
 *          - Five line buffers are used for the intermediate image.
 *            Four are filled with each 4xLI row operation; the fifth
 *            is needed because the upscale and dithering ops are
 *            out of sync.
 *      (3) If a full 4x expanded grayscale image can be kept in memory,
 *          this function is only about 5% faster than separately doing
 *          a linear interpolation to a large grayscale image, followed
 *          by error-diffusion dithering to binary.
 */
PIX *
pixScaleGray4xLIDither(PIX  *pixs)
{
l_int32    i, j, ws, hs, hsm, wd, hd, wpls, wplb, wpld;
l_uint32  *datas, *datad;
l_uint32  *lined;
l_uint32  *lineb;   /* 4 intermediate buffer lines */
l_uint32  *linebp;  /* 1 intermediate buffer line */
l_uint32  *bufs;    /* 2 source buffer lines */
PIX       *pixd;

    PROCNAME("pixScaleGray4xLIDither");

    if (!pixs)
	return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (pixGetColormap(pixs))
        L_WARNING("pixs has colormap", procName);

    ws = pixGetWidth(pixs);
    hs = pixGetHeight(pixs);
    wd = 4 * ws;
    hd = 4 * hs;
    hsm = hs - 1;
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);

        /* make line buffers for 2 lines of src image */
    if ((bufs = (l_uint32 *)CALLOC(2 * wpls, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("bufs not made", procName, NULL);

        /* make line buffer for 4 lines of virtual intermediate image */
    wplb = (wd + 3) / 4;
    if ((lineb = (l_uint32 *)CALLOC(4 * wplb, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("lineb not made", procName, NULL);

        /* make line buffer for 1 line of virtual intermediate image */
    if ((linebp = (l_uint32 *)CALLOC(wplb, sizeof(l_uint32))) == NULL)
        return (PIX *)ERROR_PTR("linebp not made", procName, NULL);

	/* make dest binary image */
    if ((pixd = pixCreate(wd, hd, 1)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    pixCopyResolution(pixd, pixs);
    pixScaleResolution(pixd, 4.0, 4.0);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);

        /* start with the first src and the first 3 dest lines */
    memcpy(bufs, datas, 4 * wpls);   /* first src line */
    memcpy(bufs + wpls, datas + wpls, 4 * wpls);  /* 2nd src line */
    scaleGray4xLILineLow(lineb, wplb, bufs, ws, wpls, 0);  /* 4 b lines */
    lined = datad;
    for (j = 0; j < 3; j++) {  /* first 3 d lines of Q */
	ditherToBinaryLineLow(lined + j * wpld, wd, lineb + j * wplb,
	                      lineb + (j + 1) * wplb,
                              DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
    }

        /* do all but last src line */
    for (i = 1; i < hsm; i++) {
	memcpy(bufs, datas + i * wpls, 4 * wpls);  /* i-th src line */
	memcpy(bufs + wpls, datas + (i + 1) * wpls, 4 * wpls);
	memcpy(linebp, lineb + 3 * wplb, 4 * wplb);
	scaleGray4xLILineLow(lineb, wplb, bufs, ws, wpls, 0);  /* 4 b lines */
	lined = datad + 4 * i * wpld;
	ditherToBinaryLineLow(lined - wpld, wd, linebp, lineb,
                              DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
	                                             /* 4th dest line of Q */
	for (j = 0; j < 3; j++) {  /* next 3 d lines of Quad */
            ditherToBinaryLineLow(lined + j * wpld, wd, lineb + j * wplb,
	                          lineb + (j + 1) * wplb,
                                 DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
        }			     
    }

       /* do the last src line and the last 5 dest lines */
    memcpy(bufs, datas + hsm * wpls, 4 * wpls);  /* hsm-th src line */
    memcpy(linebp, lineb + 3 * wplb, 4 * wplb);   /* 1 b line */
    scaleGray4xLILineLow(lineb, wplb, bufs, ws, wpls, 1);  /* 4 b lines */
    lined = datad + 4 * hsm * wpld;
    ditherToBinaryLineLow(lined - wpld, wd, linebp, lineb,
                          DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
                                                   /* 4th dest line of Q */
    for (j = 0; j < 3; j++) {  /* next 3 d lines of Quad */
	ditherToBinaryLineLow(lined + j * wpld, wd, lineb + j * wplb,
		              lineb + (j + 1) * wplb,
                              DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 0);
    }
       /* and finally, the last dest line */
    ditherToBinaryLineLow(lined + 3 * wpld, wd, lineb + 3 * wplb, NULL,
                              DEFAULT_CLIP_LOWER_1, DEFAULT_CLIP_UPPER_1, 1);

    FREE((void *)bufs);
    FREE((void *)lineb);
    FREE((void *)linebp);

    return pixd;
}

