; file: findrmscubefile.pro ; init: Jan 13 2018 Rob Rutten Deil from minmaxcubefile.pro ; last: Aug 13 2022 Rob Rutten Deil ;+ pro findrmscubefile,infile,maxrms,itbest,$ rmsfile=rmsfile,threshold=threshold,pedestal=pedestal,$ trange=trange,itskip=itskip,$ croptriangles=croptriangles,trimbox=trimbox,$ flatten=flatten,itdelay=itdelay,$ itthresholdfile=itthresholdfile,rmsrms=rmsrms,rmsplot=rmsplot,$ show=show ; get rms stuff from a cubefile or a segment of it using assoc ; cubefile may be fits or La Palma ; ; INPUT: ; infile: string with path/filename ; ; OUTPUTS: ; maxrms = maximum rms value ; itbest = it of maximum rms value ; ; OPTIONAL KEYWORD INPUTS: ; rmsfile = file to write rms (read with readcol) ; threshold: threshold to remove below (0.3 = worst 30% out) ; pedestal: increase zero level to this value (larger "rms" values) ; trange = [itstart,itend]: only sample this subrange ; itskip: skip samplings to speed up (default 0) ; croptriangles: crop flat borders and triangles (default 0) ; trimbox: partial cutout specification (default = none = -1) ; flatten: flatten extent, default 0 (may undo offlimb, umbral black) ; itdelay: print best it pairs with delay itdelay between them ; ; OPTIONAL KEYWORD OUTPUTs: ; itthresholdfile = file to write itthreshold into (read with readcol) ; rmsrms = rms of the non-zero (non-skipped) rms values ; rmsplot=2/1/0: 2 = make ps plot, 1 = show plot on screen (default 0) ; show = 1/0: show best image (cropped) ; ; HISTORY: ; Jan 14 2018 RR: start ; Apr 10 2018 RR: itthresholdfile for time-averaging sequenced ; May 5 2018 RR: itdelay ; Mar 15 2019 RR: finetuned after haphazard results, pedestal ; Jan 14 2020 RR: croptriangles ; Dec 27 2020 RR: ps plot ;- ; answer no-parameter query if (n_params() lt 3) then begin sp,findrmscubefile return endif ; keyword defaults if (n_elements(rmsfile) eq 0) then rmsfile='' if (n_elements(threshold) eq 0) then threshold=0 if (n_elements(pedestal) eq 0) then pedestal=0 if (n_elements(itthresholdfile) eq 0) then itthresholdfile='' if (n_elements(trange) eq 0) then trange=[0,-1] if (n_elements(itskip) eq 0) then itskip=0 if (n_elements(croptriangles) eq 0) then croptriangles=0 if (n_elements(trimbox) eq 0) then trimbox=-1 if (n_elements(flatten) eq 0) then flatten=0 if (n_elements(itdelay) eq 0) then itdelay=0 if (n_elements(rmsplot) eq 0) then rmsplot=0 if (n_elements(show) eq 0) then show=0 ; get filetype dummy=cgrootname(infile,extension=ext) ;RR needs coyotelib ; infile is a FITS file if (ext eq 'fits') then begin ; get cube geometry and file datatype from the fits header header_in=headfits_rr(infile) header_insize=(1+fix(n_elements(header_in)/36.))*2880 bitpix_in=fxpar(header_in,'bitpix') nx_in=fxpar(header_in,'naxis1') ny_in=fxpar(header_in,'naxis2') nt_in=fxpar(header_in,'naxis3') ; check if (nx_in eq 0 or ny_in eq 0 or nt_in eq 0) then begin print,' ###### no proper cubefile since no nx, ny, nt in header '+infile return endif ; set endian bigendian=1 endif ; infile is a La Palma file if (ext eq 'icube' or ext eq 'fcube' or ext eq 'bcube') then begin ; get cube geometry and file datatype from the LP header lp_readheader,infile,header=header,datatype=datatype_in,$ nx=nx_in,ny=ny_in,nt=nt_in,$ endian=endian_file ; check if (nx_in eq 0 or ny_in eq 0 or nt_in eq 0) then begin print,$ ' ##### findrmscubefile abort: no nx, ny, nt in header '+infile return endif ; set bitpix if (datatype_in eq 1) then bitpix_in=8 if (datatype_in eq 2) then bitpix_in=16 if (datatype_in eq 4) then bitpix_in=-32 ; set header size header_insize=512 ; set endian bigendian=0 endif ; open input file for assoc get_lun, unit_in if (bigendian) then openr,unit_in,infile,/swap_if_little_endian $ else openr,unit_in,infile if (bitpix_in eq -32) then $ inassoc=assoc(unit_in,fltarr(nx_in,ny_in),header_insize) if (bitpix_in eq 16) then $ inassoc=assoc(unit_in,intarr(nx_in,ny_in),header_insize) if (bitpix_in eq 8) then $ inassoc=assoc(unit_in,bytarr(nx_in,ny_in),header_insize) ; check input image size im0=inassoc[0] sizeim0=size(im0) nx_im0=sizeim0[1] ny_im0=sizeim0[2] if (nx_im0 ne nx_in or ny_im0 ne ny_in) then begin print, ' ##### nx or ny in header inequal to first image [nx,ny]' return endif ; trange itstart=trange[0] itend=trange[1] if (itend eq -1) then itend=nt_in-1 ;RR using sequence mean taken out because too slow? ;; ; get mean value selected cubefile part ;; meanval=0. ;; for it=itstart,itend do begin ;; im=inassoc[it] ;; if (trimbox[0] eq -1) then $ ;; reformimage,float(im),imcrop,flatten=flatten else $ ;; reformimage,float(im),imcrop,trimbox=trimbox,flatten=flatten ;; ;; imcrop=histo_opt(imcrop) ; makes no difference ;; ;; imcrop=imcrop>(0.8*avg(imcrop)<1.2*avg(imcrop)) ;; meanval=meanval+mean(imcrop) ;; endfor ;; meanval=meanval/(itend-itstart+1) ; find rms for all selected it's rms=fltarr(nt_in) if (itskip lt 1) then itskip=1 for it=itstart,itend,itskip do begin im=inassoc[it] if (trimbox[0] eq -1) then $ reformimage,float(im),imcrop,croptriangles=croptriangles,$ flatten=flatten else $ reformimage,float(im),imcrop,trimbox=trimbox,flatten=flatten imcrop=(imcrop>pedestal)-pedestal momim=moment(imcrop) rms[it]=sqrt(momim[1])/avg(imcrop) endfor CHECKAGAIN: ; find max and it its it maxrms=max(rms) itbest=where(rms eq maxrms) ; get rms of the rms rmsposind=where(rms gt 0) rmspos=rms[rmsposind] momrms=moment(rmspos) rmsrms=sqrt(momrms[1]) ; eliminate suspect spikes (only in full run with every sample) if (maxrms-avg(rms) gt rmsrms*5 and itskip lt 2) then begin rms[itbest]=0 print,' ----- findrmscubefile.pro undid suspect itbest ='+$ trimd(itbest)+'; check with showex,'+infile goto,CHECKAGAIN endif ; option: itdelay to get best/worst pairs with itdelay time difference mean2=(-1)*rms if (itdelay ne 0) then begin for it=itstart,itend-itdelay,itskip do $ mean2[it]=(rms[it]+rms[it+itdelay])/2. for ipair=0,10 do begin maxmean2=max(mean2) itfirst=where(mean2 eq maxmean2) print,' ----- best delay pairs ordered by mean: '+$ trimd(itfirst)+trimd(itfirst+itdelay)+$ trimd(rms[itfirst])+trimd(rms[itfirst+itdelay]) mean2[itfirst]=-1 endfor mean2=100*rms for it=itstart,itend-itdelay,itskip do $ mean2[it]=(rms[it]+rms[it+itdelay])/2. for ipair=0,10 do begin minmean2=min(mean2) itfirst=where(mean2 eq minmean2) print,' ----- worst delay pairs ordered by mean: '+$ trimd(itfirst)+trimd(itfirst+itdelay)+$ trimd(rms[itfirst])+trimd(rms[itfirst+itdelay]) mean2[itfirst]=100 endfor endif ; optional: plot on screen if (rmsplot eq 1) then begin window,/free plot,itstart+indgen(itend-itstart),rms[itstart:itend],/ynozero,$ xrange=[itstart,itend],xstyle=1,title=' selected it = '+trimd(itbest) plots,itbest,maxrms,psym=6,symsize=3 endif ; optional: ps plot in /tmp if (rmsplot eq 2) then begin psfilename='/tmp/rmsplot.ps' axrat=1.62 ; 1.62 ; golden ratio openpsplot,psfilename,thick=2,fontsize=7,xsize=8.8,ysize=8.8/axrat plot,itstart+indgen(itend-itstart),rms[itstart:itend],/ynozero,$ position=[0.2,0.2,0.95,0.9],$ ; margins xticklen=0.03,yticklen=0.03/axrat,$ ; same-length ticks xrange=[itstart,itend],xstyle=1,$ xtitle='time step',ytitle='rms',title=infile closepsplot,psfilename,opengv=1 endif ; optional: write rmsfile if (rmsfile ne '') then writecol,rmsfile,rms ; optional: write itthresholdfile if (itthresholdfile ne '') then begin rms_sel=histo_opt_rr(rms[itstart:itend],threshold,/bottom_only) ; 0.3 indeed cuts at 30% lowest if (rmsplot ne 0) then oplot,rms_sel,thick=3 rms_ind=where(rms gt min(rms_sel)) it_all=indgen(itend-itstart+1) it_sel=it_all[rms_ind] writecol,itthresholdfile,it_sel,fmt='(i5)' endif ; optional: show best image and optional flatten version if (show ne 0) then begin print,' ----- itbest ='+trimd(itbest)+' show best image, also flattened' im=inassoc[itbest] if (trimbox[0] eq -1) then $ reformimage,float(im),imcrop,croptriangles=croptriangles else $ reformimage,float(im),imcrop,trimbox=trimbox sv,imcrop if (flatten ne 0) then begin reformimage,imcrop,imcropflat,flatten=flatten sv,imcropflat endif endif ; free the input file free_lun,unit_in end ; =================== main for testing per IDLWAVE H-c ======================= ;; ; SST demo ;; path='/home/rutten/data/SST/2016-09-05-demo/' ;; infile=path+'crispex/wb.6563.09:48:31.corrected.aligned.icube' ;; ;; infile=path+'sst/ha_wb.fits' ;; ;; trimbox=[15,15,900,900] ; black border wb file upsets ;; croptriangles=1 ;; itskip=0 ;; show=1 ;; ; Gregal data ;; cd,'/media/rutten/SSTDATA/alldata/SST/2014-08-29-gregal' ;; infile='sst/crispex.6563.15:40:09.time_corrected.aligned_lc.icube' ;; infile='crispex/wb.6563.15:40:09.corrected.fits' ;; itskip=20 ;; ;; trange=[230,280] ;; ;; itskip=0 ;; ;; trimbox=[614,277,1000,595] ; disk segment works; with limb nothing works ;; ; DST/IBIS data from Tom Schad ;; cd,'/media/rutten/RRDATA/alldata/DST/2017-06-14-schad/' ;; infile='ibisA/ibis_wb_6563_20170614_151541.fits' ;; ; oops: itbest=152 is a spike! totally weird offset image ;; trimbox=[200,200,800,800] ; absolutely needed for bad circle aperture ;; ; rbecontrails data B = contrail pub data ;; cd,'~/data/SST/2014-06-21-quiet/sst_iris_sdo' ;; infile='wb.6563.08:02:31.corrected.aligned.fits' ;; rmsfile='rms.txt' ;; threshold=0.33 ; discard worst third; worse seeing than data A below ;; pedestal=345 ; above pore darkness, to get values roughly like data A below ;; itskip=0 ;; itdelay=21 ;; flatten=50 ;; itthresholdfile='rms_itsel.txt' ;; rmsplot=1 ;; show=1 ;; ; rbeconrails data A = original RBE pub data ;; cd,'~/data/SST/2008-06-15-RBEs/sst' ;; infile='wb_15Jun2008_strt.bcube' ;; rmsfile='rms.txt' ;; threshold=0.25 ; discard worst quarter ;; pedestal=0 ; bytecube, probably bytescaled in production, no zero ;; itskip=0 ;; itdelay=18 ;; itdelay=36 ;; flatten=50 ;; ;; trimbox=[100,100,800,800] ;; itthresholdfile='rms_itsel.txt' ;; rmsplot=1 ;; show=1 cd,'/home/rutten/data/ALMA/2018-04-12-ca/gong' infile='gongall.fits' trimbox=[1024-512,1024-512,1024+512,1024+512] rmsfile='gongallrms.dat' rmsplot=2 findrmscubefile,infile,maxrms,itbest,$ rmsfile=rmsfile,threshold=threshold,pedestal=pedestal,$ trange=trange,itskip=itskip,trimbox=trimbox,croptriangles=croptriangles,$ flatten=flatten,itdelay=itdelay,$ itthresholdfile=itthresholdfile,rmsrms=rmsrms,rmsplot=rmsplot,show=show ; print values print,' ----- maxrms ='+trimd(maxrms)+' itbest ='+trimd(itbest)+$ ' rmsrms ='+trimd(rmsrms) ; showex to check whether itbest is indeed the best showex,infile,itthis=itbest end