Viewing contents of file '../idllib/astron/contrib/varosi/code/allpro/de_gulch.pro'
function de_gulch, image_data, fit_coef, DEGREE=deg, INTERACTIVE=int, PLOT=plot
;+
; NAME:
;	de_gulch
; PURPOSE:
;	Gulch effect caused by intense source in images
;	will be removed by LLsq fitting of row thru source
;	to rows away from source and subtracting fit
;	to return flat fielded image.
; CALLING EXAMPLE:
;	image_degulched = de_gulch( image, fit_coef, /INTER )
; INPUTS:
;	image = 2D array of data
; KEYWORDS:
;	DEGREE = polynomial degree of fit, default = 1 for Linear fit.
;	/INTERACTIVE : allows specification of gulch fit percentage.
;	/PLOT : plots fit result against data.
; OUTPUTS:
;	fit_coef = coefficients of LLsq fit.
; RESULT:
;	Function returns the flat fielded image (gulch supressed).
; EXTERNAL CALLS:
;	function filter_image
;	function conv_vartype
;	function poly_fit
;	function poly
; PROCEDURE:
; MODIFICATION HISTORY:
;	Written, Frank Varosi NASA/GSFC 1991.
;-
	image = filter_image( image_data, /SMOOTH,/ALL )
	sim = size( image )
	xsiz = sim(1)
	ysiz = sim(2)
	if N_elements( deg ) NE 1 then deg=1

	maxim = max( image, imax )
	xpeak = imax MOD xsiz
	ypeak = imax / xsiz
	rowC = image(*,ypeak)


	if (ypeak LT ysiz/3) then begin
		y = ypeak + ysiz/2
		rowE = ( image(*,y-3) + image(*,y-4) )/2
	  endif else if (ypeak GT 2*ysiz/3) then begin
		y = ypeak - ysiz/2
		rowE = ( image(*,y+2) +  image(*,y+3) )/2
	   endif else begin
		rowT = ( image(*,ysiz-3) + image(*,ysiz-4) )/2
		rowB = ( image(*,2) +  image(*,3) )/2
		rowE = (rowT+rowB)/2
	    endelse

	coef = transpose( poly_fit( rowC, rowE, deg ) )

	if keyword_set( int ) then begin
		print,coef
		print," default recommended % for de-gulch =", -coef(1)*100
		text = ""
		read," enter new %  (return for default) :", text
		if (text NE "") then begin
			text = get_words( text )
			if N_elements( text ) GT 1 then coef = float( text ) $
				else coef(1)= -float( text )/100.
		   endif
	   endif

	image = image_data - poly( rowC, coef ) # replicate( 1, sim(2) )

	if keyword_set( plot ) then begin
		if N_elements( rowB ) GT 1 then begin
			plot,rowB
			oplot,rowT,LINE=2
		 endif else plot, rowE
		oplot,poly( rowC, coef ),LINE=1
	   endif

return, conv_vartype( image, TYPE=sim(sim(0)+1) )
end