affine.c

Go to the documentation of this file.

/*====================================================================*
 -  Copyright (C) 2001 Leptonica.  All rights reserved.
 -
 -  Redistribution and use in source and binary forms, with or without
 -  modification, are permitted provided that the following conditions
 -  are met:
 -  1. Redistributions of source code must retain the above copyright
 -     notice, this list of conditions and the following disclaimer.
 -  2. Redistributions in binary form must reproduce the above
 -     copyright notice, this list of conditions and the following
 -     disclaimer in the documentation and/or other materials
 -     provided with the distribution.
 -
 -  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 -  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 -  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 -  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL ANY
 -  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 -  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 -  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 -  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
 -  OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 -  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 -  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *====================================================================*/


#include <string.h>
#include <math.h>
#include "allheaders.h"

extern l_float32  AlphaMaskBorderVals[2];

#ifndef  NO_CONSOLE_IO
#define  DEBUG     0
#endif  /* ~NO_CONSOLE_IO */


/*-------------------------------------------------------------*
 *               Sampled affine image transformation           *
 *-------------------------------------------------------------*/
PIX *
pixAffineSampledPta(PIX     *pixs,
                    PTA     *ptad,
                    PTA     *ptas,
                    l_int32  incolor)
{
l_float32  *vc;
PIX        *pixd;

    PROCNAME("pixAffineSampledPta");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!ptas)
        return (PIX *)ERROR_PTR("ptas not defined", procName, NULL);
    if (!ptad)
        return (PIX *)ERROR_PTR("ptad not defined", procName, NULL);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return (PIX *)ERROR_PTR("invalid incolor", procName, NULL);
    if (ptaGetCount(ptas) != 3)
        return (PIX *)ERROR_PTR("ptas count not 3", procName, NULL);
    if (ptaGetCount(ptad) != 3)
        return (PIX *)ERROR_PTR("ptad count not 3", procName, NULL);

        /* Get backwards transform from dest to src, and apply it */
    getAffineXformCoeffs(ptad, ptas, &vc);
    pixd = pixAffineSampled(pixs, vc, incolor);
    LEPT_FREE(vc);

    return pixd;
}


PIX *
pixAffineSampled(PIX        *pixs,
                 l_float32  *vc,
                 l_int32     incolor)
{
l_int32     i, j, w, h, d, x, y, wpls, wpld, color, cmapindex;
l_uint32    val;
l_uint32   *datas, *datad, *lines, *lined;
PIX        *pixd;
PIXCMAP    *cmap;

    PROCNAME("pixAffineSampled");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!vc)
        return (PIX *)ERROR_PTR("vc not defined", procName, NULL);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return (PIX *)ERROR_PTR("invalid incolor", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 1 && d != 2 && d != 4 && d != 8 && d != 32)
        return (PIX *)ERROR_PTR("depth not 1, 2, 4, 8 or 16", procName, NULL);

        /* Init all dest pixels to color to be brought in from outside */
    pixd = pixCreateTemplate(pixs);
    if ((cmap = pixGetColormap(pixs)) != NULL) {
        if (incolor == L_BRING_IN_WHITE)
            color = 1;
        else
            color = 0;
        pixcmapAddBlackOrWhite(cmap, color, &cmapindex);
        pixSetAllArbitrary(pixd, cmapindex);
    } else {
        if ((d == 1 && incolor == L_BRING_IN_WHITE) ||
            (d > 1 && incolor == L_BRING_IN_BLACK)) {
            pixClearAll(pixd);
        } else {
            pixSetAll(pixd);
        }
    }

        /* Scan over the dest pixels */
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            affineXformSampledPt(vc, j, i, &x, &y);
            if (x < 0 || y < 0 || x >=w || y >= h)
                continue;
            lines = datas + y * wpls;
            if (d == 1) {
                val = GET_DATA_BIT(lines, x);
                SET_DATA_BIT_VAL(lined, j, val);
            } else if (d == 8) {
                val = GET_DATA_BYTE(lines, x);
                SET_DATA_BYTE(lined, j, val);
            } else if (d == 32) {
                lined[j] = lines[x];
            } else if (d == 2) {
                val = GET_DATA_DIBIT(lines, x);
                SET_DATA_DIBIT(lined, j, val);
            } else if (d == 4) {
                val = GET_DATA_QBIT(lines, x);
                SET_DATA_QBIT(lined, j, val);
            }
        }
    }

    return pixd;
}


/*---------------------------------------------------------------------*
 *               Interpolated affine image transformation              *
 *---------------------------------------------------------------------*/
