;+ ; ; Project : s3drs reconstruction software ; ; Name : (various test routines) ; ; Purpose : part of the s3drs test kit ; ; Explanation: (archived here for completeness, not intended for users use) ; ; Use : (self-documented, see notes in code) ; ; Written : Sandy Antunes, NRC, 2005-2009 ; ;- ; Given input data and noise, solve for best underlying density. ; Also accepts an initial starting image, or creates one if needed. ; ; default is even kernel redistribution to all larger voxels. ; using /gradient makes it redistribute only to larger voxels that ; are in the direction of the global center of mass. Requires the ; cmdelta be provided, typically cmdelta=( ([i,j,k]-masscenter) lt 0) ) ; FUNCTION simple3dkernel,voxel, ksize, factor,gradient=gradient,cmdelta=cmdelta ; optional 'cmdelta' is (current_kernel_center - global_center_of_mass) lt 0 ; i.e. True for x left of CM, y below CM, z below CM ; distribute horizonally ; distribute vertically ; distribute depth gradient=setyesno(gradient) if (gradient and n_elements(cmdelta) eq 0) then gradient=0 voxadd=voxel[ksize,ksize,ksize]*factor/(float(ksize)^3) answer=voxel*0.0; initial 'delta' cube of 0s for i=0,ksize*2 do for j=0,ksize*2 do for k=0,ksize*2 do begin if (voxel[ksize,ksize,ksize] lt voxel[i,j,k]) then begin ; redistribute density to neighbor ; note we optionally apply a 'gradient' towards CM voxaddadjusted=voxadd if (gradient) then begin ; compare cmdelta with range, example for x: ; cmdelta < 0 means pt is left of CM ; cmdelta > 0 means pt is right of CM ; i < ksize means region is to left of voxel center ; i > ksize means region is to right of voxel center ; so cmdelta=true and i > ksize = good, also ; cmdelta=false and i < ksize = good tilt= [(i gt ksize),(j gt ksize),(k gt ksize)] if (total(cmdelta ne tilt) ne 3) then voxaddadjusted=0 end answer[i,j,k] += voxaddadjusted answer[ksize,ksize,ksize] -= voxaddadjusted end end return,answer END ; Experimental settings: ; /gradient uses a kernel towards gradient rather than flat ; /renorm re-normalizes the solution after each solution ; FUNCTION simple3dsolve,pxwrapper,datum,noise,cube,normalize=normalize,$ ksize=ksize,gradient=gradient,renorm=renorm if (n_elements(ksize) eq 0) then ksize=1; kernel size if (n_elements(cube) eq 0) then begin cube=simple3drv(datum,pxwrapper=pxwrapper,/usemax) normalize=1; force appropriate data values as simple3drv is always 0-1 end normalize=setyesno(normalize) gradient=setyesno(gradient) renorm=setyesno(renorm) if (renorm) then print,'Will renorm each step' if (normalize) then begin ; normalize our input solution as a first guess cube=cube*max(datum)/max(pr_render(pxwrapper,cube)) end ; get our initial criteria render=pr_render(pxwrapper,cube) ; chisq=TOTAL(ratio(abs(datum-render),noise)) chisq=MEAN(ratio(abs(datum-render),noise)) masscenter=centroid3d(cube) print,'Initial CM is ',masscenter N=pxwrapper.sr.ncube numloops=0 numbetter=0 ; NOW ITERATE TO SOLVE ; try gross sharpening then reduce amount for later 'better' solutions for ifactor=5,1,-2 do begin factor=ifactor/100.0; convert to a percent print,'Sharpening at ',ifactor,'% (',factor,')' ; initially do each pixel, then start skipping for sparser focusing in improved=1 while (improved) do begin ; *** Solver newcube=cube ; make a temporary storage for our guesses ; FOR NOW, I just skip the 'ksize' border region, bad but easy :) ++numloops ; now go voxel-by-voxel for i=ksize, N-1-ksize, 1 do $ for j=ksize, N-1-ksize, 1 do $ for k=ksize, N-1-ksize, 1 do $ newcube[i-ksize:i+ksize,j-ksize:j+ksize,k-ksize:k+ksize] += $ simple3dkernel($ cube[i-ksize:i+ksize,j-ksize:j+ksize,k-ksize:k+ksize],$ ksize,factor,gradient=gradient,$ cmdelta=( ([i,j,k]-masscenter) lt 0) ) ; *** Check render=pr_render(pxwrapper,newcube) if (renorm) then begin newcube=newcube*max(datum)/max(render) render=pr_render(pxwrapper,newcube) end ; newchisq=TOTAL(ratio(abs(datum-render),noise)) newchisq=MEAN(ratio(abs(datum-render),noise)) print,'A',ksize,newchisq,chisq,max(render),max(datum),max(noise) if (newchisq le chisq) then begin cube=newcube chisq=newchisq ++numbetter end else begin improved=0 ; time to end this loop and try a lower amt end end ; end of 'do while improved' end ; end of ifactor loop print,'Out of ',numloops,' iterations, ',numbetter,' improved chisq' print,'Final CM is ',centroid3d(cube) return,cube END ; compare this against ; quickcme,64,'20071231_012220','back',/checksn,/nozero,sncutoff=1,$ ; /occultright,/occulttop pro test_simple3dsolve,ksize=ksize,gradient=gradient,renorm=renorm,N=N on_error,2 if (n_elements(N) eq 0) then N=128 fname1='20071231_025220_'+int2str(N)+'_inputdata.sav' fname2='20071231_012220_'+int2str(N)+'_inputdata.sav' if (file_exist(fname1)) then restore,/v,fname1 else $ if (file_exist(fname2)) then restore,/v,fname2 else return diff=diff[*,*,0:1]; ndiff=ndiff[*,*,0:1]; mset=mset[*,*,0:1]; pxwrapper.sr.Nd=2 sol=simple3dsolve(pxwrapper,diff*mset,ndiff,ksize=ksize,gradient=gradient,$ renorm=renorm) save,sol,file='sol.sav' renders=pr_render(pxwrapper,sol) sn=ratio(diff*mset,ndiff) lsetd=make_array(/string,pxwrapper.sr.Nd,value='Input') lsetn=make_array(/string,pxwrapper.sr.Nd,value='Noise') lsets=make_array(/string,pxwrapper.sr.Nd,value='S/N') lsetr=make_array(/string,pxwrapper.sr.Nd,value='Renders') dset=[[[renders]],[[sn]],[[ndiff]],[[diff*mset]]] lset=[lsetr,lsets,lsetn,lsetd] tv_multi,fixed=128,dset,/log,labels=lset,oom=2,line=pxwrapper.sr.Nd,$ /showrange,/unique,/cm threeview,fixed=128,sol,/log,/pixon,title='Solution Image Cube',oom=2,$ /cm,threed m=p_getmask(sr=pxwrapper.sr,iele=0,ratio=3.0) m3=[[[m]],[[m]],[[m]]] tv_multi,threed*m3,fixed=128,/raw,/log,oom=1,/local end