// Compute the distance between (x, y1) and (x, y2) using the rule that
// a decrease in textline density is weighted more heavily than an increase.
// The coordinates are in source image space, ie processed by any denorm
// already, but not yet scaled by scale_factor_.
// Going from the outside of a textline to the inside should measure much
// less distance than going from the inside of a textline to the outside.
// How it works:
// An increase is cheap (getting closer to a textline).
// Constant costs unity.
// A decrease is expensive (getting further from a textline).
// Pixels in projection map Counted distance
//              2
//              3              1/x
//              3               1
//              2               x
//              5              1/x
//              7              1/x
// Total: 1 + x + 3/x where x = kWrongWayPenalty.
int TextlineProjection::VerticalDistance(bool debug, int x,
                                         int y1, int y2) const {
  x = ImageXToProjectionX(x);
  y1 = ImageYToProjectionY(y1);
  y2 = ImageYToProjectionY(y2);
  if (y1 == y2) return 0;
  int wpl = pixGetWpl(pix_);
  int step = y1 < y2 ? 1 : -1;
  uint32_t* data = pixGetData(pix_) + y1 * wpl;
  wpl *= step;
  int prev_pixel = GET_DATA_BYTE(data, x);
  int distance = 0;
  int right_way_steps = 0;
  for (int y = y1; y != y2; y += step) {
    data += wpl;
    int pixel = GET_DATA_BYTE(data, x);
    if (debug)
      tprintf("At (%d,%d), pix = %d, prev=%d\n",
              x, y + step, pixel, prev_pixel);
    if (pixel < prev_pixel)
      distance += kWrongWayPenalty;
    else if (pixel > prev_pixel)
      ++right_way_steps;
    else
      ++distance;
    prev_pixel = pixel;
  }
  return distance * scale_factor_ +
      right_way_steps * scale_factor_ / kWrongWayPenalty;
}

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

/*!
 *  pixAddConstantGray()
 *
 *      Input:  pixs (8, 16 or 32 bpp)
 *              val  (amount to add to each pixel)
 *      Return: 0 if OK, 1 on error
 *
 *  Notes:
 *      (1) In-place operation.
 *      (2) No clipping for 32 bpp.
 *      (3) For 8 and 16 bpp, if val > 0 the result is clipped
 *          to 0xff and 0xffff, rsp.
 *      (4) For 8 and 16 bpp, if val < 0 the result is clipped to 0.
 */
l_int32
pixAddConstantGray(PIX      *pixs,
                   l_int32   val)
{
l_int32    i, j, w, h, d, wpl, pval;
l_uint32  *data, *line;

    PROCNAME("pixAddConstantGray");

    if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8 && d != 16 && d != 32)
        return ERROR_INT("pixs not 8, 16 or 32 bpp", procName, 1);

    data = pixGetData(pixs);
    wpl = pixGetWpl(pixs);
    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        if (d == 8) {
            if (val < 0) {
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_BYTE(line, j);
                    pval = L_MAX(0, pval + val);
                    SET_DATA_BYTE(line, j, pval);
                }
            } else {  /* val >= 0 */
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_BYTE(line, j);
                    pval = L_MIN(255, pval + val);
                    SET_DATA_BYTE(line, j, pval);
                }
            }
        } else if (d == 16) {
            if (val < 0) {
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_TWO_BYTES(line, j);
                    pval = L_MAX(0, pval + val);
                    SET_DATA_TWO_BYTES(line, j, pval);
                }
            } else {  /* val >= 0 */
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_TWO_BYTES(line, j);
                    pval = L_MIN(0xffff, pval + val);
                    SET_DATA_TWO_BYTES(line, j, pval);
                }
            }
        } else {  /* d == 32; no check for overflow (< 0 or > 0xffffffff) */
            for (j = 0; j < w; j++)
                *(line + j) += val;
        }
    }

    return 0;
}

/*!
 *  bilateralApply()
 *
 *      Input:  bil
 *      Return: pixd
 */
