/*
 * castel.c		Jim Piper		14 July 1987
 *
 * (C)  M E D I C A L   R E S E A R C H   C O U N C I L   1987
 *
 * Routines to (i) compute castellation correction vs. grey value
 * tables, (ii) to modify these tables to perform simultaneous
 * contrast enhancement, (iii) to apply the tables to correct and
 * enhance an image.
 */
#include <stdio.h>
#include <wstruct.h>

/*
 * Compute decastellation tables t0 (even address pixels) and t1
 * (odd address pixels) for input object obj.  Also contrast stretch
 * so that least grey value is gbot, and top pc% of grey values become
 * gtop, with linear transformation between these ranges.
 *
 * The algorithm is as follows.  Suppose that an even-address pixel
 * has grey value ge and its neighbour (either side) has odd-address
 * (obviously) and grey value go.  Part of the difference (ge-go) is
 * of course due to image variation, but part may be due to a constant
 * castellation bias.  This bias appears to depend on the local image
 * intensity.
 *
 * So the algorithm attempts to compute the necessary correction to
 * be applied to odd-valued pixels to compensate for castellation,
 * indexed by the original grey value of the pixel; thus we require
 * a table Co(go).  To compute this, we find
 * the mean of the difference (go-ge) over all neighbour pairs
 * for which the odd-address pixel has grey value go.  It is assumed
 * that over the image as a whole, image variation will be smoothed
 * out leaving only the castellation factor.
 *
 * However, this is not sufficient because we are usually working
 * with THRESHOLDED images rather than RAW DIGITISATIONS.  Consider
 * the neighbours of an odd-address pixel with a grey value go
 * rather close to the threshold value.  Then consider its neighbours.
 * If they are less in value than go, they will mostly have been
 * removed by thresholding.  Thus the set of neighbours of odd-address
 * pixels with grey value go will be
 * biased in favour of grey values higher than go; the nearer go is
 * to the threshold value, the more pronounced the effect.  A similar
 * effect occurs at pixels of near maximum grey-value; they are at
 * the "top" of density hills and will inevitably tend to have
 * neighbours of lower grey-value.
 *
 * A solution is also to compute the castellation correction factor to
 * be applied to even-address pixels, Ce(ge).  So long as the castellation
 * correction does not vary rapidly with grey value, Co(g) should be equal
 * to -Ce(g) for any grey value g, but the bias effects will still be
 * in the same direction - towards higher value neighbours for low
 * grey values, towards lower value neighbours for high grey values.
 * Thus (Co(go)-Ce(go))/2 should be an almost unbiased estimate of the
 * correction required for odd-address pixels.
 *
 * We make the assumption that the correction function is smooth, and
 * so during its construction take care to fill in missing values or
 * values where there are only very few neighbour pairs by using a
 * hysteresis filter on the set of correction values.  The function
 * is then smoothed heavily with by a linear convolution.
 *
 * The castellation correction is compounded with the required contrast
 * stretch into two look-up tables, one for odd-address pixels, the
 * other for even-address pixels.
 */
