Blob Blame Raw


/*
guardare DAFARE
write_lut_image ??? quando ?

*/

#include "autoadjust.h"
#include "cleanupcommon.h"

#define TO_ABS(x)       \
	{                   \
		if ((x) < 0)    \
			(x) = -(x); \
	}

typedef struct big {
	UINT lo, hi;
} BIG;
#define CLEAR_BIG(B) ((B).lo = 0, (B).hi = 0, (B))
#define ADD_BIG(B, X) ((B).lo += (UINT)(X), \
					   (B).hi += (B).lo >> 30, (B).lo &= 0x3fffffff, (B))
#define ADD_BIG_BIG(B1, B2) ((B1).lo += (B2).lo, (B1).hi += (B2).hi, \
							 (B1).hi += (B1).lo >> 30, (B1).lo &= 0x3fffffff, (B1))
#define BIG_TO_DOUBLE(B) ((double)(B).hi * (double)0x40000000 + (double)(B).lo)

static int Black = 0;

static int Ref_cum[256];
static int Ref_edgelen;

/*---------------------------------------------------------------------------*/

void autoadj_set_black_value(int black)
{
	Black = black;
}

/*===========================================================================*/

static void smooth_func256(float func[256], int rad)
{
	int i, j, k;
	float smooth[256];

	for (i = 0; i < 256; i++) {
		k = i - rad;
		notLessThan(0, k);
		smooth[i] = func[k] / 2;
		k = i + rad;
		notMoreThan(255, k);
		smooth[i] += func[k] / 2;
		for (j = i - rad + 1; j <= i + rad - 1; j++) {
			k = j;
			notLessThan(0, k);
			notMoreThan(255, k);
			smooth[i] += func[k];
		}
		smooth[i] /= 2 * rad;
	}
	for (i = 0; i < 256; i++)
		func[i] = smooth[i];
}

/*---------------------------------------------------------------------------*/

static void build_cum(int histo[256], int cum[256])
{
	int i;

	cum[0] = histo[0];
	for (i = 1; i < 256; i++)
		cum[i] = cum[i - 1] + histo[i];
}

/*===========================================================================*/

static int Window_x0, Window_y0, Window_x1, Window_y1;

void set_autoadjust_window(int x0, int y0, int x1, int y1)
{
	Window_x0 = x0;
	Window_y0 = y0;
	Window_x1 = x1;
	Window_y1 = y1;
	if (Window_x0 > Window_x1)
		tswap(Window_x0, Window_x1);
	if (Window_y0 > Window_y1)
		tswap(Window_y0, Window_y1);
}

/*---------------------------------------------------------------------------*/

static void get_virtual_buffer(const TRasterImageP &image,
							   int *p_lx, int *p_ly, int *p_wrap,
							   UCHAR **p_buffer)
{
	int x0, y0, x1, y1;
	int x_margin, y_margin;
	int lx, ly, wrap;
	UCHAR *buffer;

	TRasterGR8P ras8(image->getRaster());

	assert(ras8);
	double xdpi, ydpi;
	image->getDpi(xdpi, ydpi);

	/* BORDO DI MEZZO CENTIMETRO */
	x_margin = troundp(mmToPixel(5.0, xdpi));
	y_margin = troundp(mmToPixel(5.0, ydpi));
	x0 = Window_x0 + x_margin;
	y0 = Window_y0 + y_margin;
	x1 = Window_x1 - x_margin;
	y1 = Window_y1 - y_margin;
	notLessThan(x0 + 9, x1);
	notLessThan(y0 + 9, y1);
	notLessThan(0, x0);
	notMoreThan(ras8->getLx() - 1, x0);
	notLessThan(0, y0);
	notMoreThan(ras8->getLy() - 1, y0);
	notLessThan(0, x1);
	notMoreThan(ras8->getLx() - 1, x1);
	notLessThan(0, y1);
	notMoreThan(ras8->getLy() - 1, y1);

	lx = x1 - x0 + 1;
	ly = y1 - y0 + 1;
	wrap = ras8->getWrap();
	buffer = (UCHAR *)ras8->getRawData() + x0 + y0 * wrap;

	*p_lx = lx;
	*p_ly = ly;
	*p_wrap = wrap;
	*p_buffer = buffer;
}

/*===========================================================================*/

