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 | // |
---|
13 | Function 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) |
---|
38 | End |
---|
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 |
---|
44 | Function 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) |
---|
84 | End |
---|
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 | // |
---|
94 | Function 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) |
---|
143 | End |
---|
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 | // |
---|
149 | Function 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) |
---|
178 | End |
---|
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!!! |
---|
235 | Function 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 | |
---|
256 | End |
---|
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 |
---|
264 | Function 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) |
---|
276 | End |
---|
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 |
---|
280 | Function 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 | 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 |
---|
291 | if(dik<=rad) |
---|
292 | mat[ii][jj][zloc] = fill |
---|
293 | endif |
---|
294 | endfor |
---|
295 | endfor |
---|
296 | |
---|
297 | return(0) |
---|
298 | End |
---|
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) |
---|
304 | // rad is the desired radius |
---|
305 | // len is the desired length, in the X-direction |
---|
306 | Function FillXCylinder(mat,grid,rad,xc,yc,zc,len,fill) |
---|
307 | WAVE mat |
---|
308 | Variable grid,rad,xc,yc,zc,len,fill |
---|
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) |
---|
314 | FillYZCircle(mat,grid,rad,yc,zc,ii,fill) |
---|
315 | endfor |
---|
316 | |
---|
317 | return(0) |
---|
318 | End |
---|
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 |
---|
322 | Function FillYZCircle(mat,grid,rad,yc,zc,xloc,fill) |
---|
323 | WAVE mat |
---|
324 | Variable grid,rad,yc,zc,xloc,fill |
---|
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 |
---|
333 | if(dik<=rad) |
---|
334 | mat[xloc][ii][jj] = fill |
---|
335 | endif |
---|
336 | endfor |
---|
337 | endfor |
---|
338 | |
---|
339 | return(0) |
---|
340 | End |
---|
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) |
---|
347 | // rad is the desired radius |
---|
348 | // len is the desired length, in the Y-direction |
---|
349 | Function FillYCylinder(mat,grid,rad,xc,yc,zc,len,fill) |
---|
350 | WAVE mat |
---|
351 | Variable grid,rad,xc,yc,zc,len,fill |
---|
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) |
---|
357 | FillXZCircle(mat,grid,rad,xc,zc,ii,fill) |
---|
358 | endfor |
---|
359 | |
---|
360 | return(0) |
---|
361 | End |
---|
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 |
---|
365 | Function FillXZCircle(mat,grid,rad,xc,zc,yloc,fill) |
---|
366 | WAVE mat |
---|
367 | Variable grid,rad,xc,zc,yloc,fill |
---|
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 |
---|
376 | if(dik<=rad) |
---|
377 | mat[ii][yloc][jj] = fill |
---|
378 | endif |
---|
379 | endfor |
---|
380 | endfor |
---|
381 | |
---|
382 | return(0) |
---|
383 | End |
---|
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? |
---|
388 | Function SobolFill3DMat(mat,num,fill) |
---|
389 | Wave mat |
---|
390 | variable num,fill //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] = fill |
---|
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) |
---|
421 | End |
---|
422 | |
---|
423 | |
---|
424 | Function RandomFill3DMat(mat,num,fill) |
---|
425 | Wave mat |
---|
426 | variable num,fill //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] = fill |
---|
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) |
---|
458 | End |
---|
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 |
---|
464 | Function 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 | NVAR solventSLD = root:FFT_SolventSLD |
---|
471 | |
---|
472 | if(mat[xt][yt][zt] != solventSLD) //already occupied, trial is invalid |
---|
473 | return (invalid) |
---|
474 | Endif |
---|
475 | |
---|
476 | xMax=DimSize(mat,0) - 1 //max index of the array |
---|
477 | yMax=DimSize(mat,1) - 1 |
---|
478 | zMax=DimSize(mat,2) - 1 |
---|
479 | //connected? |
---|
480 | //any of the surrounding 26 points occupied? |
---|
481 | Variable xi,yi,zi,ii,jj,kk |
---|
482 | for(ii=-1;ii<=1;ii+=1) |
---|
483 | xi = xt+ii |
---|
484 | xi = (xi< 0) ? 0 : xi //keep in bounds of matrix index |
---|
485 | xi = (xi>xMax) ? xMax : xi |
---|
486 | for(jj=-1;jj<=1;jj+=1) |
---|
487 | yi = yt + jj |
---|
488 | yi = (yi<0) ? 0 : yi |
---|
489 | yi = (yi>yMax) ? yMax : yi |
---|
490 | for(kk=-1;kk<=1;kk+=1) |
---|
491 | zi = zt+kk |
---|
492 | zi = (zi<0) ? 0 : zi |
---|
493 | zi = (zi>zMax) ? zMax : zi |
---|
494 | if( (xi==xt) && (yi==yt) && (zi==zt) ) |
---|
495 | //at central point, do nothing |
---|
496 | else |
---|
497 | if(mat[xi][yi][zi] != solventSLD) |
---|
498 | ncont+=1 |
---|
499 | endif |
---|
500 | endif |
---|
501 | endfor |
---|
502 | endfor |
---|
503 | endfor |
---|
504 | |
---|
505 | //if((ncont > 0) && (ncont <2) ) // only one neighbor |
---|
506 | if(ncont==0) // no nearest neighbor contact points |
---|
507 | return(valid) |
---|
508 | else |
---|
509 | return(invalid) |
---|
510 | endif |
---|
511 | End |
---|
512 | |
---|
513 | Proc ParseMatrix3DToXYZ(matStr) |
---|
514 | String matStr="mat" |
---|
515 | ParseMatrix3D_rho($matStr) |
---|
516 | End |
---|
517 | |
---|
518 | |
---|
519 | /// |
---|
520 | /// Depricated - use ParseMatrix3D_rho() instead, to get the SLD information |
---|
521 | /// |
---|
522 | /// |
---|
523 | // parses the 3d matrix to get xyz coordinates |
---|
524 | // CHANGE - ? Make as byte waves to save space |
---|
525 | // |
---|
526 | // XYZ waves are hard-wired to be "x3d,y3d,z3d" |
---|
527 | // overwrites previous waves |
---|
528 | // |
---|
529 | //Function ParseMatrix3D(mat) |
---|
530 | // Wave mat |
---|
531 | // |
---|
532 | // Variable nptx,npty,nptz,ii,jj,kk,num |
---|
533 | // |
---|
534 | // nptx=DimSize(mat,0) |
---|
535 | // npty=DimSize(mat,1) |
---|
536 | // nptz=DimSize(mat,2) |
---|
537 | // Make/O/D/N=0 x3d,y3d,z3d //ensure that the waves are SP |
---|
538 | // num=0 |
---|
539 | // |
---|
540 | // for(ii=0;ii<nptx;ii+=1) |
---|
541 | // for(jj=0;jj<npty;jj+=1) |
---|
542 | // for(kk=0;kk<nptz;kk+=1) |
---|
543 | // if(mat[ii][jj][kk]) |
---|
544 | // num+=1 |
---|
545 | // InsertPoints num,1,x3d,y3d,z3d |
---|
546 | // x3d[num-1]=ii |
---|
547 | // y3d[num-1]=jj |
---|
548 | // z3d[num-1]=kk |
---|
549 | // endif |
---|
550 | // endfor |
---|
551 | // endfor |
---|
552 | // endfor |
---|
553 | // |
---|
554 | // return(0) |
---|
555 | //End |
---|
556 | |
---|
557 | |
---|
558 | // parses the 3d matrix to get xyz coordinates |
---|
559 | // CHANGE - ? Make as byte waves to save space |
---|
560 | // |
---|
561 | // XYZ waves are hard-wired to be "x3d,y3d,z3d" |
---|
562 | // overwrites previous waves |
---|
563 | // |
---|
564 | Function ParseMatrix3D_rho(mat) |
---|
565 | Wave mat |
---|
566 | |
---|
567 | Variable nptx,npty,nptz,ii,jj,kk,num |
---|
568 | |
---|
569 | NVAR solventSLD = root:FFT_SolventSLD |
---|
570 | |
---|
571 | nptx=DimSize(mat,0) |
---|
572 | npty=DimSize(mat,1) |
---|
573 | nptz=DimSize(mat,2) |
---|
574 | Make/O/D/N=0 x3d,y3d,z3d,rho3d //ensure that the waves are DP |
---|
575 | num=0 |
---|
576 | |
---|
577 | for(ii=0;ii<nptx;ii+=1) |
---|
578 | for(jj=0;jj<npty;jj+=1) |
---|
579 | for(kk=0;kk<nptz;kk+=1) |
---|
580 | if(mat[ii][jj][kk]!=solventSLD) |
---|
581 | num+=1 |
---|
582 | InsertPoints num,1,x3d,y3d,z3d,rho3d |
---|
583 | x3d[num-1]=ii |
---|
584 | y3d[num-1]=jj |
---|
585 | z3d[num-1]=kk |
---|
586 | rho3d[num-1]=mat[ii][jj][kk] |
---|
587 | endif |
---|
588 | endfor |
---|
589 | endfor |
---|
590 | endfor |
---|
591 | |
---|
592 | return(0) |
---|
593 | End |
---|
594 | |
---|
595 | Function ConvertXYZto3N(xw,yw,zw) |
---|
596 | Wave xw,yw,zw |
---|
597 | |
---|
598 | Variable num=numpnts(xw) |
---|
599 | Make/O/D/N=(num,3) matGiz |
---|
600 | Concatenate/O {x3d,y3d,z3d},matGiz //Igor 5+ only |
---|
601 | // |
---|
602 | // or |
---|
603 | // matGiz[][0] = xw[p] |
---|
604 | // matGiz[][1] = yw[p] |
---|
605 | // matGiz[][2] = zw[p] |
---|
606 | |
---|
607 | return(0) |
---|
608 | End |
---|
609 | |
---|
610 | ////////////////// |
---|
611 | |
---|
612 | //pad the matrix with zeros to make it bigger, ?maybe avoid errors |
---|
613 | // of box size |
---|
614 | // (this does not yet work.....) |
---|
615 | // |
---|
616 | Function PadMatrix(mat,num) |
---|
617 | Wave mat |
---|
618 | Variable num |
---|
619 | |
---|
620 | Variable numold=DimSize(mat,0) |
---|
621 | Make/O/B/U/N=(numold+2*num,numold+2*num,numold+2*num) matNew=0 |
---|
622 | matNew[num,numold+num-1][num,numold+num-1][num,numold+num-1]=mat[p-num][q-num][r-num] |
---|
623 | end |
---|
624 | |
---|
625 | // takes values at XYZ and converts to a 3D byte volume |
---|
626 | // - each voxel may have different values |
---|
627 | // |
---|
628 | // xyz values are to start from 0-> |
---|
629 | // voxel values must be integer, or they will end up byte |
---|
630 | // |
---|
631 | Function XYZV_toByteVoxels(xx,yy,zz,val) |
---|
632 | Wave xx,yy,zz,val |
---|
633 | |
---|
634 | Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt |
---|
635 | |
---|
636 | WaveStats/Q xx |
---|
637 | x1 = V_max |
---|
638 | |
---|
639 | WaveStats/Q yy |
---|
640 | y1 = V_max |
---|
641 | |
---|
642 | WaveStats/Q zz |
---|
643 | z1 = V_max |
---|
644 | |
---|
645 | //get the maximum dimension |
---|
646 | npt = max(x1,y1) |
---|
647 | npt = max(npt,z1) |
---|
648 | |
---|
649 | Make/O/B/U/N=(npt,npt,npt) VoxelMat=0 |
---|
650 | |
---|
651 | num=numpnts(xx) |
---|
652 | |
---|
653 | for(ii=0;ii<num;ii+=1) |
---|
654 | VoxelMat[xx[ii]][yy[ii]][zz[ii]] = val[ii] |
---|
655 | endfor |
---|
656 | |
---|
657 | return(0) |
---|
658 | End |
---|
659 | |
---|
660 | // takes XYZV values such as output from Ken Rubinson's converter |
---|
661 | // where xyz may take negative values, and tries to put it into |
---|
662 | // mat, shifting XYZ as needed |
---|
663 | // |
---|
664 | Function XYZV_FillMat(xx,yy,zz,val,erase) |
---|
665 | Wave xx,yy,zz,val |
---|
666 | Variable erase |
---|
667 | |
---|
668 | Variable x1,x2,y1,y2,z1,z2,ii,jj,kk,num,npt,minp,maxp |
---|
669 | |
---|
670 | WAVE mat = root:mat |
---|
671 | NVAR solventSLD = root:FFT_SolventSLD |
---|
672 | |
---|
673 | if(erase) |
---|
674 | mat = solventSLD |
---|
675 | endif |
---|
676 | |
---|
677 | WaveStats/Q xx |
---|
678 | x1 = V_min |
---|
679 | |
---|
680 | WaveStats/Q yy |
---|
681 | y1 = V_min |
---|
682 | |
---|
683 | WaveStats/Q zz |
---|
684 | z1 = V_min |
---|
685 | |
---|
686 | |
---|
687 | // find the minimum xyz |
---|
688 | minp = min(x1,y1) |
---|
689 | minp = min(minp,z1) |
---|
690 | // print "minp = ",minp |
---|
691 | |
---|
692 | if(minp < 0) |
---|
693 | minp = abs(minp) |
---|
694 | else |
---|
695 | minp = 0 // if all of the values are positive, don't shift anything |
---|
696 | endif |
---|
697 | |
---|
698 | |
---|
699 | // will the adjusted xyz values fit the mat? |
---|
700 | Variable matp = DimSize(mat,0) // assume all 3 dimensions are the same |
---|
701 | |
---|
702 | WaveStats/Q xx |
---|
703 | x1 = V_max |
---|
704 | |
---|
705 | WaveStats/Q yy |
---|
706 | y1 = V_max |
---|
707 | |
---|
708 | WaveStats/Q zz |
---|
709 | z1 = V_max |
---|
710 | |
---|
711 | // find the minimum xyz |
---|
712 | maxp = max(x1,y1) |
---|
713 | maxp = max(maxp,z1) |
---|
714 | // Print "maxp = ",maxp |
---|
715 | |
---|
716 | if(x1+minp > matp-1) |
---|
717 | Print "Mat is not big enough. x1+minp = ",x1+minp |
---|
718 | return(0) |
---|
719 | Endif |
---|
720 | if(y1+minp > matp-1) |
---|
721 | Print "Mat is not big enough. y1+minp = ",y1+minp |
---|
722 | return(0) |
---|
723 | Endif |
---|
724 | if(z1+minp > matp-1) |
---|
725 | Print "Mat is not big enough. z1+minp = ",z1+minp |
---|
726 | return(0) |
---|
727 | Endif |
---|
728 | |
---|
729 | num=numpnts(xx) |
---|
730 | |
---|
731 | for(ii=0;ii<num;ii+=1) |
---|
732 | mat[xx[ii]+minp][yy[ii]+minp][zz[ii]+minp] = val[ii] |
---|
733 | endfor |
---|
734 | |
---|
735 | return(0) |
---|
736 | End |
---|