static PIX *
bilateralApply(L_BILATERAL  *bil)
{
l_int32      i, j, k, ired, jred, w, h, wpls, wpld, ncomps, reduction;
l_int32      vals, vald, lowval, hival;
l_int32     *kindex;
l_float32    fract;
l_float32   *kfract;
l_uint32    *lines, *lined, *datas, *datad;
l_uint32  ***lineset = NULL;  /* for set of PBC */
PIX         *pixs, *pixd;
PIXA        *pixac;

    PROCNAME("bilateralApply");

    if (!bil)
        return (PIX *)ERROR_PTR("bil not defined", procName, NULL);
    pixs = bil->pixs;
    ncomps = bil->ncomps;
    kindex = bil->kindex;
    kfract = bil->kfract;
    reduction = bil->reduction;
    pixac = bil->pixac;
    lineset = bil->lineset;
    if (pixaGetCount(pixac) != ncomps)
        return (PIX *)ERROR_PTR("PBC images do not exist", procName, NULL);

    if ((pixd = pixCreateTemplate(pixs)) == NULL)
        return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    datad = pixGetData(pixd);
    wpld = pixGetWpl(pixd);
    pixGetDimensions(pixs, &w, &h, NULL);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        lined = datad + i * wpld;
        ired = i / reduction;
        for (j = 0; j < w; j++) {
            jred = j / reduction;
            vals = GET_DATA_BYTE(lines, j);
            k = kindex[vals];
            lowval = GET_DATA_BYTE(lineset[k][ired], jred);
            hival = GET_DATA_BYTE(lineset[k + 1][ired], jred);
            fract = kfract[vals];
            vald = (l_int32)((1.0 - fract) * lowval + fract * hival + 0.5);
            SET_DATA_BYTE(lined, j, vald);
        }
    }

    return pixd;
}

/*!
 *  addConstantGrayLow()
 */
void
addConstantGrayLow(l_uint32  *data,
                   l_int32    w,
                   l_int32    h,
                   l_int32    d,
                   l_int32    wpl,
                   l_int32    val)
{
    l_int32    i, j, pval;
    l_uint32  *line;

    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        if (d == 8) {
            if (val < 0) {
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_BYTE(line, j);
                    pval = L_MAX(0, pval + val);
                    SET_DATA_BYTE(line, j, pval);
                }
            }
            else {  /* val >= 0 */
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_BYTE(line, j);
                    pval = L_MIN(255, pval + val);
                    SET_DATA_BYTE(line, j, pval);
                }
            }
        }
        else if (d == 16) {
            if (val < 0) {
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_TWO_BYTES(line, j);
                    pval = L_MAX(0, pval + val);
                    SET_DATA_TWO_BYTES(line, j, pval);
                }
            }
            else {  /* val >= 0 */
                for (j = 0; j < w; j++) {
                    pval = GET_DATA_TWO_BYTES(line, j);
                    pval = L_MIN(0xffff, pval + val);
                    SET_DATA_TWO_BYTES(line, j, pval);
                }
            }
        }
        else {  /* d == 32; no check for overflow (< 0 or > 0xffffffff) */
            for (j = 0; j < w; j++)
                *(line + j) += val;
        }
    }
    return;
}