PIX *
pixAffinePta(PIX     *pixs,
             PTA     *ptad,
             PTA     *ptas,
             l_int32  incolor)
{
l_int32   d;
l_uint32  colorval;
PIX      *pixt1, *pixt2, *pixd;

    PROCNAME("pixAffinePta");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!ptas)
        return (PIX *)ERROR_PTR("ptas not defined", procName, NULL);
    if (!ptad)
        return (PIX *)ERROR_PTR("ptad not defined", procName, NULL);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return (PIX *)ERROR_PTR("invalid incolor", procName, NULL);
    if (ptaGetCount(ptas) != 3)
        return (PIX *)ERROR_PTR("ptas count not 3", procName, NULL);
    if (ptaGetCount(ptad) != 3)
        return (PIX *)ERROR_PTR("ptad count not 3", procName, NULL);

    if (pixGetDepth(pixs) == 1)
        return pixAffineSampledPta(pixs, ptad, ptas, incolor);

        /* Remove cmap if it exists, and unpack to 8 bpp if necessary */
    pixt1 = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
    d = pixGetDepth(pixt1);
    if (d < 8)
        pixt2 = pixConvertTo8(pixt1, FALSE);
    else
        pixt2 = pixClone(pixt1);
    d = pixGetDepth(pixt2);

        /* Compute actual color to bring in from edges */
    colorval = 0;
    if (incolor == L_BRING_IN_WHITE) {
        if (d == 8)
            colorval = 255;
        else  /* d == 32 */
            colorval = 0xffffff00;
    }

    if (d == 8)
        pixd = pixAffinePtaGray(pixt2, ptad, ptas, colorval);
    else  /* d == 32 */
        pixd = pixAffinePtaColor(pixt2, ptad, ptas, colorval);
    pixDestroy(&pixt1);
    pixDestroy(&pixt2);
    return pixd;
}


PIX *
pixAffine(PIX        *pixs,
          l_float32  *vc,
          l_int32     incolor)
{
l_int32   d;
l_uint32  colorval;
PIX      *pixt1, *pixt2, *pixd;

    PROCNAME("pixAffine");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!vc)
        return (PIX *)ERROR_PTR("vc not defined", procName, NULL);

    if (pixGetDepth(pixs) == 1)
        return pixAffineSampled(pixs, vc, incolor);

        /* Remove cmap if it exists, and unpack to 8 bpp if necessary */
    pixt1 = pixRemoveColormap(pixs, REMOVE_CMAP_BASED_ON_SRC);
    d = pixGetDepth(pixt1);
    if (d < 8)
        pixt2 = pixConvertTo8(pixt1, FALSE);
    else
        pixt2 = pixClone(pixt1);
    d = pixGetDepth(pixt2);

        /* Compute actual color to bring in from edges */
    colorval = 0;
    if (incolor == L_BRING_IN_WHITE) {
        if (d == 8)
            colorval = 255;
        else  /* d == 32 */
            colorval = 0xffffff00;
    }

    if (d == 8)
        pixd = pixAffineGray(pixt2, vc, colorval);
    else  /* d == 32 */
        pixd = pixAffineColor(pixt2, vc, colorval);
    pixDestroy(&pixt1);
    pixDestroy(&pixt2);
    return pixd;
}


PIX *
pixAffinePtaColor(PIX      *pixs,
                  PTA      *ptad,
                  PTA      *ptas,
                  l_uint32  colorval)
{
l_float32  *vc;
PIX        *pixd;

    PROCNAME("pixAffinePtaColor");

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

        /* Get backwards transform from dest to src, and apply it */
    getAffineXformCoeffs(ptad, ptas, &vc);
    pixd = pixAffineColor(pixs, vc, colorval);
    LEPT_FREE(vc);

    return pixd;
}


