; Copyright (c) 2011 A.P. Hitchcock All rights reserved
;+
;NAME:
;	REMOVE_LINES.PRO
;
;LAST CHANGED: ----------------------------------- 22-Apr-11
;
;PURPOSE:
;	This function identifies scan lines with values far from mean
; and replaces that line with either the average of adjacent 2 lines (if it is single)
; or the preceding line
; developed for image and stack cleanup in STXM5321 commissioning phase
;
;CATEGORY:
;	STAND ALONE: utilities
;
;CALLING SEQUENCE:
;	Result = REMOVE_LINES(2d_array,  silent=silent)
;
; INPUT
; 2d array (t.d in aXis2000 structure)
;
; KEYWORDS
;	silent			if set, do not display feedback
;
;CALLED FROM AXIS:
;	stack_analyze
;
;COMMON BLOCKS:
; axis_com
;
;MODIFICATION HISTORY:
; (22-Apr-11 aph) first version
;-
;
; *******************************************************************************
function remove_lines, image, silent=silent
@axis_com
on_Error,2

tmp = image.d

t = size(tmp)
num_row = t(2)
num_col = t(1)
; add vertical lines to make a 1-d array

spec = fltarr(num_row)
for j = 0, num_row-1 do begin
	for i = 0, num_col-1 do begin
		spec(j) = spec(j) + tmp(i,j)
	endfor
endfor
;plot, spec
;t=dialog_message('continue? ', /question)
;if t EQ 'No' then return, image

; differentiate
dspec = fltarr(num_row)
for i = 1, num_row-2 do dspec(i) = spec(i+1) - spec(i)
dspec(num_row-1) = spec(num_row-2)
dspec = shift(dspec,1)
dspec(0) = dspec(1)
; plot, dspec
; return, image

; evaluate 80% of min & % of max
min_dspec = 0.8*min(dspec)
max_dspec = 0.8*max(dspec)
print, 'lines below ',min_dspec, ' and above ', max_dspec, ' are considered bad'
; plot, dspec

; find min_index - index of lines that are below 0.2*min
min_index=where(dspec LT min_dspec, count_lo)
print, count_lo
; find index of lines above  0.2*max
; max_index=dspec(where(dspec GT max_dspec), count_hi)

; for each entry in min_index (min_ind(i)
if count_lo GT 0 then begin
; check value for the following entry. if it exceeds max_dspec and nothing for the subsequent value
; then replace min_ind(i) line in image with average of {min_ind(i)-1, min_ind(i)+1)
;if there is a value for min_ind(i)+2 then replace min_ind(i) line in image with min_ind(i)-1
; and the min_ind(i)+1) line with the min_ind(i)+2) line   (Handles 2 adjacent dropped lines)
	for i = 0, count_lo - 1 do begin
		if NOT keyword_set(silent)  then $
		    print, 'Processing bad line (min)', min_index(i),' sum of line derivative = ', dspec(min_index(i))
		if min_index(i) EQ 0 then begin
			tmp(*,0) = tmp(*,1)
			goto, continue_min
		endif
		if min_index(i) EQ num_row-2 then begin
			tmp(*,num_row-2) = tmp(*,num_row-3)
			goto, continue_min
		endif
		if min_index(i) EQ num_row-1 then begin
			tmp(*,num_row-1) = tmp(*,num_row-2)
		endif else begin
			if dspec(min_index(i)+1) GT max_dspec then begin	; MUST HAVE BOTH to define bad line
				if dspec(min_index(i)+2) LT max_dspec then begin
					tmp(*,min_index(i))= (tmp(*,min_index(i)-1) + tmp(*,min_index(i)+1))/2.
				endif else begin						; situation of 2 bad lines together
					tmp(*,min_index(i))  = tmp(*,min_index(i)-1)
	           		tmp(*,min_index(i)+1)= tmp(*,min_index(i)+2)
	           	endelse
	        endif
	    endelse
	    continue_min:
	endfor
endif

; regenerate spec, differentiate the modified array again and search for maxima
spec = fltarr(num_row)  ; reset to zero
for j = 0, num_row-1 do begin
	for i = 0, num_col-1 do begin
		spec(j) = spec(j) + tmp(i,j)
	endfor
endfor

for i = 0, num_row-2 do dspec(i) = spec(i+1) - spec(i)
dspec(num_row-1) = spec(num_row-2)
dspec = shift(dspec,1)
dspec(0) = dspec(1)
;return, dspec

; find index of lines above  0.2*max
 max_index=where(dspec GT max_dspec, count_hi)

; for each entry in max_index (max_ind(i)
if count_hi GT 0 then begin
; check the preceding value
; and nothing for 2 more then replace min_ind(i) line in image with average of {min_ind(i)-1, min_ind(i)+1)
;if there is a value for min_ind(i)+2 then replace min_ind(i) line in image with min_ind(i)-1
; and the min_ind(i)+1) line with the min_ind(i)+2) line   (Handles 2 adjacent dropped lines)
	for i = 0, count_hi - 1 do begin
		if NOT keyword_set(silent)  then $
		    print, 'Processing bad line (max)', max_index(i),' sum of line derivative = ', dspec(max_index(i))
		if max_index(i) EQ 0 then begin
			tmp(*,0) = tmp(*,1)
			goto, continue_max
		endif
		if max_index(i) EQ num_row-2 then begin
			tmp(*,num_row-2) = tmp(*,num_row-3)
			goto, continue_max
		endif
		if max_index(i) GE num_row-2 then begin
			tmp(*,max_index(i))= tmp(*,max_index(i)-1)		; handles case where bad line is last or last but one
		endif else begin
			if dspec(max_index(i)-1) LT min_dspec then begin	; MUST have up/down
				if dspec(max_index(i)-2) GT min_dspec then begin
					tmp(*,max_index(i))= (tmp(*,max_index(i)-1) + tmp(*,max_index(i)+1))/2.
				endif else begin					; situation of 2 bad lines together
					tmp(*,max_index(i))  = tmp(*,max_index(i)-1)
	           		tmp(*,max_index(i)+1)= tmp(*,max_index(i)+2)
	        	endelse
	        endif
	    endelse
	    continue_max:
	endfor
endif

if (count_lo + count_hi) GT 0 then begin
	image.dl = image.dl + ' RL'
	if NOT keyword_set(silent)  then $
   		axis_log, strtrim(string(count_lo + count_hi),2) +' lines replaced'
endif
; return original (if not bad lines)  or the corrected image
image.d = tmp
return, image
end