void black_eq_algo(const TRasterImageP &image)
{
	int lx, ly, wrap;
	int x, y, grey /*, width*/;
	int d_histo[256][256];
	int d, dd, m;
	UCHAR *buffer, *pix, *north, *south;
	UCHAR prev, darkest;
	long n;
	float mean_d[256];
	BIG s;
	int max_d_grey;
	float max_d;
	float fac;
	int val;
	UCHAR lut[256];
	image->getRaster()->lock();

	get_virtual_buffer(image, &lx, &ly, &wrap, &buffer);

	for (grey = 0; grey < 256; grey++) {
		for (d = 0; d < 256; d++)
			d_histo[grey][d] = 0;
	}
	for (y = 1; y < ly - 1; y++) {
		pix = buffer + y * wrap + 1;
		north = pix + wrap;
		south = pix - wrap;
		for (x = 1; x < lx - 1; x++, pix++, north++, south++) {
			m = 2 * ((int)north[0] + (int)pix[-1] + 2 * (int)pix[0] + (int)pix[1] + (int)south[0]) + (int)north[-1] + (int)north[1] + (int)south[-1] + (int)south[1];
			m = (m + 8) >> 4;
			d = (int)north[-1] + (int)north[0] * 2 + (int)north[1] - (int)south[-1] - (int)south[0] * 2 - (int)south[1];
			TO_ABS(d)
			dd = (int)north[-1] + (int)pix[-1] * 2 + (int)south[-1] - (int)north[1] - (int)pix[1] * 2 - (int)south[1];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			dd = (int)pix[-1] + (int)north[-1] * 2 + (int)north[0] - (int)pix[1] - (int)south[1] * 2 - (int)south[0];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			dd = (int)pix[1] + (int)north[1] * 2 + (int)north[0] - (int)pix[-1] - (int)south[-1] * 2 - (int)south[0];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			d = (d + 2) >> 2;
			d_histo[m][d]++;
		}
	}
	for (grey = 0; grey < 256; grey++) {
		n = 0;
		CLEAR_BIG(s);
		for (d = 0; d < 256; d++) {
			ADD_BIG(s, d_histo[grey][d] * d);
			n += d_histo[grey][d];
		}
		if (n)
			mean_d[grey] = (float)(BIG_TO_DOUBLE(s) / n);
		else
			mean_d[grey] = 0;
	}
	smooth_func256(mean_d, 5);
	max_d_grey = 0;
	max_d = 0.0;
	for (grey = 0; grey < 256; grey++) {
		if (max_d < mean_d[grey]) {
			max_d = mean_d[grey];
			max_d_grey = grey;
		}
	}

	n = 0;
	CLEAR_BIG(s);
	for (y = 0; y < ly; y++) {
		pix = buffer + y * wrap;
		prev = 255;
		darkest = 255;
		for (x = 0; x < lx; x++, pix++) {
			if (*pix < max_d_grey) {
				if (prev < max_d_grey)
					notMoreThan(*pix, darkest);
				else
					darkest = *pix;
			} else {
				if (prev < max_d_grey) {
					ADD_BIG(s, darkest);
					n++;
				}
			}
			prev = *pix;
		}
	}
	darkest = troundp(BIG_TO_DOUBLE(s) / n);

	/*
for (grey = 0; grey < darkest; grey++)
  lut[grey] = 0;
fac = (float)255 / (float)(255 - darkest);
for (grey = darkest; grey < 256; grey++)
  lut[grey] = 255 - troundp ((255 - grey) * fac);
*/
	notLessThan(0, Black);
	;
	notMoreThan(255, Black);
	fac = (float)(255 - Black) / (float)(255 - darkest);
	for (grey = 0; grey < 256; grey++) {
		val = 255 - troundp((255 - grey) * fac);
		notLessThan(0, val);
		notMoreThan(255, val);
		lut[grey] = val;
	}

	apply_lut(image, lut);
	image->getRaster()->unlock();
}

/*===========================================================================*/

#define MAX_WIDTH 100

/*---------------------------------------------------------------------------*/

void build_lw(const TRasterImageP &image, float lw[256])
{
	int lx, ly, wrap;
	int x, y, grey, width;
	int d_histo[256][256];
	int d, dd;
	int lw_histo[256][MAX_WIDTH + 1];
	int x_start[256];
	static int *y_start[256];
	static int y_start_alloc = 0;
	UCHAR *buffer, *pix, *north, *south;
	UCHAR cur_grey, x_prev_grey, x_next_grey, y_prev_grey, y_next_grey;
	int x_grad, y_grad;
	long n;
	float mean_d[256];
	BIG s;
	int max_d_grey;
	float max_d;
	long max_d_n;

	get_virtual_buffer(image, &lx, &ly, &wrap, &buffer);

	if (y_start_alloc < lx) {
		for (grey = 0; grey < 256; grey++) {
			delete[] y_start[grey];
			y_start[grey] = new int[lx];
			//TREALLOC (y_start[grey], lx);
		}
		y_start_alloc = lx;
	}
	for (grey = 0; grey < 256; grey++) {
		for (d = 0; d < 256; d++)
			d_histo[grey][d] = 0;
		for (x = 0; x < lx; x++)
			y_start[grey][x] = -MAX_WIDTH - 1;
		for (width = 0; width <= MAX_WIDTH; width++)
			lw_histo[grey][width] = 0;
	}

	for (y = 1; y < ly - 1; y++) {
		for (grey = 0; grey < 256; grey++)
			x_start[grey] = -MAX_WIDTH - 1;
		pix = buffer + y * wrap + 1;
		north = pix + wrap;
		south = pix - wrap;
		for (x = 1; x < lx - 1; x++, pix++, north++, south++) {
			d = (int)north[-1] + (int)north[0] * 2 + (int)north[1] - (int)south[-1] - (int)south[0] * 2 - (int)south[1];
			TO_ABS(d)
			dd = (int)north[-1] + (int)pix[-1] * 2 + (int)south[-1] - (int)north[1] - (int)pix[1] * 2 - (int)south[1];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			dd = (int)pix[-1] + (int)north[-1] * 2 + (int)north[0] - (int)pix[1] - (int)south[1] * 2 - (int)south[0];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			dd = (int)pix[1] + (int)north[1] * 2 + (int)north[0] - (int)pix[-1] - (int)south[-1] * 2 - (int)south[0];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			d >>= 3;
			d_histo[*pix][d]++;

			cur_grey = pix[0];
			x_prev_grey = pix[-1];
			x_next_grey = pix[1];
			y_prev_grey = south[0];
			y_next_grey = north[0];
			x_grad = (int)x_next_grey - (int)x_prev_grey;
			y_grad = (int)y_next_grey - (int)y_prev_grey;
			if (cur_grey < x_prev_grey) {
				if (x_grad < 0 && x_grad < y_grad && x_grad < -y_grad)
					for (grey = cur_grey; grey < x_prev_grey; grey++)
						x_start[grey] = x;
				else
					for (grey = cur_grey; grey < x_prev_grey; grey++)
						x_start[grey] = -MAX_WIDTH - 1;
			} else if (cur_grey > x_prev_grey) {
				if (x_grad > 0 && x_grad > y_grad && x_grad > -y_grad)
					for (grey = x_prev_grey; grey < cur_grey; grey++) {
						width = x - x_start[grey];
						if (width <= MAX_WIDTH)
							lw_histo[grey][width]++;
					}
			}
			if (cur_grey < y_prev_grey) {
				if (y_grad < 0 && y_grad < x_grad && y_grad < -x_grad)
					for (grey = cur_grey; grey < y_prev_grey; grey++)
						y_start[grey][x] = y;
				else
					for (grey = cur_grey; grey < y_prev_grey; grey++)
						y_start[grey][x] = -MAX_WIDTH - 1;
			} else if (cur_grey > y_prev_grey) {
				if (y_grad > 0 && y_grad > x_grad && y_grad > -x_grad)
					for (grey = y_prev_grey; grey < cur_grey; grey++) {
						width = y - y_start[grey][x];
						if (width <= MAX_WIDTH)
							lw_histo[grey][width]++;
					}
			}
		}
	}
	for (grey = 0; grey < 256; grey++) {
		n = 0;
		CLEAR_BIG(s);
		for (width = 0; width <= MAX_WIDTH; width++) {
			ADD_BIG(s, lw_histo[grey][width] * width);
			n += lw_histo[grey][width];
		}
		if (n)
			lw[grey] = (float)(BIG_TO_DOUBLE(s) / n);
		else
			lw[grey] = 0.0;
	}
	max_d_grey = 0;
	max_d = 0.0;
	for (grey = 0; grey < 256; grey++) {
		n = 0;
		CLEAR_BIG(s);
		for (d = 0; d < 256; d++) {
			ADD_BIG(s, d_histo[grey][d] * d);
			n += d_histo[grey][d];
		}
		if (n)
			mean_d[grey] = (float)(BIG_TO_DOUBLE(s) / n);
		else
			mean_d[grey] = 0.0;
		if (max_d < mean_d[grey]) {
			max_d = mean_d[grey];
			max_d_grey = grey;
		}
	}

	max_d_n = 0;
	for (width = 0; width <= MAX_WIDTH; width++)
		max_d_n += lw_histo[max_d_grey][width];
	for (grey = max_d_grey - 1; grey >= 0; grey--) {
		if (!lw[grey])
			break;
		n = 0;
		for (width = 0; width <= MAX_WIDTH; width++)
			n += lw_histo[grey][width];
		if (n < max_d_n / 10)
			break;
	}
	for (; grey >= 0; grey--)
		lw[grey] = 0.0;
	for (grey = max_d_grey + 1; grey < 256; grey++) {
		if (!lw[grey])
			break;
		if (255 - grey < (255 - max_d_grey) / 2)
			break;
		if (mean_d[grey] < max_d / 2.0)
			break;
	}
	for (; grey < 256; grey++)
		lw[grey] = 0.0;
}

