/*!
 *  pixRotateAMGrayCorner()
 *
 *      Input:  pixs
 *              angle (radians; clockwise is positive)
 *              grayval (0 to bring in BLACK, 255 for WHITE)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) Rotates the image about the UL corner.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) Specify the grayvalue to be brought in from outside the image.
 */
PIX *
pixRotateAMGrayCorner(PIX       *pixs,
                      l_float32  angle,
                      l_uint8    grayval)
{
l_int32    w, h, wpls, wpld;
l_uint32  *datas, *datad;
PIX       *pixd;

    PROCNAME("pixRotateAMGrayCorner");

    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 (L_ABS(angle) < VERY_SMALL_ANGLE)
        return pixClone(pixs);

    pixGetDimensions(pixs, &w, &h, NULL);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    pixd = pixCreateTemplate(pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    rotateAMGrayCornerLow(datad, w, h, wpld, datas, wpls, angle, grayval);

    return pixd;
}

/*!
 *  pixRotate3Shear()
 *
 *      Input:  pixs
 *              xcen, ycen (center of rotation)
 *              angle (radians)
 *              incolor (L_BRING_IN_WHITE, L_BRING_IN_BLACK);
 *      Return: pixd, or null on error.
 *
 *  Notes:
 *      (1) This rotates the image about the image center,
 *          using the 3-shear method.  It can be used for any angle, and
 *          should be used for angles larger than MAX_2_SHEAR_ANGLE.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) 3-shear rotation by a specified angle is equivalent
 *          to the sequential transformations
 *            y' = y + tan(angle/2) * (x - xcen)     for first y-shear
 *            x' = x + sin(angle) * (y - ycen)       for x-shear
 *            y' = y + tan(angle/2) * (x - xcen)     for second y-shear
 *      (4) Computation of tan(angle) is performed in the shear operations.
 *      (5) This brings in 'incolor' pixels from outside the image.
 *      (6) The algorithm was published by Alan Paeth: "A Fast Algorithm
 *          for General Raster Rotation," Graphics Interface '86,
 *          pp. 77-81, May 1986.  A description of the method, along with
 *          an implementation, can be found in Graphics Gems, p. 179,
 *          edited by Andrew Glassner, published by Academic Press, 1990.
 */
PIX *
pixRotate3Shear(PIX       *pixs,
                l_int32    xcen,
                l_int32    ycen,
                l_float32  angle,
                l_int32    incolor)
{
l_float32  hangle;
PIX              *pixt, *pixd;

    PROCNAME("pixRotate3Shear");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return (PIX *)(PIX *)ERROR_PTR("invalid incolor value", procName, NULL);

    if (L_ABS(angle) < VERY_SMALL_ANGLE)
        return pixClone(pixs);

    hangle = atan(sin(angle));
    if ((pixd = pixVShear(NULL, pixs, xcen, angle / 2., incolor)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    if ((pixt = pixHShear(NULL, pixd, ycen, hangle, incolor)) == NULL)
        return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
    pixVShear(pixd, pixt, xcen, angle / 2., incolor);
    pixDestroy(&pixt);

    return pixd;
}

/*!
 *  kernelNormalize()
 *
 *      Input:  kels (source kel, to be normalized)
 *              normsum (desired sum of elements in keld)
 *      Return: keld (normalized version of kels), or null on error
 *                   or if sum of elements is very close to 0)
 *
 *  Notes:
 *      (1) If the sum of kernel elements is close to 0, do not
 *          try to calculate the normalized kernel.  Instead,
 *          return a copy of the input kernel, with an error message.
 */
L_KERNEL *
kernelNormalize(L_KERNEL  *kels,
                l_float32  normsum)
{
l_int32    i, j, sx, sy, cx, cy;
l_float32  sum, factor;
L_KERNEL  *keld;

    PROCNAME("kernelNormalize");

    if (!kels)
        return (L_KERNEL *)ERROR_PTR("kels not defined", procName, NULL);

    kernelGetSum(kels, &sum);
    if (L_ABS(sum) < 0.01) {
        L_ERROR("null sum; not normalizing; returning a copy", procName);
        return kernelCopy(kels);
    }

    kernelGetParameters(kels, &sy, &sx, &cy, &cx);
    if ((keld = kernelCreate(sy, sx)) == NULL)
        return (L_KERNEL *)ERROR_PTR("keld not made", procName, NULL);
    keld->cy = cy;
    keld->cx = cx;

    factor = normsum / sum;
    for (i = 0; i < sy; i++)
        for (j = 0; j < sx; j++)
            keld->data[i][j] = factor * kels->data[i][j];

    return keld;
}

/*!
 *  pixBilateralGrayExact()
 *
 *      Input:  pixs (8 bpp gray)
 *              spatial_kel  (gaussian kernel)
 *              range_kel (<optional> 256 x 1, monotonically decreasing)
 *      Return: pixd (8 bpp bilateral filtered image)
 *
 *  Notes:
 *      (1) See pixBilateralExact().
 */
PIX *
pixBilateralGrayExact(PIX *pixs,
                      L_KERNEL *spatial_kel,
                      L_KERNEL *range_kel) {
    l_int32 i, j, id, jd, k, m, w, h, d, sx, sy, cx, cy, wplt, wpld;
    l_int32 val, center_val;
    l_uint32 *datat, *datad, *linet, *lined;
    l_float32 sum, weight_sum, weight;
    L_KERNEL *keli;
    PIX *pixt, *pixd;

    PROCNAME("pixBilateralGrayExact");

    if (!pixs)
        return (PIX *) ERROR_PTR("pixs not defined", procName, NULL);
    if (pixGetDepth(pixs) != 8)
        return (PIX *) ERROR_PTR("pixs must be gray", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (!spatial_kel)
        return (PIX *) ERROR_PTR("spatial kel not defined", procName, NULL);

    if (!range_kel)
        return pixConvolve(pixs, spatial_kel, 8, 1);
    if (range_kel->sx != 256 || range_kel->sy != 1)
        return (PIX *) ERROR_PTR("range kel not {256 x 1", procName, NULL);

    keli = kernelInvert(spatial_kel);
    kernelGetParameters(keli, &sy, &sx, &cy, &cx);
    if ((pixt = pixAddMirroredBorder(pixs, cx, sx - cx, cy, sy - cy)) == NULL)
        return (PIX *) ERROR_PTR("pixt not made", procName, NULL);

    pixd = pixCreate(w, h, 8);
    datat = pixGetData(pixt);
    datad = pixGetData(pixd);
    wplt = pixGetWpl(pixt);
    wpld = pixGetWpl(pixd);
    for (i = 0, id = 0; id < h; i++, id++) {
        lined = datad + id * wpld;
        for (j = 0, jd = 0; jd < w; j++, jd++) {
            center_val = GET_DATA_BYTE(datat + (i + cy) * wplt, j + cx);
            weight_sum = 0.0;
            sum = 0.0;
            for (k = 0; k < sy; k++) {
                linet = datat + (i + k) * wplt;
                for (m = 0; m < sx; m++) {
                    val = GET_DATA_BYTE(linet, j + m);
                    weight = keli->data[k][m] *
                             range_kel->data[0][L_ABS(center_val - val)];
                    weight_sum += weight;
                    sum += val * weight;
                }
            }
            SET_DATA_BYTE(lined, jd, (l_int32)(sum / weight_sum + 0.5));
        }
    }

    kernelDestroy(&keli);
    pixDestroy(&pixt);
    return pixd;
}

/*!
 *  pixRotateAMColorFast()
 *
 *      Input:  pixs
 *              angle (radians; clockwise is positive)
 *              colorval (e.g., 0 to bring in BLACK, 0xffffff00 for WHITE)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) This rotates a color image about the image center.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) It uses area mapping, dividing each pixel into
 *          16 subpixels.
 *      (4) It is about 10% to 20% faster than the more accurate linear
 *          interpolation function pixRotateAMColor(),
 *          which uses 256 subpixels.
 *      (5) For some reason it shifts the image center.
 *          No attempt is made to rotate the alpha component.
 *
 *  *** Warning: implicit assumption about RGB component ordering ***
 */
PIX *
pixRotateAMColorFast(PIX       *pixs,
                     l_float32  angle,
                     l_uint32   colorval)
{
l_int32    w, h, wpls, wpld;
l_uint32  *datas, *datad;
PIX       *pixd;

    PROCNAME("pixRotateAMColorFast");

    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);

    if (L_ABS(angle) < MIN_ANGLE_TO_ROTATE)
        return pixClone(pixs);

    pixGetDimensions(pixs, &w, &h, NULL);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    pixd = pixCreateTemplate(pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    rotateAMColorFastLow(datad, w, h, wpld, datas, wpls, angle, colorval);
    return pixd;
}

/*!
 * \brief   regTestCompareValues()
 *
 * \param[in]    rp regtest parameters
 * \param[in]    val1 typ. the golden value
 * \param[in]    val2 typ. the value computed
 * \param[in]    delta allowed max absolute difference
 * \return  0 if OK, 1 on error a failure in comparison is not an error
 */
l_int32
regTestCompareValues(L_REGPARAMS  *rp,
                     l_float32     val1,
                     l_float32     val2,
                     l_float32     delta)
{
l_float32  diff;

    PROCNAME("regTestCompareValues");

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

    rp->index++;
    diff = L_ABS(val2 - val1);

        /* Record on failure */
    if (diff > delta) {
        if (rp->fp) {
            fprintf(rp->fp,
                    "Failure in %s_reg: value comparison for index %d\n"
                    "difference = %f but allowed delta = %f\n",
                    rp->testname, rp->index, diff, delta);
        }
        fprintf(stderr,
                    "Failure in %s_reg: value comparison for index %d\n"
                    "difference = %f but allowed delta = %f\n",
                    rp->testname, rp->index, diff, delta);
        rp->success = FALSE;
    }
    return 0;
}

/*!
 * \brief   pixRotate3Shear()
 *
 * \param[in]    pixs
 * \param[in]    xcen, ycen center of rotation
 * \param[in]    angle radians
 * \param[in]    incolor L_BRING_IN_WHITE, L_BRING_IN_BLACK;
 * \return  pixd, or NULL on error.
 *
 * <pre>
 * Notes:
 *      (1) This rotates the image about the given point, using the 3-shear
 *          method.  It should only be used for angles smaller than
 *          LIMIT_SHEAR_ANGLE.  For larger angles, a warning is issued.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) 3-shear rotation by a specified angle is equivalent
 *          to the sequential transformations
 *            y' = y + tan(angle/2) * (x - xcen)     for first y-shear
 *            x' = x + sin(angle) * (y - ycen)       for x-shear
 *            y' = y + tan(angle/2) * (x - xcen)     for second y-shear
 *      (4) Computation of tan(angle) is performed in the shear operations.
 *      (5) This brings in 'incolor' pixels from outside the image.
 *      (6) If the image has an alpha layer, it is rotated separately by
 *          two shears.
 *      (7) The algorithm was published by Alan Paeth: "A Fast Algorithm
 *          for General Raster Rotation," Graphics Interface '86,
 *          pp. 77-81, May 1986.  A description of the method, along with
 *          an implementation, can be found in Graphics Gems, p. 179,
 *          edited by Andrew Glassner, published by Academic Press, 1990.
 * </pre>
 */
PIX *
pixRotate3Shear(PIX       *pixs,
                l_int32    xcen,
                l_int32    ycen,
                l_float32  angle,
                l_int32    incolor)
{
l_float32  hangle;
PIX       *pix1, *pix2, *pixd;

    PROCNAME("pixRotate3Shear");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return (PIX *)(PIX *)ERROR_PTR("invalid incolor value", procName, NULL);

    if (L_ABS(angle) < MIN_ANGLE_TO_ROTATE)
        return pixClone(pixs);
    if (L_ABS(angle) > LIMIT_SHEAR_ANGLE) {
        L_WARNING("%6.2f radians; large angle for 3-shear rotation\n",
                  procName, L_ABS(angle));
    }

    hangle = atan(sin(angle));
    if ((pixd = pixVShear(NULL, pixs, xcen, angle / 2., incolor)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    if ((pix1 = pixHShear(NULL, pixd, ycen, hangle, incolor)) == NULL) {
        pixDestroy(&pixd);
        return (PIX *)ERROR_PTR("pix1 not made", procName, NULL);
    }
    pixVShear(pixd, pix1, xcen, angle / 2., incolor);
    pixDestroy(&pix1);

    if (pixGetDepth(pixs) == 32 && pixGetSpp(pixs) == 4) {
        pix1 = pixGetRGBComponent(pixs, L_ALPHA_CHANNEL);
            /* L_BRING_IN_WHITE brings in opaque for the alpha component */
        pix2 = pixRotate3Shear(pix1, xcen, ycen, angle, L_BRING_IN_WHITE);
        pixSetRGBComponent(pixd, pix2, L_ALPHA_CHANNEL);
        pixDestroy(&pix1);
        pixDestroy(&pix2);
    }
    return pixd;
}

/*!
 *  pixRasteropHip()
 *
 *      Input:  pixd (in-place operation)
 *              by  (top of horizontal band)
 *              bh  (height of horizontal band)
 *              hshift (horizontal shift of band; hshift > 0 is to right)
 *              incolor (L_BRING_IN_WHITE, L_BRING_IN_BLACK)
 *      Return: 0 if OK; 1 on error
 *
 *  Notes:
 *      (1) This rasterop translates a horizontal band of the
 *          image either left or right, bringing in either white
 *          or black pixels from outside the image.
 *      (2) The horizontal band extends the full width of pixd.
 *      (3) If a colormap exists, the nearest color to white or black
 *          is brought in.
 */
l_int32
pixRasteropHip(PIX     *pixd,
               l_int32  by,
               l_int32  bh,
               l_int32  hshift,
               l_int32  incolor)
{
l_int32   w, h, d, index, op;
PIX      *pixt;
PIXCMAP  *cmap;

    PROCNAME("pixRasteropHip");

    if (!pixd)
        return ERROR_INT("pixd not defined", procName, 1);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return ERROR_INT("invalid value for incolor", procName, 1);
    if (bh <= 0)
        return ERROR_INT("bh must be > 0", procName, 1);

    if (hshift == 0)
        return 0;

    pixGetDimensions(pixd, &w, &h, &d);
    rasteropHipLow(pixGetData(pixd), h, d, pixGetWpl(pixd), by, bh, hshift);

    cmap = pixGetColormap(pixd);
    if (!cmap) {
        if ((d == 1 && incolor == L_BRING_IN_BLACK) ||
            (d > 1 && incolor == L_BRING_IN_WHITE))
            op = PIX_SET;
        else
            op = PIX_CLR;

            /* Set the pixels brought in at left or right */
        if (hshift > 0)
            pixRasterop(pixd, 0, by, hshift, bh, op, NULL, 0, 0);
        else  /* hshift < 0 */
            pixRasterop(pixd, w + hshift, by, -hshift, bh, op, NULL, 0, 0);
        return 0;
    }

        /* Get the nearest index and fill with that */
    if (incolor == L_BRING_IN_BLACK)
        pixcmapGetRankIntensity(cmap, 0.0, &index);
    else  /* white */
        pixcmapGetRankIntensity(cmap, 1.0, &index);
    pixt = pixCreate(L_ABS(hshift), bh, d);
    pixSetAllArbitrary(pixt, index);
    if (hshift > 0)
        pixRasterop(pixd, 0, by, hshift, bh, PIX_SRC, pixt, 0, 0);
    else  /* hshift < 0 */
        pixRasterop(pixd, w + hshift, by, -hshift, bh, PIX_SRC, pixt, 0, 0);
    pixDestroy(&pixt);
    return 0;
}

/*!
 * \brief   makeBarrelshiftString()
 */
static char *
makeBarrelshiftString(l_int32  delx,    /* j - cx */
                      l_int32  dely)    /* i - cy */
{
l_int32  absx, absy;
char     bigbuf[L_BUF_SIZE];

    PROCNAME("makeBarrelshiftString");

    if (delx < -31 || delx > 31)
        return (char *)ERROR_PTR("delx out of bounds", procName, NULL);
    if (dely < -31 || dely > 31)
        return (char *)ERROR_PTR("dely out of bounds", procName, NULL);
    absx = L_ABS(delx);
    absy = L_ABS(dely);

    if ((delx == 0) && (dely == 0))
        sprintf(bigbuf, "(*sptr)");
    else if ((delx == 0) && (dely < 0))
        sprintf(bigbuf, "(*(sptr %s))", wplstrm[absy - 1]);
    else if ((delx == 0) && (dely > 0))
        sprintf(bigbuf, "(*(sptr %s))", wplstrp[absy - 1]);
    else if ((delx < 0) && (dely == 0))
        sprintf(bigbuf, "((*(sptr) >> %d) | (*(sptr - 1) << %d))",
              absx, 32 - absx);
    else if ((delx > 0) && (dely == 0))
        sprintf(bigbuf, "((*(sptr) << %d) | (*(sptr + 1) >> %d))",
              absx, 32 - absx);
    else if ((delx < 0) && (dely < 0))
        sprintf(bigbuf, "((*(sptr %s) >> %d) | (*(sptr %s - 1) << %d))",
              wplstrm[absy - 1], absx, wplstrm[absy - 1], 32 - absx);
    else if ((delx > 0) && (dely < 0))
        sprintf(bigbuf, "((*(sptr %s) << %d) | (*(sptr %s + 1) >> %d))",
              wplstrm[absy - 1], absx, wplstrm[absy - 1], 32 - absx);
    else if ((delx < 0) && (dely > 0))
        sprintf(bigbuf, "((*(sptr %s) >> %d) | (*(sptr %s - 1) << %d))",
              wplstrp[absy - 1], absx, wplstrp[absy - 1], 32 - absx);
    else  /*  ((delx > 0) && (dely > 0))  */
        sprintf(bigbuf, "((*(sptr %s) << %d) | (*(sptr %s + 1) >> %d))",
              wplstrp[absy - 1], absx, wplstrp[absy - 1], 32 - absx);

    return stringNew(bigbuf);
}

/*!
 *  dewarpaInfo()
 *
 *      Input:  fp
 *              dewa
 *      Return: 0 if OK, 1 on error
 */
l_int32
dewarpaInfo(FILE       *fp,
            L_DEWARPA  *dewa)
{
l_int32    i, n, pageno, nnone, nvsuccess, nvvalid, nhsuccess, nhvalid, nref;
L_DEWARP  *dew;

    PROCNAME("dewarpaInfo");

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

    fprintf(fp, "\nDewarpaInfo: %p\n", dewa);
    fprintf(fp, "nalloc = %d, maxpage = %d\n", dewa->nalloc, dewa->maxpage);
    fprintf(fp, "sampling = %d, redfactor = %d, minlines = %d\n",
            dewa->sampling, dewa->redfactor, dewa->minlines);
    fprintf(fp, "maxdist = %d, useboth = %d\n",
            dewa->maxdist, dewa->useboth);

    dewarpaModelStats(dewa, &nnone, &nvsuccess, &nvvalid,
                      &nhsuccess, &nhvalid, &nref);
    n = numaGetCount(dewa->napages);
    fprintf(stderr, "Total number of pages with a dew = %d\n", n);
    fprintf(stderr, "Number of pages without any models = %d\n", nnone);
    fprintf(stderr, "Number of pages with a vert model = %d\n", nvsuccess);
    fprintf(stderr, "Number of pages with a valid vert model = %d\n", nvvalid);
    fprintf(stderr, "Number of pages with both models = %d\n", nhsuccess);
    fprintf(stderr, "Number of pages with both models valid = %d\n", nhvalid);
    fprintf(stderr, "Number of pages with a ref model = %d\n", nref);

    for (i = 0; i < n; i++) {
        numaGetIValue(dewa->napages, i, &pageno);
        if ((dew = dewarpaGetDewarp(dewa, pageno)) == NULL)
            continue;
        fprintf(stderr, "Page: %d\n", dew->pageno);
        fprintf(stderr, "  hasref = %d, refpage = %d\n",
                dew->hasref, dew->refpage);
        fprintf(stderr, "  nlines = %d\n", dew->nlines);
        fprintf(stderr, "  w = %d, h = %d, nx = %d, ny = %d\n",
                dew->w, dew->h, dew->nx, dew->ny);
        if (dew->sampvdispar)
            fprintf(stderr, "  Vertical disparity builds:\n"
                    "    (min,max,abs-diff) line curvature = (%d,%d,%d)\n",
                    dew->mincurv, dew->maxcurv, dew->maxcurv - dew->mincurv);
        if (dew->samphdispar)
            fprintf(stderr, "  Horizontal disparity builds:\n"
                    "    left edge slope = %d, right edge slope = %d\n"
                    "    (left,right,abs-diff) edge curvature = (%d,%d,%d)\n",
                    dew->leftslope, dew->rightslope, dew->leftcurv,
                    dew->rightcurv, L_ABS(dew->leftcurv - dew->rightcurv));
    }
    return 0;
}

/*!
 *  localSearchForBackground()
 *
 *      Input:  &x, &y (starting position for search; return found position)
 *              maxrad (max distance to search from starting location)
 *      Return: 0 if bg pixel found; 1 if not found
 */
static l_int32
localSearchForBackground(PIX  *pix,
                         l_int32  *px,
                         l_int32  *py,
                         l_int32  maxrad)
{
l_int32   x, y, w, h, r, i, j;
l_uint32  val;

    x = *px;
    y = *py;
    pixGetPixel(pix, x, y, &val);
    if (val == 0) return 0;

        /* For each value of r, restrict the search to the boundary
         * pixels in a square centered on (x,y), clipping to the
         * image boundaries if necessary.  */
    pixGetDimensions(pix, &w, &h, NULL);
    for (r = 1; r < maxrad; r++) {
        for (i = -r; i <= r; i++) {
            if (y + i < 0 || y + i >= h)
                continue;
            for (j = -r; j <= r; j++) {
                if (x + j < 0 || x + j >= w)
                    continue;
                if (L_ABS(i) != r && L_ABS(j) != r)  /* not on "r ring" */
                    continue;
                pixGetPixel(pix, x + j, y + i, &val);
                if (val == 0) {
                    *px = x + j;
                    *py = y + i;
                    return 0;
                }
            }
        }
    }
    return 1;
}

/*!
 *  pixRotateShear()
 *
 *      Input:  pixs
 *              xcen (x value for which there is no horizontal shear)
 *              ycen (y value for which there is no vertical shear)
 *              angle (radians)
 *              incolor (L_BRING_IN_WHITE, L_BRING_IN_BLACK);
 *      Return: pixd, or null on error.
 *
 *  Notes:
 *      (1) This rotates an image about the given point, using
 *          either 2 or 3 shears.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) This brings in 'incolor' pixels from outside the image.
 */
PIX *
pixRotateShear(PIX       *pixs,
               l_int32    xcen,
               l_int32    ycen,
               l_float32  angle,
               l_int32    incolor)
{
    PROCNAME("pixRotateShear");

    if (!pixs)
        return (PIX *)(PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    if (incolor != L_BRING_IN_WHITE && incolor != L_BRING_IN_BLACK)
        return (PIX *)(PIX *)ERROR_PTR("invalid incolor value", procName, NULL);

    if (L_ABS(angle) < VERY_SMALL_ANGLE)
        return pixClone(pixs);

    if (L_ABS(angle) <= MAX_2_SHEAR_ANGLE)
        return pixRotate2Shear(pixs, xcen, ycen, angle, incolor);
    else
        return pixRotate3Shear(pixs, xcen, ycen, angle, incolor);

}

/**
 * Test whether b1 and b2 are close enough to be a character pair.
 */
bool ValidatePairOld(BOX *b1, BOX *b2) {
  l_int32 max_w = L_MAX(b1->w, b2->w);
  l_int32 centerx1 = b1->x + b1->w / 2;
  l_int32 centerx2 = b2->x + b2->w / 2;
  l_int32 h_dist = L_ABS(centerx1 - centerx2);

  /* Horizontal distance between centers is
   * less than twice the wider character */
  if (h_dist > max_w * OLDPAIR_MAX_HDIST_RATIO)
    return false;

  l_int32 max_h = L_MAX(b1->h, b2->h);
  l_int32 centery1 = b1->y + b1->h / 2;
  l_int32 centery2 = b2->y + b2->h / 2;
  l_int32 v_dist = L_ABS(centery1 - centery2);

  /* Vertical distance between centers is
   less than 50% of the taller character */
  if (v_dist > max_h * OLDPAIR_MAX_VDIST_RATIO)
    return false;

  l_int32 min_h = L_MIN(b1->h, b2->h);
  l_float32 h_ratio = min_h / (max_h + 1.0);

  /* Height ratio is between 0.5 and 2 */
  if (h_ratio < OLDPAIR_MIN_HPAIR_RATIO)
    return false;

  l_int32 min_w = L_MIN(b1->w, b2->w);
  l_float32 w_ratio = min_w / (max_w + 1.0);

  /* Width ratio is between 0.1 and 10 */
  if (w_ratio < OLDPAIR_MIN_WPAIR_RATIO)
    return false;

  return true;
}

/*!
 *  pixRotateAM()
 *
 *      Input:  pixs (2, 4, 8 bpp gray or colormapped, or 32 bpp RGB)
 *              angle (radians; clockwise is positive)
 *              incolor (L_BRING_IN_WHITE, L_BRING_IN_BLACK)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) Rotates about image center.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) Brings in either black or white pixels from the boundary.
 */
PIX *
pixRotateAM(PIX       *pixs,
            l_float32  angle,
            l_int32    incolor)
{
l_int32   d;
l_uint32  fillval;
PIX      *pixt1, *pixt2, *pixd;

    PROCNAME("pixRotateAM");

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

    if (L_ABS(angle) < MIN_ANGLE_TO_ROTATE)
        return pixClone(pixs);

        /* 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 incoming color */
    fillval = 0;
    if (incolor == L_BRING_IN_WHITE) {
        if (d == 8)
            fillval = 255;
        else  /* d == 32 */
            fillval = 0xffffff00;
    }

    if (d == 8)
        pixd = pixRotateAMGray(pixt2, angle, fillval);
    else   /* d == 32 */
        pixd = pixRotateAMColor(pixt2, angle, fillval);

    pixDestroy(&pixt1);
    pixDestroy(&pixt2);
    return pixd;
}

/*!
 *  pixcmapColorToGray()
 *
 *      Input:  cmap
 *              rwt, gwt, bwt  (non-negative; these should add to 1.0)
 *      Return: cmap (gray), or null on error
 *
 *  Notes:
 *      (1) This creates a gray colormap from an arbitrary colormap.
 *      (2) In use, attach the output gray colormap to the pix
 *          (or a copy of it) that provided the input colormap.
 */
PIXCMAP *
pixcmapColorToGray(PIXCMAP   *cmaps,
                   l_float32  rwt,
                   l_float32  gwt,
                   l_float32  bwt)
{
l_int32    i, n, rval, gval, bval, val;
l_float32  sum;
PIXCMAP   *cmapd;

    PROCNAME("pixcmapColorToGray");

    if (!cmaps)
        return (PIXCMAP *)ERROR_PTR("cmaps not defined", procName, NULL);
    if (rwt < 0.0 || gwt < 0.0 || bwt < 0.0)
        return (PIXCMAP *)ERROR_PTR("weights not all >= 0.0", procName, NULL);

        /* Make sure the sum of weights is 1.0; otherwise, you can get
         * overflow in the gray value. */
    sum = rwt + gwt + bwt;
    if (sum == 0.0) {
        L_WARNING("all weights zero; setting equal to 1/3", procName);
        rwt = gwt = bwt = 0.33333;
        sum = 1.0;
    }
    if (L_ABS(sum - 1.0) > 0.0001) {  /* maintain ratios with sum == 1.0 */
        L_WARNING("weights don't sum to 1; maintaining ratios", procName);
        rwt = rwt / sum;
        gwt = gwt / sum;
        bwt = bwt / sum;
    }

    cmapd = pixcmapCopy(cmaps);
    n = pixcmapGetCount(cmapd);
    for (i = 0; i < n; i++) {
        pixcmapGetColor(cmapd, i, &rval, &gval, &bval);
        val = (l_int32)(rwt * rval + gwt * gval + bwt * bval + 0.5);
        pixcmapResetColor(cmapd, i, val, val, val);
    }
        
    return cmapd;
}

/*!
 *  pixRotateAMColorCorner()
 *
 *      Input:  pixs
 *              angle (radians; clockwise is positive)
 *              colorval (e.g., 0 to bring in BLACK, 0xffffff00 for WHITE)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) Rotates the image about the UL corner.
 *      (2) A positive angle gives a clockwise rotation.
 *      (3) Specify the color to be brought in from outside the image.
 */
PIX *
pixRotateAMColorCorner(PIX       *pixs,
                       l_float32  angle,
                       l_uint32   fillval)
{
l_int32    w, h, wpls, wpld;
l_uint32  *datas, *datad;
PIX       *pix1, *pix2, *pixd;

    PROCNAME("pixRotateAMColorCorner");

    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);

    if (L_ABS(angle) < MIN_ANGLE_TO_ROTATE)
        return pixClone(pixs);

    pixGetDimensions(pixs, &w, &h, NULL);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    pixd = pixCreateTemplate(pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);

    rotateAMColorCornerLow(datad, w, h, wpld, datas, wpls, angle, fillval);
    if (pixGetSpp(pixs) == 4) {
        pix1 = pixGetRGBComponent(pixs, L_ALPHA_CHANNEL);
        pix2 = pixRotateAMGrayCorner(pix1, angle, 255);  /* bring in opaque */
        pixSetRGBComponent(pixd, pix2, L_ALPHA_CHANNEL);
        pixDestroy(&pix1);
        pixDestroy(&pix2);
    }

    return pixd;
}

/*!
 *  pixcmapGetNearestGrayIndex()
 *
 *      Input:  cmap
 *              val (gray value to search for; in range [0, ... 255])
 *              &index (<return> the index of the nearest color)
 *      Return: 0 if OK, 1 on error (caller must check)
 *
 *  Notes:
 *      (1) This should be used on gray colormaps.  It uses only the
 *          green value of the colormap.
 *      (2) Returns the index of the exact color if possible, otherwise the
 *          index of the color closest to the target color.
 */
l_int32
pixcmapGetNearestGrayIndex(PIXCMAP  *cmap,
                           l_int32   val,
                           l_int32  *pindex)
{
l_int32     i, n, dist, mindist;
RGBA_QUAD  *cta;

    PROCNAME("pixcmapGetNearestGrayIndex");

    if (!pindex)
        return ERROR_INT("&index not defined", procName, 1);
    *pindex = 0;
    if (!cmap)
        return ERROR_INT("cmap not defined", procName, 1);
    if (val < 0 || val > 255)
        return ERROR_INT("val not in [0 ... 255]", procName, 1);

    if ((cta = (RGBA_QUAD *)cmap->array) == NULL)
        return ERROR_INT("cta not defined(!)", procName, 1);
    n = pixcmapGetCount(cmap);

    mindist = 256;
    for (i = 0; i < n; i++) {
        dist = cta[i].green - val;
        dist = L_ABS(dist);
        if (dist < mindist) {
            *pindex = i;
            if (dist == 0)
                break;
            mindist = dist;
        }
    }

    return 0;
}

//.........这里部分代码省略.........
    if (!naa1)
        return ERROR_INT("naa1 not defined", procName, 1);
    if (!naa2)
        return ERROR_INT("naa2 not defined", procName, 1);
    if (nperline < 1)
        return ERROR_INT("nperline < 1", procName, 1);
    if (nreq < 1)
        return ERROR_INT("nreq < 1", procName, 1);

    n1 = numaaGetCount(naa1);
    n2 = numaaGetCount(naa2);
    if (n1 < nreq || n2 < nreq)
        return 0;

        /* Find the lines in naa1 and naa2 with sufficient boxes.
         * Also, find the y-values for each of the lines, and the
         * LH x-values of the first box in each line. */
    line1 = (l_int32 *)LEPT_CALLOC(n1, sizeof(l_int32));
    line2 = (l_int32 *)LEPT_CALLOC(n2, sizeof(l_int32));
    yloc1 = (l_int32 *)LEPT_CALLOC(n1, sizeof(l_int32));
    yloc2 = (l_int32 *)LEPT_CALLOC(n2, sizeof(l_int32));
    xleft1 = (l_int32 *)LEPT_CALLOC(n1, sizeof(l_int32));
    xleft2 = (l_int32 *)LEPT_CALLOC(n2, sizeof(l_int32));
    for (i = 0; i < n1; i++) {
        na1 = numaaGetNuma(naa1, i, L_CLONE);
        numaGetIValue(na1, 0, yloc1 + i);
        numaGetIValue(na1, 1, xleft1 + i);
        nbox = (numaGetCount(na1) - 1) / 2;
        if (nbox >= nperline)
            line1[i] = 1;
        numaDestroy(&na1);
    }
    for (i = 0; i < n2; i++) {
        na2 = numaaGetNuma(naa2, i, L_CLONE);
        numaGetIValue(na2, 0, yloc2 + i);
        numaGetIValue(na2, 1, xleft2 + i);
        nbox = (numaGetCount(na2) - 1) / 2;
        if (nbox >= nperline)
            line2[i] = 1;
        numaDestroy(&na2);
    }

        /* Enumerate all possible line matches.  A 'possible' line
         * match is one where the x and y shifts for the first box
         * in each line are within the maxshiftx and maxshifty
         * constraints, and the left and right sides of the remaining
         * (nperline - 1) successive boxes are within delx of each other.
         * The result is a set of four numas giving parameters of
         * each set of matching lines. */
    nai1 = numaCreate(0);  /* line index 1 of match */
    nai2 = numaCreate(0);  /* line index 2 of match */
    nasx = numaCreate(0);  /* shiftx for match */
    nasy = numaCreate(0);  /* shifty for match */
    for (i = 0; i < n1; i++) {
        if (line1[i] == 0) continue;
        y1 = yloc1[i];
        xl1 = xleft1[i];
        na1 = numaaGetNuma(naa1, i, L_CLONE);
        for (j = 0; j < n2; j++) {
            if (line2[j] == 0) continue;
            y2 = yloc2[j];
            if (L_ABS(y1 - y2) > maxshifty) continue;
            xl2 = xleft2[j];
            if (L_ABS(xl1 - xl2) > maxshiftx) continue;
            shiftx = xl1 - xl2;  /* shift to add to x2 values */
            shifty = y1 - y2;  /* shift to add to y2 values */
            na2 = numaaGetNuma(naa2, j, L_CLONE);

                /* Now check if 'nperline' boxes in the two lines match */
            match = testLineAlignmentX(na1, na2, shiftx, delx, nperline);
            if (match) {
                numaAddNumber(nai1, i);
                numaAddNumber(nai2, j);
                numaAddNumber(nasx, shiftx);
                numaAddNumber(nasy, shifty);
            }
            numaDestroy(&na2);
        }
        numaDestroy(&na1);
    }

        /* Determine if there are a sufficient number of mutually
         * aligned matches.  Mutually aligned matches place an additional
         * constraint on the 'possible' matches, where the relative
         * shifts must not exceed the (delx, dely) distances. */
    countAlignedMatches(nai1, nai2, nasx, nasy, n1, n2, delx, dely,
                        nreq, psame, debugflag);

    LEPT_FREE(line1);
    LEPT_FREE(line2);
    LEPT_FREE(yloc1);
    LEPT_FREE(yloc2);
    LEPT_FREE(xleft1);
    LEPT_FREE(xleft2);
    numaDestroy(&nai1);
    numaDestroy(&nai2);
    numaDestroy(&nasx);
    numaDestroy(&nasy);
    return 0;
}

l_float32 RelativeDiff(l_int32 v1, l_int32 v2) {
  return L_ABS(v1 - v2) / (L_MIN(v1, v2) + 1.0);
}

l_float32 ComputePairNormalizedHorizontalDistance(BOX *b1, BOX *b2) {
  l_float32 dx = (b1->x - b2->x) + (b1->w - b2->w) / 2.0;
  l_float32 hdist = 2.0 * L_ABS(dx) / (b1->w + b2->w);

  return hdist;
}