// Helper returns the mean pixel value over the line between the start_pt and
// end_pt (inclusive), but shifted perpendicular to the line in the projection
// image by offset pixels. For simplicity, it is assumed that the vector is
// either nearly horizontal or nearly vertical. It works on skewed textlines!
// The end points are in external coordinates, and will be denormalized with
// the denorm if not NULL before further conversion to pix coordinates.
// After all the conversions, the offset is added to the direction
// perpendicular to the line direction. The offset is thus in projection image
// coordinates, which allows the caller to get a guaranteed displacement
// between pixels used to calculate gradients.
int TextlineProjection::MeanPixelsInLineSegment(const DENORM* denorm,
                                                int offset,
                                                TPOINT start_pt,
                                                TPOINT end_pt) const {
  TransformToPixCoords(denorm, &start_pt);
  TransformToPixCoords(denorm, &end_pt);
  TruncateToImageBounds(&start_pt);
  TruncateToImageBounds(&end_pt);
  int wpl = pixGetWpl(pix_);
  uint32_t* data = pixGetData(pix_);
  int total = 0;
  int count = 0;
  int x_delta = end_pt.x - start_pt.x;
  int y_delta = end_pt.y - start_pt.y;
  if (abs(x_delta) >= abs(y_delta)) {
    if (x_delta == 0)
      return 0;
    // Horizontal line. Add the offset vertically.
    int x_step = x_delta > 0 ? 1 : -1;
    // Correct offset for rotation, keeping it anti-clockwise of the delta.
    offset *= x_step;
    start_pt.y += offset;
    end_pt.y += offset;
    TruncateToImageBounds(&start_pt);
    TruncateToImageBounds(&end_pt);
    x_delta = end_pt.x - start_pt.x;
    y_delta = end_pt.y - start_pt.y;
    count = x_delta * x_step + 1;
    for (int x = start_pt.x; x != end_pt.x; x += x_step) {
      int y = start_pt.y + DivRounded(y_delta * (x - start_pt.x), x_delta);
      total += GET_DATA_BYTE(data + wpl * y, x);
    }
  } else {
    // Vertical line. Add the offset horizontally.
    int y_step = y_delta > 0 ? 1 : -1;
    // Correct offset for rotation, keeping it anti-clockwise of the delta.
    // Pix holds the image with y=0 at the top, so the offset is negated.
    offset *= -y_step;
    start_pt.x += offset;
    end_pt.x += offset;
    TruncateToImageBounds(&start_pt);
    TruncateToImageBounds(&end_pt);
    x_delta = end_pt.x - start_pt.x;
    y_delta = end_pt.y - start_pt.y;
    count = y_delta * y_step + 1;
    for (int y = start_pt.y; y != end_pt.y; y += y_step) {
      int x = start_pt.x + DivRounded(x_delta * (y - start_pt.y), y_delta);
      total += GET_DATA_BYTE(data + wpl * y, x);
    }
  }
  return DivRounded(total, count);
}

/*!
 *  pixGlobalStats()
 *
 *      Input:  pixs   (8 bpp grayscale)
 *              &mean  (<optional return> pixs mean)
 *              &var   (<optional return> pixs variance)
 *              &std   (<optional return> pixs standard deviation)
 *      Return: 0 if OK; 1 on error
 */
l_int32
pixGlobalStats(PIX       *pixs,
			   l_float32 *mean,
			   l_float32 *var,
			   l_float32 *std)
{
	l_int32    w, h, d, i, j;
	l_int32    wpl;
	l_uint32  *data, *line;
	l_float32  m, v;
	
	PROCNAME("pixGlobalStats");
	
	if (!mean && !var && !std)
        return ERROR_INT("nothing to do", procName, 1);
	if (!pixs)
        return ERROR_INT("pixs not defined", procName, 1);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return ERROR_INT("pixs not 8 bpp", procName, 1);
	
	wpl = pixGetWpl(pixs);
	data = pixGetData(pixs);
	
	/* Calculate the global mean */
	m = 0.;
	for (i = 0; i < h; i++) {
		line = data + i * wpl;
		for (j = 0; j < w; j++)
			m += GET_DATA_BYTE(line, j);
	}
	m /= (w * h);
	
	/* Calculate the global variance */
	v = 0.;
	for (i = 0; i < h; i++) {
		line = data + i * wpl;
		for (j = 0; j < w; j++)
			v += pow((GET_DATA_BYTE(line, j) - m), 2);
	}
	v /= (w * h);
	
	if (mean)
		*mean = m;
	if (var)
		*var = v;
	if (std)
		*std = sqrt(v);
	
	return 0;
}

/*!
 *  pixDilateGray3h()
 *
 *      Input:  pixs (8 bpp, not cmapped)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) Special case for horizontal 3x1 brick Sel;
 *          also used as the first step for the 3x3 brick Sel.
 */