/*---------------------------------------------------------------------------*/

void build_lw_lut(float ref_lw[256], float lw[256], UCHAR lut[256])
{
	/* crea una lut tale che l'immagine con il profilo di linewidths lw
   venga mappata in un'immagine con il profilo di riferimento ref_lw.
   Le lw sono non decrescenti e delimitate da 0.0 li' dove non sono valide
*/

	int i, j;
	float bot_ref_lw, top_ref_lw, bot_lw, top_lw;
	int bot_ref_gr, top_ref_gr, bot_gr, top_gr;
	float min_lw, max_lw;
	int min_ref_gr, max_ref_gr, min_gr, max_gr;
	float fac;

	for (i = 0; !ref_lw[i]; i++) {
	}
	bot_ref_lw = ref_lw[i];
	bot_ref_gr = i;
	for (i = 255; !ref_lw[i]; i--) {
	}
	top_ref_lw = ref_lw[i];
	top_ref_gr = i;
	for (i = 0; !lw[i]; i++) {
	}
	bot_lw = lw[i];
	bot_gr = i;
	for (i = 255; !lw[i]; i--) {
	}
	top_lw = lw[i];
	top_gr = i;

	min_lw = tmax(bot_ref_lw, bot_lw);
	max_lw = tmin(top_ref_lw, top_lw);

	if (min_lw >= max_lw) {
		for (i = 0; i < 256; i++)
			lut[i] = i;
		return;
	}

	for (i = bot_ref_gr; ref_lw[i] < min_lw; i++) {
	}
	min_ref_gr = i;
	for (i = top_ref_gr; ref_lw[i] > max_lw; i--) {
	}
	max_ref_gr = i;
	for (i = bot_gr; lw[i] < min_lw; i++) {
	}
	min_gr = i;
	for (i = top_gr; lw[i] > max_lw; i--) {
	}
	max_gr = i;

	j = min_ref_gr;
	for (i = min_gr; i <= max_gr; i++) {
		while (ref_lw[j] < lw[i] && j < max_ref_gr)
			j++;
		lut[i] = j;
	}
	fac = (float)min_ref_gr / (float)min_gr;
	for (i = 0; i < min_gr; i++)
		lut[i] = troundp(i * fac);
	fac = (float)(255 - max_ref_gr) / (float)(255 - max_gr);
	for (i = 255; i > max_gr; i--)
		lut[i] = 255 - troundp((255 - i) * fac);

	/***
printf ("-----\n\lut:\n\n");
for (i = 255; i >= 0; i--)
  printf ("%4d :%4u\n", i, lut[i]);
printf ("\n");
***/
}

/*===========================================================================*/

