function median_mad_iter_2d,x,niter=niter,weight=weight,nmad=nmad,plot=plot

;  Finds median and median absolute deviation for both I & V of all pixels bright in I (to exclude off-limb, maybe sunspot)
; x = 2-dim array of IV data = x(space,polarization)

if not keyword_set(niter) then niter=3
niter = fix(niter)>2
if not keyword_set(weight) then weight = 0.1
if not keyword_set(nmad) then nmad=2.

nx = n_elements(x(*,0))
tt = fltarr(4,niter)

xs = reform(float(x(*,0)),nx)
xs(1:nx-2) = weight*(xs(0:nx-3)+xs(2:nx-1)) + (1.-weight)*xs(1:nx-2)
ys = reform(float(x(*,1)),nx)
ys(1:nx-2) = weight*(ys(0:nx-3)+ys(2:nx-1)) + (1.-weight)*ys(1:nx-2)

tt(0,0) = median(xs)
tt(1,0) = median(abs(xs-tt(0,0)))
tt(2,0) = median(ys)
tt(3,0) = median(abs(ys-tt(2,0)))

if keyword_set(plot) then begin
	plot_hist,xs,xx,hh,bin=0.25*tt(1,0)
	hm = 0.3*max(hh)
	oplot,tt(0,0)+tt(1,0)*nmad*[-1,0,1],hm*[1,1,1],psym=-4
end	

for i=1,niter-1 do begin
	w = where(abs(xs-tt(0,i-1)) le nmad*tt(1,i-1))
	tt(0,i) = median(xs(w))
	tt(1,i) = median(abs(xs - tt(0,i)))
	tt(2,i) = median(ys(w))
	tt(3,i) = median(abs(ys(w) - tt(2,i)))
	if keyword_set(plot) then oplot,tt(0,i)+tt(1,i)*nmad*[-1,0,1],hm*(1.+i/10.)*[1,1,1],psym=-4	
end

return,reform(tt(*,niter-1))
end

Function fg_fix_shuttered_iv,x,thresh=thresh,xcen=xcen,entire=entire
; Fixes Shuttered IV images with CCD read anomalies on left hand side (LHS) of CCD, often seen in 2x2 or 4x4 summing
; x = IV 3-dimensional array
; thresh = threshold in DN for difference between V on LHS and RHS to apply correction
; xcen = location of last column of LHS of the IV array (in case offset ROI was used)
; entire = 1 causes it to process the entire image, not just the LHS of the CCD
; Returns new IV array with both I and V corrected (approximately) for the bad read
; 
;  Sample call:  fg_prep,file,-1,index,x  &  x = fg_fix_shuttered_iv(x,/entire)
;
;  Good test data may be found in /net/kokuten/archive/hinode/sot/level0/2007/11/09/FGIV/H2000/

; vers. 0.8 TDT  6-APR-2008

nx = n_elements(x(*,0,0))
ny = n_elements(x(0,*,0))
npol = n_elements(x(0,0,*))
if not keyword_set(xcen) then xcen = nx/2-1
lhs = fltarr(5)
rhs = lhs
for i=0,4 do begin
	xx = reform(x(xcen-1-i,3:ny-4,*),ny-6,npol)
	tt = median_mad_iter_2d(xx)
	lhs(i) = tt(2)
	xx = reform(x(xcen+2+i,3:ny-4,*),ny-6,npol)
	tt = median_mad_iter_2d(xx)
	rhs(i) = tt(2)
	end
	lhs = median(lhs)
	rhs = median(rhs)
if keyword_set(thresh) then if (abs(lhs-rhs) le thresh) then return,x
if (lhs lt rhs) then sgn = 1. else sgn = -1.
print,lhs,rhs,sgn
		
if keyword_set(entire) then ncols = nx-1 else ncols = xcen
mmx = fltarr(4,ncols+1)
for i=0,ncols do mmx(*,i) = median_mad_iter_2d(reform(x(i,3:ny-4,*),ny-6,npol))
for i=0,ncols do begin
	corr = mmx(2,i)*x(i,*,0)/(mmx(0,i)>10)
	x(i,*,1) = x(i,*,1) - corr
	x(i,*,0) = x(i,*,0) - sgn*corr
end
return,x
end