PIX *
pixAffineColor(PIX        *pixs,
               l_float32  *vc,
               l_uint32    colorval)
{
l_int32    i, j, w, h, d, wpls, wpld;
l_uint32   val;
l_uint32  *datas, *datad, *lined;
l_float32  x, y;
PIX       *pix1, *pix2, *pixd;

    PROCNAME("pixAffineColor");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32)
        return (PIX *)ERROR_PTR("pixs must be 32 bpp", procName, NULL);
    if (!vc)
        return (PIX *)ERROR_PTR("vc not defined", procName, NULL);

    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    pixd = pixCreateTemplate(pixs);
    pixSetAllArbitrary(pixd, colorval);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

        /* Iterate over destination pixels */
    for (i = 0; i < h; i++) {
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
                /* Compute float src pixel location corresponding to (i,j) */
            affineXformPt(vc, j, i, &x, &y);
            linearInterpolatePixelColor(datas, wpls, w, h, x, y, colorval,
                                        &val);
            *(lined + j) = val;
        }
    }

        /* If rgba, transform the pixs alpha channel and insert in pixd */
    if (pixGetSpp(pixs) == 4) {
        pix1 = pixGetRGBComponent(pixs, L_ALPHA_CHANNEL);
        pix2 = pixAffineGray(pix1, vc, 255);  /* bring in opaque */
        pixSetRGBComponent(pixd, pix2, L_ALPHA_CHANNEL);
        pixDestroy(&pix1);
        pixDestroy(&pix2);
    }

    return pixd;
}


PIX *
pixAffinePtaGray(PIX     *pixs,
                 PTA     *ptad,
                 PTA     *ptas,
                 l_uint8  grayval)
{
l_float32  *vc;
PIX        *pixd;

    PROCNAME("pixAffinePtaGray");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!ptas)
        return (PIX *)ERROR_PTR("ptas not defined", procName, NULL);
    if (!ptad)
        return (PIX *)ERROR_PTR("ptad not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (ptaGetCount(ptas) != 3)
        return (PIX *)ERROR_PTR("ptas count not 3", procName, NULL);
    if (ptaGetCount(ptad) != 3)
        return (PIX *)ERROR_PTR("ptad count not 3", procName, NULL);

        /* Get backwards transform from dest to src, and apply it */
    getAffineXformCoeffs(ptad, ptas, &vc);
    pixd = pixAffineGray(pixs, vc, grayval);
    LEPT_FREE(vc);

    return pixd;
}



PIX *
pixAffineGray(PIX        *pixs,
              l_float32  *vc,
              l_uint8     grayval)
{
l_int32    i, j, w, h, wpls, wpld, val;
l_uint32  *datas, *datad, *lined;
l_float32  x, y;
PIX       *pixd;

    PROCNAME("pixAffineGray");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs must be 8 bpp", procName, NULL);
    if (!vc)
        return (PIX *)ERROR_PTR("vc not defined", procName, NULL);

    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    pixd = pixCreateTemplate(pixs);
    pixSetAllArbitrary(pixd, grayval);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

        /* Iterate over destination pixels */
    for (i = 0; i < h; i++) {
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
                /* Compute float src pixel location corresponding to (i,j) */
            affineXformPt(vc, j, i, &x, &y);
            linearInterpolatePixelGray(datas, wpls, w, h, x, y, grayval, &val);
            SET_DATA_BYTE(lined, j, val);
        }
    }

    return pixd;
}


/*---------------------------------------------------------------------------*
 *            Affine transform including alpha (blend) component             *
 *---------------------------------------------------------------------------*/
PIX *
pixAffinePtaWithAlpha(PIX       *pixs,
                      PTA       *ptad,
                      PTA       *ptas,
                      PIX       *pixg,
                      l_float32  fract,
                      l_int32    border)
{
l_int32  ws, hs, d;
PIX     *pixd, *pixb1, *pixb2, *pixg2, *pixga;
PTA     *ptad2, *ptas2;

    PROCNAME("pixAffinePtaWithAlpha");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &ws, &hs, &d);
    if (d != 32 && pixGetColormap(pixs) == NULL)
        return (PIX *)ERROR_PTR("pixs not cmapped or 32 bpp", procName, NULL);
    if (pixg && pixGetDepth(pixg) != 8) {
        L_WARNING("pixg not 8 bpp; using 'fract' transparent alpha\n",
                  procName);
        pixg = NULL;
    }
    if (!pixg && (fract < 0.0 || fract > 1.0)) {
        L_WARNING("invalid fract; using 1.0 (fully transparent)\n", procName);
        fract = 1.0;
    }
    if (!pixg && fract == 0.0)
        L_WARNING("fully opaque alpha; image will not be blended\n", procName);
    if (!ptad)
        return (PIX *)ERROR_PTR("ptad not defined", procName, NULL);
    if (!ptas)
        return (PIX *)ERROR_PTR("ptas not defined", procName, NULL);

        /* Add border; the color doesn't matter */
    pixb1 = pixAddBorder(pixs, border, 0);

        /* Transform the ptr arrays to work on the bordered image */
    ptad2 = ptaTransform(ptad, border, border, 1.0, 1.0);
    ptas2 = ptaTransform(ptas, border, border, 1.0, 1.0);

        /* Do separate affine transform of rgb channels of pixs and of pixg */
    pixd = pixAffinePtaColor(pixb1, ptad2, ptas2, 0);
    if (!pixg) {
        pixg2 = pixCreate(ws, hs, 8);
        if (fract == 1.0)
            pixSetAll(pixg2);
        else
            pixSetAllArbitrary(pixg2, (l_int32)(255.0 * fract));
    } else {
        pixg2 = pixResizeToMatch(pixg, NULL, ws, hs);
    }
    if (ws > 10 && hs > 10) {  /* see note 7 */
        pixSetBorderRingVal(pixg2, 1,
                            (l_int32)(255.0 * fract * AlphaMaskBorderVals[0]));
        pixSetBorderRingVal(pixg2, 2,
                            (l_int32)(255.0 * fract * AlphaMaskBorderVals[1]));

    }
    pixb2 = pixAddBorder(pixg2, border, 0);  /* must be black border */
    pixga = pixAffinePtaGray(pixb2, ptad2, ptas2, 0);
    pixSetRGBComponent(pixd, pixga, L_ALPHA_CHANNEL);
    pixSetSpp(pixd, 4);

    pixDestroy(&pixg2);
    pixDestroy(&pixb1);
    pixDestroy(&pixb2);
    pixDestroy(&pixga);
    ptaDestroy(&ptad2);
    ptaDestroy(&ptas2);
    return pixd;
}


