; Procedure to fix cosmic ray hits in a single image

; input : 2-D input image
; output : 2-D output image
; sizex  : input & output x dimension size
; sizey  : input & output y dimension size

pro fixcrh,input,output,sizex,sizey

; calculate the local averages

thres = (9.0*smooth(input,3)-input)/8.0 

; calculate the local variances 

var = thres*0
for i=0,sizex-1 do for j=0,sizey-1 do begin
    count = 0
    for k=-1,1 do for l=-1,1 do begin
      if (k ne 0) and (l ne 0) and ((i+k) ge 0) and ((i+k) lt sizex) and $
	((j+l) ge 0) and ((j+l) lt sizey) then begin
	  count = count + 1
	  var(i,j) = var(i,j) + (input(i+k,j+l)-thres(i,j))^2
      endif
  	
      if count ne 0 then var(i,j) = var(i,j)/float(count) 
    endfor
endfor

var = sqrt(var)

;set threshold 

var = thres + 8*var
ind = where(var lt input,count)
output=input

; replace bad pixels

if (count) then output(ind) = thres(ind)
return
end


;beginning of main

name = ''
read,name,prompt = 'Enter the name of the fits file to fix > '

input=readfits(name)

name = ''
read,name,prompt = 'Enter the output file name > '

sd = size(input)

print,sd

tvscl, sqrt(abs(input))

output = input*0
fixcrh,input,output,sd(1),sd(2)

tvscl, sqrt(abs(output))

print,'Writing Results to : ',name 
writefits,name,output

end