static PIX *
pixDilateGray3h(PIX  *pixs)
{
l_uint32  *datas, *datad, *lines, *lined;
l_int32    w, h, wpl, i, j;
l_int32    val0, val1, val2, val3, val4, val5, val6, val7, val8, val9, maxval;
PIX       *pixd;

    PROCNAME("pixDilateGray3h");

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

    pixd = pixCreateTemplateNoInit(pixs);
    pixSetBorderVal(pixd, 4, 8, 2, 8, 0);  /* only to silence valgrind */
    pixGetDimensions(pixs, &w, &h, NULL);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpl = pixGetWpl(pixs);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpl;
        lined = datad + i * wpl;
        for (j = 1; j < w - 8; j += 8) {
            val0 = GET_DATA_BYTE(lines, j - 1);
            val1 = GET_DATA_BYTE(lines, j);
            val2 = GET_DATA_BYTE(lines, j + 1);
            val3 = GET_DATA_BYTE(lines, j + 2);
            val4 = GET_DATA_BYTE(lines, j + 3);
            val5 = GET_DATA_BYTE(lines, j + 4);
            val6 = GET_DATA_BYTE(lines, j + 5);
            val7 = GET_DATA_BYTE(lines, j + 6);
            val8 = GET_DATA_BYTE(lines, j + 7);
            val9 = GET_DATA_BYTE(lines, j + 8);
            maxval = L_MAX(val1, val2);
            SET_DATA_BYTE(lined, j, L_MAX(val0, maxval));
            SET_DATA_BYTE(lined, j + 1, L_MAX(maxval, val3));
            maxval = L_MAX(val3, val4);
            SET_DATA_BYTE(lined, j + 2, L_MAX(val2, maxval));
            SET_DATA_BYTE(lined, j + 3, L_MAX(maxval, val5));
            maxval = L_MAX(val5, val6);
            SET_DATA_BYTE(lined, j + 4, L_MAX(val4, maxval));
            SET_DATA_BYTE(lined, j + 5, L_MAX(maxval, val7));
            maxval = L_MAX(val7, val8);
            SET_DATA_BYTE(lined, j + 6, L_MAX(val6, maxval));
            SET_DATA_BYTE(lined, j + 7, L_MAX(maxval, val9));
        }
    }
    return pixd;
}

/*!
 *  pixDilateGray3v()
 *
 *      Input:  pixs (8 bpp, not cmapped)
 *      Return: pixd, or null on error
 *
 *  Notes:
 *      (1) Special case for vertical 1x3 brick Sel;
 *          also used as the second step for the 3x3 brick Sel.
 */
static PIX *
pixDilateGray3v(PIX  *pixs)
{
l_uint32  *datas, *datad, *linesi, *linedi;
l_int32    w, h, wpl, i, j;
l_int32    val0, val1, val2, val3, val4, val5, val6, val7, val8, val9, maxval;
PIX       *pixd;

    PROCNAME("pixDilateGray3v");

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

    pixd = pixCreateTemplateNoInit(pixs);
    pixGetDimensions(pixs, &w, &h, NULL);
    datas = pixGetData(pixs);
    datad = pixGetData(pixd);
    wpl = pixGetWpl(pixs);
    for (j = 0; j < w; j++) {
        for (i = 1; i < h - 8; i += 8) {
            linesi = datas + i * wpl;
            linedi = datad + i * wpl;
            val0 = GET_DATA_BYTE(linesi - wpl, j);
            val1 = GET_DATA_BYTE(linesi, j);
            val2 = GET_DATA_BYTE(linesi + wpl, j);
            val3 = GET_DATA_BYTE(linesi + 2 * wpl, j);
            val4 = GET_DATA_BYTE(linesi + 3 * wpl, j);
            val5 = GET_DATA_BYTE(linesi + 4 * wpl, j);
            val6 = GET_DATA_BYTE(linesi + 5 * wpl, j);
            val7 = GET_DATA_BYTE(linesi + 6 * wpl, j);
            val8 = GET_DATA_BYTE(linesi + 7 * wpl, j);
            val9 = GET_DATA_BYTE(linesi + 8 * wpl, j);
            maxval = L_MAX(val1, val2);
            SET_DATA_BYTE(linedi, j, L_MAX(val0, maxval));
            SET_DATA_BYTE(linedi + wpl, j, L_MAX(maxval, val3));
            maxval = L_MAX(val3, val4);
            SET_DATA_BYTE(linedi + 2 * wpl, j, L_MAX(val2, maxval));
            SET_DATA_BYTE(linedi + 3 * wpl, j, L_MAX(maxval, val5));
            maxval = L_MAX(val5, val6);
            SET_DATA_BYTE(linedi + 4 * wpl, j, L_MAX(val4, maxval));
            SET_DATA_BYTE(linedi + 5 * wpl, j, L_MAX(maxval, val7));
            maxval = L_MAX(val7, val8);
            SET_DATA_BYTE(linedi + 6 * wpl, j, L_MAX(val6, maxval));
            SET_DATA_BYTE(linedi + 7 * wpl, j, L_MAX(maxval, val9));
        }
    }
    return pixd;
}