/*-------------------------------------------------------------*
 *                 Affine coordinate transformation            *
 *-------------------------------------------------------------*/
l_ok
getAffineXformCoeffs(PTA         *ptas,
                     PTA         *ptad,
                     l_float32  **pvc)
{
l_int32     i;
l_float32   x1, y1, x2, y2, x3, y3;
l_float32  *b;   /* rhs vector of primed coords X'; coeffs returned in *pvc */
l_float32  *a[6];  /* 6x6 matrix A  */

    PROCNAME("getAffineXformCoeffs");

    if (!ptas)
        return ERROR_INT("ptas not defined", procName, 1);
    if (!ptad)
        return ERROR_INT("ptad not defined", procName, 1);
    if (!pvc)
        return ERROR_INT("&vc not defined", procName, 1);

    if ((b = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32))) == NULL)
        return ERROR_INT("b not made", procName, 1);
    *pvc = b;

    ptaGetPt(ptas, 0, &x1, &y1);
    ptaGetPt(ptas, 1, &x2, &y2);
    ptaGetPt(ptas, 2, &x3, &y3);
    ptaGetPt(ptad, 0, &b[0], &b[1]);
    ptaGetPt(ptad, 1, &b[2], &b[3]);
    ptaGetPt(ptad, 2, &b[4], &b[5]);

    for (i = 0; i < 6; i++)
        if ((a[i] = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32))) == NULL)
            return ERROR_INT("a[i] not made", procName, 1);

    a[0][0] = x1;
    a[0][1] = y1;
    a[0][2] = 1.;
    a[1][3] = x1;
    a[1][4] = y1;
    a[1][5] = 1.;
    a[2][0] = x2;
    a[2][1] = y2;
    a[2][2] = 1.;
    a[3][3] = x2;
    a[3][4] = y2;
    a[3][5] = 1.;
    a[4][0] = x3;
    a[4][1] = y3;
    a[4][2] = 1.;
    a[5][3] = x3;
    a[5][4] = y3;
    a[5][5] = 1.;

    gaussjordan(a, b, 6);

    for (i = 0; i < 6; i++)
        LEPT_FREE(a[i]);

    return 0;
}


l_ok
affineInvertXform(l_float32   *vc,
                  l_float32  **pvci)
{
l_int32     i;
l_float32  *vci;
l_float32  *a[3];
l_float32   b[3] = {1.0, 1.0, 1.0};   /* anything; results ignored */

    PROCNAME("affineInvertXform");

    if (!pvci)
        return ERROR_INT("&vci not defined", procName, 1);
    *pvci = NULL;
    if (!vc)
        return ERROR_INT("vc not defined", procName, 1);

#if 1
    for (i = 0; i < 3; i++)
        a[i] = (l_float32 *)LEPT_CALLOC(3, sizeof(l_float32));
    a[0][0] = vc[0];
    a[0][1] = vc[1];
    a[0][2] = vc[2];
    a[1][0] = vc[3];
    a[1][1] = vc[4];
    a[1][2] = vc[5];
    a[2][2] = 1.0;
    gaussjordan(a, b, 3);  /* this inverts matrix a */
    vci = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32));
    *pvci = vci;
    vci[0] = a[0][0];
    vci[1] = a[0][1];
    vci[2] = a[0][2];
    vci[3] = a[1][0];
    vci[4] = a[1][1];
    vci[5] = a[1][2];
    for (i = 0; i < 3; i++)
        LEPT_FREE(a[i]);