void build_gr_cum(const TRasterImageP &image, int cum[256])
{
	int lx, ly, wrap, true_lx, true_ly;
	int i, x, y;
	UCHAR *pix, *buffer;
	int histo[256], raster_is_savebox;

	get_virtual_buffer(image, &lx, &ly, &wrap, &buffer);

	for (i = 0; i < 256; i++)
		histo[i] = 0;
	for (y = 0; y < ly; y++) {
		pix = buffer + y * wrap;
		for (x = 0; x < lx; x++)
			histo[*pix++]++;
	}

	raster_is_savebox = 1;
	TRect saveBox = image->getSavebox();
	if ((saveBox.getLx() > 0 && saveBox.getLx() < image->getRaster()->getLx()) ||
		(saveBox.getLy() > 0 && saveBox.getLy() < image->getRaster()->getLy()))
		raster_is_savebox = 0;

	if (raster_is_savebox) {
		true_lx = saveBox.getLx() ? saveBox.getLx() : image->getRaster()->getLx();
		true_ly = saveBox.getLy() ? saveBox.getLy() : image->getRaster()->getLy();
	} else {
		true_lx = image->getRaster()->getLx();
		true_ly = image->getRaster()->getLy();
	}
	histo[255] += true_lx * true_ly - lx * ly;

	build_cum(histo, cum);
}

/*---------------------------------------------------------------------------*/

void build_gr_lut(int ref_cum[256], int cum[256], UCHAR lut[256])
{
	/* crea una lut tale che l'immagine con la distribuzione cumulativa cum
   venga mappata in un'immagine con la cumulativa di riferimento ref_cum
*/
	int i, j, n;

	j = 0;
	for (i = 0; i < 256; i++) {
		n = cum[i];
		while (ref_cum[j] < n)
			j++;
		lut[i] = j;
	}
}

/*===========================================================================*/

struct edge_config {
	char *str;
	int val;
};

static struct edge_config Edge_base[] =
	{
		{"     "
		 "  X  "
		 "     ",
		 0},

		{"X    "
		 "  X  "
		 "     ",
		 0},

		{"  X  "
		 "  X  "
		 "     ",
		 0},

		{"X X  "
		 "  X  "
		 "     ",
		 5},

		{"X   X"
		 "  X  "
		 "     ",
		 20},

		{"X    "
		 "  X X"
		 "     ",
		 24},

		{"X    "
		 "  X  "
		 "    X",
		 28},

		{"  X  "
		 "  X X"
		 "     ",
		 7},

		{"  X  "
		 "  X  "
		 "  X  ",
		 20},

		{"X X X"
		 "  X  "
		 "     ",
		 10},

		{"X X  "
		 "  X X"
		 "     ",
		 12},

		{"X X  "
		 "  X  "
		 "    X",
		 26},

		{"X X  "
		 "  X  "
		 "  X  ",
		 22},

		{"X X  "
		 "  X  "
		 "X    ",
		 22},

		{"X X  "
		 "X X  "
		 "     ",
		 5},

		{"X   X"
		 "  X  "
		 "    X",
		 34},

		{"X   X"
		 "  X  "
		 "  X  ",
		 34},

		{"X    "
		 "  X X"
		 "  X  ",
		 24},

		{"  X  "
		 "  X X"
		 "  X  ",
		 10},

		{"X X X"
		 "  X X"
		 "     ",
		 12},

		{"X X X"
		 "  X  "
		 "    X",
		 24},

		{"X X X"
		 "  X  "
		 "  X  ",
		 24},

		{"X X  "
		 "  X X"
		 "    X",
		 14},

		{"X X  "
		 "  X X"
		 "  X  ",
		 12},

		{"X X  "
		 "  X X"
		 "X    ",
		 22},

		{"X X  "
		 "X X X"
		 "     ",
		 10},

		{"X X  "
		 "  X  "
		 "  X X",
		 24},

		{"X X  "
		 "  X  "
		 "X   X",
		 32},

		{"X X  "
		 "X X  "
		 "    X",
		 24},

		{"X X  "
		 "  X  "
		 "X X  ",
		 20},

		{"X   X"
		 "  X  "
		 "X   X",
		 40},

		{"  X  "
		 "X X X"
		 "  X  ",
		 0},

		{"X X X"
		 "X X X"
		 "     ",
		 10},

		{"X X X"
		 "X X  "
		 "    X",
		 22},

		{"X X  "
		 "X X X"
		 "    X",
		 12},

		{"X   X"
		 "X X X"
		 "    X",
		 22},

		{"  X X"
		 "X X X"
		 "    X",
		 12},

		{"X X X"
		 "  X X"
		 "    X",
		 14},

		{"X X  "
		 "X X X"
		 "  X  ",
		 0},

		{"X   X"
		 "X X X"
		 "  X  ",
		 10},

		{"X   X"
		 "X X  "
		 "  X X",
		 20},

		{"X   X"
		 "X X  "
		 "X   X",
		 30},

		{"X X X"
		 "X X X"
		 "    X",
		 12},

		{"X X X"
		 "X X X"
		 "  X  ",
		 0},

		{"X X X"
		 "X X  "
		 "  X X",
		 10},

		{"  X X"
		 "X X X"
		 "X X  ",
		 0},

		{"X X X"
		 "X X  "
		 "X   X",
		 20},

		{"X   X"
		 "X X X"
		 "X   X",
		 20},

		{"X X X"
		 "X X X"
		 "  X X",
		 0},

		{"X X X"
		 "X X X"
		 "X   X",
		 10},

		{"X X X"
		 "X X X"
		 "X X X",
		 0},

		{0, 0},
};

static int Edge_value[256];
static int Edge_init_done = 0;

static int edge_int(char x[8]);
static void edge_rotate(char x[8]);
static void edge_mirror(char x[8]);

/*---------------------------------------------------------------------------*/