/*!
 *  dpixMeanSquareAccum()
 *
 *      Input:  pixs (1 bpp or 8 bpp grayscale)
 *      Return: dpix (64 bit array), or null on error
 *
 *  Notes:
 *      (1) This is an extension to the standard pixMeanSquareAccum()
 *          implementation provided by Leptonica, to handle 1bpp binary pix
 *          transparently.
 *      (1) Similar to pixBlockconvAccum(), this computes the
 *          sum of the squares of the pixel values in such a way
 *          that the value at (i,j) is the sum of all squares in
 *          the rectangle from the origin to (i,j).
 *      (2) The general recursion relation (v are squared pixel values) is
 *            a(i,j) = v(i,j) + a(i-1, j) + a(i, j-1) - a(i-1, j-1)
 *          For the first line, this reduces to the special case
 *            a(i,j) = v(i,j) + a(i, j-1)
 *          For the first column, the special case is
 *            a(i,j) = v(i,j) + a(i-1, j)
 */
DPIX *
dpixMeanSquareAccum(PIX  *pixs)
{
	l_int32     i, j, w, h, d, wpl, wpls, val;
	l_uint32   *datas, *lines;
	l_float64  *data, *line, *linep;
	DPIX       *dpix;
	
    PROCNAME("dpixMeanSquareAccum");
	
    if (!pixs)
        return (DPIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 1 && d != 8)
        return (DPIX *)ERROR_PTR("pixs not 1 bpp or 8 bpp", procName, NULL);
    if ((dpix = dpixCreate(w, h)) ==  NULL)
        return (DPIX *)ERROR_PTR("dpix not made", procName, NULL);
	
    datas = pixGetData(pixs);
    wpls = pixGetWpl(pixs);
    data = dpixGetData(dpix);
    wpl = dpixGetWpl(dpix);
	
    lines = datas;
    line = data;
    for (j = 0; j < w; j++) {   /* first line */
        val = d == 1 ? GET_DATA_BIT(lines, j) : GET_DATA_BYTE(lines, j);
        if (j == 0)
            line[0] = val * val;
        else
            line[j] = line[j - 1] + val * val;
    }
	
	/* Do the other lines */
    for (i = 1; i < h; i++) {
        lines = datas + i * wpls;
        line = data + i * wpl;  /* current dest line */
        linep = line - wpl;;  /* prev dest line */
        for (j = 0; j < w; j++) {
            val = d == 1 ? GET_DATA_BIT(lines, j) : GET_DATA_BYTE(lines, j);
            if (j == 0)
                line[0] = linep[0] + val * val;
            else
                line[j] = line[j - 1] + linep[j] - linep[j - 1] + val * val;
        }
    }
	
    return dpix;
}

// Helper computes the local 2-D gradient (dx, dy) from the 2x2 cell centered
// on the given (x,y). If the cell would go outside the image, it is padded
// with white.
static void ComputeGradient(const l_uint32* data, int wpl,
                            int x, int y, int width, int height,
                            ICOORD* gradient) {
  const l_uint32* line = data + y * wpl;
  int pix_x_y = x < width && y < height ?
      GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<const void *>(line)), x) : 255;
  int pix_x_prevy = x < width && y > 0 ?
      GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<const void *>(line - wpl)), x) : 255;
  int pix_prevx_prevy = x > 0 && y > 0 ?
      GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<void const*>(line - wpl)), x - 1) : 255;
  int pix_prevx_y = x > 0 && y < height ?
      GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<const void *>(line)), x - 1) : 255;
  gradient->set_x(pix_x_y + pix_x_prevy - (pix_prevx_y + pix_prevx_prevy));
  gradient->set_y(pix_x_prevy + pix_prevx_prevy - (pix_x_y + pix_prevx_y));
}

