PRO xrt_fourier_vacuum, image_in , image_out, $
clean_type=clean_type, nsigma=nsigma, $
nmed=nmed, debug=debug
; =========================================================================
;+
; PROJECT:
; Solar-B / XRT
;
; NAME:
; XRT_FOURIER_VACUUM
;
; CATEGORY:
; Data calibration
;
; PURPOSE:
; Remove Fourier signals introduced into the data by the
; read-out electronics.
;
; CALLING SEQUENCE:
; XRT_FOURIER_VACUUM, image_in, image_out
; [,clean_type=clean_type] [,nsigma=nsigma]
; [,nmed=nmed] [,/debug]
;
; INPUTS:
; IMAGE_IN - [Mandatory] (2-dim float array, [Nx,Ny])
; The image containing read-out signals.
;
; KEYWORDS:
; CLEAN_TYPE - [Optional] Type of Fourier cleaning to use:
; 0 = prefilter on semi-fixed streaks, then remove
; Fourier streaks, stars, & smudges [default]
; 1 = remove semi-fixed streaks only
; 2 = remove Fourier streaks, stars, & smudges
; 3 = don't Fourier clean!
; NSIGMA - [Optional] Number of standard deviations
; beyond which an FFT pixel is
; considered bad in "Fourier star" removal.
; Default is 4.5
; (Note: much below 4 not recommended)
; NMED - [Optional] Number of standard deviations
; above smoothed central minimum background
; to begin shielding (presumed) data from correction
; Default is 3.5
; (Note: outside range of 2.0 to 4.5 not recommended)
; /DEBUG - (Boolean) if set =1, give verbose diagnostic output.
;
; OUTPUTS:
; IMAGE_OUT - [Mandatory] (2-dim float array, [Nx,Ny])
; The image with the fit to read-out signals removed.
;
; EXAMPLE:
;
; Basic usage:
; IDL> xrt_fourier_vacuum, image_in, image_out
;
;
; COMMON BLOCKS:
; None.
;
; NOTES:
; 1) Procedure: In the default mode of operation (clean_type=0)
; first search for the four strongest semi-fixed Fourier "streaks"
; and suppress these, constraining both real and imaginary
; parts to be < local background plus a fraction of local rms.
; This makes the following search more effective, and does a better
; job at streak removal.
; Next conduct a ring search for local Fourier power concentrations
; ("stars", "smudges" and remaining "streaks") and suppress these.
; In both cases the lowest frequency portion of the transform
; and the FFT edges are shielded from correction due to the likely
; significant Fourier power contribution from the actual data.
; 2) Low frequency read-out signals are not corrected to avoid
; damaging the data; thus some larger scale read-out patterns are
; only partially corrected.
; 3) In some cases, typically at lower frequencies, the applied
; correction can slightly damage the data as well. Usually this
; takes the form of "sinc function ringing" (often in x) emanating
; from certain small brighter features in the data. To
; reduce/eliminate this, you can experiment decreasing nmed (to
; shield more of the data) or increasing nsigma (to correct fewer
; pixels). Running (in addition) with clean_type=2 can also help.
; Eventually the code may be improved to adjust this automatically.
; 4) You may see error messages such as
; % CURVEFIT: Failed to converge
; and/or
; % Program caused arithmetic error: Floating divide by 0
; % Program caused arithmetic error: Floating underflow
; % Program caused arithmetic error: Floating overflow
; % Program caused arithmetic error: Floating illegal operand
; These do not necessarily mean there is a problem; this can occur
; when a given searched-for Fourier feature is not found. (This
; will be improved in the future.)
; 5) If there are saturated pixels, this routine works better if
; saturation_fudge is run first (see xrt_clean_ro).
;
;
;
; MODIFICATION HISTORY:
;
progver = 'v2007.Apr.23' ; --- (SS, MW) Written by SS. Cleaned up by MW.
progver = 'v2007.May.10' ; --- (SS) Various bug fixes.
progver = 'v2007.May.18' ; --- (SS) Fixed minor bug in data shielding near
; --- the middle of the array; added note 5.
progver = 'v2007.Nov.28' ; --- (SS) Fixed bug causing failure when Nx>Ny
;
;
;-
; =========================================================================
if not keyword_set(clean_type) then clean_type=0
if not keyword_set(nsigma) then nsigma=4.5 ; threshold sigma to flag spike
if not keyword_set(nmed) then nmed=3.5 ; shield data > nsig*central_median
;===== predetermined semi-fixed streaks in Fourier space (from calib darks)
xlin=[211,263,316,870]
wlin=[1,4,1,1]
nl=n_elements(xlin)
;===== setup some constants (subject to modification)
cy=.25 ; outer fraction of image in y to exclude in streak searches
nsig=20. ; size of bins for sigma determination
damp=.25 ; frac increase over back for line to be processed
wsrch=3.6 ; size in line widths of search window on each side of line
csig=1. ; frac of sigma over background used in streak suppresion
cwid=1.0 ; size in line widths to suppress
nsmoo=11 ; # of pixels to smooth over (Fourier space) for 512x512
edg=5 ; number of y edge pixels to shield in 512x512 image
offs=0.01 ; amount to temporarily increase min array value above zero
;===== determine image size
siz=size(image_in)
nx=siz(1)
ny=siz(2)
;===== scale some of the constants appropriately
ns=round(nsmoo*max([nx,ny])/512.)
edg=round(edg*ny/512.)
;===== temporarily remove negative values
immin=min(image_in)
offset= (immin<0.) - offs*(immin lt 0.)
impos = image_in - offset
;===== forward transform
fy=fft(impos,-1)
fy0=fy
;===== construct x & y positional matricies
ymat=findgen(ny)
ymat=rotate(rebin(ymat,ny,nx),1)
xmat=findgen(nx)
xmat=rebin(xmat,nx,ny)
;===== find where FFT is > nmed*(median center backg.) or too close
;===== to edges; these areas will be shielded from correction as
;===== they likely contain real data
;===== median background in transform center
bckc=median(abs(fy0(nx*(.5-cy/4):nx*(.5+cy/4),ny*(.5-cy/4):ny*(.5+cy/4))))
;===== smooth centered transform, do edges properly
z=shift(smooth(shift(abs(fy0),nx/2,ny/2),nsmoo,/edge_truncate),-nx/2,-ny/2)
zc=smooth(abs(fy0),nsmoo,/edge_truncate)
z(nx/2-nsmoo:nx/2+nsmoo,*)=zc(nx/2-nsmoo:nx/2+nsmoo,*)
z(*,ny/2-nsmoo:ny/2+nsmoo)=zc(*,ny/2-nsmoo:ny/2+nsmoo)
;===== threshold > nmed*(median center backg.), shift to center,
;===== label contiguous regions, shift back
z = shift(z gt nmed*bckc,nx/2,ny/2)
iz=shift(label_region(z),-nx/2,-ny/2)
;===== ID region containing main Fourier power (ie, at 0,0)
iz0=iz(0)
;===== fix any region splits near nx/2: if either nx/2+1 or nx/2-2 are
;===== from region iz0 then set tweeners=iz0
izm=iz(nx/2-2,*)
izp=iz(nx/2+1,*)
iz0e=where(izm eq iz0 or izp eq iz0)
if iz0e(0) ne -1 then iz(nx/2-1:nx/2,iz0e)=iz0
;===== where either nx/2+1 or nx/2-2 is from region iz0 but the other isnt,
;===== set all iz of type <>iz0 equal to iz0
iz0not=where(((izp ne iz0) and (izm eq iz0)) or ((izp eq iz0) and (izm ne iz0)))
if iz0not(0) ne -1 then begin
typnot0=reform(izm(0,iz0not))
typnot0=[typnot0,reform(izp(0,iz0not))]
h0p=histogram(typnot0)
x0p=indgen(n_elements(h0p))+min(typnot0)
ifix=where(h0p ne 0 and x0p ne iz0,nfix)
if nfix ne 0 then begin
izfix=x0p(ifix)
for j=0,nfix-1 do iz(where(iz eq izfix(j)))=iz0
endif
endif
;===== find indicies of contiguous region containing freq.=0
;===== (main FFT power), plus edge areas
inofix=where((iz eq iz0) or (abs(ymat-ny/2) gt (ny/2-edg)) or $
(abs(xmat-nx/2) gt (nx/2-edg)))
;===== the first part of this program locates and suppresses (as needed)
;===== certain semi-fixed Fourier streaks. This allows the star, streak,
;===== smudge remover which follows (star_squasher)
;===== semi-fixed streak suppression begins here (clean_type = 0,1)
if clean_type le 1 then begin
;===== more heavily smoothed transform
fys0=shift(smooth(abs(shift(fy0,nx/2,ny/2)),ns,/edge_truncate),-nx/2,-ny/2)
;===== scale streak-finding parameters
xl0=xlin*nx/2176.
widp=round(wlin/2.*nx/512.)>1<8. ; peak width
widb=round(3*nx/2176.)>2 ; background extract width
dx0b=widp+(4*nx/2176.) ; background extract offset
amps=fltarr(nl)
bck0=amps
xcs=amps
wids=amps
xx=findgen(nx)
;===== loop over nl semi-fixed Fourier streaks
nl=n_elements(xlin)
for i=0,nl-1 do begin
;===== look in center of FFT, away from main "data" part of transform
yb=round(ny/2-cy*ny)
ye=round(ny/2+cy*ny)
;===== set search widths around each estimated streak_i location
xsrch=widp(i)*wsrch<12
xpb=round(xl0(i)-xsrch)
xpe=round(xl0(i)+xsrch)
fyli=avg(abs(fy(xpb:xpe,yb:ye)),1)
;===== typically, gaussfit to find location of streak_i
if i ne 1 then begin
ng=xpe-xpb+1
xg=findgen(ng)+xpb
g=gaussfit(xg,fyli,aa,nterms=4)
xcs(i)=aa(1)
xc=round(aa(1))
amps(i)=aa(0)
bck0(i)=aa(3)
wids(i)=aa(2)
if (aa(0) le 0) or (aa(1) le min(xg)) or (aa(1) ge max(xg)) or $
(aa(2) le 0) or (stdev(fyli) le stdev(fyli-g)) then skip=1 $
else skip=0
;===== second (i=1) streak is multicomponent, centroid it instead
endif else begin
slbase=(avg(fyli([0,1]))-avg(fyli([xpe-xpb-1,xpe-xpb])))/(xpb-xpe)
base=slbase*(xx(xpb:xpe)-xx(xpb)) + avg(fyli([0,1]))
xci=total((fyli-base)*xx(xpb:xpe))/total(fyli-base)
xc=round(xci)
xcs(i)=xci
bck0(i)=avg(base)
amps(i)=max(fyli-base)
if (amps(i) le 0) or (xc le xpb) or (xc ge xpe) then skip=1 $
else skip=0
end
;===== if gaussfit/centroid successful...
if skip ne 1 then begin
;===== use averages on either side of center to define background
fybck=(avg(abs(fy(xc-(dx0b(i)+widb):xc-dx0b(i),*)),0) + $
avg(abs(fy(xc+dx0b(i):xc+dx0b(i)+widb,*)),0))/2.
;===== use lightly smoothed (3 pix) median (7 pix) to define local background
mbck=median(fybck,7)
fyb0=smooth(mbck,3)
;===== select lower of local median background or a smoothed average one
fyb1=min([[fyb0],[reform(fys0(xc,*))]],dim=2)
;===== if streak/background > (1 + damp) continue with correction
fystr=median(abs(reform(fy(xc,*))),7)
if avg(fystr(cy*ny:(1-cy)*ny)/fyb1(cy*ny:(1-cy)*ny)) $
gt (1+damp) then begin
;===== calculate average sigma of background @ streak_i in nsig bins
sigb=fltarr(ny/nsig)
for j=0,ny/nsig-1 do begin
sigbj=moment(fyb1(nsig*j:nsig*(j+1)-1))
sigb(j)=sqrt(sigbj(1))
endfor
ysig=nsig/2+findgen(ny/nsig)*nsig
;===== cap endpoints to stabilize spline, and
;===== spline a continuous sigma_background vector
ysig=[0.,ysig,ny-1]
sigb=[sigb(0),sigb,sigb(ny/nsig-1)]
sigbk=spline(ysig,sigb,findgen(ny),3)
;===== set limiting streak column value at median(back)+sigma_back*csig,
;===== make 2D, rotate into y
fyb=rotate(rebin(fyb1+csig*sigbk,ny,nx),1)
;===== tweak x range of area to be extracted
xb=round(xc-widp(i)*cwid)
xe=round(xc+widp(i)*cwid)
;===== suppress real and imaginary parts to background level,
;===== retaining sign information
del=1d-16
sgnr=float(fy0(xb:xe,*))/(abs(float(fy0(xb:xe,*)))+del)
sgni=imaginary(fy0(xb:xe,*))/(abs(imaginary(fy0(xb:xe,*)))+del)
fyr=(abs(float(fy0(xb:xe,*)))<fyb(xb:xe,*))*sgnr
fyi=(abs(imaginary(fy0(xb:xe,*)))<fyb(xb:xe,*))*sgni
fy(xb:xe,*)=complex(fyr,fyi)
;===== do same mirror image across nyquist freq., using mirrored (reversed)
;===== background vector (note reversed imaginary sign and 1 pixel shift!)
fyr=shift(rotate(fyr,2),0,1)
fyi=shift(rotate(fyi,2),0,1)
fy(nx-xe:nx-xb,*)=complex(fyr,-fyi)
if keyword_set(debug) then begin
print,i,xc,alog10(avg(fyb(*,yb:ye)))
plot,alog10(avg(abs(fy0(*,yb:ye)),1)),xr=[xc-.05*nx,xc+.05*nx]
afy=alog10(avg(abs(fy(*,yb:ye)),1))
oplot,afy,lines=1
oplot,[1,1]*xc,[-4,1],lines=1
xbc=findgen(widb+1)
oplot,xbc+xc-(dx0b(i)+widb),afy(xbc+xc-(dx0b(i)+widb)),$
psym=1,syms=.5
oplot,xbc+xc+dx0b(i),afy(xbc+xc+dx0b(i)),psym=1,syms=.5
read,'ready?',ok
endif
endif
endif
endfor
;===== undo correction too close to y edges or data
fy(inofix) = fy0(inofix)
;==== if clean_type=0, do star/smudge correction on de-streaked FFT
if clean_type eq 0 then fyc=star_squasher(fy,ic3,i03,nsigma=nsigma,debug=0)
;===== otherwise, do star/smudge correction on raw transform (clean_type=2)
endif else fyc=star_squasher(fy0,ic1,i01,nsigma=nsigma,debug=0)
;===== undo corrections too close to data or edges
fy(inofix) = fy0(inofix)
if clean_type ne 1 then fyc(inofix) = fy0(inofix)
;===== transform back, and rescale to conserve average(image_in)
if clean_type ne 1 then image_out=float(fft(fyc,1)) else $
image_out=float(fft(fy,1))
image_out=image_out/total(image_out)*total(impos) + offset
image_out=image_out/total(image_out)*total(image_in)
if keyword_set(debug) then begin
mat=ymat*0.
mat(inofix)=1.
xtv,alog10(abs(fy0)>bckc*2)*(mat +1)
read,'ready?',ok
xtv,alog10(abs(fyc)>bckc*2)*(mat +2)
read,'ready?',ok
fixed=(abs(fyc) ne abs(fy0))
xtv,alog10(abs(fyc)>bckc*2)*(mat +1)*(fixed-.75)
print,'minmax input: ',minmax(image_in)
print,'minmax output: ',minmax(image_out)
mnmx=minmax(image_in-image_out)
print,'minmax difference: ',mnmx,' median diff: ', $
median(image_in-image_out)
read,'ready?',ok
xtv,(image_in-image_out)
print,'FINAL stop'
stop
endif
return
END ;======================================================================