1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | |
---|
3 | // |
---|
4 | // utility functions and procedures for displaying information |
---|
5 | // setting the matrix, and doing the calculations |
---|
6 | // |
---|
7 | |
---|
8 | |
---|
9 | // |
---|
10 | // TO DO: |
---|
11 | // x- I need to change a lot of routines (most notably Gizmo) to take "10" as the default SLD |
---|
12 | // rather than "1" |
---|
13 | // |
---|
14 | // |
---|
15 | // -- incorporate utility function that reads the wave note from a binary file |
---|
16 | // Function/S LoadNoteFunc(PName,FName[,FileRef]) - use this as an "inspector?" |
---|
17 | //////////// |
---|
18 | |
---|
19 | // |
---|
20 | |
---|
21 | //*********** |
---|
22 | //For a 3D slicer view of dmat: |
---|
23 | // |
---|
24 | //Duplicate/O dmat dmatView |
---|
25 | // |
---|
26 | //Open a New 3D voxelgram. Don't set any of the levels |
---|
27 | //Open the 3D Slicer under the Gizmo menu |
---|
28 | //Make the X, Y, and Z slices |
---|
29 | //(Yellow Hot seems to be the best) |
---|
30 | // |
---|
31 | //// for a better view |
---|
32 | //dmatView = log(dmat) |
---|
33 | // |
---|
34 | //// if it all goes blank after the log transform get rid of the INF |
---|
35 | //dmatView = numtype(dmatView)== 0 ? dmatView : 0 |
---|
36 | // |
---|
37 | //************* |
---|
38 | |
---|
39 | |
---|
40 | ///PANEL |
---|
41 | ///// panel procedures |
---|
42 | |
---|
43 | Proc Init_FFT() |
---|
44 | DoWindow/F FFT_Panel |
---|
45 | if(V_flag==0) |
---|
46 | SetDataFolder root: |
---|
47 | //init the globals |
---|
48 | Variable/G root:FFT_N=128 |
---|
49 | Variable/G root:FFT_T=5 |
---|
50 | Variable/G root:FFT_Qmax = 0 |
---|
51 | Variable/G root:FFT_QmaxReal = 0 |
---|
52 | Variable/G root:FFT_DQ=0 |
---|
53 | Variable/G root:FFT_estTime = 0 |
---|
54 | |
---|
55 | Variable/G root:FFT_SolventSLD = 0 |
---|
56 | Variable/G root:FFT_delRho = 1e-7 //multiplier for SLD (other value is 1e-7) |
---|
57 | |
---|
58 | FFT_Qmax :=2*pi/FFT_T |
---|
59 | FFT_QmaxReal := FFT_Qmax/2 |
---|
60 | FFT_DQ := pi/(FFT_N*FFT_T) |
---|
61 | //empirical fit (cubic) of time vs N=50 to N=256 |
---|
62 | FFT_estTime := 0.56 - 0.0156*FFT_N + 0.000116*FFT_N^2 + 8e-7*FFT_N^3 |
---|
63 | // FFT_estTime := FFT_N/128 |
---|
64 | |
---|
65 | FFT_Panel() |
---|
66 | endif |
---|
67 | End |
---|
68 | |
---|
69 | |
---|
70 | Proc FFT_Panel() |
---|
71 | PauseUpdate; Silent 1 // building window... |
---|
72 | NewPanel /W=(1452,44,1768,531)/K=1 as "FFT_Panel" |
---|
73 | DoWindow/C FFT_Panel |
---|
74 | SetDrawLayer UserBack |
---|
75 | DrawLine 5,68,311,68 |
---|
76 | DrawLine 5,142,311,142 |
---|
77 | DrawLine 5,250,311,250 |
---|
78 | SetVariable FFTSetVar_0,pos={7,7},size={150,15},title="Cells per edge (N)" |
---|
79 | SetVariable FFTSetVar_0,limits={50,512,2},value= FFT_N,live= 1 |
---|
80 | SetVariable FFTSetVar_1,pos={7,33},size={150,15},title="Length per Cell (T)" |
---|
81 | SetVariable FFTSetVar_1,limits={1,500,0.2},value= FFT_T,live= 1 |
---|
82 | SetVariable FFTSetVar_2,pos={183,26},size={120,15},title="Real Qmax" |
---|
83 | SetVariable FFTSetVar_2,limits={0,0,0},value= FFT_QmaxReal,noedit= 1,live= 1 |
---|
84 | SetVariable FFTSetVar_3,pos={183,48},size={120,15},title="delta Q (A)" |
---|
85 | SetVariable FFTSetVar_3,limits={0,0,0},value= FFT_DQ,noedit= 1,live= 1 |
---|
86 | Button FFTButton_0,pos={15,79},size={90,20},proc=FFT_MakeMatrixButtonProc,title="Make Matrix" |
---|
87 | Button FFTButton_1,pos={14,157},size={90,20},proc=FFTMakeGizmoButtonProc,title="Make Gizmo" |
---|
88 | Button FFTButton_2,pos={14,187},size={100,20},proc=FFTDrawSphereButtonProc,title="Draw Sphere" |
---|
89 | Button FFTButton_3,pos={14,265},size={70,20},proc=DoTheFFT_ButtonProc,title="Do FFT" |
---|
90 | Button FFTButton_4,pos={180,264},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot FFT Results" |
---|
91 | Button FFTButton_5,pos={13,218},size={120,20},proc=FFTDrawZCylinderButtonProc,title="Draw Cylinder" |
---|
92 | // Button FFTButton_6,pos={134,79},size={90,20},proc=FFTEraseMatrixButtonProc,title="Erase Matrix" |
---|
93 | Button FFTButton_6a,pos={160,79},size={60,20},proc=FFTSaveMatrixButtonProc,title="Save" |
---|
94 | Button FFTButton_6b,pos={240,79},size={60,20},proc=FFTLoadMatrixButtonProc,title="Load" |
---|
95 | Button FFTButton_7,pos={13,329},size={130,20},proc=FFT_BinnedSpheresButtonProc,title="Do Binned Debye" |
---|
96 | Button FFTButton_7a,pos={180,329},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Binned Results" |
---|
97 | |
---|
98 | Button FFTButton_8,pos={13,297},size={130,20},proc=FFT_AltiSpheresButtonProc,title="Do Debye Spheres" |
---|
99 | Button FFTButton_8a,pos={180,297},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Debye Results" |
---|
100 | |
---|
101 | Button FFTButton_14,pos={13,360},size={130,20},proc=FFT_BinnedSLDButtonProc,title="Do Binned SLD" |
---|
102 | Button FFTButton_14a,pos={180,360},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot SLD Results" |
---|
103 | |
---|
104 | SetVariable FFTSetVar_4,pos={201,4},size={100,15},title="FFT time(s)" |
---|
105 | SetVariable FFTSetVar_4,limits={0,0,0},value= FFT_estTime,noedit= 1,live= 1,format="%d" |
---|
106 | Button FFTButton_9,pos={200,400},size={100,20},proc=FFT_Get2DSlice,title="Get 2D Slice" |
---|
107 | Button FFTButton_10,pos={169,156},size={130,20},proc=FFT_TransposeMat,title="Transpose Matrix" |
---|
108 | Button FFTButton_11,pos={169,189},size={130,20},proc=FFT_RotateMat,title="Rotate Matrix" |
---|
109 | Button FFTButton_12,pos={168,219},size={130,20},proc=FFT_AddRotatedObject,title="Add Rotated Obj" |
---|
110 | Button FFTButton_12,disable=2 // hide this button |
---|
111 | Button FFTButton_13,pos={14,109},size={120,20},proc=FFTFillSolventMatrixProc,title="Solvent Matrix" |
---|
112 | SetVariable FFTSetVar_5,pos={155,111},size={150,15},title="Solvent SLD (10^-7)" |
---|
113 | SetVariable FFTSetVar_5,limits={-99,99,1},value= FFT_SolventSLD,live= 1 |
---|
114 | Button FFTButton_15,pos={209,430},size={90,20},proc=Interp2DSliceButton,title="Interp 2D" |
---|
115 | Button FFTButton_16,pos={14,460},size={70,20},proc=FFTHelpButton,title="Help" |
---|
116 | EndMacro |
---|
117 | |
---|
118 | // Save a matrix wave, plus the N, T, and solvent values in the wave note for reloading |
---|
119 | Function FFTSaveMatrixButtonProc(ba) : ButtonControl |
---|
120 | STRUCT WMButtonAction &ba |
---|
121 | |
---|
122 | String win = ba.win |
---|
123 | |
---|
124 | switch (ba.eventCode) |
---|
125 | case 2: |
---|
126 | // click code here |
---|
127 | String fileStr="" |
---|
128 | SaveMyMatrix(fileStr) |
---|
129 | |
---|
130 | break |
---|
131 | endswitch |
---|
132 | |
---|
133 | return 0 |
---|
134 | End |
---|
135 | |
---|
136 | // this will wave as Igor Binary, so be sure to use the ".ibw extension. |
---|
137 | // - this could possibly be enforced, but that's maybe not necessary at this stage. |
---|
138 | // |
---|
139 | Function SaveMyMatrix(fileStr) |
---|
140 | String fileStr |
---|
141 | |
---|
142 | WAVE mat=root:mat |
---|
143 | NVAR FFT_T = root:FFT_T |
---|
144 | NVAR FFT_N = root:FFT_N |
---|
145 | NVAR FFT_SolventSLD = root:FFT_SolventSLD |
---|
146 | String str="" |
---|
147 | sprintf str,"FFT_T=%g;FFT_N=%d;FFT_SolventSLD=%d;",FFT_T,FFT_N,FFT_SolventSLD |
---|
148 | Note mat,str |
---|
149 | Save/C/P=home mat as fileStr //will ask for a file name, save as Igor Binary |
---|
150 | Note/K mat //kill wave note on exiting since I don't properly update this anywhere else |
---|
151 | |
---|
152 | return(0) |
---|
153 | end |
---|
154 | |
---|
155 | |
---|
156 | // load in a previously saved matrix, and reset FFT_N, FFT_T and solvent |
---|
157 | // from the wave note when saved |
---|
158 | Function FFTLoadMatrixButtonProc(ba) : ButtonControl |
---|
159 | STRUCT WMButtonAction &ba |
---|
160 | |
---|
161 | String win = ba.win |
---|
162 | |
---|
163 | switch (ba.eventCode) |
---|
164 | case 2: |
---|
165 | // click code here |
---|
166 | String fileStr="" |
---|
167 | ReloadMatrix(fileStr) |
---|
168 | |
---|
169 | break |
---|
170 | endswitch |
---|
171 | |
---|
172 | return 0 |
---|
173 | End |
---|
174 | |
---|
175 | |
---|
176 | Function ReloadMatrix(fileStr) |
---|
177 | String fileStr |
---|
178 | |
---|
179 | LoadWave/M/O/W/P=home fileStr //will ask for a file, Igor Binary format is assumed here |
---|
180 | String str |
---|
181 | str=note(mat) |
---|
182 | NVAR FFT_T = root:FFT_T |
---|
183 | NVAR FFT_N = root:FFT_N |
---|
184 | NVAR FFT_SolventSLD = root:FFT_SolventSLD |
---|
185 | |
---|
186 | FFT_T = NumberByKey("FFT_T", str, "=" ,";") |
---|
187 | FFT_N = NumberByKey("FFT_N", str, "=" ,";") |
---|
188 | FFT_SolventSLD = NumberByKey("FFT_SolventSLD", str, "=" ,";") |
---|
189 | |
---|
190 | // if I got bad values, put in default values |
---|
191 | if(numtype(FFT_T) != 0 ) |
---|
192 | FFT_T = 5 |
---|
193 | endif |
---|
194 | if(numtype(FFT_N) != 0 ) |
---|
195 | FFT_N = DimSize(mat,0) |
---|
196 | endif |
---|
197 | if(numtype(FFT_SolventSLD) != 0 ) |
---|
198 | FFT_SolventSLD = 0 |
---|
199 | endif |
---|
200 | |
---|
201 | Print "Loaded matrix parameters = ",str |
---|
202 | Execute "NumberOfPoints()" |
---|
203 | |
---|
204 | return(0) |
---|
205 | end |
---|
206 | |
---|
207 | // utility function that reads the wave note from a binary file |
---|
208 | // |
---|
209 | Function/S LoadNoteFunc(PName,FName[,FileRef]) |
---|
210 | String PName, FName |
---|
211 | Variable FileRef |
---|
212 | |
---|
213 | Variable noteStart, noteLength, version, dependLength |
---|
214 | String noteStr, typeStr = ".ibw" |
---|
215 | if (ParamIsDefault(FileRef)) |
---|
216 | Open/R/P=$PName/T=typeStr fileRef, as FName //open the file |
---|
217 | endif |
---|
218 | FSetPos fileRef, 0 |
---|
219 | FBinRead/F=2 fileRef, version |
---|
220 | |
---|
221 | |
---|
222 | if (version == 5) |
---|
223 | |
---|
224 | FSetPos fileRef, 4 |
---|
225 | Make/N=(3)/I/Free SizeWave |
---|
226 | FBinRead FileRef,SizeWave |
---|
227 | noteStart = SizeWave[0] |
---|
228 | DependLength = SizeWave[1] |
---|
229 | NoteLength = SizeWave[2] |
---|
230 | noteStart += dependLength+64 |
---|
231 | |
---|
232 | elseif (version == 2) |
---|
233 | |
---|
234 | FBinRead/F=3 fileRef, noteStart |
---|
235 | // FBinRead/F=4 fileRef, dependLength |
---|
236 | FBinRead/F=3 fileRef, noteLength |
---|
237 | noteStart += 16 |
---|
238 | |
---|
239 | else |
---|
240 | |
---|
241 | if (ParamIsDefault(FileRef)) |
---|
242 | Close(FileRef) //close the file |
---|
243 | endif |
---|
244 | return "" |
---|
245 | |
---|
246 | endif |
---|
247 | if (!NoteLength) |
---|
248 | if (ParamIsDefault(FileRef)) |
---|
249 | Close(FileRef) //close the file |
---|
250 | endif |
---|
251 | return("") |
---|
252 | endif |
---|
253 | FSetPos fileRef, noteStart |
---|
254 | NoteStr = PadString("",NoteLength,0) |
---|
255 | FBinRead FileRef,NoteStr |
---|
256 | |
---|
257 | if (ParamIsDefault(FileRef)) |
---|
258 | Close(FileRef) //close the file |
---|
259 | endif |
---|
260 | return noteStr |
---|
261 | |
---|
262 | End //LoadNoteFunc |
---|
263 | |
---|
264 | |
---|
265 | Function FFTHelpButton(ba) : ButtonControl |
---|
266 | STRUCT WMButtonAction &ba |
---|
267 | |
---|
268 | String win = ba.win |
---|
269 | |
---|
270 | switch (ba.eventCode) |
---|
271 | case 2: |
---|
272 | // click code here |
---|
273 | DisplayHelpTopic/Z/K=1 "Real-Space Modeling of SANS Data" |
---|
274 | if(V_flag !=0) |
---|
275 | DoAlert 0,"The Real-Space Modeling Help file could not be found" |
---|
276 | endif |
---|
277 | break |
---|
278 | endswitch |
---|
279 | |
---|
280 | return 0 |
---|
281 | End |
---|
282 | |
---|
283 | |
---|
284 | Function Interp2DSliceButton(ba) : ButtonControl |
---|
285 | STRUCT WMButtonAction &ba |
---|
286 | |
---|
287 | switch( ba.eventCode ) |
---|
288 | case 2: // mouse up |
---|
289 | // click code here |
---|
290 | String folderStr="" |
---|
291 | Prompt folderStr, "Pick a 2D data folder" // Set prompt for x param |
---|
292 | DoPrompt "Enter data folder", folderStr |
---|
293 | if (V_Flag || strlen(folderStr) == 0) |
---|
294 | return -1 // User canceled, null string entered |
---|
295 | endif |
---|
296 | |
---|
297 | Interpolate2DSliceToData(folderStr) |
---|
298 | |
---|
299 | //display the 2D plane |
---|
300 | DoWindow FFT_Interp2D_log |
---|
301 | if(V_flag == 0) |
---|
302 | Execute "Display2DInterpSlice_log()" |
---|
303 | endif |
---|
304 | |
---|
305 | break |
---|
306 | endswitch |
---|
307 | |
---|
308 | return 0 |
---|
309 | End |
---|
310 | |
---|
311 | |
---|
312 | |
---|
313 | Function FFT_Get2DSlice(ctrlName) : ButtonControl |
---|
314 | String ctrlName |
---|
315 | |
---|
316 | WAVE/Z dmat = root:dmat |
---|
317 | if(WaveExists(dmat)==1) |
---|
318 | get2DSlice(dmat) |
---|
319 | endif |
---|
320 | |
---|
321 | //display the 2D plane |
---|
322 | DoWindow FFT_Slice2D |
---|
323 | if(V_flag == 0) |
---|
324 | Execute "Display2DSlice()" |
---|
325 | endif |
---|
326 | |
---|
327 | //display the 2D plane |
---|
328 | DoWindow FFT_Slice2D_log |
---|
329 | if(V_flag == 0) |
---|
330 | Execute "Display2DSlice_log()" |
---|
331 | endif |
---|
332 | |
---|
333 | //display the 1D binning of the 2D plane |
---|
334 | DoWindow FFT_Slice1D |
---|
335 | if(V_flag == 0) |
---|
336 | Execute "Slice2_1D()" |
---|
337 | endif |
---|
338 | |
---|
339 | return(0) |
---|
340 | End |
---|
341 | |
---|
342 | Proc Display2DSlice() |
---|
343 | PauseUpdate; Silent 1 // building window... |
---|
344 | Display /W=(1038,44,1404,403) |
---|
345 | DoWindow/C FFT_Slice2D |
---|
346 | AppendImage/T detPlane |
---|
347 | ModifyImage detPlane ctab= {*,*,YellowHot,0} |
---|
348 | ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 |
---|
349 | ModifyGraph mirror=2 |
---|
350 | ModifyGraph nticks=4 |
---|
351 | ModifyGraph minor=1 |
---|
352 | ModifyGraph fSize=9 |
---|
353 | ModifyGraph standoff=0 |
---|
354 | ModifyGraph tkLblRot(left)=90 |
---|
355 | ModifyGraph btLen=3 |
---|
356 | ModifyGraph tlOffset=-2 |
---|
357 | SetAxis/A/R left |
---|
358 | End |
---|
359 | |
---|
360 | Proc Display2DInterpSlice_log() |
---|
361 | PauseUpdate; Silent 1 // building window... |
---|
362 | Display /W=(1038,44,1404,403) |
---|
363 | DoWindow/C FFT_Interp2D_log |
---|
364 | AppendImage/T interp2DSlice_log |
---|
365 | ModifyImage interp2DSlice_log ctab= {*,*,YellowHot,0} |
---|
366 | ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 |
---|
367 | ModifyGraph mirror=2 |
---|
368 | ModifyGraph nticks=4 |
---|
369 | ModifyGraph minor=1 |
---|
370 | ModifyGraph fSize=9 |
---|
371 | ModifyGraph standoff=0 |
---|
372 | ModifyGraph tkLblRot(left)=90 |
---|
373 | ModifyGraph btLen=3 |
---|
374 | ModifyGraph tlOffset=-2 |
---|
375 | SetAxis/A/R left |
---|
376 | End |
---|
377 | |
---|
378 | Proc Display2DSlice_log() |
---|
379 | PauseUpdate; Silent 1 // building window... |
---|
380 | Display /W=(1038,44,1404,403) |
---|
381 | DoWindow/C FFT_Slice2D_log |
---|
382 | AppendImage/T logP |
---|
383 | ModifyImage logP ctab= {*,*,YellowHot,0} |
---|
384 | ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 |
---|
385 | ModifyGraph mirror=2 |
---|
386 | ModifyGraph nticks=4 |
---|
387 | ModifyGraph minor=1 |
---|
388 | ModifyGraph fSize=9 |
---|
389 | ModifyGraph standoff=0 |
---|
390 | ModifyGraph tkLblRot(left)=90 |
---|
391 | ModifyGraph btLen=3 |
---|
392 | ModifyGraph tlOffset=-2 |
---|
393 | SetAxis/A/R left |
---|
394 | End |
---|
395 | Proc Slice2_1D() |
---|
396 | PauseUpdate; Silent 1 // building window... |
---|
397 | Display /W=(1034,425,1406,763) iBin_2d vs qBin_2d |
---|
398 | DoWindow/C FFT_Slice1D |
---|
399 | ModifyGraph mode=4 |
---|
400 | ModifyGraph marker=19 |
---|
401 | ModifyGraph msize=2 |
---|
402 | ModifyGraph grid=1 |
---|
403 | ModifyGraph log=1 |
---|
404 | ModifyGraph mirror=2 |
---|
405 | Legend |
---|
406 | EndMacro |
---|
407 | |
---|
408 | |
---|
409 | Function FFT_TransposeMat(ctrlName) : ButtonControl |
---|
410 | String ctrlName |
---|
411 | |
---|
412 | Variable mode=1 |
---|
413 | Prompt mode,"Transform XYZ to:",popup,"XZY;ZXY;ZYX;YZX;YXZ;" |
---|
414 | DoPrompt "Transform mode",mode |
---|
415 | if (V_Flag) |
---|
416 | return 0 // user canceled |
---|
417 | endif |
---|
418 | |
---|
419 | fFFT_TransposeMat(mode) |
---|
420 | |
---|
421 | return (0) |
---|
422 | End |
---|
423 | |
---|
424 | // starting in XYZ |
---|
425 | // mode 1;2;3;4;5; |
---|
426 | //correspond to |
---|
427 | // "XZY;ZXY;ZYX;YZX;YXZ;" |
---|
428 | // |
---|
429 | Function fFFT_TransposeMat(mode) |
---|
430 | Variable mode |
---|
431 | |
---|
432 | WAVE/Z mat=mat |
---|
433 | ImageTransform /G=(mode) transposeVol mat |
---|
434 | WAVE M_VolumeTranspose=M_VolumeTranspose |
---|
435 | Duplicate/O M_VolumeTranspose mat |
---|
436 | |
---|
437 | return(0) |
---|
438 | End |
---|
439 | |
---|
440 | Function FFT_RotateMat(ctrlName) : ButtonControl |
---|
441 | String ctrlName |
---|
442 | |
---|
443 | Variable degree=45,sense=1 |
---|
444 | Prompt degree, "Degrees of rotation around Z-axis:" |
---|
445 | // Prompt sense, "Direction of rotation:",popup,"CW;CCW;" |
---|
446 | DoPrompt "Enter paramters for rotation", degree |
---|
447 | |
---|
448 | if (V_Flag) |
---|
449 | return 0 // user canceled |
---|
450 | endif |
---|
451 | |
---|
452 | fFFT_RotateMat(degree) |
---|
453 | return(0) |
---|
454 | End |
---|
455 | // note that the rotation is not perfect. if the rotation produces an |
---|
456 | // odd number of pixels, then the object will "walk" one pixel. Also, some |
---|
457 | // small artifacts are introduced, simply due to the voxelization of the object |
---|
458 | // as it is rotated. rotating a cylinder 10 then 350 shows a few extra "bumps" on |
---|
459 | // the surface, but the calculated scattering is virtually identical. |
---|
460 | // |
---|
461 | // these issues, may be correctable, if needed |
---|
462 | // |
---|
463 | Function fFFT_RotateMat(degree) |
---|
464 | Variable degree |
---|
465 | |
---|
466 | Variable fill=0 |
---|
467 | |
---|
468 | WAVE mat = mat |
---|
469 | Variable dx,dy,dz,nx,ny,nz |
---|
470 | dx = DimSize(mat,0) |
---|
471 | dy = DimSize(mat,1) |
---|
472 | dz = DimSize(mat,2) |
---|
473 | |
---|
474 | //? the /W and /C flags seem to be special cases for 90 deg rotations |
---|
475 | ImageRotate /A=(degree)/E=(fill) mat //mat is not overwritten, unless /O |
---|
476 | |
---|
477 | nx = DimSize(M_RotatedImage,0) |
---|
478 | ny = DimSize(M_RotatedImage,1) |
---|
479 | nz = DimSize(M_RotatedImage,2) |
---|
480 | // Print "rotated dims = ",dx,dy,dz,nx,ny,nz |
---|
481 | Variable delx,dely,odd=0 |
---|
482 | delx = nx-dx |
---|
483 | dely = ny-dy |
---|
484 | |
---|
485 | if(mod(delx,2) != 0) //sometimes the new x,y dims are odd! |
---|
486 | odd = 1 |
---|
487 | endif |
---|
488 | // delx = (delx + odd)/2 |
---|
489 | // dely = (dely + odd)/2 |
---|
490 | delx = trunc(delx/2) |
---|
491 | dely = trunc(dely/2) |
---|
492 | |
---|
493 | //+odd removes an extra point if there is one |
---|
494 | Duplicate/O/R=[delx+odd,nx-delx-1][dely+odd,ny-dely-1][0,nz-1] M_RotatedImage mat |
---|
495 | |
---|
496 | // - not sure why the duplicate operation changes the start and delta of mat, but I |
---|
497 | // need to reset the start and delta to be sure that the display is correct, and that the scaling |
---|
498 | // is read correctly later |
---|
499 | SetScale/P x 0,1,"", mat |
---|
500 | SetScale/P y 0,1,"", mat |
---|
501 | |
---|
502 | nx = DimSize(mat,0) |
---|
503 | ny = DimSize(mat,1) |
---|
504 | nz = DimSize(mat,2) |
---|
505 | // Print "redim = ",dx,dy,dz,nx,ny,nz |
---|
506 | |
---|
507 | KillWaves/Z M_RotatedImage |
---|
508 | |
---|
509 | End |
---|
510 | |
---|
511 | // also look for ImageTransform operations that will allow translation of the objects. |
---|
512 | // then simple objects could be oriented @0,0,0, and translated to the correct position (and added) |
---|
513 | // |
---|
514 | Function FFT_AddRotatedObject(ctrlName) : ButtonControl |
---|
515 | String ctrlName |
---|
516 | |
---|
517 | Print "Not yet implemented" |
---|
518 | End |
---|
519 | |
---|
520 | |
---|
521 | Function FFT_MakeMatrixButtonProc(ctrlName) : ButtonControl |
---|
522 | String ctrlName |
---|
523 | |
---|
524 | NVAR nn=root:FFT_N |
---|
525 | if(mod(nn, 2 ) ==1) //force an even number for FFT |
---|
526 | nn +=1 |
---|
527 | endif |
---|
528 | MakeMatrix("mat",nn,nn,nn) |
---|
529 | End |
---|
530 | |
---|
531 | Function FFTMakeGizmoButtonProc(ctrlName) : ButtonControl |
---|
532 | String ctrlName |
---|
533 | |
---|
534 | DoWindow/F Gizmo_VoxelMat |
---|
535 | if(V_flag==0) |
---|
536 | Execute "Gizmo_VoxelMat()" |
---|
537 | endif |
---|
538 | End |
---|
539 | |
---|
540 | Function FFTDrawSphereButtonProc(ctrlName) : ButtonControl |
---|
541 | String ctrlName |
---|
542 | |
---|
543 | Execute "FFTDrawSphereProc()" |
---|
544 | End |
---|
545 | |
---|
546 | Proc FFTDrawSphereProc(matStr,rad,xc,yc,zc,fill) |
---|
547 | String matStr="mat" |
---|
548 | Variable rad=25,xc=50,yc=50,zc=50,fill=10 |
---|
549 | Prompt matStr,"the wave" //,popup,WaveList("*",";","") |
---|
550 | Prompt rad,"enter real radius (A)" |
---|
551 | Prompt xc,"enter the X-center" |
---|
552 | Prompt yc,"enter the Y-center" |
---|
553 | Prompt zc,"enter the Z-center" |
---|
554 | Prompt fill,"fill SLD value" |
---|
555 | |
---|
556 | Variable grid=root:FFT_T |
---|
557 | |
---|
558 | FillSphereRadius($matStr,grid,rad,xc,yc,zc,fill) |
---|
559 | |
---|
560 | End |
---|
561 | |
---|
562 | Function FFTDrawZCylinderButtonProc(ctrlName) : ButtonControl |
---|
563 | String ctrlName |
---|
564 | |
---|
565 | Execute "FFTDrawCylinder()" |
---|
566 | End |
---|
567 | |
---|
568 | Proc FFTDrawCylinder(direction,matStr,rad,len,xc,yc,zc,fill) |
---|
569 | String direction |
---|
570 | String matStr="mat" |
---|
571 | Variable rad=25,len=300,xc=50,yc=50,zc=50,fill=10 |
---|
572 | Prompt direction, "Direction", popup "X;Y;Z;" |
---|
573 | Prompt matStr,"the wave" //,popup,WaveList("*",";","") |
---|
574 | Prompt rad,"enter real radius (A)" |
---|
575 | Prompt len,"enter length (A)" |
---|
576 | Prompt xc,"enter the X-center" |
---|
577 | Prompt yc,"enter the Y-center" |
---|
578 | Prompt zc,"enter the Z-center" |
---|
579 | Prompt fill,"fill SLD value" |
---|
580 | |
---|
581 | |
---|
582 | Variable grid=root:FFT_T |
---|
583 | |
---|
584 | if(cmpstr(direction,"X")==0) |
---|
585 | FillXCylinder($matStr,grid,rad,xc,yc,zc,len,fill) |
---|
586 | endif |
---|
587 | if(cmpstr(direction,"Y")==0) |
---|
588 | FillYCylinder($matStr,grid,rad,xc,yc,zc,len,fill) |
---|
589 | endif |
---|
590 | if(cmpstr(direction,"Z")==0) |
---|
591 | FillZCylinder($matStr,grid,rad,xc,yc,zc,len,fill) |
---|
592 | endif |
---|
593 | |
---|
594 | End |
---|
595 | |
---|
596 | |
---|
597 | Function DoTheFFT_ButtonProc(ctrlName) : ButtonControl |
---|
598 | String ctrlName |
---|
599 | |
---|
600 | Execute "DoFFT()" |
---|
601 | End |
---|
602 | |
---|
603 | Function FFT_PlotResultsButtonProc(ctrlName) : ButtonControl |
---|
604 | String ctrlName |
---|
605 | |
---|
606 | Variable first=0 |
---|
607 | DoWindow/F FFT_IQ |
---|
608 | if(V_flag==0) |
---|
609 | first = 1 |
---|
610 | Display /W=(295,44,627,302) |
---|
611 | DoWindow/C FFT_IQ |
---|
612 | Endif |
---|
613 | |
---|
614 | // append the desired data, if it's not already there |
---|
615 | // FFTButton_4 = FFT = iBin |
---|
616 | // FFTButton_7a = binned = _XOP |
---|
617 | // FFTButton_8a = Debye = _full |
---|
618 | // FFTButton_14a = SLD = _SLD |
---|
619 | strswitch(ctrlName) |
---|
620 | case "FFTButton_4": |
---|
621 | if(!isTraceOnGraph("iBin","FFT_IQ") && exists("iBin")==1) //only append if it's not already there |
---|
622 | AppendToGraph /W=FFT_IQ iBin vs qBin |
---|
623 | ModifyGraph mode=4,marker=19,msize=2,rgb(iBin)=(65535,0,0) |
---|
624 | endif |
---|
625 | break |
---|
626 | case "FFTButton_7a": |
---|
627 | if(!isTraceOnGraph("ival_XOP","FFT_IQ") && exists("ival_XOP")==1) //only append if it's not already there |
---|
628 | AppendToGraph /W=FFT_IQ ival_XOP vs qval_XOP |
---|
629 | ModifyGraph mode=4,marker=19,msize=2,rgb(ival_XOP)=(1,12815,52428) |
---|
630 | endif |
---|
631 | break |
---|
632 | case "FFTButton_8a": |
---|
633 | if(!isTraceOnGraph("ival_full","FFT_IQ") && exists("ival_full")==1) //only append if it's not already there |
---|
634 | AppendToGraph /W=FFT_IQ ival_full vs qval_full |
---|
635 | ModifyGraph mode=4,marker=19,msize=2,rgb(ival_full)=(0,0,0) |
---|
636 | endif |
---|
637 | break |
---|
638 | case "FFTButton_14a": |
---|
639 | if(!isTraceOnGraph("ival_SLD","FFT_IQ") && exists("ival_SLD")==1) //only append if it's not already there |
---|
640 | AppendToGraph /W=FFT_IQ ival_SLD vs qval_SLD |
---|
641 | ModifyGraph mode=4,marker=19,msize=2,rgb(ival_SLD)=(2,39321,1) |
---|
642 | endif |
---|
643 | break |
---|
644 | endswitch |
---|
645 | |
---|
646 | if(first) |
---|
647 | ModifyGraph mode=4 |
---|
648 | ModifyGraph marker=19 |
---|
649 | ModifyGraph msize=2 |
---|
650 | ModifyGraph gaps=0 |
---|
651 | ModifyGraph grid=1 |
---|
652 | ModifyGraph log=1 |
---|
653 | ModifyGraph mirror=2 |
---|
654 | Legend |
---|
655 | endif |
---|
656 | |
---|
657 | return(0) |
---|
658 | End |
---|
659 | |
---|
660 | Function isTraceOnGraph(traceStr,winStr) |
---|
661 | String traceStr,winStr |
---|
662 | |
---|
663 | Variable isOn=0 |
---|
664 | String str |
---|
665 | str = TraceNameList("winStr",";",1) //only normal traces |
---|
666 | isOn = stringmatch(str,traceStr) //must be an exact match |
---|
667 | |
---|
668 | return(isOn) |
---|
669 | End |
---|
670 | |
---|
671 | Function FFTEraseMatrixButtonProc(ctrlName) : ButtonControl |
---|
672 | String ctrlName |
---|
673 | |
---|
674 | Wave mat=root:mat |
---|
675 | FastOp mat=0 |
---|
676 | return(0) |
---|
677 | End |
---|
678 | |
---|
679 | Function FFTFillSolventMatrixProc(ctrlName) : ButtonControl |
---|
680 | String ctrlName |
---|
681 | |
---|
682 | Wave mat=root:mat |
---|
683 | NVAR val=root:FFT_SolventSLD |
---|
684 | mat=val |
---|
685 | return(0) |
---|
686 | End |
---|
687 | |
---|
688 | Function FFT_AltiSpheresButtonProc(ctrlName) : ButtonControl |
---|
689 | String ctrlName |
---|
690 | |
---|
691 | Execute "DoSpheresCalcFFTPanel()" |
---|
692 | End |
---|
693 | |
---|
694 | Function FFT_BinnedSpheresButtonProc(ctrlName) : ButtonControl |
---|
695 | String ctrlName |
---|
696 | |
---|
697 | Execute "DoBinnedSpheresCalcFFTPanel()" |
---|
698 | End |
---|
699 | |
---|
700 | Function FFT_BinnedSLDButtonProc(ctrlName) : ButtonControl |
---|
701 | String ctrlName |
---|
702 | |
---|
703 | Execute "DoBinnedSLDCalcFFTPanel()" |
---|
704 | End |
---|
705 | |
---|
706 | |
---|
707 | |
---|
708 | |
---|
709 | |
---|
710 | /////UTILITIES |
---|
711 | |
---|
712 | // inverts the values in the matrix 0<->1 |
---|
713 | Function InvertMatrixFill(mat) |
---|
714 | Wave mat |
---|
715 | |
---|
716 | mat = (mat==1) ? 0 : 1 |
---|
717 | End |
---|
718 | |
---|
719 | //overwrites any existing matrix |
---|
720 | // matrix is byte, to save space |
---|
721 | //forces the name to be "mat" |
---|
722 | // |
---|
723 | // switched to signed byte |
---|
724 | // |
---|
725 | Function MakeMatrix(nam,xd,yd,zd) |
---|
726 | String nam |
---|
727 | Variable xd,yd,zd |
---|
728 | |
---|
729 | nam="mat" |
---|
730 | Print "Matrix has been created and named \"mat\"" |
---|
731 | // Make/O/U/B/N=(xd,yd,zd) $nam |
---|
732 | Make/O/B/N=(xd,yd,zd) $nam |
---|
733 | End |
---|
734 | |
---|
735 | |
---|
736 | // calculate the average center-to-center distance between points |
---|
737 | // assumes evenly spaced on a cubic grid |
---|
738 | // |
---|
739 | Proc Center_to_Center(np) |
---|
740 | Variable np = 100 |
---|
741 | |
---|
742 | Variable Nedge = root:FFT_N |
---|
743 | Variable Tscale = root:FFT_T |
---|
744 | |
---|
745 | Variable davg |
---|
746 | |
---|
747 | davg = (Nedge*Tscale)^3 / np |
---|
748 | davg = davg^(1/3) |
---|
749 | Print "Average separation (A) = ",davg |
---|
750 | |
---|
751 | End |
---|
752 | |
---|
753 | // calculate the number of points required on a grid |
---|
754 | // to yield a given average center-to-center distance between points |
---|
755 | // assumes evenly spaced on a cubic grid |
---|
756 | // |
---|
757 | // davg and Tscale are the same units, Angstroms |
---|
758 | // |
---|
759 | Proc Davg_to_Np(davg) |
---|
760 | Variable davg = 400 |
---|
761 | |
---|
762 | Variable Nedge = root:FFT_N |
---|
763 | Variable Tscale = root:FFT_T |
---|
764 | |
---|
765 | Variable np |
---|
766 | |
---|
767 | np = (Nedge*Tscale)^3 / davg^3 |
---|
768 | Print "Number of points required = ",np |
---|
769 | |
---|
770 | End |
---|
771 | |
---|
772 | |
---|
773 | |
---|
774 | |
---|
775 | |
---|
776 | |
---|
777 | // The matrix is not necessarily 0|1, this reports the number of filled voxels |
---|
778 | // - needed to estimate the time required for the AltiVec_Spheres calculation |
---|
779 | Proc NumberOfPoints() |
---|
780 | |
---|
781 | Print "Number of points = ",NumberOfPoints_Occ(root:mat) |
---|
782 | Print "Fraction occupied = ",VolumeFraction_Occ(root:mat) |
---|
783 | Print "Overall Cube Edge [A] = ",root:FFT_T * root:FFT_N |
---|
784 | Print "Found values in matrix = ",ListOfValues(root:mat) |
---|
785 | |
---|
786 | End |
---|
787 | |
---|
788 | Function NumberOfPoints_Occ(m) |
---|
789 | Wave m |
---|
790 | |
---|
791 | Variable num = NonZeroValues(m) |
---|
792 | |
---|
793 | return(num) |
---|
794 | End |
---|
795 | |
---|
796 | |
---|
797 | Function VolumeFraction_Occ(m) |
---|
798 | Wave m |
---|
799 | |
---|
800 | Variable num = NonZeroValues(m) |
---|
801 | Variable dim = DimSize(m,0) |
---|
802 | |
---|
803 | return(num/dim^3) |
---|
804 | End |
---|
805 | |
---|
806 | //it's a byte wave, so I can't use the trick of setting the zero values to NaN |
---|
807 | // since NaN can't be expressed as byte. So make it binary 0|1 |
---|
808 | // |
---|
809 | // - the assumption here is that the defined solvent value is not part of the |
---|
810 | // particle volume. |
---|
811 | // |
---|
812 | Function NonZeroValues(m) |
---|
813 | Wave m |
---|
814 | |
---|
815 | Variable num |
---|
816 | NVAR val = root:FFT_SolventSLD |
---|
817 | |
---|
818 | // Variable t1=ticks,dim |
---|
819 | // dim = DimSize(m, 0 ) //assume NxNxN |
---|
820 | Duplicate/O m,mz |
---|
821 | |
---|
822 | MultiThread mz = (m[p][q] != val) ? 1 : 0 |
---|
823 | |
---|
824 | WaveStats/Q/M=1 mz // NaN and Inf are not reported in V_npnts |
---|
825 | |
---|
826 | num = V_npnts*V_avg |
---|
827 | |
---|
828 | KillWaves/Z mz |
---|
829 | // Print "Number of points = ",num |
---|
830 | // Print "Fraction occupied = ",num/V_npnts |
---|
831 | |
---|
832 | // Print "time = ",(ticks-t1)/60.15 |
---|
833 | return(num) |
---|
834 | End |
---|
835 | |
---|
836 | // |
---|
837 | // return a list of the different values of the voxels in the matrix |
---|
838 | // |
---|
839 | Function/S ListOfValues(m) |
---|
840 | Wave m |
---|
841 | |
---|
842 | String list="" |
---|
843 | Variable done |
---|
844 | |
---|
845 | Duplicate/O m,mz |
---|
846 | |
---|
847 | done=0 |
---|
848 | do |
---|
849 | WaveStats/Q/M=1 mz // NaN and Inf are not reported in V_npnts |
---|
850 | if(V_max == V_min) |
---|
851 | list += num2str(V_min) + ";" |
---|
852 | done = 1 |
---|
853 | else |
---|
854 | list += num2str(V_max) + ";" |
---|
855 | MultiThread mz = mz[p][q] == V_max ? V_min : mz[p][q] // replace the max with min |
---|
856 | endif |
---|
857 | while(!done) |
---|
858 | |
---|
859 | // Print "Found values in matrix = ",list |
---|
860 | KillWaves/Z mz |
---|
861 | |
---|
862 | return(list) |
---|
863 | End |
---|
864 | |
---|
865 | |
---|
866 | |
---|
867 | // returns estimate in seconds |
---|
868 | // |
---|
869 | // updated for quad-core iMac, 2010 |
---|
870 | // |
---|
871 | // type 3 = binned distances and SLD |
---|
872 | // type 2 = binned distances |
---|
873 | // type 1 is rather meaningless |
---|
874 | // type 0 = is the Debye, double sum (multithreaded) |
---|
875 | // |
---|
876 | // |
---|
877 | // - types 2 and 3, the binned methods are inherently AAO, so there is no significant |
---|
878 | // dependence on the number of q-values (unless it's >> 1000). So values reported are |
---|
879 | // for 100 q-values. |
---|
880 | // |
---|
881 | // There is a function Timing_Method(type) in FFT_ConcentratedSpheres.ipf to automate some of this |
---|
882 | // |
---|
883 | Function EstimatedTime(nx,nq,type) |
---|
884 | Variable nx,nq,type |
---|
885 | |
---|
886 | Variable est=0 |
---|
887 | |
---|
888 | if(type==0) // "full" XOP |
---|
889 | est = 0.4 + 1.622e-5*nx+4.86e-7*nx^2 |
---|
890 | est /= 100 |
---|
891 | est *= nq |
---|
892 | endif |
---|
893 | |
---|
894 | if(type==1) //Igor function (much slower, 9x) |
---|
895 | est = (nx^2)*0.0000517 //empirical values (seconds) for igor code, 100 q-values |
---|
896 | est /= 100 // per q now... |
---|
897 | est *= nq |
---|
898 | endif |
---|
899 | |
---|
900 | if(type==2) //binned distances, slow parts XOPed |
---|
901 | est = 0.0680 + 2.94e-6*nx+4.51e-9*nx^2 //with threading |
---|
902 | |
---|
903 | // est /= 100 // per q now... |
---|
904 | // est *= nq |
---|
905 | endif |
---|
906 | |
---|
907 | if(type==3) //binned distances AND sld, slow parts XOPed |
---|
908 | est = 0.576 + 4.22e-7*nx+1.76e-8*nx^2 //with threading |
---|
909 | |
---|
910 | // est /= 100 // per q now... |
---|
911 | // est *= nq |
---|
912 | endif |
---|
913 | |
---|
914 | return(est) |
---|
915 | End |
---|