// Helper evaluates a horizontal difference, (x,y) - (x-1,y), where y is implied
// by the input image line, returning true if the difference matches diff_sign
// and updating the best_diff, best_sum, best_x if a new max.
static bool EvaluateHorizontalDiff(const l_uint32* line, int diff_sign,
                                   int x, int width,
                                   int* best_diff, int* best_sum, int* best_x) {
  if (x <= 0 || x >= width)
    return false;
  int pixel1 = GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<const void *>(line)), x - 1);
  int pixel2 = GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<const void *>(line)), x);
  int diff = (pixel2 - pixel1) * diff_sign;
  if (diff > *best_diff) {
    *best_diff = diff;
    *best_sum = pixel1 + pixel2;
    *best_x = x;
  }
  return diff > 0;
}

struct corners* RunFastDetector9(PIX *pix,
				 unsigned int w, 
				 unsigned int h) {

  xy *rawcorners = (xy *) malloc(sizeof(xy));
  unsigned char *im = (unsigned char*) malloc(sizeof(unsigned char) * (w*h));
  void **pix_lines = pixGetLinePtrs(pix, NULL);
  unsigned int x,y, num_corners;
  unsigned int i = 0;
  for(y=0;y<h;y++) 
    for(x=0;x<w;x++){
      im[i] = (unsigned char) GET_DATA_BYTE(pix_lines[y],x);
      i++;
    }

  free(pix_lines);
  rawcorners = fast9_detect(im, w, h, w, 88, &num_corners);
  free(im);

  unsigned int mx,my;
  mx = (unsigned int)w/2;
  my = (unsigned int)h/2;
  
  float skew_angle = 0.0;

  struct corners *corners = ParseRawCorners(rawcorners, 
					    num_corners,
					    mx,my,
					    skew_angle);
  free(rawcorners);
  return corners;
}

// Create a window and display the projection in it.
void TextlineProjection::DisplayProjection() const {
#ifndef GRAPHICS_DISABLED
  int width = pixGetWidth(pix_);
  int height = pixGetHeight(pix_);
  Pix* pixc = pixCreate(width, height, 32);
  int src_wpl = pixGetWpl(pix_);
  int col_wpl = pixGetWpl(pixc);
  uint32_t* src_data = pixGetData(pix_);
  uint32_t* col_data = pixGetData(pixc);
  for (int y = 0; y < height; ++y, src_data += src_wpl, col_data += col_wpl) {
    for (int x = 0; x < width; ++x) {
      int pixel = GET_DATA_BYTE(src_data, x);
      l_uint32 result;
      if (pixel <= 17)
        composeRGBPixel(0, 0, pixel * 15, &result);
      else if (pixel <= 145)
        composeRGBPixel(0, (pixel - 17) * 2, 255, &result);
      else
        composeRGBPixel((pixel - 145) * 2, 255, 255, &result);
      col_data[x] = result;
    }
  }
  ScrollView* win = new ScrollView("Projection", 0, 0,
                                   width, height, width, height);
  win->Image(pixc, 0, 0);
  win->Update();
  pixDestroy(&pixc);
#endif  // GRAPHICS_DISABLED
}

// Helper evaluates a vertical difference, (x,y) - (x,y-1), returning true if
// the difference, matches diff_sign and updating the best_diff, best_sum,
// best_y if a new max.
static bool EvaluateVerticalDiff(const l_uint32* data, int wpl, int diff_sign,
                                 int x, int y, int height,
                                 int* best_diff, int* best_sum, int* best_y) {
  if (y <= 0 || y >= height)
    return false;
  const l_uint32* line = data + y * wpl;
  int pixel1 = GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<const void *>(line - wpl)), x);
  int pixel2 = GET_DATA_BYTE(const_cast<void*> (reinterpret_cast<const void *>(line)), x);
  int diff = (pixel2 - pixel1) * diff_sign;
  if (diff > *best_diff) {
    *best_diff = diff;
    *best_sum = pixel1 + pixel2;
    *best_y = y;
  }
  return diff > 0;
}

