source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_FillMatrixShapes.ipf @ 953

Last change on this file since 953 was 953, checked in by srkline, 8 years ago

some minor changes to the FFT filling procedures to help out Abhiram

Replaced cansasXML.ipf v1.11 with v1.12 (current version)

some restructuring of the VSANS file names.

File size: 19.5 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// random point fills are also here, either random or Sobol
8//
9
10//fills specified sphere
11//***Zeros the rest of the matrix iff yesZero==1
12//
13Function FillSphere(mat,rad,xc,yc,zc,yesZero)
14        WAVE mat
15        Variable rad,xc,yc,zc,yesZero
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 )
26                                if(dik<=rad)
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
44Function FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
45        WAVE mat
46        Variable grid,rad,xc,yc,zc,fill
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
72                                if(dik<=rad)
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// Now enforces periodic boundary conditions
90Function FillSphereRadiusPeriodic(mat,grid,rad,xc,yc,zc,fill)
91        WAVE mat
92        Variable grid,rad,xc,yc,zc,fill
93       
94        Variable ii,jj,kk,dik,numVox
95        Variable x1,y1,z1
96        Variable nDimx,nDimy,nDimz
97
98        numVox= round(rad/grid*3)               //*3 for extra room, *2 is diam
99       
100        nDimx=DimSize(mat,0)-1
101        nDimy=DimSize(mat,1)-1
102        nDimz=DimSize(mat,2)-1
103       
104// keep the bounding box within bounds 
105//      x1 = (xc >= numVox) ? (xc-numVox) : 0
106//      x2 = (xc + numVox) <= numx ? (xc + numVox) : numx
107//      y1 = (yc >= numVox) ? (yc-numVox) : 0
108//      y2 = (yc + numVox) <= numy ? (yc + numVox) : numy
109//      z1 = (zc >= numVox) ? (zc-numVox) : 0
110//      z2 = (zc + numVox) <= numz ? (zc + numVox) : numz
111
112
113        for(ii=xc-numVox;ii<xc+numVox;ii+=1)
114                for(jj=yc-numVox;jj<yc+numVox;jj+=1)
115                        for(kk=zc-numVox;kk<zc+numVox;kk+=1)
116                                //calculate the distance, even if it's off the matrix
117                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
118                                if(dik<=rad)
119                                        //then if it's part of the sphere, get the voxel in bounds before filling
120                                        x1=ii
121                                        y1=jj
122                                        z1=kk
123                                        if( (ii>nDimx) || (ii<0) )
124                                                x1 = abs(nDimx+1 - abs(x1))
125                                        endif
126                                        if( (jj>nDimy) || (jj<0) )
127                                                y1 = abs(nDimy+1 - abs(y1))
128                                        endif
129                                        if( (kk>nDimz) || (kk<0) )
130                                                z1 = abs(nDimz+1 - abs(z1))
131                                        endif
132                               
133                                        mat[x1][y1][z1] = fill
134
135
136                                endif           //if dist < rad
137                        endfor
138                endfor
139        endfor
140       
141        return(0)       
142End
143
144//grid is the size of each voxel in the matrix in angstroms (grid^3 is the volume)
145//radius is the real-space radius that you want the sphere to be
146//
147// fills sphere with fill value if the space is clear, that is there is nothing
148// there but solvent.
149//
150// actually runs 2x to make sure that there is no overlap
151//
152Function FillSphereRadiusNoOverlap(mat,grid,rad,xc,yc,zc,fill)
153        WAVE mat
154        Variable grid,rad,xc,yc,zc,fill
155       
156        Variable ii,jj,kk,dik,numx,numy,numz,numVox
157        Variable x1,x2,y1,y2,z1,z2
158       
159        NVAR FFT_SolventSLD = root:FFT_SolventSLD
160
161        numVox= rad/grid*3              //*3 for extra room, *2 is diam
162       
163        numx=DimSize(mat,0)
164        numy=DimSize(mat,1)
165        numz=DimSize(mat,2)
166
167// keep the bounding box within bounds 
168        x1 = (xc >= numVox) ? (xc-numVox) : 0
169        x2 = (xc + numVox) <= numx ? (xc + numVox) : numx
170        y1 = (yc >= numVox) ? (yc-numVox) : 0
171        y2 = (yc + numVox) <= numy ? (yc + numVox) : numy
172        z1 = (zc >= numVox) ? (zc-numVox) : 0
173        z2 = (zc + numVox) <= numz ? (zc + numVox) : numz
174
175        Variable err=0
176// check for an overlapping point in the sphere to fill
177        for(ii=x1;ii<x2;ii+=1)
178                for(jj=y1;jj<y2;jj+=1)
179                        for(kk=z1;kk<z2;kk+=1)
180                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
181                                if( (dik<=rad) && (mat[ii][jj][kk] != FFT_SolventSLD) )                 //if in sphere to fill AND already occupied
182                                        err += 1
183                                        return err              //get out now
184                                endif
185                        endfor
186                endfor
187        endfor
188       
189        // OK, go ahead and fill
190        for(ii=x1;ii<x2;ii+=1)
191                for(jj=y1;jj<y2;jj+=1)
192                        for(kk=z1;kk<z2;kk+=1)
193                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
194                                if(dik<=rad)
195                                        mat[ii][jj][kk] = fill
196                                endif
197                        endfor
198                endfor
199        endfor
200        return(0)       
201End
202
203// fills the sphere at the specified center and radius
204// the input matrix is Nx3 and will be extended by this routine
205// as needed to add the sphere
206//
207Function FillSphere3(mat3,rad,xc,yc,zc,yesZero)
208        WAVE mat3
209        Variable rad,xc,yc,zc,yesZero
210       
211        Variable ii,jj,kk,dik,np
212       
213        np = DimSize(mat3,0)
214        for(ii=xc-rad;ii<xc+rad;ii+=1)
215                for(jj=yc-rad;jj<yc+rad;jj+=1)
216                        for(kk=zc-rad;kk<zc+rad;kk+=1)
217                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )
218                                if(dik<=rad)
219                                        InsertPoints/M=0 np, 1, mat3
220//                                      InsertPoints/M=1 np, 1, mat3
221//                                      InsertPoints/M=2 np, 1, mat3
222                                        mat3[np][0] = ii
223                                        mat3[np][1] = jj
224                                        mat3[np][2] = kk
225                                        np += 1                         //one point has been added
226                                else
227                                        if(yesZero==1)
228                                                //mat[ii][jj][kk] = 0
229                                        endif
230                                endif
231                        endfor
232                endfor
233        endfor
234       
235        return(0)       
236End
237
238//erases specified sphere
239//
240//obsolete
241//Function EraseSphere(mat,rad,xc,yc,zc)
242//      WAVE mat
243//      Variable rad,xc,yc,zc
244//     
245//      Variable ii,jj,kk,dik,numx,numy,numz
246//     
247//      numx=DimSize(mat,0)
248//      numy=DimSize(mat,1)
249//      numz=DimSize(mat,2)
250//      for(ii=0;ii<numx;ii+=1)
251//              for(jj=0;jj<numy;jj+=1)
252//                      for(kk=0;kk<numz;kk+=1)
253//                              dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )
254//                              if(dik<=rad)
255//                                      mat[ii][jj][kk] = 0
256//                              endif
257//                      endfor
258//              endfor
259//      endfor
260//     
261//      return(0)       
262//End
263
264//erases specified sphere
265//
266//obsolete
267//Function EraseSphereRadius(mat,grid,rad,xc,yc,zc)
268//      WAVE mat
269//      Variable grid,rad,xc,yc,zc
270//     
271//      Variable ii,jj,kk,dik,numx,numy,numz
272//     
273//      numx=DimSize(mat,0)
274//      numy=DimSize(mat,1)
275//      numz=DimSize(mat,2)
276//      for(ii=0;ii<numx;ii+=1)
277//              for(jj=0;jj<numy;jj+=1)
278//                      for(kk=0;kk<numz;kk+=1)
279//                              dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )*grid
280//                              if(dik<=rad)
281//                                      mat[ii][jj][kk] = 0
282//                              endif
283//                      endfor
284//              endfor
285//      endfor
286//     
287//      return(0)       
288//End
289
290//fill = 1, erase = 0
291//
292// !! this does not check for overlap!!!
293Function SphereAtEachPoint(mat,rad,pd,fill)
294        Wave mat
295        Variable rad,pd,fill
296       
297        Wave x3d=x3d
298        Wave y3d=y3d
299        Wave z3d=z3d
300       
301        Variable ii=0,num=numpnts(x3d),meanRad
302        meanRad=rad
303        NVAR  grid=root:FFT_T
304       
305        for(ii=0;ii<num;ii+=1)
306                //print "sphere ",ii
307                //FillSphere(mat,rad,x3d[ii],y3d[ii],z3d[ii],0)
308                if(pd !=0 )
309                        rad = meanRad + gnoise(pd*meanRad)
310                endif
311                FillSphereRadiusPeriodic(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],fill)
312        endfor 
313       
314End
315
316//fills cylinder (fill == 1) or erases it (fill == 0)
317//
318// xc,yc,zc is the center
319// grid is the grid spacing (A)
320// rad is the desired radius
321// len is the desired length, in the z-direction
322Function FillZCylinder(mat,grid,rad,xc,yc,zc,len,fill)
323        WAVE mat
324        Variable grid,rad,xc,yc,zc,len,fill
325       
326        Variable ii,pts
327        //put half the length above - half below
328        pts = trunc(len/2/grid)
329        for(ii=(zc-pts);ii<(zc+pts);ii+=1)
330                FillXYCircle(mat,grid,rad,xc,yc,ii,fill)
331        endfor
332       
333        return(0)       
334End
335
336//fills the specified circle (fill==1), or erases it (fill == 0)
337//zloc assumed to be a valid z-index of the matrix
338Function FillXYCircle(mat,grid,rad,xc,yc,zloc,fill)
339        WAVE mat
340        Variable grid,rad,xc,yc,zloc,fill
341       
342        Variable ii,jj,dik,numx,numy,numz
343       
344        numx=DimSize(mat,0)
345        numy=DimSize(mat,1)
346       
347        Variable x1,x2,y1,y2
348        x1 = xc - trunc(rad/grid) - 2
349        x2 = xc + trunc(rad/grid) + 2
350        y1 = yc - trunc(rad/grid) - 2
351        y2 = yc + trunc(rad/grid) + 2
352       
353        y1 = y1 < 0 ? 0 : y1
354        x1 = x1 < 0 ? 0 : x1
355       
356        y2 = y2 > numy ? numy : y2
357        x2 = x2 > numx ? numx : x2
358       
359        for(ii=x1;ii<x2;ii+=1)
360                for(jj=y1;jj<y2;jj+=1)
361                                dik=sqrt( (ii-xc)^2+(jj-yc)^2)*grid
362                                if(dik<=rad)
363                                        mat[ii][jj][zloc] = fill
364                                endif
365                endfor
366        endfor
367       
368//      for(ii=0;ii<numx;ii+=1)
369//              for(jj=0;jj<numy;jj+=1)
370//                              dik=sqrt( (ii-xc)^2+(jj-yc)^2)*grid
371//                              if(dik<=rad)
372//                                      mat[ii][jj][zloc] = fill
373//                              endif
374//              endfor
375//      endfor
376       
377        return(0)       
378End
379
380//fills cylinder (fill == 1) or erases it (fill == 0)
381//
382// xc,yc,zc is the center
383// grid is the grid spacing (A)
384// rad is the desired radius
385// len is the desired length, in the X-direction
386Function FillXCylinder(mat,grid,rad,xc,yc,zc,len,fill)
387        WAVE mat
388        Variable grid,rad,xc,yc,zc,len,fill
389
390//      Variable tref = startMSTimer
391               
392        Variable ii,pts
393        //put half the length above - half below
394        pts = trunc(len/2/grid)
395        for(ii=(xc-pts);ii<(xc+pts);ii+=1)
396                FillYZCircle(mat,grid,rad,yc,zc,ii,fill)
397        endfor
398
399//      Variable ms = stopMSTimer(tref)
400//      print "Time elapsed = ", ms/1e6, "s"
401               
402        return(0)       
403End
404
405//fills the specified circle (fill==1), or erases it (fill == 0)
406//xloc assumed to be a valid x-index of the matrix
407Function FillYZCircle(mat,grid,rad,yc,zc,xloc,fill)
408        WAVE mat
409        Variable grid,rad,yc,zc,xloc,fill
410       
411        Variable ii,jj,dik,numx,numy,numz
412       
413        numy=DimSize(mat,0)
414        numz=DimSize(mat,1)
415       
416        Variable y1,y2,z1,z2
417        y1 = yc - trunc(rad/grid) - 2
418        y2 = yc + trunc(rad/grid) + 2
419        z1 = zc - trunc(rad/grid) - 2
420        z2 = zc + trunc(rad/grid) + 2
421       
422        y1 = y1 < 0 ? 0 : y1
423        z1 = z1 < 0 ? 0 : z1
424       
425        y2 = y2 > numy ? numy : y2
426        z2 = z2 > numz ? numz : z2
427       
428
429        for(ii=y1;ii<y2;ii+=1)
430                for(jj=z1;jj<z2;jj+=1)
431                                dik=sqrt( (ii-yc)^2+(jj-zc)^2)*grid
432                                if(dik<=rad)
433                                        mat[xloc][ii][jj] = fill
434                                endif
435                endfor
436        endfor
437
438
439////slow way           
440//      for(ii=0;ii<numy;ii+=1)
441//              for(jj=0;jj<numz;jj+=1)
442//                              dik=sqrt( (ii-yc)^2+(jj-zc)^2)*grid
443//                              if(dik<=rad)
444//                                      mat[xloc][ii][jj] = fill
445//                              endif
446//              endfor
447//      endfor
448
449        return(0)       
450End
451
452
453//fills cylinder (fill == 1) or erases it (fill == 0)
454//
455// xc,yc,zc is the center
456// grid is the grid spacing (A)
457// rad is the desired radius
458// len is the desired length, in the Y-direction
459Function FillYCylinder(mat,grid,rad,xc,yc,zc,len,fill)
460        WAVE mat
461        Variable grid,rad,xc,yc,zc,len,fill
462       
463        Variable ii,pts
464        //put half the length above - half below
465        pts = trunc(len/2/grid)
466        for(ii=(yc-pts);ii<(yc+pts);ii+=1)
467                FillXZCircle(mat,grid,rad,xc,zc,ii,fill)
468        endfor
469       
470        return(0)       
471End
472
473//fills the specified circle (fill==1), or erases it (fill == 0)
474//yloc assumed to be a valid y-index of the matrix
475Function FillXZCircle(mat,grid,rad,xc,zc,yloc,fill)
476        WAVE mat
477        Variable grid,rad,xc,zc,yloc,fill
478       
479        Variable ii,jj,dik,numx,numy,numz
480       
481        numx=DimSize(mat,0)
482        numz=DimSize(mat,1)
483       
484        Variable x1,x2,z1,z2
485        x1 = xc - trunc(rad/grid) - 2
486        x2 = xc + trunc(rad/grid) + 2
487        z1 = zc - trunc(rad/grid) - 2
488        z2 = zc + trunc(rad/grid) + 2
489
490       
491        z1 = z1 < 0 ? 0 : z1
492        x1 = x1 < 0 ? 0 : x1
493       
494        z2 = z2 > numz ? numz : z2
495        x2 = x2 > numx ? numx : x2
496       
497        for(ii=x1;ii<x2;ii+=1)
498                for(jj=z1;jj<z2;jj+=1)
499                                dik=sqrt( (ii-xc)^2+(jj-zc)^2)*grid
500                                if(dik<=rad)
501                                        mat[ii][yloc][jj] = fill
502                                endif
503                endfor
504        endfor
505       
506//      for(ii=0;ii<numx;ii+=1)
507//              for(jj=0;jj<numz;jj+=1)
508//                              dik=sqrt( (ii-xc)^2+(jj-zc)^2)*grid
509//                              if(dik<=rad)
510//                                      mat[ii][yloc][jj] = fill
511//                              endif
512//              endfor
513//      endfor
514       
515        return(0)       
516End
517
518// uses the Sobol' sequence to fill the points
519// - likely more "evenly" spread out than a random pick
520// -- but is it too regular?
521Function SobolFill3DMat(mat,num,fill)
522        Wave mat
523        variable num,fill               //number of spheres to add
524       
525        Variable row,col,lay,ii,xt,yt,zt,fail=0
526       
527        Make/O/D/N=3 Sobol3D
528       
529        // initialize. XOP should prevent re-initialization
530        SobolX(-1,Sobol3D)
531
532        row=DimSize(mat,0)             
533        col=DimSize(mat,1)
534        lay=DimSize(mat,2)
535        ii=0
536        do
537                SobolX(3,Sobol3D)
538                xt = trunc(row*Sobol3D[0])
539                yt = trunc(col*Sobol3D[1])
540                zt = trunc(lay*Sobol3D[2])
541                if( NoNeighbors(xt,yt,zt,mat) )
542                        mat[xt][yt][zt] = fill
543                        ii+=1           //increment number of spheres actually added
544                        //Print "point ",ii
545                else
546                        fail +=1
547                endif
548        while(ii<num)
549//      Print "failures = ",fail
550       
551//      ParseMatrix3D_rho(mat)          //convert to XYZ
552       
553        return(0)
554End
555
556
557Function RandomFill3DMat(mat,num,fill)
558        Wave mat
559        variable num,fill               //number of spheres to add
560       
561        Variable row,col,lay,ii,xt,yt,zt,fail=0
562       
563        row=DimSize(mat,0)             
564        col=DimSize(mat,1)
565        lay=DimSize(mat,2)     
566        //place first 2 points
567        //xt=trunc(abs(enoise(row)))            //random integer distr betw (0,row-1)
568        //yt=trunc(abs(enoise(col)))
569        //zt=trunc(abs(enoise(lay)))
570        //mat[xt][yt][1] = 1
571        //mat[xt][yt][2] = 1
572        //loop over remaining spheres to add
573        ii=0
574        do
575                xt=trunc(abs(enoise(row)))              //distr betw (0,npt)
576                yt=trunc(abs(enoise(col)))
577                zt=trunc(abs(enoise(lay)))
578                if( NoNeighbors(xt,yt,zt,mat) )
579                        mat[xt][yt][zt] = fill
580                        ii+=1           //increment number of spheres actually added
581                        //Print "point ",ii
582                else
583                        fail +=1
584                endif
585        while(ii<num)
586        Print "failures = ",fail
587       
588//      ParseMatrix3D_rho(mat)          //convert to XYZ
589       
590        return(0)
591End
592
593//point is considered valid if (x,y,z) is unoccupied
594// AND there is a NO neighbor (no connectivity)
595//
596// could speed this up by immediately bouncing out when the first neighbor is found
597Function NoNeighbors(xt,yt,zt,mat)
598        Variable xt,yt,zt
599        WAVE mat
600        //returns 1 if valid, 0 if not valid
601        Variable valid=1,invalid=0,xMax,yMax,zMax,ncont=0
602       
603        NVAR    solventSLD = root:FFT_SolventSLD
604
605        if(mat[xt][yt][zt] != solventSLD)               //already occupied, trial is invalid
606                return (invalid)
607        Endif
608       
609        xMax=DimSize(mat,0) - 1                 //max index of the array
610        yMax=DimSize(mat,1) - 1
611        zMax=DimSize(mat,2) - 1
612        //connected?
613        //any of the surrounding 26 points occupied?
614        Variable xi,yi,zi,ii,jj,kk
615        for(ii=-1;ii<=1;ii+=1)
616                xi = xt+ii
617                xi = (xi< 0) ? 0 : xi           //keep in bounds of matrix index
618                xi = (xi>xMax) ? xMax : xi
619                for(jj=-1;jj<=1;jj+=1)
620                        yi = yt + jj
621                        yi = (yi<0) ? 0 : yi
622                        yi = (yi>yMax) ? yMax : yi
623                        for(kk=-1;kk<=1;kk+=1)
624                                zi = zt+kk
625                                zi = (zi<0) ? 0 : zi
626                                zi = (zi>zMax) ? zMax : zi
627                                if( (xi==xt) && (yi==yt) && (zi==zt) )
628                                        //at central point, do nothing
629                                else
630                                        if(mat[xi][yi][zi] != solventSLD)
631                                                ncont+=1
632                                        endif
633                                endif
634                        endfor
635                endfor
636        endfor
637       
638        //if((ncont > 0) && (ncont <2) )        // only one neighbor
639        if(ncont==0)                    // no nearest neighbor contact points
640                return(valid)
641        else
642                return(invalid)
643        endif
644End
645
646Proc ParseMatrix3DToXYZ(matStr)
647        String matStr="mat"
648        ParseMatrix3D_rho($matStr)
649End
650
651
652///
653/// Depricated - use ParseMatrix3D_rho() instead, to get the SLD information
654///
655///
656// parses the 3d matrix to get xyz coordinates
657// CHANGE - ? Make as byte waves to save space
658//
659// XYZ waves are hard-wired to be "x3d,y3d,z3d"
660// overwrites previous waves
661//
662//Function ParseMatrix3D(mat)
663//      Wave mat
664//     
665//      Variable nptx,npty,nptz,ii,jj,kk,num
666//     
667//      nptx=DimSize(mat,0)
668//      npty=DimSize(mat,1)
669//      nptz=DimSize(mat,2)     
670//      Make/O/D/N=0 x3d,y3d,z3d                //ensure that the waves are SP
671//      num=0
672//     
673//      for(ii=0;ii<nptx;ii+=1)
674//              for(jj=0;jj<npty;jj+=1)
675//                      for(kk=0;kk<nptz;kk+=1)
676//                              if(mat[ii][jj][kk])
677//                                      num+=1
678//                                      InsertPoints num,1,x3d,y3d,z3d
679//                                      x3d[num-1]=ii
680//                                      y3d[num-1]=jj
681//                                      z3d[num-1]=kk
682//                              endif
683//                      endfor
684//              endfor
685//      endfor
686//     
687//      return(0)
688//End
689
690
691// parses the 3d matrix to get xyz coordinates
692// CHANGE - ? Make as byte waves to save space
693//
694// XYZ waves are hard-wired to be "x3d,y3d,z3d"
695// overwrites previous waves
696//
697Function ParseMatrix3D_rho(mat)
698        Wave mat
699
700//      Variable tref = startMSTimer
701               
702        Variable nptx,npty,nptz,ii,jj,kk,num
703       
704        NVAR    solventSLD = root:FFT_SolventSLD
705
706        nptx=DimSize(mat,0)
707        npty=DimSize(mat,1)
708        nptz=DimSize(mat,2)     
709        Make/O/D/N=0 x3d,y3d,z3d,rho3d          //ensure that the waves are DP
710        num=0
711       
712        for(ii=0;ii<nptx;ii+=1)
713                for(jj=0;jj<npty;jj+=1)
714                        for(kk=0;kk<nptz;kk+=1)
715                                if(mat[ii][jj][kk]!=solventSLD)
716                                        num+=1
717                                        InsertPoints num,1,x3d,y3d,z3d,rho3d
718                                        x3d[num-1]=ii
719                                        y3d[num-1]=jj
720                                        z3d[num-1]=kk
721                                        rho3d[num-1]=mat[ii][jj][kk]
722                                endif
723                        endfor
724                endfor
725        endfor
726
727
728//      Variable ms = stopMSTimer(tref)
729//      print "Time elapsed = ", ms/1e6, "s"
730               
731        return(0)
732End
733
734Function ConvertXYZto3N(xw,yw,zw)
735        Wave xw,yw,zw
736       
737        Variable num=numpnts(xw)
738        Make/O/D/N=(num,3) matGiz
739        Concatenate/O {x3d,y3d,z3d},matGiz              //Igor 5+ only
740//
741// or
742//      matGiz[][0] = xw[p]
743//      matGiz[][1] = yw[p]
744//      matGiz[][2] = zw[p]
745       
746        return(0)
747End
748
749//////////////////
750
751//pad the matrix with zeros to make it bigger, ?maybe avoid errors
752// of box size
753// (this does not yet work.....)
754//
755Function PadMatrix(mat,num)
756        Wave mat
757        Variable num
758       
759        Variable numold=DimSize(mat,0)
760        Make/O/B/U/N=(numold+2*num,numold+2*num,numold+2*num) matNew=0
761        matNew[num,numold+num-1][num,numold+num-1][num,numold+num-1]=mat[p-num][q-num][r-num]
762end
763
764// takes values at XYZ and converts to a 3D byte volume
765// - each voxel may have different values
766//
767// xyz values are to start from 0->
768// voxel values must be integer, or they will end up byte
769//
770Function XYZV_toByteVoxels(xx,yy,zz,val)
771        Wave xx,yy,zz,val
772       
773        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt
774       
775        WaveStats/Q xx
776        x1 = V_max
777
778        WaveStats/Q yy
779        y1 = V_max
780       
781        WaveStats/Q zz
782        z1 = V_max
783       
784        //get the maximum dimension
785        npt = max(x1,y1)
786        npt = max(npt,z1)
787       
788        Make/O/B/U/N=(npt,npt,npt) VoxelMat=0
789       
790        num=numpnts(xx)
791       
792        for(ii=0;ii<num;ii+=1)
793                VoxelMat[xx[ii]][yy[ii]][zz[ii]] = val[ii]
794        endfor
795       
796        return(0)
797End
798
799
800// takes XYZV values such as output from  Ken Rubinson's converter
801// where xyz may take negative values, and tries to put it into
802// mat, shifting XYZ as needed
803//
804// now the xyz center is passed in (for a single object) that is treated as the center of that object
805//
806Function XYZV_FillMat_Centered(xx,yy,zz,xc,yc,zc,rad,len,val,erase)
807        Wave xx,yy,zz
808        Variable xc,yc,zc,rad,len,val,erase
809       
810        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt,minp,maxp
811       
812        WAVE mat = root:mat
813        NVAR solventSLD = root:FFT_SolventSLD
814        NVAR FFT_T=root:FFT_T
815
816        if(erase == 1)
817                mat = solventSLD
818        endif
819       
820        num=numpnts(xx)
821       
822        for(ii=0;ii<num;ii+=1)
823                mat[xx[ii]][yy[ii]][zz[ii]] = val
824        endfor
825       
826        return(0)
827End
828
829
830
831
832///// replaced by XYZV_FillMat_Centered() after using Ken's routines, but may have some use in the future
833//
834// takes XYZV values such as output from  Ken Rubinson's converter
835// where xyz may take negative values, and tries to put it into
836// mat, shifting XYZ as needed
837//
838Function XYZV_FillMat(xx,yy,zz,val,erase)
839        Wave xx,yy,zz,val
840        Variable erase
841       
842        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt,minp,maxp
843       
844        WAVE mat = root:mat
845        NVAR solventSLD = root:FFT_SolventSLD
846
847        if(erase)
848                mat = solventSLD
849        endif
850       
851        WaveStats/Q xx
852        x1 = V_min
853
854        WaveStats/Q yy
855        y1 = V_min
856       
857        WaveStats/Q zz
858        z1 = V_min
859       
860       
861        // find the minimum xyz
862        minp = min(x1,y1)
863        minp = min(minp,z1)
864//      print "minp = ",minp
865
866        if(minp < 0)
867                minp = abs(minp)
868        else
869                minp = 0                // if all of the values are positive, don't shift anything
870        endif
871       
872       
873        // will the adjusted xyz values fit the mat?
874        Variable matp = DimSize(mat,0)          // assume all 3 dimensions are the same
875       
876        WaveStats/Q xx
877        x1 = V_max
878
879        WaveStats/Q yy
880        y1 = V_max
881       
882        WaveStats/Q zz
883        z1 = V_max
884       
885        // find the minimum xyz
886        maxp = max(x1,y1)
887        maxp = max(maxp,z1)
888//      Print "maxp = ",maxp
889       
890        if(x1+minp > matp-1)
891                Print "Mat is not big enough. x1+minp = ",x1+minp
892                return(0)
893        Endif
894        if(y1+minp > matp-1)
895                Print "Mat is not big enough. y1+minp = ",y1+minp
896                return(0)
897        Endif
898        if(z1+minp > matp-1)
899                Print "Mat is not big enough. z1+minp = ",z1+minp
900                return(0)
901        Endif
902       
903        num=numpnts(xx)
904       
905        for(ii=0;ii<num;ii+=1)
906                mat[xx[ii]+minp][yy[ii]+minp][zz[ii]+minp] = val[ii]
907        endfor
908       
909        return(0)
910End
Note: See TracBrowser for help on using the repository browser.