#else

        /* Alternative version, inverting a 2x2 matrix */
    { l_float32 *a2[2];
    for (i = 0; i < 2; i++)
        a2[i] = (l_float32 *)LEPT_CALLOC(2, sizeof(l_float32));
    a2[0][0] = vc[0];
    a2[0][1] = vc[1];
    a2[1][0] = vc[3];
    a2[1][1] = vc[4];
    b[0] = vc[2];
    b[1] = vc[5];
    gaussjordan(a2, b, 2);  /* this inverts matrix a2 */
    vci = (l_float32 *)LEPT_CALLOC(6, sizeof(l_float32));
    *pvci = vci;
    vci[0] = a2[0][0];
    vci[1] = a2[0][1];
    vci[2] = -b[0];   /* note sign */
    vci[3] = a2[1][0];
    vci[4] = a2[1][1];
    vci[5] = -b[1];   /* note sign */
    for (i = 0; i < 2; i++)
        LEPT_FREE(a2[i]);
    }
#endif

    return 0;
}


l_ok
affineXformSampledPt(l_float32  *vc,
                     l_int32     x,
                     l_int32     y,
                     l_int32    *pxp,
                     l_int32    *pyp)
{
    PROCNAME("affineXformSampledPt");

    if (!vc)
        return ERROR_INT("vc not defined", procName, 1);

    *pxp = (l_int32)(vc[0] * x + vc[1] * y + vc[2] + 0.5);
    *pyp = (l_int32)(vc[3] * x + vc[4] * y + vc[5] + 0.5);
    return 0;
}


l_ok
affineXformPt(l_float32  *vc,
              l_int32     x,
              l_int32     y,
              l_float32  *pxp,
              l_float32  *pyp)
{
    PROCNAME("affineXformPt");

    if (!vc)
        return ERROR_INT("vc not defined", procName, 1);

    *pxp = vc[0] * x + vc[1] * y + vc[2];
    *pyp = vc[3] * x + vc[4] * y + vc[5];
    return 0;
}


/*-------------------------------------------------------------*
 *                 Interpolation helper functions              *
 *-------------------------------------------------------------*/
l_ok
linearInterpolatePixelColor(l_uint32  *datas,
                            l_int32    wpls,
                            l_int32    w,
                            l_int32    h,
                            l_float32  x,
                            l_float32  y,
                            l_uint32   colorval,
                            l_uint32  *pval)
{
l_int32    valid, xpm, ypm, xp, xp2, yp, xf, yf;
l_int32    rval, gval, bval;
l_uint32   word00, word01, word10, word11;
l_uint32  *lines;

    PROCNAME("linearInterpolatePixelColor");

    if (!pval)
        return ERROR_INT("&val not defined", procName, 1);
    *pval = colorval;
    if (!datas)
        return ERROR_INT("datas not defined", procName, 1);

        /* Skip if x or y are invalid. (x,y) must be in the source image.
         * Failure to detect an invalid point will cause a mem address fault.
         * Occasionally, x or y will be a nan, and relational checks always
         * fail for nans.  Therefore we check if the point is inside the pix */
    valid = (x >= 0.0 && y >= 0.0 && x < w && y < h);
    if (!valid) return 0;

    xpm = (l_int32)(16.0 * x);
    ypm = (l_int32)(16.0 * y);
    xp = xpm >> 4;
    xp2 = xp + 1 < w ? xp + 1 : xp;
    yp = ypm >> 4;
    if (yp + 1 >= h) wpls = 0;
    xf = xpm & 0x0f;
    yf = ypm & 0x0f;

#if  DEBUG
    if (xf < 0 || yf < 0)
        fprintf(stderr, "xp = %d, yp = %d, xf = %d, yf = %d\n", xp, yp, xf, yf);
#endif  /* DEBUG */

        /* Do area weighting (eqiv. to linear interpolation) */
    lines = datas + yp * wpls;
    word00 = *(lines + xp);
    word10 = *(lines + xp2);
    word01 = *(lines + wpls + xp);
    word11 = *(lines + wpls + xp2);
    rval = ((16 - xf) * (16 - yf) * ((word00 >> L_RED_SHIFT) & 0xff) +
        xf * (16 - yf) * ((word10 >> L_RED_SHIFT) & 0xff) +
        (16 - xf) * yf * ((word01 >> L_RED_SHIFT) & 0xff) +
        xf * yf * ((word11 >> L_RED_SHIFT) & 0xff)) / 256;
    gval = ((16 - xf) * (16 - yf) * ((word00 >> L_GREEN_SHIFT) & 0xff) +
        xf * (16 - yf) * ((word10 >> L_GREEN_SHIFT) & 0xff) +
        (16 - xf) * yf * ((word01 >> L_GREEN_SHIFT) & 0xff) +
        xf * yf * ((word11 >> L_GREEN_SHIFT) & 0xff)) / 256;
    bval = ((16 - xf) * (16 - yf) * ((word00 >> L_BLUE_SHIFT) & 0xff) +
        xf * (16 - yf) * ((word10 >> L_BLUE_SHIFT) & 0xff) +
        (16 - xf) * yf * ((word01 >> L_BLUE_SHIFT) & 0xff) +
        xf * yf * ((word11 >> L_BLUE_SHIFT) & 0xff)) / 256;
    composeRGBPixel(rval, gval, bval, pval);
    return 0;
}


