source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Panel.ipf @ 798

Last change on this file since 798 was 798, checked in by srkline, 11 years ago

I'm not sure this is a great idea, but I'm putting the FFT / Debye sphere work that I have completed into SVN. It's all really rough - the math, I believe is correct, but the interface if really, really rough. But it's not going to develop without help.
File size: 16.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3// utility functions and procedures for displaying information
4// setting the matrix, and doing the calculations
5//
6
7
8
9///////////
10//
11// put the SLD multiplier on the panel somewhere - only 2 values are allowed - so use
12// radio buttons, or just display the value of the global, make the change at the command
13// line -- or keep the multiplier at 10^-7 and always use 2 digits
14//
15////////////
16
17//
18
19//***********
20//For a 3D slicer view of dmat:
21//
22//Duplicate/O dmat dmatView
23//
24//Open a New 3D voxelgram. Don't set any of the levels
25//Open the 3D Slicer under the Gizmo menu
26//Make the X, Y, and Z slices
27//(Yellow Hot seems to be the best)
28//
29//// for a better view
30//dmatView = log(dmat)
31//
32//// if it all goes blank after the log transform get rid of the INF
33//dmatView = numtype(dmatView)== 0 ? dmatView : 0
34//*************
35
36
37///PANEL
38///// panel procedures
39
40Proc Init_FFT()
41        DoWindow/F FFT_Panel
42        if(V_flag==0)
43                SetDataFolder root:
44                //init the globals
45                Variable/G root:FFT_N=128
46                Variable/G root:FFT_T=5
47                Variable/G root:FFT_Qmax = 0
48                Variable/G root:FFT_QmaxReal = 0
49                Variable/G root:FFT_DQ=0
50                Variable/G root:FFT_estTime = 0
51               
52                Variable/G root:FFT_SolventSLD = 0
53                Variable/G root:FFT_delRho = 1e-6                       //multiplier for SLD (other value is 1e-7)
54               
55                FFT_Qmax :=2*pi/FFT_T
56                FFT_QmaxReal := FFT_Qmax/2
57                FFT_DQ := pi/(FFT_N*FFT_T)
58                //empirical fit (cubic) of time vs N=50 to N=256
59                FFT_estTime := 0.56 - 0.0156*FFT_N + 0.000116*FFT_N^2 + 8e-7*FFT_N^3
60//              FFT_estTime := FFT_N/128
61               
62                FFT_Panel()
63        endif
64End
65
66Proc FFT_Panel()
67        PauseUpdate; Silent 1           // building window...
68        NewPanel /W=(1452,44,1768,531)/K=1 as "FFT_Panel"
69        DoWindow/C FFT_Panel
70        SetDrawLayer UserBack
71        DrawLine 5,68,311,68
72        DrawLine 5,142,311,142
73        DrawLine 5,250,311,250
74        SetVariable FFTSetVar_0,pos={7,7},size={150,15},title="Cells per edge (N)"
75        SetVariable FFTSetVar_0,limits={50,512,2},value= FFT_N,live= 1
76        SetVariable FFTSetVar_1,pos={7,33},size={150,15},title="Length per Cell (T)"
77        SetVariable FFTSetVar_1,limits={1,500,0.2},value= FFT_T,live= 1
78        SetVariable FFTSetVar_2,pos={183,26},size={120,15},title="Real Qmax"
79        SetVariable FFTSetVar_2,limits={0,0,0},value= FFT_QmaxReal,noedit= 1,live= 1
80        SetVariable FFTSetVar_3,pos={183,48},size={120,15},title="delta Q (A)"
81        SetVariable FFTSetVar_3,limits={0,0,0},value= FFT_DQ,noedit= 1,live= 1
82        Button FFTButton_0,pos={15,79},size={90,20},proc=FFT_MakeMatrixButtonProc,title="Make Matrix"
83        Button FFTButton_1,pos={14,157},size={90,20},proc=FFTMakeGizmoButtonProc,title="Make Gizmo"
84        Button FFTButton_2,pos={14,187},size={100,20},proc=FFTDrawSphereButtonProc,title="Draw Sphere"
85        Button FFTButton_3,pos={14,265},size={70,20},proc=DoTheFFT_ButtonProc,title="Do FFT"
86        Button FFTButton_4,pos={191,264},size={110,20},proc=FFT_PlotResultsButtonProc,title="Plot 1D Results"
87        Button FFTButton_5,pos={13,218},size={120,20},proc=FFTDrawZCylinderButtonProc,title="Draw Cylinder"
88        Button FFTButton_6,pos={134,79},size={90,20},proc=FFTEraseMatrixButtonProc,title="Erase Matrix"
89        Button FFTButton_7,pos={13,329},size={130,20},proc=FFT_BinnedSpheresButtonProc,title="Do Binned Debye"
90        Button FFTButton_8,pos={13,297},size={130,20},proc=FFT_AltiSpheresButtonProc,title="Do Debye Spheres"
91        Button FFTButton_14,pos={13,360},size={130,20},proc=FFT_BinnedSLDButtonProc,title="Do Binned SLD"
92        SetVariable FFTSetVar_4,pos={201,4},size={100,15},title="FFT time(s)"
93        SetVariable FFTSetVar_4,limits={0,0,0},value= FFT_estTime,noedit= 1,live= 1,format="%d"
94        Button FFTButton_9,pos={200,295},size={100,20},proc=FFT_Get2DSlice,title="Get Slice"
95        Button FFTButton_10,pos={169,156},size={130,20},proc=FFT_TransposeMat,title="Transpose Matrix"
96        Button FFTButton_11,pos={169,189},size={130,20},proc=FFT_RotateMat,title="Rotate Matrix"
97        Button FFTButton_12,pos={168,219},size={130,20},proc=FFT_AddRotatedObject,title="Add Rotated Obj"
98        Button FFTButton_13,pos={14,109},size={120,20},proc=FFTFillSolventMatrixProc,title="Solvent Matrix"
99        SetVariable FFTSetVar_5,pos={155,111},size={100,15},title="Solvent SLD"
100        SetVariable FFTSetVar_5,limits={-99,99,1},value= FFT_SolventSLD,live= 1
101        Button FFTButton_15,pos={209,327},size={90,20},proc=Interp2DSliceButton,title="Interp 2D"
102EndMacro
103
104Function Interp2DSliceButton(ba) : ButtonControl
105        STRUCT WMButtonAction &ba
106
107        switch( ba.eventCode )
108                case 2: // mouse up
109                        // click code here
110                        String folderStr=""
111                        Prompt folderStr, "Pick a 2D data folder"               // Set prompt for x param
112                        DoPrompt "Enter data folder", folderStr
113                        if (V_Flag || strlen(folderStr) == 0)
114                                return -1                                                               // User canceled, null string entered
115                        endif                   
116       
117                        Interpolate2DSliceToData(folderStr)
118                       
119                        //display the 2D plane
120                        DoWindow FFT_Interp2D_log
121                        if(V_flag == 0)
122                                Execute "Display2DInterpSlice_log()"
123                        endif
124                       
125                        break
126        endswitch
127
128        return 0
129End
130
131
132
133Function FFT_Get2DSlice(ctrlName) : ButtonControl
134        String ctrlName
135               
136        WAVE/Z dmat = root:dmat
137        if(WaveExists(dmat)==1)
138                get2DSlice(dmat)
139        endif
140       
141        //display the 2D plane
142        DoWindow FFT_Slice2D
143        if(V_flag == 0)
144                Execute "Display2DSlice()"
145        endif
146       
147        //display the 2D plane
148        DoWindow FFT_Slice2D_log
149        if(V_flag == 0)
150                Execute "Display2DSlice_log()"
151        endif
152       
153        //display the 1D binning of the 2D plane
154        DoWindow FFT_Slice1D
155        if(V_flag == 0)
156                Execute "Slice2_1D()"
157        endif
158       
159        return(0)       
160End
161
162Proc Display2DSlice()
163        PauseUpdate; Silent 1           // building window...
164        Display /W=(1038,44,1404,403)
165        DoWindow/C FFT_Slice2D
166        AppendImage/T detPlane
167        ModifyImage detPlane ctab= {*,*,YellowHot,0}
168        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
169        ModifyGraph mirror=2
170        ModifyGraph nticks=4
171        ModifyGraph minor=1
172        ModifyGraph fSize=9
173        ModifyGraph standoff=0
174        ModifyGraph tkLblRot(left)=90
175        ModifyGraph btLen=3
176        ModifyGraph tlOffset=-2
177        SetAxis/A/R left
178End
179
180Proc Display2DInterpSlice_log()
181        PauseUpdate; Silent 1           // building window...
182        Display /W=(1038,44,1404,403)
183        DoWindow/C FFT_Interp2D_log
184        AppendImage/T interp2DSlice_log
185        ModifyImage interp2DSlice_log ctab= {*,*,YellowHot,0}
186        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
187        ModifyGraph mirror=2
188        ModifyGraph nticks=4
189        ModifyGraph minor=1
190        ModifyGraph fSize=9
191        ModifyGraph standoff=0
192        ModifyGraph tkLblRot(left)=90
193        ModifyGraph btLen=3
194        ModifyGraph tlOffset=-2
195        SetAxis/A/R left
196End
197
198Proc Display2DSlice_log()
199        PauseUpdate; Silent 1           // building window...
200        Display /W=(1038,44,1404,403)
201        DoWindow/C FFT_Slice2D_log
202        AppendImage/T logP
203        ModifyImage logP ctab= {*,*,YellowHot,0}
204        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
205        ModifyGraph mirror=2
206        ModifyGraph nticks=4
207        ModifyGraph minor=1
208        ModifyGraph fSize=9
209        ModifyGraph standoff=0
210        ModifyGraph tkLblRot(left)=90
211        ModifyGraph btLen=3
212        ModifyGraph tlOffset=-2
213        SetAxis/A/R left
214End
215Proc Slice2_1D()
216        PauseUpdate; Silent 1           // building window...
217        Display /W=(1034,425,1406,763) iBin_2d vs qBin_2d
218        DoWindow/C FFT_Slice1D
219        ModifyGraph mode=4
220        ModifyGraph marker=19
221        ModifyGraph msize=2
222        ModifyGraph grid=1
223        ModifyGraph log=1
224        ModifyGraph mirror=2
225        Legend
226EndMacro
227
228
229Function FFT_TransposeMat(ctrlName) : ButtonControl
230        String ctrlName
231
232        Variable mode=1
233        Prompt mode,"Transform XYZ to:",popup,"XZY;ZXY;ZYX;YZX;YXZ;"
234        DoPrompt "Transform mode",mode
235        if (V_Flag)
236                return 0                                                                        // user canceled
237        endif
238       
239        fFFT_TransposeMat(mode)
240       
241        return (0)
242End
243
244// starting in XYZ
245// mode 1;2;3;4;5;
246//correspond to
247// "XZY;ZXY;ZYX;YZX;YXZ;"
248//
249Function fFFT_TransposeMat(mode)
250        Variable mode
251
252        WAVE/Z mat=mat
253        ImageTransform /G=(mode) transposeVol mat
254        WAVE M_VolumeTranspose=M_VolumeTranspose
255        Duplicate/O M_VolumeTranspose mat
256       
257        return(0)       
258End
259
260Function FFT_RotateMat(ctrlName) : ButtonControl
261        String ctrlName
262       
263        Variable degree=45,sense=1
264        Prompt degree, "Degrees of rotation around Z-axis:"
265//      Prompt sense, "Direction of rotation:",popup,"CW;CCW;"
266        DoPrompt "Enter paramters for rotation", degree
267       
268        if (V_Flag)
269                return 0                                                                        // user canceled
270        endif
271       
272        fFFT_RotateMat(degree)
273        return(0)
274End
275// note that the rotation is not perfect. if the rotation produces an
276// odd number of pixels, then the object will "walk" one pixel. Also, some
277// small artifacts are introduced, simply due to the voxelization of the object
278// as it is rotated. rotating a cylinder 10 then 350 shows a few extra "bumps" on
279// the surface, but the calculated scattering is virtually identical.
280//
281// these issues, may be correctable, if needed
282//
283Function fFFT_RotateMat(degree)
284        Variable degree
285       
286        Variable fill=0
287       
288        WAVE mat = mat
289        Variable dx,dy,dz,nx,ny,nz
290        dx = DimSize(mat,0)
291        dy = DimSize(mat,1)
292        dz = DimSize(mat,2)
293               
294        //? the /W and /C flags seem to be special cases for 90 deg rotations
295        ImageRotate /A=(degree)/E=(fill) mat            //mat is not overwritten, unless /O
296       
297        nx = DimSize(M_RotatedImage,0)
298        ny = DimSize(M_RotatedImage,1)
299        nz = DimSize(M_RotatedImage,2)
300//      Print "rotated dims = ",dx,dy,dz,nx,ny,nz
301        Variable delx,dely,odd=0
302        delx = nx-dx
303        dely = ny-dy
304       
305        if(mod(delx,2) != 0)            //sometimes the new x,y dims are odd!
306                odd = 1
307        endif
308//      delx = (delx + odd)/2
309//      dely = (dely + odd)/2
310        delx = trunc(delx/2)
311        dely = trunc(dely/2)
312
313        //+odd removes an extra point if there is one
314        Duplicate/O/R=[delx+odd,nx-delx-1][dely+odd,ny-dely-1][0,nz-1] M_RotatedImage mat
315
316// - not sure why the duplicate operation changes the start and delta of mat, but I
317// need to reset the start and delta to be sure that the display is correct, and that the scaling
318// is read correctly later
319        SetScale/P x 0,1,"", mat
320        SetScale/P y 0,1,"", mat
321       
322        nx = DimSize(mat,0)
323        ny = DimSize(mat,1)
324        nz = DimSize(mat,2)
325//      Print "redim = ",dx,dy,dz,nx,ny,nz
326       
327        KillWaves/Z M_RotatedImage
328       
329End
330
331// also look for ImageTransform operations that will allow translation of the objects.
332// then simple objects could be oriented @0,0,0, and translated to the correct position (and added)
333//
334Function FFT_AddRotatedObject(ctrlName) : ButtonControl
335        String ctrlName
336
337        Print "Not yet implemented"
338End
339
340
341Function FFT_MakeMatrixButtonProc(ctrlName) : ButtonControl
342        String ctrlName
343
344        NVAR nn=root:FFT_N
345        if(mod(nn, 2 ) ==1)             //force an even number for FFT
346                nn +=1
347        endif
348        MakeMatrix("mat",nn,nn,nn)
349End
350
351Function FFTMakeGizmoButtonProc(ctrlName) : ButtonControl
352        String ctrlName
353       
354        DoWindow/F Gizmo_VoxelMat
355        if(V_flag==0)
356                Execute "Gizmo_VoxelMat()"
357        endif
358End
359
360Function FFTDrawSphereButtonProc(ctrlName)  : ButtonControl
361        String ctrlName
362       
363        Execute "FFTDrawSphereProc()"
364End
365
366Proc FFTDrawSphereProc(matStr,rad,xc,yc,zc,fill)
367        String matStr="mat"
368        Variable rad=25,xc=50,yc=50,zc=50,fill=1
369        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
370        Prompt rad,"enter real radius (A)"
371        Prompt xc,"enter the X-center"
372        Prompt yc,"enter the Y-center"
373        Prompt zc,"enter the Z-center"
374        Prompt fill,"1= fill, 0 = erase"
375       
376        Variable grid=root:FFT_T
377       
378        FillSphereRadius($matStr,grid,rad,xc,yc,zc,fill)
379       
380End
381
382Function FFTDrawZCylinderButtonProc(ctrlName)  : ButtonControl
383        String ctrlName
384       
385        Execute "FFTDrawCylinder()"
386End
387
388Proc FFTDrawCylinder(direction,matStr,rad,len,xc,yc,zc,fill)
389        String direction
390        String matStr="mat"
391        Variable rad=25,len=300,xc=50,yc=50,zc=50,fill=1
392        Prompt direction, "Direction", popup "X;Y;Z;"
393        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
394        Prompt rad,"enter real radius (A)"
395        Prompt len,"enter length (A)"
396        Prompt xc,"enter the X-center"
397        Prompt yc,"enter the Y-center"
398        Prompt zc,"enter the Z-center"
399        Prompt fill,"1= fill, 0 = erase"
400       
401       
402        Variable grid=root:FFT_T
403       
404        if(cmpstr(direction,"X")==0)
405                FillXCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
406        endif
407        if(cmpstr(direction,"Y")==0)
408                FillYCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
409        endif
410        if(cmpstr(direction,"Z")==0)
411                FillZCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
412        endif
413       
414End
415
416
417Function DoTheFFT_ButtonProc(ctrlName) : ButtonControl
418        String ctrlName
419
420        Execute "DoFFT()"
421End
422
423Function FFT_PlotResultsButtonProc(ctrlName) : ButtonControl
424        String ctrlName
425
426        DoWindow/F FFT_IQ
427        if(V_flag==0)
428                Execute "FFT_IQ()"
429        Endif
430End
431
432
433Function FFTEraseMatrixButtonProc(ctrlName) : ButtonControl
434        String ctrlName
435       
436        Wave mat=root:mat
437        FastOp mat=0
438        return(0)
439End
440
441Function FFTFillSolventMatrixProc(ctrlName) : ButtonControl
442        String ctrlName
443       
444        Wave mat=root:mat
445        NVAR val=root:FFT_SolventSLD
446        mat=val
447        return(0)
448End
449
450Function FFT_AltiSpheresButtonProc(ctrlName) : ButtonControl
451        String ctrlName
452
453        Execute "DoSpheresCalcFFTPanel()"
454End
455
456Function FFT_BinnedSpheresButtonProc(ctrlName) : ButtonControl
457        String ctrlName
458
459        Execute "DoBinnedSpheresCalcFFTPanel()"
460End
461
462Function FFT_BinnedSLDButtonProc(ctrlName) : ButtonControl
463        String ctrlName
464
465        Execute "DoBinnedSLDCalcFFTPanel()"
466End
467
468
469Proc FFT_IQ() : Graph
470        PauseUpdate; Silent 1           // building window...
471        Display /W=(295,44,627,302) iBin vs qBin
472        DoWindow/C FFT_IQ
473        ModifyGraph mode=4
474        ModifyGraph marker=19
475        ModifyGraph msize=2
476        ModifyGraph gaps=0
477        ModifyGraph grid=1
478        ModifyGraph log=1
479        ModifyGraph mirror=2
480        Legend/N=text0/J "\\s(iBin) iBin"
481EndMacro
482
483
484
485
486/////UTILITIES
487
488// inverts the values in the matrix 0<->1
489Function InvertMatrixFill(mat)
490        Wave mat
491       
492        mat = (mat==1) ? 0 : 1
493End
494
495//overwrites any existing matrix
496// matrix is byte, to save space
497//forces the name to be "mat"
498//
499// switched to signed byte
500//
501Function MakeMatrix(nam,xd,yd,zd)
502        String nam
503        Variable xd,yd,zd
504       
505        nam="mat"
506        Print "Matrix has been created and named \"mat\""
507//      Make/O/U/B/N=(xd,yd,zd) $nam
508        Make/O/B/N=(xd,yd,zd) $nam
509End
510
511
512// calculate the average center-to-center distance between points
513// assumes evenly spaced on a cubic grid
514//
515Proc Center_to_Center(np)
516        Variable np = 100
517       
518        Variable Nedge = root:FFT_N
519        Variable Tscale = root:FFT_T
520       
521        Variable davg
522       
523        davg = (Nedge*Tscale)^3 / np
524        davg = davg^(1/3)
525        Print "Average separation (A) = ",davg
526       
527End
528
529// calculate the number of points required on a grid
530// to yield a given average center-to-center distance between points
531// assumes evenly spaced on a cubic grid
532//
533// davg and Tscale are the same units, Angstroms
534//
535Proc Davg_to_Np(davg)
536        Variable davg = 400
537       
538        Variable Nedge = root:FFT_N
539        Variable Tscale = root:FFT_T
540       
541        Variable np
542       
543        np = (Nedge*Tscale)^3 / davg^3
544        Print "Number of points required = ",np
545       
546End
547
548// The matrix is not necessarily 0|1, this reports the number of filled voxels
549// - needed to estimate the time required for the AltiVec_Spheres calculation
550Proc NumberOfPoints()
551       
552        Print "Number of points = ",NumberOfPoints_Occ(root:mat)
553        Print "Fraction occupied = ",VolumeFraction_Occ(root:mat)
554        Print "Overall Cube Edge [A] = ",root:FFT_T * root:FFT_N
555       
556End
557
558Function NumberOfPoints_Occ(m)
559        Wave m
560       
561        Variable num = NonZeroValues(m)
562
563        return(num)
564End
565
566
567Function VolumeFraction_Occ(m)
568        Wave m
569       
570        Variable num = NonZeroValues(m)
571        Variable dim = DimSize(m,0)
572
573        return(num/dim^3)
574End
575
576//it's a byte wave, so I can't use the trick of setting the zero values to NaN
577// since NaN can't be expressed as byte. So make it binary 0|1
578//
579// - the assumption here is that the defined solvent value is not part of the
580// particle volume.
581//
582Function NonZeroValues(m)
583        Wave m
584       
585        Variable num
586        NVAR val = root:FFT_SolventSLD
587       
588//      Variable t1=ticks,dim
589//      dim = DimSize(m, 0 )            //assume NxNxN
590        Duplicate/O m,mz
591       
592        MultiThread mz = (m[p][q] != val) ? 1 : 0
593       
594        WaveStats/Q/M=1 mz              // NaN and Inf are not reported in V_npnts
595       
596        num = V_npnts*V_avg
597       
598        KillWaves/Z mz
599//      Print "Number of points = ",num
600//      Print "Fraction occupied = ",num/V_npnts
601       
602//      Print "time = ",(ticks-t1)/60.15
603        return(num)
604End
605
606// returns estimate in seconds
607//
608//              updated for quad-core iMac, 2010
609//
610// type 3 = binned distances and SLD
611// type 2 = binned distances
612// type 1 is rather meaningless
613// type 0 = is the Debye, double sum (multithreaded)
614//
615//
616// - types 2 and 3, the binned methods are inherently AAO, so there is no significant
617// dependence on the number of q-values (unless it's  >> 1000). So values reported are
618// for 100 q-values.
619//
620// There is a function Timing_Method(type) in FFT_ConcentratedSpheres.ipf to automate some of this
621//
622Function EstimatedTime(nx,nq,type)
623        Variable nx,nq,type
624       
625        Variable est=0
626       
627        if(type==0)             // "full" XOP
628                est = 0.4 + 1.622e-5*nx+4.86e-7*nx^2
629                est /= 100
630                est *= nq
631        endif
632       
633        if(type==1)             //Igor function (much slower, 9x)
634                est = (nx^2)*0.0000517          //empirical values (seconds) for igor code, 100 q-values
635                est /= 100                                      // per q now...
636                est *= nq
637        endif
638       
639        if(type==2)             //binned distances, slow parts XOPed
640                est = 0.0680 + 2.94e-6*nx+4.51e-9*nx^2                  //with threading
641               
642//              est /= 100                                      // per q now...
643//              est *= nq
644        endif
645       
646        if(type==3)             //binned distances AND sld, slow parts XOPed
647                est = 0.576 + 4.22e-7*nx+1.76e-8*nx^2           //with threading
648               
649//              est /= 100                                      // per q now...
650//              est *= nq
651        endif
652       
653        return(est)
654End
Note: See TracBrowser for help on using the repository browser.