source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Cubes.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: 14.5 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4// FFT and DebyeSphere calculations are now properly normalized to the conditions of:
5// deltaRho = 1e-6 (set as a global)
6// volume fraction = occupied fraction
7//
8
9// another thing to add - for the FFT when taking a slice. To test the rotational average, set up a loop
10// that does the FFT, rotates, sums the intensity (the q's are the same), and reports the average. May need to
11// rotate a large number of times to get a "real" average
12//
13//
14// -- Feb 2011
15//      You can use N=512. The FFT takes about 140 s (first time), then 75s after that.
16//      The size of the whole experiment must be enormous.
17//      512^2 = 134 MB for the byte matrix.
18//      (512^3)/2*4 = 268 MB for the FFT result
19//
20// Function Interpolate2DSliceToData(folderStr)
21//              -- this function takes the 2D slice and interpolates the q-values to the experimental data for comparison (or a fit!)
22//
23Proc DoFFT()
24       
25        Variable t0=ticks,t1=ticks
26       
27        fDoFFT()
28//      Print "FFT time (s) = ",(ticks - t1)/60.15
29       
30        t1=ticks
31        fDoCorrection(dmat)     // dmat will exist, and be complex on input
32//      Print "Cube correction time (s) = ",(ticks - t1)/60.15
33       
34        t1 = ticks
35        fGetMagnitude(tmp)              //complex on input, returns as real-valued
36//      Print "Magnitude Calculation time (s) = ",(ticks - t1)/60.15
37       
38        t1=ticks
39        fDoBinning(dmat,0)       // dmat is purely real when the binning is done, do not normalize result
40//      Print "Binning time (s) = ",(ticks - t1)/60.15
41       
42        Print "Total FFT time (s) = ",(ticks - t0)/60.15
43        Print " "
44       
45        //normalize the data
46        Variable Tscale=root:FFT_T
47
48        Variable delRho,vol,nx,phi
49       
50        delRho = root:FFT_delRho                                        //simply the units of SLD, not the absolute value
51       
52//      WaveStats/Q mat
53// must do this to get the number of non-zero voxels
54        Variable FFT_N = root:FFT_N
55//      nx = NonZeroValues(mat)
56//      phi = nx/FFT_N^3
57//      vol = nx*(Tscale)^3
58       
59//      iBin *= phi
60//      iBin /= vol
61       
62//      iBin *= delRho*delRho
63//      iBin *= 1e8
64        //
65//      iBin *= (Tscale)^6                      //this puts units on the FFT(magnitude), Vol^2 = T^6
66       
67       
68        // iBin *= phi/vol*(1e-6)^2*(1e8)*T^6
69        // is equivalent to
70        // iBin *= (T/N)^3 * (1e-6)^2*(1e8)
71       
72        iBin *= delRho*delRho
73        iBin *= 1e8
74        iBin *= (Tscale/FFT_N)^3
75       
76end
77
78// for future speedups...
79// 1/2 the time is spent on the Duplicate and redimension step
80// other 1/2 is spent on the FFT
81Function fDoFFT()
82       
83        Wave mat=mat
84        Variable t0,t1
85       
86//      t1 = ticks
87        Duplicate/O mat dmat
88        Redimension/S dmat                      //need single or double precision for FFT
89        //use single precision to use 1/2 memory space - since I'm averaging over the FFT result
90        // I don't need double precision.- the results are identical
91//      print (ticks-t1)/60.15
92       
93        //Set the proper real-space dimensions for the cubes
94        //
95        NVAR Tscale=root:FFT_T
96//      Variable Tscale=5
97        SetScale/P x 0,(Tscale),"", dmat
98        SetScale/P y 0,(Tscale),"", dmat
99        SetScale/P z 0,(Tscale),"", dmat
100        //
101//      t1=ticks
102        FFT dmat                // dmat is now COMPLEX
103//      print (ticks-t1)/60.15
104        // need to multiply the scaling of the FFT result by 2Pi to get proper q-units
105       
106        SetScale/P x 0,(2*pi*DimDelta(dmat, 0)),"", dmat
107        // the FFT rotates zero to the center of the y and z dimensions, so use the DimOffset
108        SetScale/P y (2*pi*DimOffset(dmat, 1)),(2*pi*DimDelta(dmat, 1)),"", dmat
109        SetScale/P z (2*pi*DimOffset(dmat, 2)),(2*pi*DimDelta(dmat, 2)),"", dmat
110       
111End
112
113
114Function fDoCorrection_old(dmat)
115        Wave/C dmat
116
117        //  do the convolution with the cube here
118        // - a lengthy triple for loop
119        // Note that x-dimension is N/2+1 now, since it was a real input
120        // y,z dims are unchanged
121        Variable xDim=DimSize(dmat,0),yDim=DimSize(dmat,1),zDim=DimSize(dmat,2)
122        Variable ii,jj,kk
123        Variable xVal,yVal,zVal
124        Variable xfac,yfac,zfac
125        Variable dimOff0,dimOff1,dimOff2,delta0,delta1,delta2
126       
127        NVAR Tscale=root:FFT_T
128        NVAR Nedge = root:FFT_N
129        if(WaveExists(SincWave) && SameDimensions(note(SincWave)) )
130                //just do the matrix multiplication
131                Wave SincWave=SincWave
132                MultiThread dmat *= SincWave
133                return(0)
134        else
135                //create the cube correction matrix
136                Make/O/N=(xDim,yDim,zDim) SincWave
137               
138                dimOff0 = DimOffset(dmat, 0)
139                dimOff1 = DimOffset(dmat, 1)
140                dimOff2 = DimOffset(dmat, 2)
141                delta0 = DimDelta(dmat,0)
142                delta1 = DimDelta(dmat,1)
143                delta2 = DimDelta(dmat,2)       
144               
145                for(ii=0;ii<xDim;ii+=1)
146                        xVal = dimOff0 + ii *delta0
147                        xVal *= Tscale/2
148                        if(ii==0)
149                                xFac = 1
150                        else
151                                xFac = sinc(xVal)
152                        endif
153                        for(jj=0;jj<yDim;jj+=1)
154                                yVal = dimOff1 + jj *delta1
155                                yVal *= Tscale/2
156                                if(jj==0)
157                                        yFac = 1
158                                else
159                                        yFac = sinc(yVal)
160                                endif
161                                for(kk=0;kk<zDim;kk+=1)
162                                        zVal = dimOff2 + kk *delta2
163                                        zval *= Tscale/2
164                                        if(kk==0)
165                                                zFac = 1
166                                        else
167                                                zFac = sinc(zVal)
168                                        endif
169                                        SincWave[ii][jj][kk] = xfac*yfac*zfac
170                                        dmat[ii][jj][kk] *= (xfac*yfac*zfac)//*(xfac*yfac*zfac)
171                                endfor
172                        endfor
173                endfor
174                Note/K SincWave
175                Note SincWave, "T="+num2str(Tscale)+";N="+num2str(Nedge)+";"
176        endif
177       
178        return(0)
179End
180
181//returns 1 if SincWave (the note) is the same as the matrix dimensions
182//  !! takes the values from the PANEL
183Function SameDimensions(str)
184        String str
185
186        NVAR Tscale=root:FFT_T
187        NVAR Nedge = root:FFT_N
188        if(NumberByKey("T", str ,"=",";") == Tscale && NumberByKey("N", str ,"=",";") == Nedge)
189                return(1)
190        else
191                return(0)
192        endif
193end
194
195Function fGetMagnitude_old(dmat)
196        Wave/C dmat
197        //now get the magnitude, I(q) = |F(q)|^2
198        // do I want magnitude, or magnitude squared? = magsqr(z)
199       
200        MultiThread dmat = r2polar(dmat)
201        Redimension/R dmat                      //just the real part is the magnitude, throw away the complex part
202        MultiThread dmat *= dmat                                //square it
203        //dmat is now a purely real matrix
204               
205End
206
207//dmat is a real matrix at this point
208// (1) use channel sharing (not yet implemented) - then will need to normalize by q^2 rather than nBins
209// --- channel sharing method does not seem to be as good...
210// (2) don't bin points beyond Qmax (done)
211//
212Function fDoBinning(dmat,normalize)
213        Wave dmat
214        Variable normalize
215       
216        Variable xDim=DimSize(dmat,0),yDim=DimSize(dmat,1),zDim=DimSize(dmat,2)
217        Variable ii,jj,kk
218        Variable qX,qY,qZ,qTot
219        Variable binIndex,val
220        NVAR FFT_QmaxReal= root:FFT_QmaxReal
221       
222        Make/O/D/N=(yDim*2) iBin,qBin,nBin
223        qBin[] = 0 + p*DimDelta(dmat, 1)/2              //use bins of 1/2 the width of the reciprocal lattice
224        SetScale/P x,0,DimDelta(dmat,1)/2,"",qBin
225       
226        iBin = 0
227        nBin = 0        //number of intensities added to each bin
228       
229        //loop through everything, and add it all up
230        Variable dimOff0,dimOff1,dimOff2,delta0,delta1,delta2
231        Variable pt,pt1,fra
232        dimOff0 = DimOffset(dmat, 0)
233        dimOff1 = DimOffset(dmat, 1)
234        dimOff2 = DimOffset(dmat, 2)
235        delta0 = DimDelta(dmat,0)
236        delta1 = DimDelta(dmat,1)
237        delta2 = DimDelta(dmat,2)       
238       
239        for(ii=0;ii<xDim;ii+=1)
240                qX = dimOff0 + ii *delta0
241                 
242                for(jj=0;jj<yDim;jj+=1)
243                        qY = dimOff1 + jj *delta1
244                       
245                        for(kk=0;kk<zDim;kk+=1)
246                                qZ = dimOff2 + kk*delta2
247                                qTot = sqrt(qX^2 + qY^2 + qZ^2)
248                                if(qTot < FFT_QmaxReal)         //only take the time to use the good values
249                                // binning method
250                                        binIndex = trunc(x2pnt(qBin, qTot))
251                                        val = dmat[ii][jj][kk]
252                                        if (numType(val)==0)            //count only the good points, ignore Nan or Inf
253                                                iBin[binIndex] += val
254                                                nBin[binIndex] += 1
255                                        endif
256                                //
257                                // channel sharing method (not as good)
258//                                      pt = (qTot - DimOffset(qBin,0))/DimDelta(qBin,0)                //fractional point
259//                                      val = dmat[ii][jj][kk]
260//                                      if (numType(val)==0)            //count only the good points, ignore Nan or Inf
261//                                              pt1 = trunc(pt)
262//                                              fra = pt - pt1
263//                                              iBin[pt1] += (1-fra)*val
264//                                              iBin[pt1+1] += fra*val
265//                                      endif
266                                endif
267                        endfor         
268                endfor
269        endfor
270       
271        iBin /= nBin
272        //normalize, iBin[0] will be NaN so use point[1], or the first "real" value
273        ii=0
274        do
275                if(nBin[ii] == 0 || nbin[ii] == 1)
276                        DeletePoints 0,1, iBin,qBin,nBin                //keep deleting the first point if there were zero or one bins
277                else
278                        val = ibin[ii]
279                        if(numtype(val)==0)
280                                break
281                        endif
282                        ii+=1
283                endif
284        while(ii<numpnts(ibin))
285
286// by channel sharing....(not as good)
287//      iBin /= qBin^2 
288        //delete first two points (q=0 and q=Qmin)
289//      DeletePoints 0,2, iBin,qBin,nBin
290//      val = iBin[0]
291
292        if(normalize)
293                iBin /= val     //then normalize by the first point
294        endif
295               
296        Duplicate/O iBin, iBinAll
297        Duplicate/O qBin, qBinAll
298        Duplicate/O nBin, nBinAll
299        //now truncate qBin at FFT_QmaxReal
300        FindLevel/Q/P qBin, FFT_QmaxReal                //level is reported as a point number, V_LevelX
301        DeletePoints  trunc(V_LevelX)+1, (numpnts(qBin) - trunc(V_LevelX)-1) , iBin,qBin,nBin
302       
303        return(0)
304End
305
306//look for NaN in w1, delete point in w1 and w2
307Function DeleteNaNInf_XY(w1,w2)
308        Wave w1,w2
309       
310        Variable num=numpnts(w1),ii
311       
312        ii=0
313        do
314                if(numtype(w1[ii]) !=0)
315                        //bad point, delete
316                        DeletePoints ii, 1, w1,w2
317                        //don't increment ii, but update (decrement) num
318                        num = numpnts(w1)
319                else
320                        //increment ii
321                        ii += 1
322                endif           
323        while(ii<num)
324End
325
326
327
328Function fDoCorrection(dmat)
329        Wave/C dmat
330
331        //  do the convolution with the cube here
332        // - a lengthy triple for loop
333        // Note that x-dimension is N/2+1 now, since it was a real input
334        // y,z dims are unchanged
335        Variable xDim=DimSize(dmat,0),yDim=DimSize(dmat,1),zDim=DimSize(dmat,2)
336        Variable ii,jj,kk
337        Variable xVal,yVal,zVal
338        Variable xfac,yfac,zfac
339        Variable dimOff0,dimOff1,dimOff2,delta0,delta1,delta2
340       
341        NVAR Tscale=root:FFT_T
342        NVAR Nedge = root:FFT_N
343        if(WaveExists(SincWave) && SameDimensions(note(SincWave)) )
344                //just do the matrix multiplication
345                Wave SincWave=SincWave
346                MatrixOP/C tmp = dmat*SincWave
347               
348                // preserve the scaling
349                SetScale/P x 0,DimDelta(dmat, 0),"", tmp
350                SetScale/P y DimOffset(dmat, 1),DimDelta(dmat, 1),"", tmp
351                SetScale/P z DimOffset(dmat, 2),DimDelta(dmat, 2),"", tmp
352               
353                //dmat=tmp
354                return(0)
355        else
356                //create the cube correction matrix
357                Make/O/N=(xDim,yDim,zDim) SincWave
358               
359                dimOff0 = DimOffset(dmat, 0)
360                dimOff1 = DimOffset(dmat, 1)
361                dimOff2 = DimOffset(dmat, 2)
362                delta0 = DimDelta(dmat,0)
363                delta1 = DimDelta(dmat,1)
364                delta2 = DimDelta(dmat,2)       
365               
366                for(ii=0;ii<xDim;ii+=1)
367                        xVal = dimOff0 + ii *delta0
368                        xVal *= Tscale/2
369                        if(ii==0)
370                                xFac = 1
371                        else
372                                xFac = sinc(xVal)
373                        endif
374                        for(jj=0;jj<yDim;jj+=1)
375                                yVal = dimOff1 + jj *delta1
376                                yVal *= Tscale/2
377                                if(jj==0)
378                                        yFac = 1
379                                else
380                                        yFac = sinc(yVal)
381                                endif
382                                for(kk=0;kk<zDim;kk+=1)
383                                        zVal = dimOff2 + kk *delta2
384                                        zval *= Tscale/2
385                                        if(kk==0)
386                                                zFac = 1
387                                        else
388                                                zFac = sinc(zVal)
389                                        endif
390                                        SincWave[ii][jj][kk] = xfac*yfac*zfac
391                                        dmat[ii][jj][kk] *= (xfac*yfac*zfac)//*(xfac*yfac*zfac)
392                                endfor
393                        endfor
394                endfor
395                Duplicate/O/C dmat,tmp          //slow the first time
396                Note/K SincWave
397                Note SincWave, "T="+num2str(Tscale)+";N="+num2str(Nedge)+";"
398        endif
399       
400        return(0)
401End
402
403Function fGetMagnitude(tmp)
404        Wave/C tmp
405        //now get the magnitude, I(q) = |F(q)|^2
406        // do I want magnitude, or magnitude squared? = magsqr(z)
407       
408        MatrixOP tmp2 = real(r2polar(tmp))
409       
410        MatrixOP/O dmat=tmp2*tmp2
411
412        // preserve the scaling - MatrixOP does NOT
413        SetScale/P x 0,DimDelta(tmp, 0),"", dmat
414        SetScale/P y DimOffset(tmp, 1),DimDelta(tmp, 1),"", dmat
415        SetScale/P z DimOffset(tmp, 2),DimDelta(tmp, 2),"", dmat
416        Killwaves/Z tmp,tmp2
417        //
418        // clean up before exiting
419        //dmat is now a purely real matrix
420        return(0)
421End
422
423
424//      w is a real (2D) matrix at this point
425// -- it comes from a slice of dmat, so it has proper wave scaling
426// (2) don't bin points beyond Qmax (done)
427//
428Function fDoBinning_Scaled2D(w,normalize)
429        Wave w
430        Variable normalize
431       
432        Variable xDim=DimSize(w,0),yDim=DimSize(w,1)
433        Variable ii,jj
434        Variable qX,qY,qTot
435        Variable binIndex,val
436        NVAR FFT_QmaxReal= root:FFT_QmaxReal
437       
438        Make/O/D/N=(yDim*2) iBin_2d,qBin_2d,nBin_2d
439        qBin_2d[] = 0 + p*DimDelta(w, 1)/2              //use bins of 1/2 the width of the reciprocal lattice
440        SetScale/P x,0,DimDelta(w,1)/2,"",qBin_2d
441       
442        iBin_2d = 0
443        nBin_2d = 0     //number of intensities added to each bin
444       
445        //loop through everything, and add it all up
446        Variable dimOff0,dimOff1,delta0,delta1
447        Variable pt,pt1,fra
448        dimOff0 = DimOffset(w, 0)
449        dimOff1 = DimOffset(w, 1)
450        delta0 = DimDelta(w,0)
451        delta1 = DimDelta(w,1)
452       
453        for(ii=0;ii<xDim;ii+=1)
454                qX = dimOff0 + ii *delta0
455                 
456                for(jj=0;jj<yDim;jj+=1)
457                        qY = dimOff1 + jj *delta1
458                       
459//                      qZ = dimOff2 + kk*delta2
460                        qTot = sqrt(qX^2 + qY^2)
461                        if(qTot < FFT_QmaxReal)         //only take the time to use the good values
462                        // binning method
463                                binIndex = trunc(x2pnt(qBin_2d, qTot))
464                                val = w[ii][jj]
465                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
466                                        iBin_2d[binIndex] += val
467                                        nBin_2d[binIndex] += 1
468                                endif
469                        //
470
471                                endif
472                endfor
473        endfor
474       
475        iBin_2d /= nBin_2d
476        //normalize, iBin[0] will be NaN so use point[1], or the first "real" value
477        ii=0
478        do
479                if(nBin_2d[ii] == 0 || nBin_2d[ii] == 1)
480                        DeletePoints 0,1, iBin_2d,qBin_2d,nBin_2d               //keep deleting the first point if there were zero or one bins
481                else
482                        val = iBin_2d[ii]
483                        if(numtype(val)==0)
484                                break
485                        endif
486                        ii+=1
487                endif
488        while(ii<numpnts(iBin_2d))
489
490// by channel sharing....(not as good)
491//      iBin /= qBin^2 
492        //delete first two points (q=0 and q=Qmin)
493//      DeletePoints 0,2, iBin,qBin,nBin
494//      val = iBin[0]
495
496        if(normalize)
497                iBin_2d /= val  //then normalize by the first point
498        endif
499               
500        Duplicate/O iBin_2d, iBinAll_2d
501        Duplicate/O qBin_2d, qBinAll_2d
502        Duplicate/O nBin_2d, nBinAll_2d
503        //now truncate qBin at FFT_QmaxReal
504        FindLevel/Q/P qBin_2d, FFT_QmaxReal             //level is reported as a point number, V_LevelX
505        DeletePoints  trunc(V_LevelX)+1, (numpnts(qBin_2d) - trunc(V_LevelX)-1) , iBin_2d,qBin_2d,nBin_2d
506       
507        return(0)
508End
509
510
511// this is how the 3D viewer in the image Processing macros works
512// to get the correct plane from the 3D matrix
513Function get2DSlice(dmat)
514        wave dmat
515       
516        ImageTransform/G=4 transposeVol dmat
517        WAVE M_VolumeTranspose
518        ImageTransform/P=0 getPlane M_VolumeTranspose
519        WAVE M_ImagePlane
520        Duplicate/O M_ImagePlane detPlane
521        Duplicate/O M_ImagePlane logP
522       
523        NVAR FFT_N=root:FFT_N
524        detPlane[FFT_N/2][FFT_N/2] = NaN                        //hopefully this trims out the singularity at the center
525        logP = log(detPlane)
526       
527        fDoBinning_Scaled2D(detPlane,0)         //last param = normalize y/n
528       
529        KillWaves/Z M_VolumeTranspose,M_ImagePlane
530       
531        return(0)
532end
533
534
535// using a 2D FFT slice, and a folder of real data, interpolate the FFT result to
536// match the q values of the data
537//
538// qz is assumed to be zero.
539//
540Function Interpolate2DSliceToData(folderStr)
541        String folderStr
542
543        Wave calc = root:detPlane
544        Wave data = $("root:"+folderStr+":"+folderStr+"_lin")
545       
546        Duplicate/O data,interp2DSlice                  //keeps the scaling of the data
547       
548        Variable rowOff,colOff,rowDel,colDel
549        rowOff = DimOffset(data,0)
550        colOff = DimOffset(data,1)
551        rowDel = DimDelta(data,0)
552        colDel = DimDelta(data,1)
553
554// see pnt2x for explanation
555//      DimOffset(data, 0) + p *DimDelta(data,0)
556//      DimOffset(data, 1) + q *DimDelta(data,1)
557
558        interp2DSlice = Interp2D (calc, rowOff + p*rowDel, colOff + q*colDel )
559
560        Duplicate/O interp2DSlice interp2DSlice_log
561        interp2DSlice_log = log(interp2DSlice)
562
563        return(0)
564end
Note: See TracBrowser for help on using the repository browser.