/*!
 *  pixAdaptiveMeanFilter()
 *
 *      Input:  pixs   (8 bpp grayscale)
 *              wc, hc (half width/height of convolution kernel)
 *              varn   (value of overall noise variance)
 *      Return: pixd (8 bpp, filtered image)
 *
 *  Notes:
 *      (1) The filter reduces gaussian noise, achieving results similar
 *          to the arithmetic and geometric mean filters but avoiding the
 *          considerable image blurring effect introduced by those filters.
 *      (2) The filter can be expressed mathematically by:
 *            f'(x, y) = g(x, y) - varN / varL * [ g(x, y) - meanL ]
 *          where:
 *            -- g(x, y) is the pixel at the center of local region S of
 *               width (2 * wc + 1) and height (2 * wh + 1)
 *            -- varN and varL are the overall noise variance (given in input)
 *               and local variance of S, respectively
 *            -- meanL is the local mean of S
 *      (3) Typically @varn is estimated by studying the PDFs produced by
 *          the camera or equipment sensors.
 */
PIX *
pixAdaptiveMeanFilter(PIX       *pixs,
					  l_int32    wc,
					  l_int32    hc,
					  l_float32  varn)
{
	l_int32    i, j, w, h, d, wplt, wpld, wincr, hincr;
	l_uint32   val;
	l_uint32  *datat, *datad, *linet, *lined;
	l_float32  norm, meanl, varl, ratio;
	PIX       *pixt, *pixd;
	
    PROCNAME("pixAdaptiveMeanFilter");
    
    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 8)
        return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL);
    if (wc < 1 || hc < 1)
        return (PIX *)ERROR_PTR("wc and hc not >= 1", procName, NULL);
	
	/* Add wc to each side, and hc to top and bottom of the image,
	 * mirroring for accuracy and to avoid special-casing the boundary. */
    if ((pixt = pixAddMirroredBorder(pixs, wc, wc, hc, hc)) == NULL)
        return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
	
	/* Place the filter center at (0, 0).  This is just a
	 * convenient location, because it allows us to perform
	 * the filtering over x:(0 ... w - 1) and y:(0 ... h - 1). */
	pixd = pixCreateTemplate(pixs);
	wplt = pixGetWpl(pixt);
    wpld = pixGetWpl(pixd);
	datat = pixGetData(pixt);
    datad = pixGetData(pixd);
	
    wincr = 2 * wc + 1;
    hincr = 2 * hc + 1;
	norm = 1.0 / (wincr * hincr);
	for (i = 0; i < h; i++) {
        linet = datat + (i + hc) * wplt;
		lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            /* Calculate mean intensity value */
			meanl = calculateLocalMeanLow(datat, wplt, wincr, hincr, i, j);
			/* Calculate local variance */
			varl = calculateLocalVarianceLow(datat, wplt, wincr, hincr, i, j, meanl);
			/* Account for special case in which varN is more than varL */
			ratio = (varn > varl) ? 1 : varn / varl;
			val = GET_DATA_BYTE(linet, j + wc);
			SET_DATA_BYTE(lined, j, (l_uint8) (val - ratio * (val - meanl)));
        } 
    }
	
	pixDestroy(&pixt);	
    return pixd;
}

void SetFourierData(void **lines,              
                    struct fouriercomponents *fc,
                    unsigned int signal_size,
                    unsigned int signal_count,
                    int scan_mode,
                    int start,
                    int end) {
  
  fc->scan_mode = scan_mode;
  fc->signal_count = signal_count;
  fc->signal_size = signal_size;
  
  fc->signals   = (double **) malloc(sizeof(double)*signal_count);  
  fc->real      = (double **) malloc(sizeof(double)*signal_count);
  fc->imag      = (double **) malloc(sizeof(double)*signal_count);
  fc->freq      = (double **) malloc(sizeof(double)*signal_count);  
  fc->magnitude = (double **) malloc(sizeof(double)*signal_count);  
  fc->phase     = (double **) malloc(sizeof(double)*signal_count);    
  
  unsigned int x,y;
  if (scan_mode==0) {
    for (y=0;y<signal_count;y++) {
      fc->signals[y] = (double*)malloc(sizeof(double)*signal_size);
      for (x=0;x<signal_size;x++) {
        fc->signals[y][x] = GET_DATA_BYTE(lines[y],x);
      }
    }
  }

  else if (scan_mode==1) {
    unsigned int i,k;
    i=k=0;
    for (x=0;x<signal_count;x++) {
      fc->signals[i] = (double*)malloc(sizeof(double)*signal_size);
      k = 0;
      for (y=start;y<end;y++) {
        fc->signals[i][k] = GET_DATA_BYTE(lines[y],x);
        k++;
      }
      i++;
    }
  }
}