static void edge_init(void)
{
	int b;
	char *str;
	int val;
	char x[8];

	for (b = 0; Edge_base[b].str; b++) {
		str = Edge_base[b].str;
		val = Edge_base[b].val;
		x[0] = str[0];
		x[1] = str[2];
		x[2] = str[4];
		x[3] = str[5];
		x[4] = str[9];
		x[5] = str[10];
		x[6] = str[12];
		x[7] = str[14];

		Edge_value[edge_int(x)] = val;
		edge_rotate(x);
		Edge_value[edge_int(x)] = val;
		edge_rotate(x);
		Edge_value[edge_int(x)] = val;
		edge_rotate(x);
		Edge_value[edge_int(x)] = val;
		edge_rotate(x);
		edge_mirror(x);
		Edge_value[edge_int(x)] = val;
		edge_rotate(x);
		Edge_value[edge_int(x)] = val;
		edge_rotate(x);
		Edge_value[edge_int(x)] = val;
		edge_rotate(x);
		Edge_value[edge_int(x)] = val;
	}
	Edge_init_done = 1;
}

/*---------------------------------------------------------------------------*/

static int edge_int(char x[8])
{
	return (x[0] != ' ') << 7 |
		   (x[1] != ' ') << 6 |
		   (x[2] != ' ') << 5 |
		   (x[3] != ' ') << 4 |
		   (x[4] != ' ') << 3 |
		   (x[5] != ' ') << 2 |
		   (x[6] != ' ') << 1 |
		   (x[7] != ' ');
}

/*---------------------------------------------------------------------------*/

static void edge_rotate(char x[8])
{
	char tmp;

	tmp = x[0];
	x[0] = x[2];
	x[2] = x[7];
	x[7] = x[5];
	x[5] = tmp;
	tmp = x[1];
	x[1] = x[4];
	x[4] = x[6];
	x[6] = x[3];
	x[3] = tmp;
}

/*---------------------------------------------------------------------------*/

static void edge_mirror(char x[8])
{
	char tmp;

	tmp = x[0];
	x[0] = x[2];
	x[2] = tmp;
	tmp = x[3];
	x[3] = x[4];
	x[4] = tmp;
	tmp = x[5];
	x[5] = x[7];
	x[7] = tmp;
}

/*---------------------------------------------------------------------------*/

void histo_l_algo(const TRasterImageP &image, int reference)
{
	int lx, ly, wrap;
	int x, y, grey, ref_grey, m, d, dd;
	UCHAR *buffer, *pix, *north, *south;
	int g_histo[256];
	int d_histo[256][256];
	int cum[256];
	long n;
	BIG s;
	float mean_d[256];
	int max_d_grey;
	float max_d;
	int conf;
	int edgelen;
	float fac;
	UCHAR lut[256];

	if (!Edge_init_done)
		edge_init();

	get_virtual_buffer(image, &lx, &ly, &wrap, &buffer);

	for (grey = 0; grey < 256; grey++) {
		g_histo[grey] = 0;
		for (d = 0; d < 256; d++)
			d_histo[grey][d] = 0;
	}

	/* istogramma dei grigi e istogramma delle derivate */
	for (y = 1; y < ly - 1; y++) {
		pix = buffer + y * wrap + 1;
		north = pix + wrap;
		south = pix - wrap;
		for (x = 1; x < lx - 1; x++, pix++, north++, south++) {
			m = 2 * ((int)north[0] + (int)pix[-1] + 2 * (int)pix[0] + (int)pix[1] + (int)south[0]) + (int)north[-1] + (int)north[1] + (int)south[-1] + (int)south[1];
			m = (m + 8) >> 4;
			d = (int)north[-1] + (int)north[0] * 2 + (int)north[1] - (int)south[-1] - (int)south[0] * 2 - (int)south[1];
			TO_ABS(d)
			dd = (int)north[-1] + (int)pix[-1] * 2 + (int)south[-1] - (int)north[1] - (int)pix[1] * 2 - (int)south[1];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			dd = (int)pix[-1] + (int)north[-1] * 2 + (int)north[0] - (int)pix[1] - (int)south[1] * 2 - (int)south[0];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			dd = (int)pix[1] + (int)north[1] * 2 + (int)north[0] - (int)pix[-1] - (int)south[-1] * 2 - (int)south[0];
			TO_ABS(dd)
			if (dd > d)
				d = dd;
			d = (d + 2) >> 2;
			d_histo[m][d]++;
			g_histo[*pix]++;
		}
	}
	build_cum(g_histo, cum);

	/* costruzione pendenze medie */
	for (grey = 0; grey < 256; grey++) {
		n = 0;
		CLEAR_BIG(s);
		for (d = 0; d < 256; d++) {
			ADD_BIG(s, d_histo[grey][d] * d);
			n += d_histo[grey][d];
		}
		if (n)
			mean_d[grey] = (float)(BIG_TO_DOUBLE(s) / n);
		else
			mean_d[grey] = 0.0;
	}
	smooth_func256(mean_d, 5);

	/* determinazione grigio di massima pendenza */
	max_d_grey = 0;
	max_d = 0.0;
	for (grey = 0; grey < 256; grey++) {
		if (max_d < mean_d[grey]) {
			max_d = mean_d[grey];
			max_d_grey = grey;
		}
	}

	/* stima della lunghezza dei bordi */
	edgelen = 0;
	for (y = 1; y < ly - 1; y++) {
		pix = buffer + y * wrap + 1;
		north = pix + wrap;
		south = pix - wrap;
		conf = -1;
		for (x = 1; x < lx - 1; x++, pix++, north++, south++) {
			if (pix[0] <= max_d_grey) {
				if (conf >= 0) {
					conf = ((conf << 1) & ((1 << 7) | (1 << 6) | (1 << 2) | (1 << 1))) |
						   ((north[1] <= max_d_grey) << 5) |
						   ((1) << 4) |
						   ((pix[1] <= max_d_grey) << 3) |
						   ((south[1] <= max_d_grey));
				} else {
					conf = ((north[-1] <= max_d_grey) << 7) |
						   ((north[0] <= max_d_grey) << 6) |
						   ((north[1] <= max_d_grey) << 5) |
						   ((pix[-1] <= max_d_grey) << 4) |
						   ((pix[1] <= max_d_grey) << 3) |
						   ((south[-1] <= max_d_grey) << 2) |
						   ((south[0] <= max_d_grey) << 1) |
						   ((south[1] <= max_d_grey));
				}
				edgelen += Edge_value[conf];
			} else
				conf = -1;
		}
	}

	if (reference) {
		for (grey = 0; grey < 256; grey++)
			Ref_cum[grey] = cum[grey];
		Ref_edgelen = edgelen;

		return;
	}

	/* normalizza la cumulativa per il numero di linee */
	fac = (float)Ref_edgelen / edgelen;

	/***
printf ("fac: %f\n", fac);
printf ("cum prima:\n");
for (grey = 0; grey < 256; grey++)
  printf ("%9d\n", cum[grey]);
***/

	for (grey = 0; grey < 255; grey++) {
		cum[grey] = (int)(cum[grey] * fac);
		notMoreThan(cum[255], cum[grey]);
	}

	/***
printf ("cum dopo:\n");
for (grey = 0; grey < 256; grey++)
  printf ("%9d\n", cum[grey]);
printf ("Ref_cum:\n");
for (grey = 0; grey < 256; grey++)
  printf ("%9d\n", Ref_cum[grey]);
***/

	/* equalizza l'istogramma */
	ref_grey = 0;
	for (grey = 0; grey < 255; grey++) {
		while (Ref_cum[ref_grey] < cum[grey])
			ref_grey++;
		lut[grey] = ref_grey;
	}
	lut[255] = 255;
	/* DAFARE
if (Wl_flag)
  write_lut_image (lut);
*/
	apply_lut(image, lut);
}

