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,*)))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 ;======================================================================