// Sends each pixel as hex value like html, e.g. #00FF00 for green.
void ScrollView::Transfer32bppImage(PIX* image) {
  int ppL = pixGetWidth(image);
  int h = pixGetHeight(image);
  int wpl = pixGetWpl(image);
  int transfer_size= ppL * 7 + 2;
  char* pixel_data = new char[transfer_size];
  for (int y = 0; y < h; ++y) {
    l_uint32* data = pixGetData(image) + y*wpl;
    for (int x = 0; x < ppL; ++x, ++data) {
      snprintf(&pixel_data[x*7], 7, "#%.2x%.2x%.2x",
               GET_DATA_BYTE(data, COLOR_RED),
               GET_DATA_BYTE(data, COLOR_GREEN),
               GET_DATA_BYTE(data, COLOR_BLUE));
    }
    pixel_data[transfer_size - 2] = '\n';
    pixel_data[transfer_size - 1] = '\0';
    SendRawMessage(pixel_data);
  }
  delete[] pixel_data;
}

/*!
 *  pixApplyLocalThreshold()
 *
 *      Input:  pixs (8 bpp grayscale; not colormapped)
 *              pixth (8 bpp array of local thresholds)
 *              redfactor ( ... )
 *      Return: pixd (1 bpp, thresholded image), or null on error
 */
PIX *
pixApplyLocalThreshold(PIX     *pixs,
                       PIX     *pixth,
                       l_int32  redfactor)
{
l_int32    i, j, w, h, wpls, wplt, wpld, vals, valt;
l_uint32  *datas, *datat, *datad, *lines, *linet, *lined;
PIX       *pixd;

    PROCNAME("pixApplyLocalThreshold");

    if (!pixs || pixGetDepth(pixs) != 8)
        return (PIX *)ERROR_PTR("pixs undefined or not 8 bpp", procName, NULL);
    if (pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("pixs is colormapped", procName, NULL);
    if (!pixth || pixGetDepth(pixth) != 8)
        return (PIX *)ERROR_PTR("pixth undefined or not 8 bpp", procName, NULL);

    pixGetDimensions(pixs, &w, &h, NULL);
    pixd = pixCreate(w, h, 1);
    datas = pixGetData(pixs);
    datat = pixGetData(pixth);
    datad = pixGetData(pixd);
    wpls = pixGetWpl(pixs);
    wplt = pixGetWpl(pixth);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        linet = datat + i * wplt;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            vals = GET_DATA_BYTE(lines, j);
            valt = GET_DATA_BYTE(linet, j);
            if (vals < valt)
                SET_DATA_BIT(lined, j);
        }
    }

    return pixd;
}

/*!
 *  addGrayLow()
 */
void
addGrayLow(l_uint32  *datad,
           l_int32    w,
           l_int32    h,
           l_int32    d,
           l_int32    wpld,
           l_uint32  *datas,
           l_int32    wpls)
{
    l_int32    i, j, val, sum;
    l_uint32  *lines, *lined;


    for (i = 0; i < h; i++) {
        lined = datad + i * wpld;
        lines = datas + i * wpls;
        if (d == 8) {
            for (j = 0; j < w; j++) {
                sum = GET_DATA_BYTE(lines, j) + GET_DATA_BYTE(lined, j);
                val = L_MIN(sum, 255);
                SET_DATA_BYTE(lined, j, val);
            }
        }
        else if (d == 16) {
            for (j = 0; j < w; j++) {
                sum = GET_DATA_TWO_BYTES(lines, j)
                      + GET_DATA_TWO_BYTES(lined, j);
                val = L_MIN(sum, 0xffff);
                SET_DATA_TWO_BYTES(lined, j, val);
            }
        }
        else {   /* d == 32; no clipping */
            for (j = 0; j < w; j++)
                *(lined + j) += *(lines + j);
        }
    }

    return;
}