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 | // Now enforces periodic boundary conditions |
---|
90 | Function 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) |
---|
142 | End |
---|
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 | // |
---|
152 | Function 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) |
---|
201 | End |
---|
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 | // |
---|
207 | Function 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) |
---|
236 | End |
---|
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!!! |
---|
293 | Function 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 | |
---|
314 | End |
---|
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 |
---|
322 | Function 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) |
---|
334 | End |
---|
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 |
---|
338 | Function 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 | |
---|
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 |
---|
363 | if(dik<=rad) |
---|
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 |
---|
372 | // if(dik<=rad) |
---|
373 | // mat[ii][jj][zloc] = fill |
---|
374 | // endif |
---|
375 | // endfor |
---|
376 | // endfor |
---|
377 | |
---|
378 | return(0) |
---|
379 | End |
---|
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) |
---|
385 | // rad is the desired radius |
---|
386 | // len is the desired length, in the X-direction |
---|
387 | Function FillXCylinder(mat,grid,rad,xc,yc,zc,len,fill) |
---|
388 | WAVE mat |
---|
389 | Variable grid,rad,xc,yc,zc,len,fill |
---|
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) |
---|
397 | FillYZCircle(mat,grid,rad,yc,zc,ii,fill) |
---|
398 | endfor |
---|
399 | |
---|
400 | // Variable ms = stopMSTimer(tref) |
---|
401 | // print "Time elapsed = ", ms/1e6, "s" |
---|
402 | |
---|
403 | return(0) |
---|
404 | End |
---|
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 |
---|
408 | Function FillYZCircle(mat,grid,rad,yc,zc,xloc,fill) |
---|
409 | WAVE mat |
---|
410 | Variable grid,rad,yc,zc,xloc,fill |
---|
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 |
---|
433 | if(dik<=rad) |
---|
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 |
---|
444 | // if(dik<=rad) |
---|
445 | // mat[xloc][ii][jj] = fill |
---|
446 | // endif |
---|
447 | // endfor |
---|
448 | // endfor |
---|
449 | |
---|
450 | return(0) |
---|
451 | End |
---|
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) |
---|
458 | // rad is the desired radius |
---|
459 | // len is the desired length, in the Y-direction |
---|
460 | Function FillYCylinder(mat,grid,rad,xc,yc,zc,len,fill) |
---|
461 | WAVE mat |
---|
462 | Variable grid,rad,xc,yc,zc,len,fill |
---|
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) |
---|
468 | FillXZCircle(mat,grid,rad,xc,zc,ii,fill) |
---|
469 | endfor |
---|
470 | |
---|
471 | return(0) |
---|
472 | End |
---|
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 |
---|
476 | Function FillXZCircle(mat,grid,rad,xc,zc,yloc,fill) |
---|
477 | WAVE mat |
---|
478 | Variable grid,rad,xc,zc,yloc,fill |
---|
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 |
---|
501 | if(dik<=rad) |
---|
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 |
---|
510 | // if(dik<=rad) |
---|
511 | // mat[ii][yloc][jj] = fill |
---|
512 | // endif |
---|
513 | // endfor |
---|
514 | // endfor |
---|
515 | |
---|
516 | return(0) |
---|
517 | End |
---|
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? |
---|
522 | Function 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) |
---|
555 | End |
---|
556 | |
---|
557 | |
---|
558 | Function 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) |
---|
592 | End |
---|
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 |
---|
598 | Function 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 |
---|
645 | End |
---|
646 | |
---|
647 | Proc ParseMatrix3DToXYZ(matStr) |
---|
648 | String matStr="mat" |
---|
649 | ParseMatrix3D_rho($matStr) |
---|
650 | End |
---|
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 | // |
---|
698 | Function 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) |
---|
733 | End |
---|
734 | |
---|
735 | Function 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) |
---|
748 | End |
---|
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 | // |
---|
756 | Function PadMatrix(mat,num) |
---|
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] |
---|
763 | end |
---|
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 | // |
---|
771 | Function 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) |
---|
798 | End |
---|
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 | Function 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) |
---|
876 | End |
---|