l_ok
linearInterpolatePixelGray(l_uint32  *datas,
                           l_int32    wpls,
                           l_int32    w,
                           l_int32    h,
                           l_float32  x,
                           l_float32  y,
                           l_int32    grayval,
                           l_int32   *pval)
{
l_int32    valid, xpm, ypm, xp, xp2, yp, xf, yf, v00, v10, v01, v11;
l_uint32  *lines;

    PROCNAME("linearInterpolatePixelGray");

    if (!pval)
        return ERROR_INT("&val not defined", procName, 1);
    *pval = grayval;
    if (!datas)
        return ERROR_INT("datas not defined", procName, 1);

        /* Skip if x or y is invalid. (x,y) must be in the source image.
         * Failure to detect an invalid point will cause a mem address fault.
         * Occasionally, x or y will be a nan, and relational checks always
         * fail for nans.  Therefore we check if the point is inside the pix */
    valid = (x >= 0.0 && y >= 0.0 && x < w && y < h);
    if (!valid) return 0;

    xpm = (l_int32)(16.0 * x);
    ypm = (l_int32)(16.0 * y);
    xp = xpm >> 4;
    xp2 = xp + 1 < w ? xp + 1 : xp;
    yp = ypm >> 4;
    if (yp + 1 >= h) wpls = 0;
    xf = xpm & 0x0f;
    yf = ypm & 0x0f;

#if  DEBUG
    if (xf < 0 || yf < 0)
        fprintf(stderr, "xp = %d, yp = %d, xf = %d, yf = %d\n", xp, yp, xf, yf);
#endif  /* DEBUG */

        /* Interpolate by area weighting. */
    lines = datas + yp * wpls;
    v00 = (16 - xf) * (16 - yf) * GET_DATA_BYTE(lines, xp);
    v10 = xf * (16 - yf) * GET_DATA_BYTE(lines, xp2);
    v01 = (16 - xf) * yf * GET_DATA_BYTE(lines + wpls, xp);
    v11 = xf * yf * GET_DATA_BYTE(lines + wpls, xp2);
    *pval = (v00 + v01 + v10 + v11) / 256;
    return 0;
}



/*-------------------------------------------------------------*
 *               Gauss-jordan linear equation solver           *
 *-------------------------------------------------------------*/
#define  SWAP(a,b)   {temp = (a); (a) = (b); (b) = temp;}

