include	<imhdr.h>
include	<foc.h>

# OVERFLOWX -- correct 8 bit overflow in FOC images as far as possible
# Two modes:
# give 2 images, and check ratio for each pixel. If ratio of overflowed
# to normal image exceeds low_thresh then we add 255 to 1st image, unless
# this would give too high a ratio. Negative ratios are interpreted as
# multiple of rms statistics ( sqrt(count) )
#
# If only one image is given then we proceed along each line and
# starting at a given pixel, and if we detect a downward shift by more than
# low_thresh we add 255, unless we would end up with a positive shift
# of greater than high_thresh.
#
# D. Giaretta, ?  Original.
# Phil Hodge,  12-Aug-1993  Call imunmap.
 
procedure t_overflowx()

char	input1[SZ_FNAME]		# overflowed image
char	input2[SZ_FNAME]		# non-overflowed image
char	output[SZ_FNAME]		# output image
real	ratio				# expected ratio of counts
real	low_thresh			# lower threshold for acceptable ratio
real	high_thresh			# upper acceptable measured ratio
int	start_pix			# starting pixel to process single lines
#--
pointer	x_foc_immap()
pointer	in1, in2, out, foc_in1, foc_in2, foc_out
int	nowhite(), clgeti(), junk
real	clgetr()

begin

	call clgstr ( "input1", input1, SZ_FNAME)
	if ( nowhite( input1, input1, SZ_FNAME) == 0 )
	    call error ( 0, " must supply input image")
	call salloc( foc_in1, SZ_FOC_COORD, TY_STRUCT)
	in1  = x_foc_immap( input1, foc_in1, READ_ONLY, 0)

	call clgstr ( "input2", input2, SZ_FNAME)

	call clgstr ( "output", output, SZ_FNAME)
	call salloc( foc_out, SZ_FOC_COORD, TY_STRUCT)
	out = x_foc_immap( output, foc_out, NEW_COPY, foc_in1)
	low_thresh  = clgetr( "lowthresh")
	high_thresh = clgetr( "highthresh")
	start_pix = clgeti( "startpix")

	
	junk = nowhite( input2, input2, SZ_FNAME)
	if ( input2[1] == EOS ) { 
	    call single_im( in1, out, start_pix, low_thresh, high_thresh)
	} else {
	    ratio = clgetr( "ratio" )
	    call salloc( foc_in2, SZ_FOC_COORD, TY_STRUCT)
	    in2  = x_foc_immap( input2, foc_in2, READ_ONLY, foc_in1 )
	    call pair_im ( in1, in2, out, low_thresh, high_thresh, 
					start_pix, ratio)
	    call imunmap (in2)
	}
	call imunmap (out)
	call imunmap (in1)
end
# SINGLE_IM -- single image for overflow correction

procedure single_im( in, out, start_pix, low_thresh, high_thresh)

pointer	in		# i: overflowed image
pointer	out		# i: output image
int	start_pix	# i: start pix in line
real	low_thresh	# i: lower thresh for acceptable ratio
real	high_thresh	# i: upper acceptable measured ratio
#--

pointer	i_buf, o_buf

long	i_v[IM_MAXDIM]
long	o_v[IM_MAXDIM]
real	low, high		# threshold values actually used
int	imgnlr(), impnlr()
int	naxis1, pix

begin

	call amovkl( long(1), i_v, IM_MAXDIM)
	call amovkl( long(1), o_v, IM_MAXDIM)

	naxis1 = IM_LEN( in, 1)
	low  = low_thresh
	high = high_thresh

	while ( imgnlr( in,  i_buf, i_v ) != EOF &&
		impnlr( out, o_buf, o_v ) != EOF   ) {

	    call amovr( Memr[i_buf], Memr[o_buf], naxis1 )
	    do pix = start_pix-1, naxis1-2 {

	        if ( low_thresh  < 0 ) 
		    low  = abs(low_thresh)*sqrt( max(0.0,Memr[ o_buf + pix]) )
		if ( high_thresh < 0 )
		    high = abs(high_thresh)*sqrt(max(0.0, Memr[ o_buf + pix]) )

		if (  (Memr[ o_buf + pix + 1] < Memr[ o_buf + pix] - 
						low ) &&
		      (Memr[ o_buf + pix + 1] + 255 < Memr[ o_buf + pix] +
						high ) ) 
		    Memr[ o_buf + pix + 1] = Memr[ o_buf + pix + 1] + 255

	    }
	}

end
# PAIR_IM -- correct for overflow given two images, the first overflowed
# the second not overflowed
procedure pair_im ( in1, in2, out, low_thresh, high_thresh, start_pix, ratio)

pointer	in1		# i: overflowed image
pointer	in2		# i: non-overflowed image
pointer	out		# i: output image
real	low_thresh	# i: lower threshold for acceptable ratio
real	high_thresh	# i: upper acceptable measured ratio
int	start_pix	# i: starting pixel to process single lines
real	ratio		# i: expected ratio of counts
#--

pointer	i1_buf, i2_buf, o_buf

long	i1_v[IM_MAXDIM]
long	i2_v[IM_MAXDIM]
long	o_v[IM_MAXDIM]
real	low, high		# threshold values actually used
int	imgnlr(), impnlr()
int	naxis1, pix

begin

	call amovkl( long(1), i1_v, IM_MAXDIM)
	call amovkl( long(1), i2_v, IM_MAXDIM)
	call amovkl( long(1), o_v,  IM_MAXDIM)

	naxis1 = IM_LEN( in1, 1)
	low  = low_thresh
	high = high_thresh

	while ( imgnlr( in1, i1_buf, i1_v ) != EOF &&
	        imgnlr( in2, i2_buf, i2_v ) != EOF &&
		impnlr( out, o_buf,   o_v ) != EOF   ) {

	    call amovr( Memr[i1_buf], Memr[o_buf], naxis1 )
	    do pix = start_pix-1, naxis1-1 {

	        if ( low_thresh  < 0 ) 
		    low  = abs(low_thresh)*sqrt( max(0.0, Memr[i2_buf + pix]) )
		if ( high_thresh < 0 )
		    high = abs(high_thresh)*sqrt(max(0.0, Memr[i2_buf + pix]) )

		if ( ratio*Memr[ o_buf + pix ] <
			     Memr[ i2_buf + pix ] - low &&
			 ratio*Memr[ o_buf + pix ] + 255 <
			     Memr[ i2_buf + pix ] + high    )
		    Memr[ o_buf + pix ] = Memr[ o_buf + pix ] + 255
	    }
	}

end