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

Last change on this file since 942 was 942, checked in by srkline, 9 years ago

Some cleanup of the FFT routines to be more exact in declaring mat

Rearranged the Event mode panel so that it's a little more obvious what to do, and in what order

Cleaned up the examples for the Simulation Scripting.

File size: 18.8 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//
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// Now enforces periodic boundary conditions
91        WAVE mat
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
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//
153        WAVE mat
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
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//
208        WAVE mat3
210
211        Variable ii,jj,kk,dik,np
212
213        np = DimSize(mat3,0)
217                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )
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
242//      WAVE mat
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 )
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
268//      WAVE mat
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
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!!!
294        Wave mat
296
297        Wave x3d=x3d
298        Wave y3d=y3d
299        Wave z3d=z3d
300
303        NVAR  grid=root:FFT_T
304
305        for(ii=0;ii<num;ii+=1)
306                //print "sphere ",ii
308                if(pd !=0 )
310                endif
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)
321// len is the desired length, in the z-direction
323        WAVE mat
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)
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
339        WAVE mat
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
354        y1 = y1 < 0 ? 0 : y1
355        x1 = x1 < 0 ? 0 : x1
356
357        y2 = y2 > numy ? numy : y2
358        x2 = x2 > numx ? numx : x2
359
360        for(ii=x1;ii<x2;ii+=1)
361                for(jj=y1;jj<y2;jj+=1)
362                                dik=sqrt( (ii-xc)^2+(jj-yc)^2)*grid
364                                        mat[ii][jj][zloc] = fill
365                                endif
366                endfor
367        endfor
368
369//      for(ii=0;ii<numx;ii+=1)
370//              for(jj=0;jj<numy;jj+=1)
371//                              dik=sqrt( (ii-xc)^2+(jj-yc)^2)*grid
373//                                      mat[ii][jj][zloc] = fill
374//                              endif
375//              endfor
376//      endfor
377
378        return(0)
379End
380
381//fills cylinder (fill == 1) or erases it (fill == 0)
382//
383// xc,yc,zc is the center
384// grid is the grid spacing (A)
386// len is the desired length, in the X-direction
388        WAVE mat
390
391//      Variable tref = startMSTimer
392
393        Variable ii,pts
394        //put half the length above - half below
395        pts = trunc(len/2/grid)
396        for(ii=(xc-pts);ii<(xc+pts);ii+=1)
398        endfor
399
400//      Variable ms = stopMSTimer(tref)
401//      print "Time elapsed = ", ms/1e6, "s"
402
403        return(0)
404End
405
406//fills the specified circle (fill==1), or erases it (fill == 0)
407//xloc assumed to be a valid x-index of the matrix
409        WAVE mat
411
412        Variable ii,jj,dik,numx,numy,numz
413
414        numy=DimSize(mat,0)
415        numz=DimSize(mat,1)
416
417        Variable y1,y2,z1,z2
418        y1 = yc - trunc(rad/grid) - 2
419        y2 = yc + trunc(rad/grid) + 2
420        z1 = zc - trunc(rad/grid) - 2
421        z2 = zc + trunc(rad/grid) + 2
422
423        y1 = y1 < 0 ? 0 : y1
424        z1 = z1 < 0 ? 0 : z1
425
426        y2 = y2 > numy ? numy : y2
427        z2 = z2 > numz ? numz : z2
428
429
430        for(ii=y1;ii<y2;ii+=1)
431                for(jj=z1;jj<z2;jj+=1)
432                                dik=sqrt( (ii-yc)^2+(jj-zc)^2)*grid
434                                        mat[xloc][ii][jj] = fill
435                                endif
436                endfor
437        endfor
438
439
440////slow way
441//      for(ii=0;ii<numy;ii+=1)
442//              for(jj=0;jj<numz;jj+=1)
443//                              dik=sqrt( (ii-yc)^2+(jj-zc)^2)*grid
445//                                      mat[xloc][ii][jj] = fill
446//                              endif
447//              endfor
448//      endfor
449
450        return(0)
451End
452
453
454//fills cylinder (fill == 1) or erases it (fill == 0)
455//
456// xc,yc,zc is the center
457// grid is the grid spacing (A)
459// len is the desired length, in the Y-direction
461        WAVE mat
463
464        Variable ii,pts
465        //put half the length above - half below
466        pts = trunc(len/2/grid)
467        for(ii=(yc-pts);ii<(yc+pts);ii+=1)
469        endfor
470
471        return(0)
472End
473
474//fills the specified circle (fill==1), or erases it (fill == 0)
475//yloc assumed to be a valid y-index of the matrix
477        WAVE mat
479
480        Variable ii,jj,dik,numx,numy,numz
481
482        numx=DimSize(mat,0)
483        numz=DimSize(mat,1)
484
485        Variable x1,x2,z1,z2
486        x1 = xc - trunc(rad/grid) - 2
487        x2 = xc + trunc(rad/grid) + 2
488        z1 = zc - trunc(rad/grid) - 2
489        z2 = zc + trunc(rad/grid) + 2
490
491
492        z1 = z1 < 0 ? 0 : z1
493        x1 = x1 < 0 ? 0 : x1
494
495        z2 = z2 > numz ? numz : z2
496        x2 = x2 > numx ? numx : x2
497
498        for(ii=x1;ii<x2;ii+=1)
499                for(jj=z1;jj<z2;jj+=1)
500                                dik=sqrt( (ii-xc)^2+(jj-zc)^2)*grid
502                                        mat[ii][yloc][jj] = fill
503                                endif
504                endfor
505        endfor
506
507//      for(ii=0;ii<numx;ii+=1)
508//              for(jj=0;jj<numz;jj+=1)
509//                              dik=sqrt( (ii-xc)^2+(jj-zc)^2)*grid
511//                                      mat[ii][yloc][jj] = fill
512//                              endif
513//              endfor
514//      endfor
515
516        return(0)
517End
518
519// uses the Sobol' sequence to fill the points
520// - likely more "evenly" spread out than a random pick
521// -- but is it too regular?
522Function SobolFill3DMat(mat,num,fill)
523        Wave mat
524        variable num,fill               //number of spheres to add
525
526        Variable row,col,lay,ii,xt,yt,zt,fail=0
527
528        Make/O/D/N=3 Sobol3D
529
530        // initialize. XOP should prevent re-initialization
531        SobolX(-1,Sobol3D)
532
533        row=DimSize(mat,0)
534        col=DimSize(mat,1)
535        lay=DimSize(mat,2)
536        ii=0
537        do
538                SobolX(3,Sobol3D)
539                xt = trunc(row*Sobol3D[0])
540                yt = trunc(col*Sobol3D[1])
541                zt = trunc(lay*Sobol3D[2])
542                if( NoNeighbors(xt,yt,zt,mat) )
543                        mat[xt][yt][zt] = fill
544                        ii+=1           //increment number of spheres actually added
545                        //Print "point ",ii
546                else
547                        fail +=1
548                endif
549        while(ii<num)
550//      Print "failures = ",fail
551
552//      ParseMatrix3D_rho(mat)          //convert to XYZ
553
554        return(0)
555End
556
557
558Function RandomFill3DMat(mat,num,fill)
559        Wave mat
560        variable num,fill               //number of spheres to add
561
562        Variable row,col,lay,ii,xt,yt,zt,fail=0
563
564        row=DimSize(mat,0)
565        col=DimSize(mat,1)
566        lay=DimSize(mat,2)
567        //place first 2 points
568        //xt=trunc(abs(enoise(row)))            //random integer distr betw (0,row-1)
569        //yt=trunc(abs(enoise(col)))
570        //zt=trunc(abs(enoise(lay)))
571        //mat[xt][yt][1] = 1
572        //mat[xt][yt][2] = 1
573        //loop over remaining spheres to add
574        ii=0
575        do
576                xt=trunc(abs(enoise(row)))              //distr betw (0,npt)
577                yt=trunc(abs(enoise(col)))
578                zt=trunc(abs(enoise(lay)))
579                if( NoNeighbors(xt,yt,zt,mat) )
580                        mat[xt][yt][zt] = fill
581                        ii+=1           //increment number of spheres actually added
582                        //Print "point ",ii
583                else
584                        fail +=1
585                endif
586        while(ii<num)
587        Print "failures = ",fail
588
589//      ParseMatrix3D_rho(mat)          //convert to XYZ
590
591        return(0)
592End
593
594//point is considered valid if (x,y,z) is unoccupied
595// AND there is a NO neighbor (no connectivity)
596//
597// could speed this up by immediately bouncing out when the first neighbor is found
598Function NoNeighbors(xt,yt,zt,mat)
599        Variable xt,yt,zt
600        WAVE mat
601        //returns 1 if valid, 0 if not valid
602        Variable valid=1,invalid=0,xMax,yMax,zMax,ncont=0
603
604        NVAR    solventSLD = root:FFT_SolventSLD
605
606        if(mat[xt][yt][zt] != solventSLD)               //already occupied, trial is invalid
607                return (invalid)
608        Endif
609
610        xMax=DimSize(mat,0) - 1                 //max index of the array
611        yMax=DimSize(mat,1) - 1
612        zMax=DimSize(mat,2) - 1
613        //connected?
614        //any of the surrounding 26 points occupied?
615        Variable xi,yi,zi,ii,jj,kk
616        for(ii=-1;ii<=1;ii+=1)
617                xi = xt+ii
618                xi = (xi< 0) ? 0 : xi           //keep in bounds of matrix index
619                xi = (xi>xMax) ? xMax : xi
620                for(jj=-1;jj<=1;jj+=1)
621                        yi = yt + jj
622                        yi = (yi<0) ? 0 : yi
623                        yi = (yi>yMax) ? yMax : yi
624                        for(kk=-1;kk<=1;kk+=1)
625                                zi = zt+kk
626                                zi = (zi<0) ? 0 : zi
627                                zi = (zi>zMax) ? zMax : zi
628                                if( (xi==xt) && (yi==yt) && (zi==zt) )
629                                        //at central point, do nothing
630                                else
631                                        if(mat[xi][yi][zi] != solventSLD)
632                                                ncont+=1
633                                        endif
634                                endif
635                        endfor
636                endfor
637        endfor
638
639        //if((ncont > 0) && (ncont <2) )        // only one neighbor
640        if(ncont==0)                    // no nearest neighbor contact points
641                return(valid)
642        else
643                return(invalid)
644        endif
645End
646
647Proc ParseMatrix3DToXYZ(matStr)
648        String matStr="mat"
649        ParseMatrix3D_rho(\$matStr)
650End
651
652
653///
654/// Depricated - use ParseMatrix3D_rho() instead, to get the SLD information
655///
656///
657// parses the 3d matrix to get xyz coordinates
658// CHANGE - ? Make as byte waves to save space
659//
660// XYZ waves are hard-wired to be "x3d,y3d,z3d"
661// overwrites previous waves
662//
663//Function ParseMatrix3D(mat)
664//      Wave mat
665//
666//      Variable nptx,npty,nptz,ii,jj,kk,num
667//
668//      nptx=DimSize(mat,0)
669//      npty=DimSize(mat,1)
670//      nptz=DimSize(mat,2)
671//      Make/O/D/N=0 x3d,y3d,z3d                //ensure that the waves are SP
672//      num=0
673//
674//      for(ii=0;ii<nptx;ii+=1)
675//              for(jj=0;jj<npty;jj+=1)
676//                      for(kk=0;kk<nptz;kk+=1)
677//                              if(mat[ii][jj][kk])
678//                                      num+=1
679//                                      InsertPoints num,1,x3d,y3d,z3d
680//                                      x3d[num-1]=ii
681//                                      y3d[num-1]=jj
682//                                      z3d[num-1]=kk
683//                              endif
684//                      endfor
685//              endfor
686//      endfor
687//
688//      return(0)
689//End
690
691
692// parses the 3d matrix to get xyz coordinates
693// CHANGE - ? Make as byte waves to save space
694//
695// XYZ waves are hard-wired to be "x3d,y3d,z3d"
696// overwrites previous waves
697//
698Function ParseMatrix3D_rho(mat)
699        Wave mat
700
701//      Variable tref = startMSTimer
702
703        Variable nptx,npty,nptz,ii,jj,kk,num
704
705        NVAR    solventSLD = root:FFT_SolventSLD
706
707        nptx=DimSize(mat,0)
708        npty=DimSize(mat,1)
709        nptz=DimSize(mat,2)
710        Make/O/D/N=0 x3d,y3d,z3d,rho3d          //ensure that the waves are DP
711        num=0
712
713        for(ii=0;ii<nptx;ii+=1)
714                for(jj=0;jj<npty;jj+=1)
715                        for(kk=0;kk<nptz;kk+=1)
716                                if(mat[ii][jj][kk]!=solventSLD)
717                                        num+=1
718                                        InsertPoints num,1,x3d,y3d,z3d,rho3d
719                                        x3d[num-1]=ii
720                                        y3d[num-1]=jj
721                                        z3d[num-1]=kk
722                                        rho3d[num-1]=mat[ii][jj][kk]
723                                endif
724                        endfor
725                endfor
726        endfor
727
728
729//      Variable ms = stopMSTimer(tref)
730//      print "Time elapsed = ", ms/1e6, "s"
731
732        return(0)
733End
734
735Function ConvertXYZto3N(xw,yw,zw)
736        Wave xw,yw,zw
737
738        Variable num=numpnts(xw)
739        Make/O/D/N=(num,3) matGiz
740        Concatenate/O {x3d,y3d,z3d},matGiz              //Igor 5+ only
741//
742// or
743//      matGiz[][0] = xw[p]
744//      matGiz[][1] = yw[p]
745//      matGiz[][2] = zw[p]
746
747        return(0)
748End
749
750//////////////////
751
752//pad the matrix with zeros to make it bigger, ?maybe avoid errors
753// of box size
754// (this does not yet work.....)
755//
757        Wave mat
758        Variable num
759
760        Variable numold=DimSize(mat,0)
761        Make/O/B/U/N=(numold+2*num,numold+2*num,numold+2*num) matNew=0
762        matNew[num,numold+num-1][num,numold+num-1][num,numold+num-1]=mat[p-num][q-num][r-num]
763end
764
765// takes values at XYZ and converts to a 3D byte volume
766// - each voxel may have different values
767//
768// xyz values are to start from 0->
769// voxel values must be integer, or they will end up byte
770//
771Function XYZV_toByteVoxels(xx,yy,zz,val)
772        Wave xx,yy,zz,val
773
774        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt
775
776        WaveStats/Q xx
777        x1 = V_max
778
779        WaveStats/Q yy
780        y1 = V_max
781
782        WaveStats/Q zz
783        z1 = V_max
784
785        //get the maximum dimension
786        npt = max(x1,y1)
787        npt = max(npt,z1)
788
789        Make/O/B/U/N=(npt,npt,npt) VoxelMat=0
790
791        num=numpnts(xx)
792
793        for(ii=0;ii<num;ii+=1)
794                VoxelMat[xx[ii]][yy[ii]][zz[ii]] = val[ii]
795        endfor
796
797        return(0)
798End
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//
804Function XYZV_FillMat(xx,yy,zz,val,erase)
805        Wave xx,yy,zz,val
806        Variable erase
807
808        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt,minp,maxp
809
810        WAVE mat = root:mat
811        NVAR solventSLD = root:FFT_SolventSLD
812
813        if(erase)
814                mat = solventSLD
815        endif
816
817        WaveStats/Q xx
818        x1 = V_min
819
820        WaveStats/Q yy
821        y1 = V_min
822
823        WaveStats/Q zz
824        z1 = V_min
825
826
827        // find the minimum xyz
828        minp = min(x1,y1)
829        minp = min(minp,z1)
830//      print "minp = ",minp
831
832        if(minp < 0)
833                minp = abs(minp)
834        else
835                minp = 0                // if all of the values are positive, don't shift anything
836        endif
837
838
839        // will the adjusted xyz values fit the mat?
840        Variable matp = DimSize(mat,0)          // assume all 3 dimensions are the same
841
842        WaveStats/Q xx
843        x1 = V_max
844
845        WaveStats/Q yy
846        y1 = V_max
847
848        WaveStats/Q zz
849        z1 = V_max
850
851        // find the minimum xyz
852        maxp = max(x1,y1)
853        maxp = max(maxp,z1)
854//      Print "maxp = ",maxp
855
856        if(x1+minp > matp-1)
857                Print "Mat is not big enough. x1+minp = ",x1+minp
858                return(0)
859        Endif
860        if(y1+minp > matp-1)
861                Print "Mat is not big enough. y1+minp = ",y1+minp
862                return(0)
863        Endif
864        if(z1+minp > matp-1)
865                Print "Mat is not big enough. z1+minp = ",z1+minp
866                return(0)
867        Endif
868
869        num=numpnts(xx)
870
871        for(ii=0;ii<num;ii+=1)
872                mat[xx[ii]+minp][yy[ii]+minp][zz[ii]+minp] = val[ii]
873        endfor
874
875        return(0)
876End
Note: See TracBrowser for help on using the repository browser.