l_int32
gaussjordan(l_float32  **a,
            l_float32   *b,
            l_int32      n)
{
l_int32    i, icol, irow, j, k, col, row, success;
l_int32   *indexc, *indexr, *ipiv;
l_float32  maxval, val, pivinv, temp;

    PROCNAME("gaussjordan");

    if (!a)
        return ERROR_INT("a not defined", procName, 1);
    if (!b)
        return ERROR_INT("b not defined", procName, 1);

    success = TRUE;
    indexc = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
    indexr = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
    ipiv = (l_int32 *)LEPT_CALLOC(n, sizeof(l_int32));
    if (!indexc || !indexr || !ipiv) {
        L_ERROR("array not made\n", procName);
        success = FALSE;
        goto cleanup_arrays;
    }

    icol = irow = 0;  /* silence static checker */
    for (i = 0; i < n; i++) {
        maxval = 0.0;
        for (j = 0; j < n; j++) {
            if (ipiv[j] != 1) {
                for (k = 0; k < n; k++) {
                    if (ipiv[k] == 0) {
                        if (fabs(a[j][k]) >= maxval) {
                            maxval = fabs(a[j][k]);
                            irow = j;
                            icol = k;
                        }
                    } else if (ipiv[k] > 1) {
                        L_ERROR("singular matrix\n", procName);
                        success = FALSE;
                        goto cleanup_arrays;
                    }
                }
            }
        }
        ++(ipiv[icol]);

        if (irow != icol) {
            for (col = 0; col < n; col++)
                SWAP(a[irow][col], a[icol][col]);
            SWAP(b[irow], b[icol]);
        }

        indexr[i] = irow;
        indexc[i] = icol;
        if (a[icol][icol] == 0.0) {
            L_ERROR("singular matrix\n", procName);
            success = FALSE;
            goto cleanup_arrays;
        }
        pivinv = 1.0 / a[icol][icol];
        a[icol][icol] = 1.0;
        for (col = 0; col < n; col++)
            a[icol][col] *= pivinv;
        b[icol] *= pivinv;

        for (row = 0; row < n; row++) {
            if (row != icol) {
                val = a[row][icol];
                a[row][icol] = 0.0;
                for (col = 0; col < n; col++)
                    a[row][col] -= a[icol][col] * val;
                b[row] -= b[icol] * val;
            }
        }
    }

    for (col = n - 1; col >= 0; col--) {
        if (indexr[col] != indexc[col]) {
            for (k = 0; k < n; k++)
                SWAP(a[k][indexr[col]], a[k][indexc[col]]);
        }
    }

cleanup_arrays:
    LEPT_FREE(indexr);
    LEPT_FREE(indexc);
    LEPT_FREE(ipiv);
    return (success) ? 0 : 1;
}


/*-------------------------------------------------------------*
 *              Sequential affine image transformation         *
 *-------------------------------------------------------------*/
