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

Last change on this file since 849 was 849, checked in by srkline, 11 years ago

Changes to the real-space modeling to sped up the drawing of cylinders, and to provide a few more examples of scritped, unattended calculations, and saving the 3D matrix.
File size: 17.2 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// 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//
94Function FillSphereRadiusNoOverlap(mat,grid,rad,xc,yc,zc,fill)
95        WAVE mat
96        Variable grid,rad,xc,yc,zc,fill
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
136                                if(dik<=rad)
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//
149Function FillSphere3(mat3,rad,xc,yc,zc,yesZero)
150        WAVE mat3
151        Variable rad,xc,yc,zc,yesZero
152       
153        Variable ii,jj,kk,dik,np
154       
155        np = DimSize(mat3,0)
156        for(ii=xc-rad;ii<xc+rad;ii+=1)
157                for(jj=yc-rad;jj<yc+rad;jj+=1)
158                        for(kk=zc-rad;kk<zc+rad;kk+=1)
159                                dik=sqrt( (ii-xc)^2+(jj-yc)^2+(kk-zc)^2 )
160                                if(dik<=rad)
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
183//Function EraseSphere(mat,rad,xc,yc,zc)
184//      WAVE mat
185//      Variable rad,xc,yc,zc
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 )
196//                              if(dik<=rad)
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
209//Function EraseSphereRadius(mat,grid,rad,xc,yc,zc)
210//      WAVE mat
211//      Variable grid,rad,xc,yc,zc
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
222//                              if(dik<=rad)
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!!!
235Function  SphereAtEachPoint(mat,rad,pd)
236        Wave mat
237        Variable rad,pd
238       
239        Wave x3d=x3d
240        Wave y3d=y3d
241        Wave z3d=z3d
242       
243        Variable ii=0,num=numpnts(x3d),meanRad
244        meanRad=rad
245        NVAR  grid=root:FFT_T
246       
247        for(ii=0;ii<num;ii+=1)
248                //print "sphere ",ii
249                //FillSphere(mat,rad,x3d[ii],y3d[ii],z3d[ii],0)
250                if(pd !=0 )
251                        rad = meanRad + gnoise(pd*meanRad)
252                endif
253                FillSphereRadius(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],1)
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)
262// rad is the desired radius
263// len is the desired length, in the z-direction
264Function FillZCylinder(mat,grid,rad,xc,yc,zc,len,fill)
265        WAVE mat
266        Variable grid,rad,xc,yc,zc,len,fill
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)
272                FillXYCircle(mat,grid,rad,xc,yc,ii,fill)
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
280Function FillXYCircle(mat,grid,rad,xc,yc,zloc,fill)
281        WAVE mat
282        Variable grid,rad,xc,yc,zloc,fill
283       
284        Variable ii,jj,dik,numx,numy,numz
285       
286        numx=DimSize(mat,0)
287        numy=DimSize(mat,1)
288       
289        Variable x1,x2,y1,y2
290        x1 = xc - trunc(rad/grid) - 2
291        x2 = xc + trunc(rad/grid) + 2
292        y1 = yc - trunc(rad/grid) - 2
293        y2 = yc + trunc(rad/grid) + 2
294
295       
296        y1 = y1 < 0 ? 0 : y1
297        x1 = x1 < 0 ? 0 : x1
298       
299        y2 = y2 > numy ? numy : y2
300        x2 = x2 > numx ? numx : x2
301       
302        for(ii=x1;ii<x2;ii+=1)
303                for(jj=y1;jj<y2;jj+=1)
304                                dik=sqrt( (ii-xc)^2+(jj-yc)^2)*grid
305                                if(dik<=rad)
306                                        mat[ii][jj][zloc] = fill
307                                endif
308                endfor
309        endfor
310       
311//      for(ii=0;ii<numx;ii+=1)
312//              for(jj=0;jj<numy;jj+=1)
313//                              dik=sqrt( (ii-xc)^2+(jj-yc)^2)*grid
314//                              if(dik<=rad)
315//                                      mat[ii][jj][zloc] = fill
316//                              endif
317//              endfor
318//      endfor
319       
320        return(0)       
321End
322
323//fills cylinder (fill == 1) or erases it (fill == 0)
324//
325// xc,yc,zc is the center
326// grid is the grid spacing (A)
327// rad is the desired radius
328// len is the desired length, in the X-direction
329Function FillXCylinder(mat,grid,rad,xc,yc,zc,len,fill)
330        WAVE mat
331        Variable grid,rad,xc,yc,zc,len,fill
332
333//      Variable tref = startMSTimer
334               
335        Variable ii,pts
336        //put half the length above - half below
337        pts = trunc(len/2/grid)
338        for(ii=(xc-pts);ii<(xc+pts);ii+=1)
339                FillYZCircle(mat,grid,rad,yc,zc,ii,fill)
340        endfor
341
342//      Variable ms = stopMSTimer(tref)
343//      print "Time elapsed = ", ms/1e6, "s"
344               
345        return(0)       
346End
347
348//fills the specified circle (fill==1), or erases it (fill == 0)
349//xloc assumed to be a valid x-index of the matrix
350Function FillYZCircle(mat,grid,rad,yc,zc,xloc,fill)
351        WAVE mat
352        Variable grid,rad,yc,zc,xloc,fill
353       
354        Variable ii,jj,dik,numx,numy,numz
355       
356        numy=DimSize(mat,0)
357        numz=DimSize(mat,1)
358       
359        Variable y1,y2,z1,z2
360        y1 = yc - trunc(rad/grid) - 2
361        y2 = yc + trunc(rad/grid) + 2
362        z1 = zc - trunc(rad/grid) - 2
363        z2 = zc + trunc(rad/grid) + 2
364       
365        y1 = y1 < 0 ? 0 : y1
366        z1 = z1 < 0 ? 0 : z1
367       
368        y2 = y2 > numy ? numy : y2
369        z2 = z2 > numz ? numz : z2
370       
371
372        for(ii=y1;ii<y2;ii+=1)
373                for(jj=z1;jj<z2;jj+=1)
374                                dik=sqrt( (ii-yc)^2+(jj-zc)^2)*grid
375                                if(dik<=rad)
376                                        mat[xloc][ii][jj] = fill
377                                endif
378                endfor
379        endfor
380
381
382////slow way           
383//      for(ii=0;ii<numy;ii+=1)
384//              for(jj=0;jj<numz;jj+=1)
385//                              dik=sqrt( (ii-yc)^2+(jj-zc)^2)*grid
386//                              if(dik<=rad)
387//                                      mat[xloc][ii][jj] = fill
388//                              endif
389//              endfor
390//      endfor
391
392        return(0)       
393End
394
395
396//fills cylinder (fill == 1) or erases it (fill == 0)
397//
398// xc,yc,zc is the center
399// grid is the grid spacing (A)
400// rad is the desired radius
401// len is the desired length, in the Y-direction
402Function FillYCylinder(mat,grid,rad,xc,yc,zc,len,fill)
403        WAVE mat
404        Variable grid,rad,xc,yc,zc,len,fill
405       
406        Variable ii,pts
407        //put half the length above - half below
408        pts = trunc(len/2/grid)
409        for(ii=(yc-pts);ii<(yc+pts);ii+=1)
410                FillXZCircle(mat,grid,rad,xc,zc,ii,fill)
411        endfor
412       
413        return(0)       
414End
415
416//fills the specified circle (fill==1), or erases it (fill == 0)
417//yloc assumed to be a valid y-index of the matrix
418Function FillXZCircle(mat,grid,rad,xc,zc,yloc,fill)
419        WAVE mat
420        Variable grid,rad,xc,zc,yloc,fill
421       
422        Variable ii,jj,dik,numx,numy,numz
423       
424        numx=DimSize(mat,0)
425        numz=DimSize(mat,1)
426       
427        Variable x1,x2,z1,z2
428        x1 = xc - trunc(rad/grid) - 2
429        x2 = xc + trunc(rad/grid) + 2
430        z1 = zc - trunc(rad/grid) - 2
431        z2 = zc + trunc(rad/grid) + 2
432
433       
434        z1 = z1 < 0 ? 0 : z1
435        x1 = x1 < 0 ? 0 : x1
436       
437        z2 = z2 > numz ? numz : z2
438        x2 = x2 > numx ? numx : x2
439       
440        for(ii=x1;ii<x2;ii+=1)
441                for(jj=z1;jj<z2;jj+=1)
442                                dik=sqrt( (ii-xc)^2+(jj-zc)^2)*grid
443                                if(dik<=rad)
444                                        mat[ii][yloc][jj] = fill
445                                endif
446                endfor
447        endfor
448       
449//      for(ii=0;ii<numx;ii+=1)
450//              for(jj=0;jj<numz;jj+=1)
451//                              dik=sqrt( (ii-xc)^2+(jj-zc)^2)*grid
452//                              if(dik<=rad)
453//                                      mat[ii][yloc][jj] = fill
454//                              endif
455//              endfor
456//      endfor
457       
458        return(0)       
459End
460
461// uses the Sobol' sequence to fill the points
462// - likely more "evenly" spread out than a random pick
463// -- but is it too regular?
464Function SobolFill3DMat(mat,num,fill)
465        Wave mat
466        variable num,fill               //number of spheres to add
467       
468        Variable row,col,lay,ii,xt,yt,zt,fail=0
469       
470        Make/O/D/N=3 Sobol3D
471       
472        // initialize. XOP should prevent re-initialization
473        SobolX(-1,Sobol3D)
474
475        row=DimSize(mat,0)             
476        col=DimSize(mat,1)
477        lay=DimSize(mat,2)
478        ii=0
479        do
480                SobolX(3,Sobol3D)
481                xt = trunc(row*Sobol3D[0])
482                yt = trunc(col*Sobol3D[1])
483                zt = trunc(lay*Sobol3D[2])
484                if( NoNeighbors(xt,yt,zt,mat) )
485                        mat[xt][yt][zt] = fill
486                        ii+=1           //increment number of spheres actually added
487                        //Print "point ",ii
488                else
489                        fail +=1
490                endif
491        while(ii<num)
492//      Print "failures = ",fail
493       
494//      ParseMatrix3D_rho(mat)          //convert to XYZ
495       
496        return(0)
497End
498
499
500Function RandomFill3DMat(mat,num,fill)
501        Wave mat
502        variable num,fill               //number of spheres to add
503       
504        Variable row,col,lay,ii,xt,yt,zt,fail=0
505       
506        row=DimSize(mat,0)             
507        col=DimSize(mat,1)
508        lay=DimSize(mat,2)     
509        //place first 2 points
510        //xt=trunc(abs(enoise(row)))            //random integer distr betw (0,row-1)
511        //yt=trunc(abs(enoise(col)))
512        //zt=trunc(abs(enoise(lay)))
513        //mat[xt][yt][1] = 1
514        //mat[xt][yt][2] = 1
515        //loop over remaining spheres to add
516        ii=0
517        do
518                xt=trunc(abs(enoise(row)))              //distr betw (0,npt)
519                yt=trunc(abs(enoise(col)))
520                zt=trunc(abs(enoise(lay)))
521                if( NoNeighbors(xt,yt,zt,mat) )
522                        mat[xt][yt][zt] = fill
523                        ii+=1           //increment number of spheres actually added
524                        //Print "point ",ii
525                else
526                        fail +=1
527                endif
528        while(ii<num)
529        Print "failures = ",fail
530       
531//      ParseMatrix3D_rho(mat)          //convert to XYZ
532       
533        return(0)
534End
535
536//point is considered valid if (x,y,z) is unoccupied
537// AND there is a NO neighbor (no connectivity)
538//
539// could speed this up by immediately bouncing out when the first neighbor is found
540Function NoNeighbors(xt,yt,zt,mat)
541        Variable xt,yt,zt
542        WAVE mat
543        //returns 1 if valid, 0 if not valid
544        Variable valid=1,invalid=0,xMax,yMax,zMax,ncont=0
545       
546        NVAR    solventSLD = root:FFT_SolventSLD
547
548        if(mat[xt][yt][zt] != solventSLD)               //already occupied, trial is invalid
549                return (invalid)
550        Endif
551       
552        xMax=DimSize(mat,0) - 1                 //max index of the array
553        yMax=DimSize(mat,1) - 1
554        zMax=DimSize(mat,2) - 1
555        //connected?
556        //any of the surrounding 26 points occupied?
557        Variable xi,yi,zi,ii,jj,kk
558        for(ii=-1;ii<=1;ii+=1)
559                xi = xt+ii
560                xi = (xi< 0) ? 0 : xi           //keep in bounds of matrix index
561                xi = (xi>xMax) ? xMax : xi
562                for(jj=-1;jj<=1;jj+=1)
563                        yi = yt + jj
564                        yi = (yi<0) ? 0 : yi
565                        yi = (yi>yMax) ? yMax : yi
566                        for(kk=-1;kk<=1;kk+=1)
567                                zi = zt+kk
568                                zi = (zi<0) ? 0 : zi
569                                zi = (zi>zMax) ? zMax : zi
570                                if( (xi==xt) && (yi==yt) && (zi==zt) )
571                                        //at central point, do nothing
572                                else
573                                        if(mat[xi][yi][zi] != solventSLD)
574                                                ncont+=1
575                                        endif
576                                endif
577                        endfor
578                endfor
579        endfor
580       
581        //if((ncont > 0) && (ncont <2) )        // only one neighbor
582        if(ncont==0)                    // no nearest neighbor contact points
583                return(valid)
584        else
585                return(invalid)
586        endif
587End
588
589Proc ParseMatrix3DToXYZ(matStr)
590        String matStr="mat"
591        ParseMatrix3D_rho($matStr)
592End
593
594
595///
596/// Depricated - use ParseMatrix3D_rho() instead, to get the SLD information
597///
598///
599// parses the 3d matrix to get xyz coordinates
600// CHANGE - ? Make as byte waves to save space
601//
602// XYZ waves are hard-wired to be "x3d,y3d,z3d"
603// overwrites previous waves
604//
605//Function ParseMatrix3D(mat)
606//      Wave mat
607//     
608//      Variable nptx,npty,nptz,ii,jj,kk,num
609//     
610//      nptx=DimSize(mat,0)
611//      npty=DimSize(mat,1)
612//      nptz=DimSize(mat,2)     
613//      Make/O/D/N=0 x3d,y3d,z3d                //ensure that the waves are SP
614//      num=0
615//     
616//      for(ii=0;ii<nptx;ii+=1)
617//              for(jj=0;jj<npty;jj+=1)
618//                      for(kk=0;kk<nptz;kk+=1)
619//                              if(mat[ii][jj][kk])
620//                                      num+=1
621//                                      InsertPoints num,1,x3d,y3d,z3d
622//                                      x3d[num-1]=ii
623//                                      y3d[num-1]=jj
624//                                      z3d[num-1]=kk
625//                              endif
626//                      endfor
627//              endfor
628//      endfor
629//     
630//      return(0)
631//End
632
633
634// parses the 3d matrix to get xyz coordinates
635// CHANGE - ? Make as byte waves to save space
636//
637// XYZ waves are hard-wired to be "x3d,y3d,z3d"
638// overwrites previous waves
639//
640Function ParseMatrix3D_rho(mat)
641        Wave mat
642
643//      Variable tref = startMSTimer
644               
645        Variable nptx,npty,nptz,ii,jj,kk,num
646       
647        NVAR    solventSLD = root:FFT_SolventSLD
648
649        nptx=DimSize(mat,0)
650        npty=DimSize(mat,1)
651        nptz=DimSize(mat,2)     
652        Make/O/D/N=0 x3d,y3d,z3d,rho3d          //ensure that the waves are DP
653        num=0
654       
655        for(ii=0;ii<nptx;ii+=1)
656                for(jj=0;jj<npty;jj+=1)
657                        for(kk=0;kk<nptz;kk+=1)
658                                if(mat[ii][jj][kk]!=solventSLD)
659                                        num+=1
660                                        InsertPoints num,1,x3d,y3d,z3d,rho3d
661                                        x3d[num-1]=ii
662                                        y3d[num-1]=jj
663                                        z3d[num-1]=kk
664                                        rho3d[num-1]=mat[ii][jj][kk]
665                                endif
666                        endfor
667                endfor
668        endfor
669
670
671//      Variable ms = stopMSTimer(tref)
672//      print "Time elapsed = ", ms/1e6, "s"
673               
674        return(0)
675End
676
677Function ConvertXYZto3N(xw,yw,zw)
678        Wave xw,yw,zw
679       
680        Variable num=numpnts(xw)
681        Make/O/D/N=(num,3) matGiz
682        Concatenate/O {x3d,y3d,z3d},matGiz              //Igor 5+ only
683//
684// or
685//      matGiz[][0] = xw[p]
686//      matGiz[][1] = yw[p]
687//      matGiz[][2] = zw[p]
688       
689        return(0)
690End
691
692//////////////////
693
694//pad the matrix with zeros to make it bigger, ?maybe avoid errors
695// of box size
696// (this does not yet work.....)
697//
698Function PadMatrix(mat,num)
699        Wave mat
700        Variable num
701       
702        Variable numold=DimSize(mat,0)
703        Make/O/B/U/N=(numold+2*num,numold+2*num,numold+2*num) matNew=0
704        matNew[num,numold+num-1][num,numold+num-1][num,numold+num-1]=mat[p-num][q-num][r-num]
705end
706
707// takes values at XYZ and converts to a 3D byte volume
708// - each voxel may have different values
709//
710// xyz values are to start from 0->
711// voxel values must be integer, or they will end up byte
712//
713Function XYZV_toByteVoxels(xx,yy,zz,val)
714        Wave xx,yy,zz,val
715       
716        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt
717       
718        WaveStats/Q xx
719        x1 = V_max
720
721        WaveStats/Q yy
722        y1 = V_max
723       
724        WaveStats/Q zz
725        z1 = V_max
726       
727        //get the maximum dimension
728        npt = max(x1,y1)
729        npt = max(npt,z1)
730       
731        Make/O/B/U/N=(npt,npt,npt) VoxelMat=0
732       
733        num=numpnts(xx)
734       
735        for(ii=0;ii<num;ii+=1)
736                VoxelMat[xx[ii]][yy[ii]][zz[ii]] = val[ii]
737        endfor
738       
739        return(0)
740End
741
742// takes XYZV values such as output from  Ken Rubinson's converter
743// where xyz may take negative values, and tries to put it into
744// mat, shifting XYZ as needed
745//
746Function XYZV_FillMat(xx,yy,zz,val,erase)
747        Wave xx,yy,zz,val
748        Variable erase
749       
750        Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt,minp,maxp
751       
752        WAVE mat = root:mat
753        NVAR solventSLD = root:FFT_SolventSLD
754
755        if(erase)
756                mat = solventSLD
757        endif
758       
759        WaveStats/Q xx
760        x1 = V_min
761
762        WaveStats/Q yy
763        y1 = V_min
764       
765        WaveStats/Q zz
766        z1 = V_min
767       
768       
769        // find the minimum xyz
770        minp = min(x1,y1)
771        minp = min(minp,z1)
772//      print "minp = ",minp
773
774        if(minp < 0)
775                minp = abs(minp)
776        else
777                minp = 0                // if all of the values are positive, don't shift anything
778        endif
779       
780       
781        // will the adjusted xyz values fit the mat?
782        Variable matp = DimSize(mat,0)          // assume all 3 dimensions are the same
783       
784        WaveStats/Q xx
785        x1 = V_max
786
787        WaveStats/Q yy
788        y1 = V_max
789       
790        WaveStats/Q zz
791        z1 = V_max
792       
793        // find the minimum xyz
794        maxp = max(x1,y1)
795        maxp = max(maxp,z1)
796//      Print "maxp = ",maxp
797       
798        if(x1+minp > matp-1)
799                Print "Mat is not big enough. x1+minp = ",x1+minp
800                return(0)
801        Endif
802        if(y1+minp > matp-1)
803                Print "Mat is not big enough. y1+minp = ",y1+minp
804                return(0)
805        Endif
806        if(z1+minp > matp-1)
807                Print "Mat is not big enough. z1+minp = ",z1+minp
808                return(0)
809        Endif
810       
811        num=numpnts(xx)
812       
813        for(ii=0;ii<num;ii+=1)
814                mat[xx[ii]+minp][yy[ii]+minp][zz[ii]+minp] = val[ii]
815        endfor
816       
817        return(0)
818End
Note: See TracBrowser for help on using the repository browser.