compute_decastel_tables(obj,t0,t1,gbot,gtop,pc)
struct object *obj, **t0, **t1;
double pc;
{
	struct histogramdomain *hdom1, *hdom0;
	struct object *makehistogram();
	int obot,otoppc,area,areapc,orange,grange,h2[256],h3[256];
	register int v,w,d,k,*h0,*h1;
	struct iwspace iwsp;
	struct gwspace gwsp;
	register GREY *g;

	/*
	 * make the histograms, set up the domain pointers etc..
	 */
	*t0 = makehistogram(256);
	*t1 = makehistogram(256);
	hdom0 = (struct histogramdomain *)(*t0)->idom;
	hdom1 = (struct histogramdomain *)(*t1)->idom;
	h0 = hdom0->hv;
	h1 = hdom1->hv;
	for (k=0; k<256; k++)
		h2[k] = h3[k] = 0;
	/*
	 * scan object, making histogram of even-address pixel values h2[],
	 * of odd-address pixels values h3[], and use histogram structures to
	 * build the sums of even-odd pixel
	 * value differences addressed by the even-address pixel value h0[],
	 * and the sums of odd-even pixel value differences addressed by
	 * the odd-address pixel value h1[].  Every pair of adjacent pixels
	 * contributes to both of these latter two structures, but
	 * generally at different addresses.  The first two histograms
	 * are in fact used as pixel-difference COUNTS; therefore all
	 * pixels except those at interval end points are counted twice.
	 */
	initgreyscan(obj,&iwsp,&gwsp);
	while (nextgreyinterval(&iwsp) == 0) {
		k = iwsp.lftpos;
		g = gwsp.grintptr;
		/* is first pixel address odd ? */
		if ((k&1) && k<iwsp.rgtpos) {
			v = *g++;	/* odd */
			w = *g;		/* even */
			d = v - w;	/* odd-even */
			h0[w] -= d;
			h1[v] += d;
			++h2[w];
			++h3[v];
			k++;
		}
		/* now process in pairs, even-address first */
		for (; k < iwsp.rgtpos; k++) {
			w = *g++;	/* even */
			v = *g++;	/* odd */
			d = v - w;	/* odd-even */
			h0[w] -= d;
			h1[v] += d;
			++h2[w];
			++h3[v];
			if (++k < iwsp.rgtpos) {
				w = *g;		/* even */
				d = v - w;	/* odd-even */
				h0[w] -= d;
				h1[v] += d;
				++h2[w];
				++h3[v];
			}
		}
	}
#ifdef DEBUG
	elapsedtime();
	hprint("histogram of even-odd pixel differences",h0,stdout);
	hprint("histogram of odd-even pixel differences",h1,stdout);
	hprint("histogram of even-address pixels",h2,stdout);
	hprint("histogram of odd-address pixels",h3,stdout);
#endif DEBUG
	/*
	 * Scan the histograms to find (i) obot, the bottom
	 * even-address grey-value in the image, (ii) otoppc,
	 * the even-address grey value which is exceeded by
	 * pc% of even-address pixels.
	 */
	area = obot = 0;
	for (k=0; k<=255; k++) {
		if (area == 0)
			obot = k;
		area += h2[k];
	}
	area = (int) (area * pc / 100.0);
	areapc = 0;
	for (k=255; k>=0; k--) {
		areapc += h2[k];
		if (areapc > area) {
			otoppc = k;
			break;
		}
	}
#ifdef DEBUG
	fprintf(stdout,"OTOPPC %d\n\n",otoppc);
#endif DEBUG
	orange = otoppc-obot+1;
	grange = gtop-gbot+1;
	/*
	 * Scan the histograms, computing in h1[] the unbiased estimate of
	 * odd-even castellation (Co(g)-Ce(g))/2.  Take care when
	 * either difference-count histogram is empty; and fill h2[]
	 * with the sum of contributing differences
	 * so that we know how reliable is our estimate of h1[].
	 */
	for (k=0; k<=255; k++) {
		if (h2[k] != 0 && h3[k] != 0) {
			d = h1[k]*h2[k] - h0[k]*h3[k];
			w = h2[k]*h3[k];
			/* take care with rounding */
			if (d >= 0)
				d += w;
			else
				d -= w;
			w *= 2;		/* cannot include in next line (OS9 bug) */
			h1[k] = d / w;
			h2[k] += h3[k];
		}
		else if (h2[k] != 0) {
			/* take care with rounding */
			if (h0[k] >= 0)
				h0[k] += h2[k]/2;
			else
				h0[k] -= h2[k]/2;
			h1[k] = -(h0[k]/h2[k]);
		}
		else if (h3[k] != 0) {
			/* take care with rounding */
			if (h1[k] >= 0)
				h1[k] += h3[k]/2;
			else
				h1[k] -= h3[k]/2;
			h1[k] /= h3[k];
			h2[k] = h3[k];
		}
	}
#ifdef DEBUG
	hprint("histogram of mean odd-even pixel differences",h1,stdout);
#endif DEBUG
	/*
	 * Apply a hysteresis smooth for empty mean-difference
	 * locations, and then a linear smooth.
	 */
	for (k=1; k<=255; k++)
		if (h2[k] == 0)
			h1[k] = h1[k-1];
#ifdef DEBUG
	hprint("histogram of mean differences, after hysteresis smooth",h1,stdout);
#endif DEBUG
	histothreshsmooth(*t1);
#ifdef DEBUG
	hprint("histogram of mean differences, after linear smooth",h1,stdout);
#endif DEBUG
	/*
	 * Compute the transformation histograms.
	 * Note that we overwrite the work-space histograms,
	 * so take care when modifying this section.......
	 */
	for (k=0; k<=255; k++) {
		/*
		 * odd pixels : compute the castellation-corrected value
		 * by adding the mean odd-even difference, use this to
		 * determine the transformed value
		 */
		v = k - h1[k];
		if (v < obot)
			h1[k] = gbot;
		else if (v > otoppc)
			h1[k] = gtop;
		else
			h1[k] = gbot + grange*(v-obot)/orange;
		/* even pixels */
		if (k < obot)
			h0[k] = gbot;
		else if (k > otoppc)
			h0[k] = gtop;
		else
			h0[k] = gbot + grange*(k-obot)/orange;
	}
#ifdef DEBUG
	hprint("even-pixel value conversion table",h0,stdout);
	hprint("odd-pixel value conversion table",h1,stdout);
#endif DEBUG
}

#ifdef DEBUG
hprint(s,h,f)
char *s;
int *h;
FILE *f;
{
	int i;
	fprintf(f,"%s\n",s);
	for (i=1; i<=256; i++) {
		fprintf(f,"%5d",h[i-1]);
		if (i%16 == 0)
			fprintf(f,"\n");
	}
	fprintf(f,"\n");
}
#endif DEBUG

struct object *makehistogram(size)
{
	struct object *makemain();
	return(makemain(13,makehistodmn(1,size),NULL,NULL,NULL));
}


/*
 * apply decastelation/contrast enhancement tables t0 to
 * even-address pixels and t1 to odd-address pixels
 */
apply_decastel_table(obj,t0,t1)
struct object *obj,*t0,*t1;
{
	register GREY *g,*gtop;
	struct iwspace iwsp;
	struct gwspace gwsp;
	register int *h0,*h1;
	struct histogramdomain *hdom0,*hdom1;
	hdom0 = (struct histogramdomain *)t0->idom;
	hdom1 = (struct histogramdomain *)t1->idom;
	h0 = hdom0->hv;
	h1 = hdom1->hv;
	initgreyscan(obj,&iwsp,&gwsp);
	while(nextgreyinterval(&iwsp) == 0) {
		g = gwsp.grintptr;
		gtop = g + iwsp.colrmn;
		/* has first pixel an odd address ? */
		if (iwsp.lftpos & 1) {
			*g = h1[*g];
			g++;
		}
		/* process pixels in even-odd address pairs */
		while (g < gtop) {
			*g = h0[*g];
			g++;
			if (g < gtop) {
				*g = h1[*g];
				g++;
			}
		}
	}
}