;
;+
;
; Sharpens SXT images using iterative (Jansen-Van Citters type)
;	algorithm.  Forces result to be > 0, conserves total
;	flux.

PRO sxt_sharpen, index, imagein, result, noise, gain=gain, cutoff=cutoff, maxcount=maxcount, dim=dim

; Input parameters:
;	index = Yohkoh index structure for image.  Used to get
;		image location, filter
;	imagein = SXT flux image, decompressed and dark subtracted.
;
; Output parameters:
;	result = sharpened image
;
; Optional Output Parameter:
;	noise = noise image, approximation to noise in input image
;
; Optional Input Keywords:
;	gain = gain factor for iteration.  Fraction of error that is
;		used to make next correction, 0. < gain < 1.. 
;		Default = 0.7
;	cutoff = threshold for iteration.  This is the ratio
;		 of the largest absolute difference to the largest
;		 image intensity.  Default = 0.02
;	maxcount = maximum number of iterations to compute.
;		 Default = 10
;	dim = dimension of PSF image to use.  This trades off
;		speed of comvolution with accuracy.  Default = 33.
;
; Procedure:
;	Convolve the image with the point-spread-function.
;	Difference of image and convolved image is a correction
;	to the image.  Iterate to stable result.  Forces result
;	to be positive and conserve total flux at each step.
;	Separation of noise and signal is based on the idea that
;	noise acts like spatial delta functions, while signal
;	is always convolved with PSF.  Noise estimation is best for
;	smooth images.
;
; History:
;	Written June 4, 1993  Barry LaBonte
;	New psf, better flux accounting, noise determination  April 5, 1994  BJL
;	Bug fixes, add /DIM keyword, progress reporting    May 23, 1994   BJL
;
;-
;

;-----------------------------------------------------------------
; Handle inputs
filter = GT_FILTB(index)
IF(N_ELEMENTS(filter) GT 1) THEN filter = filter(0)
pix = GT_CENTER(index)
IF(N_ELEMENTS(filter) GT 1) THEN BEGIN
	filter = filter(0)
	pix = pix(*,0)
ENDIF

; Gain factor for iteration
IF(KEYWORD_SET(gain) LE 0) THEN gain = 0.7
; Cutoff threshhold for iteration - largest difference as fraction
IF(KEYWORD_SET(cutoff) LE 0) THEN cutoff = 0.02
; Cutoff iteration count
IF(KEYWORD_SET(maxcount) LE 0) THEN maxcount = 10

IF( KEYWORD_SET(dim) LE 0) THEN dim = 33

; Noise upper limit
noisepass = 2.
; Noise spike cutoff, stardard deviations
noisecut = 8.

asize = SIZE(imagein)
dim1 = asize(1)
dim2 = asize(2)

; Get the PSF for this image location and filter
psf = SXT_PSF( filter, pix(0), pix(1), dim=dim )
psize = SIZE(psf)
p1 = psize(1)
p2 = psize(2)

;-----------------------------------------------------------------

; Flux to conserve
target = FLOAT(imagein)
totflux = TOTAL(target)

; Define circumference to do as a 1-dimensional array
circum = [LINDGEN(dim1), dim1*LINDGEN(dim2-2)+2*dim1-1, $
	 (dim1-1-LINDGEN(dim1)) + dim1*(dim2-1),	$
	 dim1*(dim2-3-LINDGEN(dim2-2))+dim1]

; Try to separate noise
m3 = MEDIAN(target, 3)
m3l = MEDIAN(target(circum), 3)
m3(circum) = m3l
noise = target - m3
tnoise = noise
; High pass filter twice to avoid structure
FOR i=0,1 DO BEGIN
	high = WHERE(ABS(tnoise) GT noisepass*STDEV(tnoise))
	tnoise(high) = 0.
ENDFOR
; Try to find spikes
t2noise = FLTARR(dim1,dim2)
pos = WHERE(m3 GT 0.)
t2noise(pos) = noise(pos)/m3(pos)
std = STDEV(t2noise,aver)
bad = WHERE( ABS((t2noise-aver)) GT noisecut*std, bcount )
IF( bcount GT 0 ) THEN tnoise(bad) = noise(bad)
; Force zero average
avnoise = REBIN(tnoise,1,1)
noise = tnoise - avnoise(0)
result = target - noise

; Use big image for convolution
current = FLTARR(dim1+2*p1, dim2+2*p2)
current(p1:p1+dim1-1,p2:p2+dim2-1) = result

; Fill out to match flux loss - avoid edge problems
PRINT, ' Doing 2 convolutions to match flux flow at the image edge. '
; This much flux is at the image curcumference
fluxcir = TOTAL(target(circum))
blur = CONVOL( current, psf )
temp = blur(p1:p1+dim1-1,p2:p2+dim2-1)
blur(p1:p1+dim1-1,p2:p2+dim2-1) = 0.
; This much flux is at the blurred circumference
fluxblur = TOTAL(temp(circum))
; This much flux flows out during convolution
fluxout = TOTAL(blur)
current = CONVOL( blur, psf )
; This much would flow back
fluxback = TOTAL(current(p1:p1+dim1-1,p2:p2+dim2-1))
backfrac = fluxback/fluxout
; We want the flow back to match the loss
fluxloss = fluxcir - fluxblur
current = blur * fluxloss/fluxback
current(p1:p1+dim1-1,p2:p2+dim2-1) = result

;--------------------------------------------------------------
criterion = 10. * cutoff
itercount = 0
PRINT, ' Beginning iteration loop. '
; Iteration loop
WHILE( criterion GT cutoff AND itercount LT maxcount) DO BEGIN
	blur = CONVOL( current, psf )
	bcur = blur(p1:p1+dim1-1,p2:p2+dim2-1)
	diff = target - noise - bcur

; Update the noise image
	oldnoise = noise
	m3 = MEDIAN(diff, 3)
	m3l = MEDIAN(diff(circum), 3)
	m3(circum) = m3l
	diffnoise = diff - m3
	FOR i=0,1 DO BEGIN
		high = WHERE( ABS(diffnoise) GT noisepass*STDEV(diffnoise) )
		diffnoise(high) = 0.
	ENDFOR
	m3s = MEDIAN(result, 3)
	pos = WHERE( m3s GT 0. )
	t2noise(*,*) = 0.
	t2noise(pos) = (result(pos) - m3s(pos))/m3s(pos)
	std = STDEV(t2noise, aver)
	bad = WHERE( ABS(t2noise-aver) GT noisecut*std, count )
	noise = noise + diffnoise * gain
	IF(count GT 0) THEN BEGIN
		noise(bad) = noise(bad) + gain * (result(bad) - m3s(bad))
	ENDIF
	avnoise = REBIN(noise,1,1)
	noise = noise - avnoise(0)
	diffnoise = noise - oldnoise

; Force result to be > 0, conserve flux
        result = result + (diff * gain - diffnoise)
        IF(MIN(result) LT 0.) THEN result = result - MIN(result)
        result = result * totflux / TOTAL(result)

; Test for convergence
	criterion = MAX(ABS(diff)) / MAX(target)
	PRINT, itercount, criterion
	itercount = itercount + 1

; Rematch the flux flow 
	fluxblur = TOTAL(bcur(circum))
	blur(p1:p1+dim1-1,p2:p2+dim2-1) = 0.
	fluxnow = TOTAL(blur)
	current = current + gain * (blur-current) * $
		 (fluxcir - fluxblur)/(backfrac * fluxnow)
	fluxcir = TOTAL(result(circum))
	current(p1:p1+dim1-1,p2:p2+dim2-1) = result

ENDWHILE

RETURN
END