PRO imgcleanf,ima,jma,header=header
;
; ima is the input image, and jma is the output image
; to get ima use              ima=readfits('file.fit',header)
; to run this code use        imgcleanf,ima,jma,header=header
; to get the final file use   fits_write,'finalfile.fit',jma,header
; to get deleted pixels use   del =ima ne jma
;                             fits_write,'del.fit',del,header
; We added letter 'f' to 'imgclean' to show that the code is not the same as imgclean.
;
; Besides imgclean.pro, also starchck.pro is used by imgcleanf.
;
; If it is needed, you can add new key words in the header here.
; We do not change the header.
;
; (1)
; Dec. 6, 2005.
; This code restores the edges of Deep Impact MRI and HRI VIS images 
; after the work of 'imgclean' (written by E. Deutsch)
; which deletes CRs (and most of pixels at the edges of a DI VIS image).
; For images which are not made by Deep Impact, one may not need
; to restore the edges, and the edges can be different.
;
; (2)
; Jan. 18, 2006
; A new replacement of brightness of objects considered as CRs was made.
; It was needed, because near the border of a comet, 'imgclean'
; considers too many pixels as CRs and replaces their brightness with a wrong
; brightness, so part of comet becomes 'eaten' and a lot of 'dust' objects
; appear near the comet.
; Now replaced pixels get the brightness such as that of background.
; If this brightness is greater that their initial brightness,
; then the old brightness is left.
; It allows to make pixels deleted near the border to be 
; practically the same as the nearest background. 
; So the deleted pixels do not spoil the final image of a comet 
; in such a way as imgclean.
;
; The revised code 'imgcleanf' is recommended for use in all cases, 
; not only for images with a comet.
; We recommend 'imgcleanf' for automatic consideration of a large
; number of MRI or HRI images (RADREV or RAD, not raw images),
; because 'imgclean' is a reliable (it always normally comes to the end for 
; all DI images, in contrast e.g. to crfind) and fast code
; [it is not recommended for ITS and HRI IR images].
; For specific images, other codes can give better results.
; You can use the replacement suggested below for other codes (e.g. crfind).
;
; All corrections were made by Sergei Ipatov 
; (ipatov@astro.umd.edu)  
;
; (3)
; Jan 24, 2006 - Mark Desnoyer (md246@cornell.edu)
;
; Changed the way overcorrected pixels are replaced with the originals.
; it is more efficient now although the end result is the same.


imgmode=sxpar(header,'imgmode')
inst=sxpar(header, 'INSTRUME')
nx=sxpar(header,'naxis1')-1
ny=sxpar(header,'naxis2')-1  ; size of the image [0:ny]
jma=ima  ; jma will be changed by imgclean
jmae=ima ; we use jmae to restore the edges
;
imgclean,jma   ; jma is changed by imgclean ; running imgclean ===
;
 cr=jma ne ima  ; pixels considered as CRs
 sky=median(ima,7)  ; the number can be taken to be another 
; 

; if new brightness is greater than old brightness
; then we take an old brightness because CRs
; are brighter than background                    
; in order to try to get a better result
w = where(cr,cnt)
if cnt GT 0 then begin
	jma[w] = sky[w]
	tooBright = where(sky[w] GT ima[w],c2)
	if c2 GT 0 then jma[w[tooBright]] = ima[w[tooBright]]
endif
 
;
; if you do not need to store the edges, you can deleted 
; the text below (before 'return')
;
; We find the sizes of the edges based on imgmode.
; Notes: The number of restored raws and lines in the edges is
; greater than that in Table IV on page 71 in Space Science 
; Reviews, v. 117, Nos. 1-2, 2005 (e.g. 8 for mode=1),
; because e.g. for mode=1 9th row and 9th line are bright,
; 10th row and 10th line can have a lot of bright pixels,
; and 11th line can have some bright pixels which are not cosmic rays.
; The values of iy presented below were suggested 
; by Dennis Wellnitz. Dennis suggested to use ix=iy-1, 
; but for mode=3 I found that the 6th column also can contain a lot of 
; deleted pixels. So below the values of ix were taken to be
; equal to iy suggested by Dennis. Probably for some modes
; the values of ix can be the same as it was suggested to Dennis,
; but as there are some examples which give larger edges, I took ix=iy.
; I will be glad to get any comments about the sizes of the edges
; which must be restored in order to find true values.
; For what modes can we take ix=iy-1?
; Depending on one's problem, one can consider other values of ix and iy.
ix0=0	;;Number of rows to delete on the left
ix1=0	;;the right
iy0=0	;;the bottom
iy1=0	;;the top
if strtrim(inst,2) EQ 'HRIIR' then begin	;;IR instrument
	switch imgmode of
	1:
	5:begin
		ix0=3
		iy0=3
		iy1=1
		break
	end
	2:
	3:begin
		ix0=3
		iy0=2
		iy1=2
		break
	end
	4:
	6:begin
		ix0=5
		iy0=5
		iy1=1
		break
	end
	endswitch
	ix1=ix0

endif else begin	;;VIS instruments
	switch imgmode of
	1:
	9: BEGIN
		ix0=11 & break
	END
	2:
	3: BEGIN
		ix0=6 & break
	END
	5:
	6: BEGIN
		ix0=4 & break
	END
	7: BEGIN 
		ix0=2 & break
	END
	8: BEGIN
		ix0=1 & break
	END
	endswitch

	ix1=ix0
	iy0=ix0
	iy1=ix0
endelse

;
; now restore the edges
jmae(ix0:nx-ix1,iy0:ny-iy1) = $
	jma(ix0:nx-ix1,iy0:ny-iy1)
jma = jmae

return
end