PIX *
pixAffineSequential(PIX     *pixs,
                    PTA     *ptad,
                    PTA     *ptas,
                    l_int32  bw,
                    l_int32  bh)
{
l_int32    x1, y1, x2, y2, x3, y3;    /* ptas */
l_int32    x1p, y1p, x2p, y2p, x3p, y3p;   /* ptad */
l_int32    x1sc, y1sc;  /* scaled origin */
l_float32  x2s, x2sp, scalex, scaley;
l_float32  th3, th3p, ph2, ph2p;
#if  DEBUG
l_float32  rad2deg;
#endif  /* DEBUG */
PIX       *pix1, *pix2, *pixd;

    PROCNAME("pixAffineSequential");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (!ptas)
        return (PIX *)ERROR_PTR("ptas not defined", procName, NULL);
    if (!ptad)
        return (PIX *)ERROR_PTR("ptad not defined", procName, NULL);

    if (ptaGetCount(ptas) != 3)
        return (PIX *)ERROR_PTR("ptas count not 3", procName, NULL);
    if (ptaGetCount(ptad) != 3)
        return (PIX *)ERROR_PTR("ptad count not 3", procName, NULL);
    ptaGetIPt(ptas, 0, &x1, &y1);
    ptaGetIPt(ptas, 1, &x2, &y2);
    ptaGetIPt(ptas, 2, &x3, &y3);
    ptaGetIPt(ptad, 0, &x1p, &y1p);
    ptaGetIPt(ptad, 1, &x2p, &y2p);
    ptaGetIPt(ptad, 2, &x3p, &y3p);

    pix1 = pix2 = pixd = NULL;

    if (y1 == y3)
        return (PIX *)ERROR_PTR("y1 == y3!", procName, NULL);
    if (y1p == y3p)
        return (PIX *)ERROR_PTR("y1p == y3p!", procName, NULL);

    if (bw != 0 || bh != 0) {
            /* resize all points and add border to pixs */
        x1 = x1 + bw;
        y1 = y1 + bh;
        x2 = x2 + bw;
        y2 = y2 + bh;
        x3 = x3 + bw;
        y3 = y3 + bh;
        x1p = x1p + bw;
        y1p = y1p + bh;
        x2p = x2p + bw;
        y2p = y2p + bh;
        x3p = x3p + bw;
        y3p = y3p + bh;

        if ((pix1 = pixAddBorderGeneral(pixs, bw, bw, bh, bh, 0)) == NULL)
            return (PIX *)ERROR_PTR("pix1 not made", procName, NULL);
    } else {
        pix1 = pixCopy(NULL, pixs);
    }

    /*-------------------------------------------------------------*
        The horizontal shear is done to move the 3rd point to the
        y axis.  This moves the 2nd point either towards or away
        from the y axis, depending on whether it is above or below
        the x axis.  That motion must be computed so that we know
        the angle of vertical shear to use to get the 2nd point
        on the x axis.  We must also know the x coordinate of the
        2nd point in order to compute how much scaling is required
        to match points on the axis.
     *-------------------------------------------------------------*/

        /* Shear angles required to put src points on x and y axes */
    th3 = atan2((l_float64)(x1 - x3), (l_float64)(y1 - y3));
    x2s = (l_float32)(x2 - ((l_float32)(y1 - y2) * (x3 - x1)) / (y1 - y3));
    if (x2s == (l_float32)x1) {
        L_ERROR("x2s == x1!\n", procName);
        goto cleanup_pix;
    }
    ph2 = atan2((l_float64)(y1 - y2), (l_float64)(x2s - x1));

        /* Shear angles required to put dest points on x and y axes.
         * Use the negative of these values to instead move the
         * src points from the axes to the actual dest position.
         * These values are also needed to scale the image. */
    th3p = atan2((l_float64)(x1p - x3p), (l_float64)(y1p - y3p));
    x2sp = (l_float32)(x2p -
                       ((l_float32)(y1p - y2p) * (x3p - x1p)) / (y1p - y3p));
    if (x2sp == (l_float32)x1p) {
        L_ERROR("x2sp == x1p!\n", procName);
        goto cleanup_pix;
    }
    ph2p = atan2((l_float64)(y1p - y2p), (l_float64)(x2sp - x1p));

        /* Shear image to first put src point 3 on the y axis,
         * and then to put src point 2 on the x axis */
    pixHShearIP(pix1, y1, th3, L_BRING_IN_WHITE);
    pixVShearIP(pix1, x1, ph2, L_BRING_IN_WHITE);

        /* Scale image to match dest scale.  The dest scale
         * is calculated above from the angles th3p and ph2p
         * that would be required to move the dest points to
         * the x and y axes. */
    scalex = (l_float32)(x2sp - x1p) / (x2s - x1);
    scaley = (l_float32)(y3p - y1p) / (y3 - y1);
    if ((pix2 = pixScale(pix1, scalex, scaley)) == NULL) {
        L_ERROR("pix2 not made\n", procName);
        goto cleanup_pix;
    }

#if  DEBUG
    rad2deg = 180. / 3.1415926535;
    fprintf(stderr, "th3 = %5.1f deg, ph2 = %5.1f deg\n",
            rad2deg * th3, rad2deg * ph2);
    fprintf(stderr, "th3' = %5.1f deg, ph2' = %5.1f deg\n",
            rad2deg * th3p, rad2deg * ph2p);
    fprintf(stderr, "scalex = %6.3f, scaley = %6.3f\n", scalex, scaley);
#endif  /* DEBUG */

    /*-------------------------------------------------------------*
        Scaling moves the 1st src point, which is the origin.
        It must now be moved again to coincide with the origin
        (1st point) of the dest.  After this is done, the 2nd
        and 3rd points must be sheared back to the original
        positions of the 2nd and 3rd dest points.  We use the
        negative of the angles that were previously computed
        for shearing those points in the dest image to x and y
        axes, and take the shears in reverse order as well.
     *-------------------------------------------------------------*/
        /* Shift image to match dest origin. */
    x1sc = (l_int32)(scalex * x1 + 0.5);   /* x comp of origin after scaling */
    y1sc = (l_int32)(scaley * y1 + 0.5);   /* y comp of origin after scaling */
    pixRasteropIP(pix2, x1p - x1sc, y1p - y1sc, L_BRING_IN_WHITE);

        /* Shear image to take points 2 and 3 off the axis and
         * put them in the original dest position */
    pixVShearIP(pix2, x1p, -ph2p, L_BRING_IN_WHITE);
    pixHShearIP(pix2, y1p, -th3p, L_BRING_IN_WHITE);

    if (bw != 0 || bh != 0) {
        if ((pixd = pixRemoveBorderGeneral(pix2, bw, bw, bh, bh)) == NULL)
            L_ERROR("pixd not made\n", procName);
    } else {
        pixd = pixClone(pix2);
    }

cleanup_pix:
    pixDestroy(&pix1);
    pixDestroy(&pix2);
    return pixd;
}