/*===========================================================================*/
/*===========================================================================*/

#define MAX_N_CHAINS 100

#define MAX_WIDTH 100

#define MIN_COUNT 100

#define MAX_HGREY (244 >> 1)

/*---------------------------------------------------------------------------*/
#ifdef DAFARE
static int build_th_histo(const TRasterImageP &image, int histo[256 >> 1][MAX_WIDTH + 1])
{
	int lx, ly, wrap, x, y, hgrey, width;
	UCHAR *buffer, *cur, *x_prev, *y_prev, *x_next, *y_next;
	UCHAR cur_hgrey, x_prev_hgrey, y_prev_hgrey, x_next_hgrey, y_next_hgrey;
	int x_grad, y_grad;
	int hgrey_histo[256 >> 1];
	int x_start[256 >> 1];
	static int *y_start[256 >> 1];
	static int y_start_alloc = 0;

	get_virtual_buffer(image, &lx, &ly, &wrap, &buffer);

	if (y_start_alloc < lx) {
		for (hgrey = 0; hgrey<256>> 1; hgrey++) {
			delete[] y_start[hgrey];
			y_start[hgrey] = new int[lx];
			//    TREALLOC (y_start[hgrey], lx);
		}
		y_start_alloc = lx;
	}
	for (hgrey = 0; hgrey<256>> 1; hgrey++) {
		for (x = 0; x < lx; x++)
			y_start[hgrey][x] = -MAX_WIDTH - 1;
		for (width = 0; width <= MAX_WIDTH; width++)
			histo[hgrey][width] = 0;
		hgrey_histo[hgrey] = 0;
	}

	for (y = 1; y < ly - 1; y++) {
		for (hgrey = 0; hgrey<256>> 1; hgrey++)
			x_start[hgrey] = -MAX_WIDTH - 1;
		cur = buffer + y * wrap + 1;
		x_prev = cur - 1;
		y_prev = cur - wrap;
		x_next = cur + 1;
		y_next = cur + wrap;
		for (x = 1; x < lx - 1; x++, cur++, x_prev++, y_prev++, x_next++, y_next++) {
			cur_hgrey = *cur >> 1;
			x_prev_hgrey = *x_prev >> 1;
			x_next_hgrey = *x_next >> 1;
			y_prev_hgrey = *y_prev >> 1;
			y_next_hgrey = *y_next >> 1;
			hgrey_histo[cur_hgrey]++;
			x_grad = (int)x_next_hgrey - (int)x_prev_hgrey;
			y_grad = (int)y_next_hgrey - (int)y_prev_hgrey;
			if (cur_hgrey < x_prev_hgrey) {
				if (x_grad < 0 && x_grad < y_grad && x_grad < -y_grad)
					for (hgrey = cur_hgrey; hgrey < x_prev_hgrey; hgrey++)
						x_start[hgrey] = x;
				else
					for (hgrey = cur_hgrey; hgrey < x_prev_hgrey; hgrey++)
						x_start[hgrey] = -MAX_WIDTH - 1;
			} else if (cur_hgrey > x_prev_hgrey) {
				if (x_grad > 0 && x_grad > y_grad && x_grad > -y_grad)
					for (hgrey = x_prev_hgrey; hgrey < cur_hgrey; hgrey++) {
						width = x - x_start[hgrey];
						if (width <= MAX_WIDTH)
							histo[hgrey][width]++;
					}
			}
			if (cur_hgrey < y_prev_hgrey) {
				if (y_grad < 0 && y_grad < x_grad && y_grad < -x_grad)
					for (hgrey = cur_hgrey; hgrey < y_prev_hgrey; hgrey++)
						y_start[hgrey][x] = y;
				else
					for (hgrey = cur_hgrey; hgrey < y_prev_hgrey; hgrey++)
						y_start[hgrey][x] = -MAX_WIDTH - 1;
			} else if (cur_hgrey > y_prev_hgrey) {
				if (y_grad > 0 && y_grad > x_grad && y_grad > -x_grad)
					for (hgrey = y_prev_hgrey; hgrey < cur_hgrey; hgrey++) {
						width = y - y_start[hgrey][x];
						if (width <= MAX_WIDTH)
							histo[hgrey][width]++;
					}
			}
		}
	}
	return 1;
}
#endif
/*---------------------------------------------------------------------------*/
#ifdef DAFARE
//SEMBRA NON ESSERE UTILIZZATO!
int eval_image_th(const TRasterImageP &image, int *threshold, float *linewidth)
{
#define MIN_COUNT 100
#define MAX_HGREY (244 >> 1)
	int histo[256 >> 1][MAX_WIDTH + 1];
	char is_max[256 >> 1][MAX_WIDTH + 1];
	int chain_map[256 >> 1][MAX_WIDTH + 1];
	int chain_n[MAX_N_CHAINS];
	int n_chains;
	int peak[256 >> 1];
	float mean_width[256 >> 1];
	int min_hgrey, max_hgrey, min_chain_hgrey, max_chain_hgrey;
	int hgrey, width, chain, c, last;
	long n, nx;
	float delta, min_delta;

	TRasterGR8P ras8(image->getRaster());

	if (!ras8) {
		TERROR("image is not RAS_GR8");
		return 0;
	}
	if (!(*threshold != 0 && *linewidth == 0 ||
		  *threshold == 0 && *linewidth != 0))
		return 0;

	build_th_histo(image, histo);

	if (*threshold != 0) {
		min_hgrey = *threshold >> 1;
		max_hgrey = *threshold >> 1;
		notLessThan(0, min_hgrey);
		notMoreThan<int>(MAX_HGREY, max_hgrey);
	} else {
		for (hgrey = 0; hgrey<256>> 1; hgrey++)
			for (width = 1; width < MAX_WIDTH; width++)
				if (histo[hgrey][width] >= MIN_COUNT)
					goto endfor;
	endfor:
		if (hgrey == 256 >> 1)
			return 0;

		min_hgrey = hgrey;
		max_hgrey = MAX_HGREY;
	}
	if (min_hgrey > max_hgrey)
		return 0;

	/* costruisci la matrice dai massimi */
	for (hgrey = min_hgrey; hgrey <= max_hgrey; hgrey++) {
		is_max[hgrey][0] = 0;
		is_max[hgrey][MAX_WIDTH] = 0;
		for (width = 1; width < MAX_WIDTH; width++) {
			is_max[hgrey][width] = histo[hgrey][width] > histo[hgrey][width - 1] &&
								   histo[hgrey][width] > histo[hgrey][width + 1] &&
								   histo[hgrey][width] >= MIN_COUNT;
		}
	}
	/* costruisci la matrice delle dorsali */
	n_chains = 0;
	for (width = 0; width <= MAX_WIDTH; width++) {
		if (is_max[min_hgrey][width] && n_chains < MAX_N_CHAINS) {
			chain_map[min_hgrey][width] = n_chains++;
			chain_n[chain_map[min_hgrey][width]] = histo[min_hgrey][width];
		} else {
			chain_map[min_hgrey][width] = -1;
		}
	}
	for (hgrey = min_hgrey + 1; hgrey <= max_hgrey; hgrey++) {
		for (width = 0; width <= MAX_WIDTH; width++) {
			if (is_max[hgrey][width]) {
				if (chain_map[hgrey - 1][width] >= 0) {
					chain_map[hgrey][width] = chain_map[hgrey - 1][width];
					chain_n[chain_map[hgrey][width]] += histo[hgrey][width];
				} else if (chain_map[hgrey - 1][width - 1] >= 0 &&
						   chain_map[hgrey - 1][width - 1] != chain_map[hgrey][width - 2]) {
					chain_map[hgrey][width] = chain_map[hgrey - 1][width - 1];
					chain_n[chain_map[hgrey][width]] += histo[hgrey][width];
				} else if (chain_map[hgrey - 1][width + 1] >= 0) {
					chain_map[hgrey][width] = chain_map[hgrey - 1][width + 1];
					chain_n[chain_map[hgrey][width]] += histo[hgrey][width];
				} else if (n_chains < MAX_N_CHAINS) {
					chain_map[hgrey][width] = n_chains++;
					chain_n[chain_map[hgrey][width]] = histo[hgrey][width];
				} else {
					chain_map[hgrey][width] = -1;
				}
			} else {
				chain_map[hgrey][width] = -1;
			}
		}
	}
	if (!n_chains)
		return 0;

	/* scegli la dorsale principale */
	chain = 0;
	for (c = 1; c < n_chains; c++)
		if (chain_n[c] > chain_n[chain])
			chain = c;

	/* costruisci il vettore dei massimi */
	for (hgrey = min_hgrey; hgrey <= max_hgrey; hgrey++)
		for (width = 1; width < MAX_WIDTH; width++)
			if (chain_map[hgrey][width] == chain)
				goto chain_start_found;
chain_start_found:
	min_chain_hgrey = hgrey;
	peak[min_chain_hgrey] = width;
	for (hgrey = min_chain_hgrey + 1; hgrey <= max_hgrey; hgrey++) {
		if (chain_map[hgrey][width] == chain) {
		} else if (chain_map[hgrey][width + 1] == chain)
			width++;
		else if (chain_map[hgrey][width - 1] == chain)
			width--;
		else
			break;
		peak[hgrey] = width;
	}
	max_chain_hgrey = hgrey - 1;

	/* costruisci il vettore delle medie */
	for (hgrey = min_chain_hgrey; hgrey <= max_chain_hgrey; hgrey++) {
		width = peak[hgrey];
		n = histo[hgrey][width];
		nx = n * width;
		last = n;
		for (width = peak[hgrey] - 1; width > 0; width--) {
			if (histo[hgrey][width] > last)
				break;
			n += histo[hgrey][width];
			nx += histo[hgrey][width] * width;
			last = histo[hgrey][width];
		}
		last = histo[hgrey][peak[hgrey]];
		for (width = peak[hgrey] + 1; width <= MAX_WIDTH; width++) {
			if (width > peak[hgrey] * 5)
				break;
			if (histo[hgrey][width] > last)
				break;
			n += histo[hgrey][width];
			nx += histo[hgrey][width] * width;
			last = histo[hgrey][width];
		}
		mean_width[hgrey] = (float)nx / (float)n;
	}

	if (*threshold != 0) {
		*linewidth = mean_width[*threshold >> 1];
	} else {
		min_delta = 123456789.0;
		for (hgrey = min_chain_hgrey; hgrey <= max_chain_hgrey; hgrey++) {
			delta = mean_width[hgrey] - *linewidth;
			if (delta < 0)
				delta = -delta;
			if (delta < min_delta) {
				min_delta = delta;
				*threshold = hgrey << 1;
			}
		}
	}

	return 1;
}
#endif
/*---------------------------------------------------------------------------*/
#ifdef DAFARE
//capire che cosa fa !?? chi la usa ?
//sembra prendere una gr8, modificare il buffer e cambiare il tipo in BW
void thresh_image(const TRasterImageP &image, int threshold, int oversample_factor)
{
	UCHAR *buffer, *cur, *xnext, *ynext, *xynext;
	UCHAR *bigline, *out;
	UCHAR tmp;
	int lx, ly, wrap, outwrap_bytes, outwrap_pix;
	int x, y, bit;
	UINT thresh, thresh_2, thresh_4;

	TRasterGR8P ras8(image->getRaster());

	assert(ras8);
	thresh = threshold;
	thresh_2 = thresh * 2;
	thresh_4 = thresh * 4;
	buffer = ras8->getRawData();
	lx = ras8->getLx();
	ly = ras8->getLy();
	wrap = ras8->getWrap();
	outwrap_pix = lx * oversample_factor;
	outwrap_bytes = (outwrap_pix + 7) / 8;
	TMALLOC(bigline, outwrap_bytes * oversample_factor)

	switch (oversample_factor) {
		CASE 1 : for (y = 0; y < ly; y++)
		{
			cur = buffer + y * wrap;
			out = y == 0 ? bigline : buffer + y * outwrap_bytes;
			tmp = 0;
			bit = 7;
			for (x = 0; x < lx; x++, cur++) {
				if (*cur <= thresh)
					tmp |= 1 << bit;
				bit--;
				if (bit < 0) {
					*out++ = tmp;
					tmp = 0;
					bit = 7;
				}
			}
			if (bit != 7)
				*out++ = tmp;
			if (y == 0)
				memcpy(buffer, bigline, outwrap_bytes);
		}

		CASE 2 : for (y = 0; y < ly; y++)
		{
			cur = buffer + y * wrap;
			xnext = cur + 1;
			out = y == 0 ? bigline : buffer + y * 2 * outwrap_bytes;
			tmp = 0;
			bit = 7;
			for (x = 0; x < lx; x++, cur++, xnext += x < lx - 1) {
				if (*cur <= thresh)
					tmp |= 1 << bit;
				bit--;
				if (*cur + *xnext <= thresh_2)
					tmp |= 1 << bit;
				bit--;
				if (bit < 0) {
					*out++ = tmp;
					tmp = 0;
					bit = 7;
				}
			}
			if (bit != 7)
				*out++ = tmp;

			cur = buffer + y * wrap;
			xnext = cur + 1;
			ynext = y == ly - 1 ? cur : cur + wrap;
			xynext = ynext + 1;
			out = y == 0 ? bigline + outwrap_bytes : buffer + (y * 2 + 1) * outwrap_bytes;
			tmp = 0;
			bit = 7;
			for (x = 0; x < lx; x++, cur++, xnext += x < lx - 1, ynext++, xynext += x < lx - 1) {
				if (*cur + *ynext <= thresh_2)
					tmp |= 1 << bit;
				bit--;
				if (*cur + *xnext + *ynext + *xynext <= thresh_4)
					tmp |= 1 << bit;
				bit--;
				if (bit < 0) {
					*out++ = tmp;
					tmp = 0;
					bit = 7;
				}
			}
			if (bit != 7)
				*out++ = tmp;
			if (y == 0)
				memcpy(buffer, bigline, outwrap_bytes * 2);
		}

	DEFAULT:
		assert(0);
	}
	image->ras.type = RAS_WB;
	image->ras.lx = lx * oversample_factor;
	image->ras.ly = ly * oversample_factor;
	image->ras.wrap = outwrap_pix;
	image->x_dpi *= oversample_factor;
	image->y_dpi *= oversample_factor;
	image->sb_x0 *= oversample_factor;
	image->sb_y0 *= oversample_factor;
	image->sb_or_img_lx *= oversample_factor;
	image->sb_or_img_ly *= oversample_factor;

	TFREE(bigline)
}
#endif
/*===========================================================================*/
/*===========================================================================*/

void apply_lut(const TRasterImageP &image, UCHAR lut[256])
{
	int x, y, lx, ly, wrap;
	UCHAR *buffer, *pix;

	TRasterGR8P ras8(image->getRaster());
	assert(ras8);

	lx = ras8->getLx();
	ly = ras8->getLy();
	wrap = ras8->getWrap();
	ras8->lock();
	buffer = (UCHAR *)ras8->getRawData();

	for (y = 0; y < ly; y++) {
		pix = buffer + y * wrap;
		for (x = 0; x < lx; x++, pix++)
			*pix = lut[*pix];
	}
	ras8->unlock();
}