function extend, mask

mask2 = mask * 0

for i = 1, 510 do begin
	for j = 1, 510 do begin
		mask2[i, j] = mask[i - 1, j - 1] * mask[i - 1, j] * mask[i - 1, j + 1] * $
			mask[i, j - 1] * mask[i, j] * mask[i, j + 1] * $
			mask[i + 1, j - 1] * mask[i + 1, j] * mask[i + 1, j + 1]
	endfor
endfor

return, mask2
end


function fixrows, data, rgain, rmask

data2 = make_array(/float, (dim(data))[0], (dim(data))[1], value=0)
for i = 0, (dim(data))[1]-1 do begin
	if (rmask[i] and rgain[i] > 0) then data2[*, i] = data[*, i] / rgain[i]
endfor

return, data2

end


function rowmed, data, mask

rowm = make_array(/float, (dim(data))[1], value=0)

for i=0, (dim(data))[1]-1 do begin
	row = data[*, i]
	rmask = mask[*, i]
	good = where(rmask)
	if dim(good) gt 0 then rowm[i] = median(row[good])
endfor

return, rowm

end


pro rowdemo

dw
device, retain=2
data = float(readfits('images/data1136.a.fits'))
mask = data gt 0
data = data > 0
window, 0, xs=512, ys=512, title='Data'
tvscl, data
window, 1, xs=512, ys=512, title='Mask'
tvscl, mask
rm = rowmed(data, mask)
window, 2, title='Row medians'
plot, rm
rmf = fft(rm)
window, 3, title='FFT row medians'
plot, abs(rmf) < 100
stop, 'Data, rowmeds, fft rowmeds displayed. Type .c to continue.'

rmf[256-150:256+150] = 0
plot, abs(rmf) < 100
rm2 = real_part(fft(rmf, /inv)) > 0
window, 4, title='Filtered rowmeds'
plot, rm2
stop, 'Filtered rowmeds displayed. Type .c to continue.'

rmask = rm2 gt 1000
rg = make_array(/float, 512, value=0)
good = where(rmask)
rg[good] = rm[good] / rm2[good]
window, 3, title='Row gain'
plot, rg
data2 = fixrows(data, rg, rmask)
window, 5, xs=512, ys=512, title='Filtered data'
tvscl, data2
rm3 = rowmed(data2, mask)
window, 6, title='Filtered data rowmeds'
plot, rm3
stop, 'Row gains, filtered data, and filtered data rowmeds displayed. Type .c to continue.'

wset, 0
tvscl, rebin(data, 1024, 1024, /sam)
wset, 5
tvscl, rebin(data2, 1024, 1024, /sam)
stop, 'Closeups displayed.  Type .c to return.'

maskm = median(mask, 5)
mask = mask * maskm
wset, 1
tvscl, mask
stop, 'Mask1 displayed. Type .c to continue.'

filt = make_array(/float, 5, 3, value=1)
filt[*, 0] = -0.5
filt[*, 2] = -0.5
resp = convolve(data2, filt)
tvscl, abs(resp)
stop, 'Filter response displayed. Type .c to continue.'

mask2 = abs(resp - mean(resp)) lt 2.5 * stddev(resp)
tvscl, mask2
stop, 'Mask2 displayed. Type .c to continue.'

mask = mask * mask2
tvscl, mask
stop, 'Combined mask displayed. Type .c to continue.'

mask2 = mask
for i=0, 25 do begin
	mask2 = median(mask2, 5)
endfor
mask = mask * mask2
tvscl, mask
wset, 5
tvscl, data2 *  mask
writefits, 'images/mask1.fits', mask
writefits, 'images/data1.fits', data2 *  mask
stop, 'Mask and masked filtered data displayed. --All done--.  Type .c to return.'

end