; Copyright (c) 1998-2007 A.P. Hitchcock All rights reserved
;+
;NAME:
;	REMOVE_ZEROS.PRO
;
;LAST CHANGED: ----------------------------------- 18-Sep-07
;
;PURPOSE:
;	This function replaces all zeros in an image file
;    with the adjacent non-zero data
;   Optionally, it will remove up to 10 hot or cold spots
;     hot spots are points at least 5 times larger than the mean
;     cold spots are points less than 1/10th of the mean
;        (CAUTION - many images have zero spaces)
;
;CATEGORY:
;	STAND ALONE: utilities
;
;CALLING SEQUENCE:
;	Result = REMOVE_ZEROS(2d_array,  silent=silent, remove_hot=remove_hot, remove_cold = remove_cold)
;
; KEYWORDS
;	silent			if set, do not display feedback
;	remove_hot 		if set, carry out hot spot removal
;	remove_cold 	if set, carry out cold spot removal
;
;CALLED FROM AXIS:
;	stack_analyze
;
;COMMON BLOCKS:
; axis_com
;
;MODIFICATION HISTORY:
; (10-Apr-04 aph) first version
; (18-Sep-07 aph) hot and cold spot remover added
;-
; *******************************************************************************
function remove_zeros, image, silent=silent, remove_hot=remove_hot, remove_cold = remove_cold
@axis_com
on_Error,2
tmp = image

; ------------------------- replace zeroes with previous points
 count_all = 0
count = 1
test = 5
 while (count NE 0) do begin
	index_zero = long64(where(tmp EQ 0, count))
	IF count NE 0 THEN begin
;		if index_zero(0) GT 0 then print, index_zero
		for i = 0, n_elements(index_zero)-1 do begin
		  	if index_zero(i) GT 0 AND index_zero(i) LT n_elements(tmp)-1 then $
		  	    tmp(index_zero(i)) = 0.5*(tmp(index_zero(i)-1)+ tmp(index_zero(i)+1))
		  	if index_zero(i) EQ 0 then tmp(index_zero(i)) =  tmp(1)
		  	if index_zero(i) EQ n_elements(tmp)-1 then $
		  	          tmp(index_zero(i)) =  tmp(index_zero(i)-1)
		endfor
	endif
	count_all = count_all + count
	test = test - 1
	if test EQ 0 then count = 0
 endwhile
if NOT keyword_set(silent) AND count_all GT 0 then $
         axis_log, strtrim(string(count_all),2) +' zero removed'

; --------- processing to remove hot  spots -----------
if keyword_Set(remove_hot) then qtest = 'Yes' else qtest = 'No'
if NOT keyword_set(silent) then qtest=dialog_message('Remove hot spots ? ', /question)
if qtest EQ 'Yes' then begin
;	print, 'median = ', median(tmp)
	index_hot  = where(tmp GT 5*median(tmp), count2)
	IF count2 NE 0 THEN begin
		if NOT keyword_set(silent) then axis_log, 'found ' + strtrim(string(count2),2) + ' hot spots ...'
; sort these in order of size (difference from the median)  and take out first 10 only
		hot = tmp(index_hot)
;		for i = 0, n_elements(index_hot)-1 do begin
;			print, index_hot(i), tmp(index_hot)
;		endfor
		nhot = n_elements(hot)-1
		mhot=9
		t = [mhot, nhot]
		nhot = min(t)
		index_hot = intarr(nhot+1)
		for i = 0, nhot do begin
			index_hot(i) = where(tmp EQ max(hot))
			hot(where(hot EQ max(hot))) = 0.		; assumes unique valued !!
		endfor
		for i = 0, n_elements(index_hot)-1 do begin
		  	if index_hot(i) NE 0 then $
		  	    tmp(index_hot(i)) = 0.5*(tmp(index_hot(i)-1)+tmp(index_hot(i)+1)) else $
		  	    tmp(index_hot(i)) =  tmp(1)
		endfor
	endif
	if NOT keyword_set(silent) AND count2 GT 0 then $
    	axis_log, strtrim(string(count2),2) +' hot spot removed'
endif

; ---------  processing to remove cold spots -----------
if keyword_Set(remove_cold) then qtest = 'Yes' else qtest = 'No'
if NOT keyword_set(silent) then qtest=dialog_message('Remove cold spots ? ', /question)
if qtest EQ 'Yes' then begin
;	print, 'median = ', median(tmp)
	index_cold  = where(tmp LT 0.1*median(tmp), count2)
	IF count2 NE 0 THEN begin
		axis_log, 'found ' + strtrim(string(count2),2) + ' cold spots ...'
		if NOT keyword_set(silent) then axis_log, 'found ' + strtrim(string(count2),2) + ' cold spots ...'
; sort these in order of size (difference from the median)  and take out first 10 only
		cold = tmp(index_cold)
;		for i = 0, n_elements(index_cold)-1 do begin
;			print, index_cold(i), tmp(index_cold)
;		endfor
		ncold = n_elements(cold)-1
		mcold=9
		t = [mcold, ncold]
		ncold = min(t)
		index_cold = intarr(ncold+1)
		for i = 0, ncold do begin
			index_cold(i) = where(tmp EQ max(cold))
			cold(where(cold EQ max(cold))) = 0.		; assumes unique valued !!
		endfor
		for i = 0, n_elements(index_cold)-1 do begin
		  	if index_cold(i) NE 0 then $
		  	    tmp(index_cold(i)) = 0.5*(tmp(index_cold(i)-1)+tmp(index_cold(i)+1)) else $
		  	    tmp(index_cold(i)) =  tmp(1)
		endfor
	endif
	if NOT keyword_set(silent) AND count2 GT 0 then $
    	axis_log, strtrim(string(count2),2) +' cold spot removed'
endif

return, tmp
end
