; fixbadpixel.pro, v1.1.1.
; Author: R Painter, BSc Project Student '10/'11, 30.03.11.
; Calling sequence: result=fixbadpixel(data, mask)

; Input mask supplied by the user must consist of pixel values of 1 and 0 only and must be of the same dimensions as the input data array.

; NB: This code has some flaws in that a very small minority of pixels are given very large values as a result of this code (e.g. 40,000 when the average value is <1,000). The 5-sigma mask in Steps 2 & 3 from v1.1 correct for this, however this method only works for dark frames and is not a feasible method to use on stellar images (where a 5 sigma cut-off will remove stellar pixels). Step 1 is therefore all that is used for stellar images. If time permits, it may be worth attempting to correct this routine so that the afore mentioned anomoly does not occur. W. Foxall, project parnter '10/'11 encountered the same issue in his "maskfix.pro" routine.


FUNCTION fixbadpixel,data,mask

data_fix=data
mask=mask

size_data_fix=size(data_fix)
size_mask=size(mask[where(mask eq 0)])

data_fix=float(data_fix)
data_fix[where(mask eq 0)]=!values.f_nan

width=size_data_fix[1]
height=size_data_fix[2]

rows=long(size_mask[3])
bad=long(size_mask[3])

a=long(0)


; Add single pixel 'border' around data_fix:

data_fixborder=fltarr(width+2,height+2)

For i=0,width-1 do begin &$
For j=0,height-1 do begin &$
data_fixborder[i+1,j+1]=data_fix[i,j] &$

Endfor &$
Endfor


; Let border values equal NaN

data_fixborder[0:width+1,0]=!values.f_nan
data_fixborder[0:width+1,height+1]=!values.f_nan
data_fixborder[0,0:height+1]=!values.f_nan
data_fixborder[width+1,0:height+1]=!values.f_nan


; Get bad x & y coordinates (add 1 when translating to info due to extra border dimension)

info=lonarr(3,rows)

For i=0,width-1 do begin &$
For j=0,height-1 do begin &$
value = mask[i,j] &$

If value eq 0 then begin &$
info(0,a) = i+1 &$
info(1,a) = j+1 &$
a=a+1 &$
Endif &$

Endfor &$
Endfor


; While bad pixels remain

While (bad ne 0) do begin &$


; Calculate number of free pixels ('availability') around bad pixels (CW from top-left of bad pixel)

For i=0L,rows-1 do begin &$

If (info[2,i] ne -1) then begin &$
x=info[0,i] &$
y=info[1,i] &$
free=0 &$

If (data_fixborder[x-1,y-1] gt 0) then begin &$
free=free+1 &$
Endif &$
If (data_fixborder[x,y-1] gt 0) then begin &$
free=free+1 &$
Endif &$
If (data_fixborder[x+1,y-1] gt 0) then begin &$
free=free+1 &$
Endif &$
If (data_fixborder[x+1,y] gt 0) then begin &$
free=free+1 &$
Endif &$
If (data_fixborder[x+1,y+1] gt 0) then begin &$
free=free+1 &$
Endif &$
If (data_fixborder[x,y+1] gt 0) then begin &$
free=free+1 &$
Endif &$
If (data_fixborder[x-1,y+1] gt 0) then begin &$
free=free+1 &$
Endif &$
If (data_fixborder[x-1,y] gt 0) then begin &$
free=free+1 &$
Endif &$

info[2,i]=free &$

Endif &$

Endfor &$


; Sort info by 3rd column (availability), decreasing

order=reverse(sort(info[2,*])) &$
info=info(*,order) &$


; For pixels with best availability, let bad pixel value equal mean of surrounding nonzero pixels

max=max(info[2,*]) &$

For i=0L,rows-1 do begin &$

If (max ne -1) then begin &$
While (info[2,i] eq max) do begin &$

x=info[0,i] &$
y=info[1,i] &$
a=[data_fixborder[x-1,y-1],data_fixborder[x,y-1],data_fixborder[x+1,y-1],data_fixborder[x+1,y],data_fixborder[x+1,y+1],data_fixborder[x,y+1],data_fixborder[x-1,y+1],data_fixborder[x-1,y]] &$
data_fixborder[x,y]=mean(a, /nan) &$
info[2,i]=-1 &$
bad=bad-1 &$

Endwhile &$
Endif &$
Endfor &$

Endwhile


; Transfer data_fix from 'result' array back to original array:

For i=0,width-1 do begin &$
For j=0,height-1 do begin &$
data_fix[i,j] = data_fixborder[i+1,j+1] &$

Endfor &$
Endfor


return,data_fix
end