# source:sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_FillMatrixShapes.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: 15.6 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2
3// Routines to fill a matrix with spheres in a desired pattern or method
4// - visualization of the 3-d object...
5// also contains routines to convert matrix to XYZ values for export to the calculations
6//
7// to add - more specialized visualization using CreateSlicer
8//
9
10//fills specified sphere
11//***Zeros the rest of the matrix iff yesZero==1
12//
14        WAVE mat
16
17        Variable ii,jj,kk,dik,numx,numy,numz
18
19        numx=DimSize(mat,0)
20        numy=DimSize(mat,1)
21        numz=DimSize(mat,2)
22        for(ii=0;ii<numx;ii+=1)
23                for(jj=0;jj<numy;jj+=1)
24                        for(kk=0;kk<numz;kk+=1)
25                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )
27                                        mat[ii][jj][kk] = 1
28                                else
29                                        if(yesZero==1)
30                                                mat[ii][jj][kk] = 0
31                                        endif
32                                endif
33                        endfor
34                endfor
35        endfor
36
37        return(0)
38End
39
40//grid is the size of each voxel in the matrix in angstroms (grid^3 is the volume)
41//radius is the real-space radius that you want the sphere to be
42//
43// fills sphere if fill=1, clears if fill=0
45        WAVE mat
47
48        Variable ii,jj,kk,dik,numx,numy,numz,numVox
49        Variable x1,x2,y1,y2,z1,z2
50        numVox= rad/grid*3              //*3 for extra room, *2 is diam
51
52        numx=DimSize(mat,0)
53        numy=DimSize(mat,1)
54        numz=DimSize(mat,2)
55
56// keep the bounding box within bounds
57        x1 = (xc >= numVox) ? (xc-numVox) : 0
58        x2 = (xc + numVox) <= numx ? (xc + numVox) : numx
59        y1 = (yc >= numVox) ? (yc-numVox) : 0
60        y2 = (yc + numVox) <= numy ? (yc + numVox) : numy
61        z1 = (zc >= numVox) ? (zc-numVox) : 0
62        z2 = (zc + numVox) <= numz ? (zc + numVox) : numz
63//      print numvox
64//      Print x1,x2
65//      Print y1,y2
66//      print z1,z2
67
68        for(ii=x1;ii<x2;ii+=1)
69                for(jj=y1;jj<y2;jj+=1)
70                        for(kk=z1;kk<z2;kk+=1)
71                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
73                                        mat[ii][jj][kk] = fill
74//                              else
75//                                      if(yesZero==1)
76//                                              mat[ii][jj][kk] = 0
77//                                      endif
78                                endif
79                        endfor
80                endfor
81        endfor
82
83        return(0)
84End
85
86//grid is the size of each voxel in the matrix in angstroms (grid^3 is the volume)
87//radius is the real-space radius that you want the sphere to be
88//
89// fills sphere with fill value if the space is clear, that is there is nothing
90// there but solvent.
91//
92// actually runs 2x to make sure that there is no overlap
93//
95        WAVE mat
97
98        Variable ii,jj,kk,dik,numx,numy,numz,numVox
99        Variable x1,x2,y1,y2,z1,z2
100
101        NVAR FFT_SolventSLD = root:FFT_SolventSLD
102
103        numVox= rad/grid*3              //*3 for extra room, *2 is diam
104
105        numx=DimSize(mat,0)
106        numy=DimSize(mat,1)
107        numz=DimSize(mat,2)
108
109// keep the bounding box within bounds
110        x1 = (xc >= numVox) ? (xc-numVox) : 0
111        x2 = (xc + numVox) <= numx ? (xc + numVox) : numx
112        y1 = (yc >= numVox) ? (yc-numVox) : 0
113        y2 = (yc + numVox) <= numy ? (yc + numVox) : numy
114        z1 = (zc >= numVox) ? (zc-numVox) : 0
115        z2 = (zc + numVox) <= numz ? (zc + numVox) : numz
116
117        Variable err=0
118// check for an overlapping point in the sphere to fill
119        for(ii=x1;ii<x2;ii+=1)
120                for(jj=y1;jj<y2;jj+=1)
121                        for(kk=z1;kk<z2;kk+=1)
122                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
123                                if( (dik<=rad) && (mat[ii][jj][kk] != FFT_SolventSLD) )                 //if in sphere to fill AND already occupied
124                                        err += 1
125                                        return err              //get out now
126                                endif
127                        endfor
128                endfor
129        endfor
130
131        // OK, go ahead and fill
132        for(ii=x1;ii<x2;ii+=1)
133                for(jj=y1;jj<y2;jj+=1)
134                        for(kk=z1;kk<z2;kk+=1)
135                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
137                                        mat[ii][jj][kk] = fill
138                                endif
139                        endfor
140                endfor
141        endfor
142        return(0)
143End
144
145// fills the sphere at the specified center and radius
146// the input matrix is Nx3 and will be extended by this routine
147// as needed to add the sphere
148//
150        WAVE mat3
152
153        Variable ii,jj,kk,dik,np
154
155        np = DimSize(mat3,0)
159                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )
161                                        InsertPoints/M=0 np, 1, mat3
162//                                      InsertPoints/M=1 np, 1, mat3
163//                                      InsertPoints/M=2 np, 1, mat3
164                                        mat3[np][0] = ii
165                                        mat3[np][1] = jj
166                                        mat3[np][2] = kk
167                                        np += 1                         //one point has been added
168                                else
169                                        if(yesZero==1)
170                                                //mat[ii][jj][kk] = 0
171                                        endif
172                                endif
173                        endfor
174                endfor
175        endfor
176
177        return(0)
178End
179
180//erases specified sphere
181//
182//obsolete
184//      WAVE mat
186//
187//      Variable ii,jj,kk,dik,numx,numy,numz
188//
189//      numx=DimSize(mat,0)
190//      numy=DimSize(mat,1)
191//      numz=DimSize(mat,2)
192//      for(ii=0;ii<numx;ii+=1)
193//              for(jj=0;jj<numy;jj+=1)
194//                      for(kk=0;kk<numz;kk+=1)
195//                              dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )
197//                                      mat[ii][jj][kk] = 0
198//                              endif
199//                      endfor
200//              endfor
201//      endfor
202//
203//      return(0)
204//End
205
206//erases specified sphere
207//
208//obsolete
210//      WAVE mat
212//
213//      Variable ii,jj,kk,dik,numx,numy,numz
214//
215//      numx=DimSize(mat,0)
216//      numy=DimSize(mat,1)
217//      numz=DimSize(mat,2)
218//      for(ii=0;ii<numx;ii+=1)
219//              for(jj=0;jj<numy;jj+=1)
220//                      for(kk=0;kk<numz;kk+=1)
221//                              dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
223//                                      mat[ii][jj][kk] = 0
224//                              endif
225//                      endfor
226//              endfor
227//      endfor
228//
229//      return(0)
230//End
231
232//fill = 1, erase = 0
233//
234// !! this does not check for overlap!!!
236        Wave mat
238
239        Wave x3d=x3d
240        Wave y3d=y3d
241        Wave z3d=z3d
242
245        NVAR  grid=root:FFT_T
246
247        for(ii=0;ii<num;ii+=1)
248                //print "sphere ",ii
250                if(pd !=0 )
252                endif
254        endfor
255
256End
257
258//fills cylinder (fill == 1) or erases it (fill == 0)
259//
260// xc,yc,zc is the center
261// grid is the grid spacing (A)
263// len is the desired length, in the z-direction
265        WAVE mat
267
268        Variable ii,pts
269        //put half the length above - half below
270        pts = trunc(len/2/grid)
271        for(ii=(zc-pts);ii<(zc+pts);ii+=1)
273        endfor
274
275        return(0)
276End
277
278//fills the specified circle (fill==1), or erases it (fill == 0)
279//zloc assumed to be a valid z-index of the matrix
281        WAVE mat
283
284        Variable ii,jj,dik,numx,numy,numz
285
286        numx=DimSize(mat,0)
287        numy=DimSize(mat,1)
288        for(ii=0;ii<numx;ii+=1)
289                for(jj=0;jj<numy;jj+=1)
290                                dik=sqrt( (ii-xc)^2+(jj-yc)^2)*grid
292                                        mat[ii][jj][zloc] = fill
293                                endif
294                endfor
295        endfor
296
297        return(0)
298End
299
300//fills cylinder (fill == 1) or erases it (fill == 0)
301//
302// xc,yc,zc is the center
303// grid is the grid spacing (A)
305// len is the desired length, in the X-direction
307        WAVE mat
309
310        Variable ii,pts
311        //put half the length above - half below
312        pts = trunc(len/2/grid)
313        for(ii=(xc-pts);ii<(xc+pts);ii+=1)
315        endfor
316
317        return(0)
318End
319
320//fills the specified circle (fill==1), or erases it (fill == 0)
321//xloc assumed to be a valid x-index of the matrix
323        WAVE mat
325
326        Variable ii,jj,dik,numx,numy,numz
327
328        numy=DimSize(mat,0)
329        numz=DimSize(mat,1)
330        for(ii=0;ii<numy;ii+=1)
331                for(jj=0;jj<numz;jj+=1)
332                                dik=sqrt( (ii-yc)^2+(jj-zc)^2)*grid
334                                        mat[xloc][ii][jj] = fill
335                                endif
336                endfor
337        endfor
338
339        return(0)
340End
341
342
343//fills cylinder (fill == 1) or erases it (fill == 0)
344//
345// xc,yc,zc is the center
346// grid is the grid spacing (A)
348// len is the desired length, in the Y-direction
350        WAVE mat
352
353        Variable ii,pts
354        //put half the length above - half below
355        pts = trunc(len/2/grid)
356        for(ii=(yc-pts);ii<(yc+pts);ii+=1)
358        endfor
359
360        return(0)
361End
362
363//fills the specified circle (fill==1), or erases it (fill == 0)
364//yloc assumed to be a valid y-index of the matrix
366        WAVE mat
368
369        Variable ii,jj,dik,numx,numy,numz
370
371        numx=DimSize(mat,0)
372        numz=DimSize(mat,1)
373        for(ii=0;ii<numx;ii+=1)
374                for(jj=0;jj<numz;jj+=1)
375                                dik=sqrt( (ii-xc)^2+(jj-zc)^2)*grid
377                                        mat[ii][yloc][jj] = fill
378                                endif
379                endfor
380        endfor
381
382        return(0)
383End
384
385// uses the Sobol' sequence to fill the points
386// - likely more "evenly" spread out than a random pick
387// -- but is it too regular?
388Function SobolFill3DMat(mat,num)
389        Wave mat
390        variable num            //number of spheres to add
391
392        Variable row,col,lay,ii,xt,yt,zt,fail=0
393
394        Make/O/D/N=3 Sobol3D
395
396        // initialize. XOP should prevent re-initialization
397        SobolX(-1,Sobol3D)
398
399        row=DimSize(mat,0)
400        col=DimSize(mat,1)
401        lay=DimSize(mat,2)
402        ii=0
403        do
404                SobolX(3,Sobol3D)
405                xt = trunc(row*Sobol3D[0])
406                yt = trunc(col*Sobol3D[1])
407                zt = trunc(lay*Sobol3D[2])
408                if( NoNeighbors(xt,yt,zt,mat) )
409                        mat[xt][yt][zt] = 1
410                        ii+=1           //increment number of spheres actually added
411                        //Print "point ",ii
412                else
413                        fail +=1
414                endif
415        while(ii<num)
416//      Print "failures = ",fail
417
418//      ParseMatrix3D_rho(mat)          //convert to XYZ
419
420        return(0)
421End
422
423
424Function RandomFill3DMat(mat,num)
425        Wave mat
426        variable num            //number of spheres to add
427
428        Variable row,col,lay,ii,xt,yt,zt,fail=0
429
430        row=DimSize(mat,0)
431        col=DimSize(mat,1)
432        lay=DimSize(mat,2)
433        //place first 2 points
434        //xt=trunc(abs(enoise(row)))            //random integer distr betw (0,row-1)
435        //yt=trunc(abs(enoise(col)))
436        //zt=trunc(abs(enoise(lay)))
437        //mat[xt][yt][1] = 1
438        //mat[xt][yt][2] = 1
439        //loop over remaining spheres to add
440        ii=0
441        do
442                xt=trunc(abs(enoise(row)))              //distr betw (0,npt)
443                yt=trunc(abs(enoise(col)))
444                zt=trunc(abs(enoise(lay)))
445                if( NoNeighbors(xt,yt,zt,mat) )
446                        mat[xt][yt][zt] = 1
447                        ii+=1           //increment number of spheres actually added
448                        //Print "point ",ii
449                else
450                        fail +=1
451                endif
452        while(ii<num)
453        Print "failures = ",fail
454
455//      ParseMatrix3D_rho(mat)          //convert to XYZ
456
457        return(0)
458End
459
460//point is considered valid if (x,y,z) is unoccupied
461// AND there is a NO neighbor (no connectivity)
462//
463// could speed this up by immediately bouncing out when the first neighbor is found
464Function NoNeighbors(xt,yt,zt,mat)
465        Variable xt,yt,zt
466        WAVE mat
467        //returns 1 if valid, 0 if not valid
468        Variable valid=1,invalid=0,xMax,yMax,zMax,ncont=0
469
470        if(mat[xt][yt][zt])             //already occupied, trial is invalid
471                return (invalid)
472        Endif
473
474        xMax=DimSize(mat,0) - 1                 //max index of the array
475        yMax=DimSize(mat,1) - 1
476        zMax=DimSize(mat,2) - 1
477        //connected?
478        //any of the surrounding 26 points occupied?
479        Variable xi,yi,zi,ii,jj,kk
480        for(ii=-1;ii<=1;ii+=1)
481                xi = xt+ii
482                xi = (xi< 0) ? 0 : xi           //keep in bounds of matrix index
483                xi = (xi>xMax) ? xMax : xi
484                for(jj=-1;jj<=1;jj+=1)
485                        yi = yt + jj
486                        yi = (yi<0) ? 0 : yi
487                        yi = (yi>yMax) ? yMax : yi
488                        for(kk=-1;kk<=1;kk+=1)
489                                zi = zt+kk
490                                zi = (zi<0) ? 0 : zi
491                                zi = (zi>zMax) ? zMax : zi
492                                if( (xi==xt) && (yi==yt) && (zi==zt) )
493                                        //at central point, do nothing
494                                else
495                                        if(mat[xi][yi][zi])
496                                                ncont+=1
497                                        endif
498                                endif
499                        endfor
500                endfor
501        endfor
502
503        //if((ncont > 0) && (ncont <2) )        // only one neighbor
504        if(ncont==0)                    // no nearest neighbor contact points
505                return(valid)
506        else
507                return(invalid)
508        endif
509End
510
511Proc ParseMatrix3DToXYZ(matStr)
512        String matStr="mat"
513        ParseMatrix3D_rho(\$matStr)
514End
515
516
517///
518/// Depricated - use ParseMatrix3D_rho() instead, to get the SLD information
519///
520///
521// parses the 3d matrix to get xyz coordinates
522// CHANGE - ? Make as byte waves to save space
523//
524// XYZ waves are hard-wired to be "x3d,y3d,z3d"
525// overwrites previous waves
526//
527//Function ParseMatrix3D(mat)
528//      Wave mat
529//
530//      Variable nptx,npty,nptz,ii,jj,kk,num
531//
532//      nptx=DimSize(mat,0)
533//      npty=DimSize(mat,1)
534//      nptz=DimSize(mat,2)
535//      Make/O/D/N=0 x3d,y3d,z3d                //ensure that the waves are SP
536//      num=0
537//
538//      for(ii=0;ii<nptx;ii+=1)
539//              for(jj=0;jj<npty;jj+=1)
540//                      for(kk=0;kk<nptz;kk+=1)
541//                              if(mat[ii][jj][kk])
542//                                      num+=1
543//                                      InsertPoints num,1,x3d,y3d,z3d
544//                                      x3d[num-1]=ii
545//                                      y3d[num-1]=jj
546//                                      z3d[num-1]=kk
547//                              endif
548//                      endfor
549//              endfor
550//      endfor
551//
552//      return(0)
553//End
554
555
556// parses the 3d matrix to get xyz coordinates
557// CHANGE - ? Make as byte waves to save space
558//
559// XYZ waves are hard-wired to be "x3d,y3d,z3d"
560// overwrites previous waves
561//
562Function ParseMatrix3D_rho(mat)
563        Wave mat
564
565        Variable nptx,npty,nptz,ii,jj,kk,num
566
567        NVAR    solventSLD = root:FFT_SolventSLD
568
569        nptx=DimSize(mat,0)
570        npty=DimSize(mat,1)
571        nptz=DimSize(mat,2)
572        Make/O/D/N=0 x3d,y3d,z3d,rho3d          //ensure that the waves are DP
573        num=0
574
575        for(ii=0;ii<nptx;ii+=1)
576                for(jj=0;jj<npty;jj+=1)
577                        for(kk=0;kk<nptz;kk+=1)
578                                if(mat[ii][jj][kk]!=solventSLD)
579                                        num+=1
580                                        InsertPoints num,1,x3d,y3d,z3d,rho3d
581                                        x3d[num-1]=ii
582                                        y3d[num-1]=jj
583                                        z3d[num-1]=kk
584                                        rho3d[num-1]=mat[ii][jj][kk]
585                                endif
586                        endfor
587                endfor
588        endfor
589
590        return(0)
591End
592
593Function ConvertXYZto3N(xw,yw,zw)
594        Wave xw,yw,zw
595
596        Variable num=numpnts(xw)
597        Make/O/D/N=(num,3) matGiz
598        Concatenate/O {x3d,y3d,z3d},matGiz              //Igor 5+ only
599//
600// or
601//      matGiz[][0] = xw[p]
602//      matGiz[][1] = yw[p]
603//      matGiz[][2] = zw[p]
604
605        return(0)
606End
607
608//////////////////
609
610//pad the matrix with zeros to make it bigger, ?maybe avoid errors
611// of box size
612// (this does not yet work.....)
613//
615        Wave mat
616        Variable num
617
618        Variable numold=DimSize(mat,0)
619        Make/O/B/U/N=(numold+2*num,numold+2*num,numold+2*num) matNew=0
620        matNew[num,numold+num-1][num,numold+num-1][num,numold+num-1]=mat[p-num][q-num][r-num]
621end
622
623// takes values at XYZ and converts to a 3D byte volume
624// - each voxel may have different values
625//
626// xyz values are to start from 0->
627// voxel values must be integer, or they will end up byte
628//
629Function XYZV_toByteVoxels(xx,yy,zz,val)
630        Wave xx,yy,zz,val
631
632        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt
633
634        WaveStats/Q xx
635        x1 = V_max
636
637        WaveStats/Q yy
638        y1 = V_max
639
640        WaveStats/Q zz
641        z1 = V_max
642
643        //get the maximum dimension
644        npt = max(x1,y1)
645        npt = max(npt,z1)
646
647        Make/O/B/U/N=(npt,npt,npt) VoxelMat=0
648
649        num=numpnts(xx)
650
651        for(ii=0;ii<num;ii+=1)
652                VoxelMat[xx[ii]][yy[ii]][zz[ii]] = val[ii]
653        endfor
654
655        return(0)
656End
657
658// takes XYZV values such as output from  Ken Rubinson's converter
659// where xyz may take negative values, and tries to put it into
660// mat, shifting XYZ as needed
661//
662Function XYZV_FillMat(xx,yy,zz,val,erase)
663        Wave xx,yy,zz,val
664        Variable erase
665
666        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt,minp,maxp
667
668        WAVE mat = root:mat
669
670        if(erase)
671                mat = 0
672        endif
673
674        WaveStats/Q xx
675        x1 = V_min
676
677        WaveStats/Q yy
678        y1 = V_min
679
680        WaveStats/Q zz
681        z1 = V_min
682
683
684        // find the minimum xyz
685        minp = min(x1,y1)
686        minp = min(minp,z1)
687//      print "minp = ",minp
688
689        if(minp < 0)
690                minp = abs(minp)
691        else
692                minp = 0                // if all of the values are positive, don't shift anything
693        endif
694
695
696        // will the adjusted xyz values fit the mat?
697        Variable matp = DimSize(mat,0)          // assume all 3 dimensions are the same
698
699        WaveStats/Q xx
700        x1 = V_max
701
702        WaveStats/Q yy
703        y1 = V_max
704
705        WaveStats/Q zz
706        z1 = V_max
707
708        // find the minimum xyz
709        maxp = max(x1,y1)
710        maxp = max(maxp,z1)
711//      Print "maxp = ",maxp
712
713        if(x1+minp > matp-1)
714                Print "Mat is not big enough. x1+minp = ",x1+minp
715                return(0)
716        Endif
717        if(y1+minp > matp-1)
718                Print "Mat is not big enough. y1+minp = ",y1+minp
719                return(0)
720        Endif
721        if(z1+minp > matp-1)
722                Print "Mat is not big enough. z1+minp = ",z1+minp
723                return(0)
724        Endif
725
726        num=numpnts(xx)
727
728        for(ii=0;ii<num;ii+=1)
729                mat[xx[ii]+minp][yy[ii]+minp][zz[ii]+minp] = val[ii]
730        endfor
731
732        return(0)
733End
Note: See TracBrowser for help on using the repository browser.