1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | |
---|
3 | // VERY incomplete methods for filling volumes with connected cylinders |
---|
4 | // |
---|
5 | // most is still being worked out and verified in 2D |
---|
6 | // |
---|
7 | |
---|
8 | |
---|
9 | |
---|
10 | //These macros may or may not work - designed to help automate the |
---|
11 | //process of calculating a series of configurations. |
---|
12 | |
---|
13 | |
---|
14 | |
---|
15 | Proc ConnectDots3D(maxNumConn) |
---|
16 | Variable maxNumConn=1 |
---|
17 | |
---|
18 | Variable num=numpnts(x3d) |
---|
19 | Make/O/N=(num) numConnection3D=0 |
---|
20 | Make/O/T/N=(num) connectedTo3D="" |
---|
21 | fConnectDots3D(maxNumConn) |
---|
22 | end |
---|
23 | |
---|
24 | // connect the dots, no periodic boundary conditions |
---|
25 | //thickness is set to 1 |
---|
26 | Function fConnectDots3D(maxNumConn) |
---|
27 | Variable maxNumConn |
---|
28 | |
---|
29 | |
---|
30 | WAVE x3d=x3d |
---|
31 | WAVE y3d=y3d |
---|
32 | WAVE z3d=z3d |
---|
33 | WAVE numConnection3D = numConnection3D |
---|
34 | WAVE/T connectedTo3D = connectedTo3D |
---|
35 | WAVE mat=mat |
---|
36 | Variable num=numpnts(x3d),ii=0 |
---|
37 | Variable matSize=DimSize(mat,0),jj |
---|
38 | |
---|
39 | Variable nnX,nnY,nnZ,nnInd,nnDist //return values of the nearest neighbor |
---|
40 | Variable thick = 1 |
---|
41 | do |
---|
42 | for(ii=0;ii<num;ii+=1) |
---|
43 | if(numConnection3D[ii] < maxNumConn) |
---|
44 | FindNearestNeighbor3D(x3d,y3d,z3d,ii,maxNumConn,nnX,nnY,nnZ,nnInd,nnDist) |
---|
45 | |
---|
46 | numConnection3D[nnInd] += 1 |
---|
47 | numConnection3D[ii] += 1 |
---|
48 | connectedTo3D[nnInd] += num2str(ii)+"," |
---|
49 | connectedTo3D[ii] += num2str(nnInd)+"," |
---|
50 | |
---|
51 | ConnectPoints3D(mat, x3d[ii],y3d[ii],z3d[ii],x3d[nnInd],y3d[nnInd],z3d[nnInd],thick,1) //always fill |
---|
52 | endif |
---|
53 | endfor |
---|
54 | //without rotation, a very closed network structure tends to form (bias?) |
---|
55 | // rotate (10?) points so that we're not starting from the same place every time |
---|
56 | // -- but the "connectedTo" information is LOST... |
---|
57 | //Rotate 10, x2d,y2d, numConnection,connectedTo |
---|
58 | //Print sum(numConnection,-inf,inf), num*maxNumConn |
---|
59 | while(sum(numConnection3D,-inf,inf) < num*maxNumConn) |
---|
60 | |
---|
61 | End |
---|
62 | |
---|
63 | Function FindNearestNeighbor3D(xw,yw,zw,startPt,maxConn,nnX,nnY,nnZ,nnInd,nnDist) |
---|
64 | Wave xw,yw,zw |
---|
65 | Variable startPt,maxConn,&nnX,&nnY,&nnZ,&nnInd,&nnDist |
---|
66 | |
---|
67 | Variable num,matSize,numToFind |
---|
68 | |
---|
69 | //make waves to hold the answers |
---|
70 | num=numpnts(xw) |
---|
71 | Make/O/N=(num) distWave3D,distIndex3D |
---|
72 | |
---|
73 | //now calculate the distances |
---|
74 | Variable ii,dist,testPt |
---|
75 | |
---|
76 | //this is probably the slowest step by far in the process...(tiny XOP?) |
---|
77 | for(ii=0;ii<num;ii+=1) |
---|
78 | distWave3D[ii] = (xw[ii]-xw[startPt])^2 + (yw[ii]-yw[startPt])^2 + (zw[ii]-zw[startPt])^2 //d^2 |
---|
79 | endfor |
---|
80 | |
---|
81 | MakeIndex distWave3D distIndex3D //then distWave[distIndex[ii]] will loop through in order of distance |
---|
82 | WAVE numConnection3D = numConnection3D |
---|
83 | WAVE/T connectedTo3D = connectedTo3D |
---|
84 | WAVE mat = mat |
---|
85 | |
---|
86 | for(ii=1;ii<num;ii+=1) //[0] point is the test point, dist == 0 |
---|
87 | testPt = distIndex3D[ii] //index |
---|
88 | if(numConnection3D[testPt] < maxConn ) //can test pt accept another connection? |
---|
89 | if(WhichListItem(num2str(testPt),connectedTo3D[startPt],",",0) == -1) // not already connected to the starting point? |
---|
90 | // Print ii,testPt |
---|
91 | // Printf "nearest point is (%d,%d), dist^2 = %g\r",xw[testPt],yw[testPt],distWave[testPt] |
---|
92 | nnX = xw[testPt] //return values as pass-by-reference |
---|
93 | nnY = yw[testPt] |
---|
94 | nnZ = zw[testPt] |
---|
95 | nnInd = testPt |
---|
96 | nnDist = distWave3D[testPt] |
---|
97 | |
---|
98 | return(0) //found a point, return |
---|
99 | endif |
---|
100 | endif |
---|
101 | endfor |
---|
102 | |
---|
103 | return (1) //did not find a point, return error |
---|
104 | End |
---|
105 | |
---|
106 | ///!!! fails if the points are along an axis (or even close) - then the box becomes infinitely thin |
---|
107 | Function ConnectPoints3D_old(w, x1,y1,z1,x2,y2,z2,thick,fill) |
---|
108 | Wave w |
---|
109 | Variable x1,y1,z1,x2,y2,z2,thick,fill |
---|
110 | |
---|
111 | //find all of the points within dist=thick of the line defined by p1 and p2 |
---|
112 | // and give these points (in w) the value 1 |
---|
113 | |
---|
114 | // search through the box defined by p1 and p2 |
---|
115 | Variable ii,jj,kk,dist |
---|
116 | |
---|
117 | for(ii=min(x1,x2);ii<=max(x1,x2);ii+=1) |
---|
118 | for(jj=min(y1,y2);jj<=max(y1,y2);jj+=1) |
---|
119 | for(kk=min(z1,z2);kk<=max(z1,z2);kk+=1) |
---|
120 | //find distance between current point and P1P2 |
---|
121 | dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2) |
---|
122 | //print dist |
---|
123 | if(dist<=thick) |
---|
124 | w[ii][jj][kk] = fill //add the point |
---|
125 | endif |
---|
126 | endfor |
---|
127 | endfor |
---|
128 | endfor |
---|
129 | |
---|
130 | End |
---|
131 | |
---|
132 | // |
---|
133 | // thickness is given in voxels |
---|
134 | // |
---|
135 | // the search box is expanded by the thickness in each direction (within the box) |
---|
136 | // to account for connections that are along an axis direction |
---|
137 | // |
---|
138 | Function ConnectPoints3D(w, x1,y1,z1,x2,y2,z2,thick,fill) |
---|
139 | Wave w |
---|
140 | Variable x1,y1,z1,x2,y2,z2,thick,fill |
---|
141 | |
---|
142 | //find all of the points within dist=thick of the line defined by p1 and p2 |
---|
143 | // and give these points (in w) the value 1 |
---|
144 | |
---|
145 | // search through the box defined by p1 and p2 |
---|
146 | Variable ii,jj,kk,dist,minX,minY,minZ,maxX,maxY,maxZ |
---|
147 | Variable boxMax,slop |
---|
148 | boxMax = DimSize(w,0) - 1 //max index of box |
---|
149 | slop = thick + 0 // to make the bounding box large enough the get it all |
---|
150 | |
---|
151 | minX=min(x1,x2) //give me the smaller x |
---|
152 | minX=max(minX-slop,0) //give me zero if it's negative |
---|
153 | minY=min(y1,y2) |
---|
154 | minY=max(minY-slop,0) |
---|
155 | minZ=min(z1,z2) |
---|
156 | minZ=max(minZ-slop,0) |
---|
157 | |
---|
158 | maxX=max(x1,x2) //give me the larger x |
---|
159 | maxX=min(maxX+slop,boxMax) //give the the box dimension if it's too big |
---|
160 | maxY=max(y1,y2) |
---|
161 | maxY=min(maxY+slop,boxMax) |
---|
162 | maxZ=max(z1,z2) |
---|
163 | maxZ=min(maxZ+slop,boxMax) |
---|
164 | |
---|
165 | for(ii=minX;ii<=maxX;ii+=1) |
---|
166 | for(jj=minY;jj<=maxY;jj+=1) |
---|
167 | for(kk=minZ;kk<=maxZ;kk+=1) |
---|
168 | //find distance between current point and P1P2 |
---|
169 | dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2) |
---|
170 | //print dist |
---|
171 | if(dist<=thick) |
---|
172 | w[ii][jj][kk] = fill //add the point |
---|
173 | endif |
---|
174 | endfor |
---|
175 | endfor |
---|
176 | endfor |
---|
177 | |
---|
178 | End |
---|
179 | |
---|
180 | //solution courtesy of the math forum @ Drexel "Ask Dr. Math" |
---|
181 | Function Pt_Line_Distance3D(xx,yy,zz, x1,y1,z1,x2,y2,z2) |
---|
182 | Variable xx,yy,zz, x1,y1,z1,x2,y2,z2 |
---|
183 | |
---|
184 | Variable dist,x3,y3,z3,x4,y4,z4,aa,ab |
---|
185 | |
---|
186 | //AB |
---|
187 | x3 = x2-x1 |
---|
188 | y3 = y2-y1 |
---|
189 | z3 = z2-z1 |
---|
190 | //AP |
---|
191 | x4 = xx-x1 |
---|
192 | y4 = yy-y1 |
---|
193 | z4 = zz-z1 |
---|
194 | |
---|
195 | //aa = length of ABxAP |
---|
196 | aa = sqrt( (y3*z4-z3*y4)^2 + (x3*z4-x4*z3)^2 + (x3*y4-x4*y3)^2 ) |
---|
197 | |
---|
198 | //ab = length of AB |
---|
199 | ab = sqrt(x3^2+y3^2+z3^2) |
---|
200 | |
---|
201 | dist = aa/ab |
---|
202 | |
---|
203 | return(dist) |
---|
204 | End |
---|
205 | |
---|
206 | ////////////////////////// |
---|
207 | //// random rod fills - just a voxel at each point, so use T=diameter |
---|
208 | |
---|
209 | // fill a volume with "rods" that are made up of a numer of spheres |
---|
210 | // lined up in a random direction |
---|
211 | // |
---|
212 | // all of the spheres/rods that are added must be connected = "adjacent" to |
---|
213 | // an occupied voxel. |
---|
214 | // |
---|
215 | // this implementation is a random walk, end-to-end connection |
---|
216 | // |
---|
217 | Function ConnectedRodFill(mat,len,num,periodic) |
---|
218 | Wave mat |
---|
219 | variable len,num //length in direction, number of spheres to add |
---|
220 | Variable periodic //==1 if periodic boundaries |
---|
221 | |
---|
222 | Variable nptx,npty,nptz,ii,fail=0,nAdd |
---|
223 | Make/O/I/N=3 stVec,dirVec,endVec |
---|
224 | |
---|
225 | nptx=DimSize(mat,0) |
---|
226 | npty=DimSize(mat,1) |
---|
227 | nptz=DimSize(mat,2) |
---|
228 | |
---|
229 | //place first rod randomly |
---|
230 | stVec[0]=trunc(abs(enoise(nptx))) //random integer distr betw (0,npt) |
---|
231 | stVec[1]=trunc(abs(enoise(npty))) |
---|
232 | stVec[2]=trunc(abs(enoise(nptz))) |
---|
233 | dirVec[0] = PickDirection() |
---|
234 | dirVec[1] = PickDirection() |
---|
235 | dirVec[2] = PickDirection() |
---|
236 | |
---|
237 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic) |
---|
238 | |
---|
239 | stVec = endVec //end point is the new starting point |
---|
240 | |
---|
241 | //loop over remaining |
---|
242 | ii=nAdd |
---|
243 | do |
---|
244 | |
---|
245 | dirVec[0] = PickDirection() |
---|
246 | dirVec[1] = PickDirection() |
---|
247 | dirVec[2] = PickDirection() |
---|
248 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic) |
---|
249 | stVec = endVec //end point is the new starting point |
---|
250 | |
---|
251 | ii+=nAdd |
---|
252 | |
---|
253 | if(mod(ii,100000) == 0) |
---|
254 | print "Point = ",ii |
---|
255 | endif |
---|
256 | |
---|
257 | while(ii<num) |
---|
258 | |
---|
259 | |
---|
260 | return(0) |
---|
261 | End |
---|
262 | |
---|
263 | // fill a volume with "rods" that are made up of a numer of spheres |
---|
264 | // lined up in a random direction |
---|
265 | // |
---|
266 | // |
---|
267 | Function UnConnectedRodFill(mat,len,num,periodic) |
---|
268 | Wave mat |
---|
269 | variable len,num //length in direction, number of spheres to add |
---|
270 | Variable periodic //==1 if periodic boundaries |
---|
271 | |
---|
272 | Variable nptx,npty,nptz,ii,fail=0,nAdd |
---|
273 | Make/O/I/N=3 stVec,dirVec,endVec |
---|
274 | |
---|
275 | nptx=DimSize(mat,0) |
---|
276 | npty=DimSize(mat,1) |
---|
277 | nptz=DimSize(mat,2) |
---|
278 | |
---|
279 | //place first rod randomly |
---|
280 | stVec[0]=trunc(abs(enoise(nptx))) //random integer distr betw (0,npt) |
---|
281 | stVec[1]=trunc(abs(enoise(npty))) |
---|
282 | stVec[2]=trunc(abs(enoise(nptz))) |
---|
283 | dirVec[0] = PickDirection() |
---|
284 | dirVec[1] = PickDirection() |
---|
285 | dirVec[2] = PickDirection() |
---|
286 | |
---|
287 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic) |
---|
288 | //loop over remaining |
---|
289 | ii=nAdd |
---|
290 | do |
---|
291 | //pick a new point |
---|
292 | stVec[0]=trunc(abs(enoise(nptx))) //random integer distr betw (0,npt) |
---|
293 | stVec[1]=trunc(abs(enoise(npty))) |
---|
294 | stVec[2]=trunc(abs(enoise(nptz))) |
---|
295 | dirVec[0] = PickDirection() |
---|
296 | dirVec[1] = PickDirection() |
---|
297 | dirVec[2] = PickDirection() |
---|
298 | //add as many spheres as possible in that direction |
---|
299 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic) |
---|
300 | ii+=nAdd |
---|
301 | //Print "point ",ii |
---|
302 | |
---|
303 | while(ii<num) |
---|
304 | Print "failures = ",fail |
---|
305 | |
---|
306 | // ParseMatrix3D_rho(mat) //convert to XYZ |
---|
307 | |
---|
308 | return(0) |
---|
309 | End |
---|
310 | |
---|
311 | |
---|
312 | Function PickDirection() |
---|
313 | return( trunc(abs(enoise(3))) - 1 ) //integer -1,0,1 |
---|
314 | End |
---|
315 | |
---|
316 | //need to add periodic boundary conditions (done) |
---|
317 | // |
---|
318 | // need to insist on connectivity |
---|
319 | // |
---|
320 | // |
---|
321 | Function AddSpheresInDirection(m,stVec,dirVec,endVec,nTry,periodic) |
---|
322 | WAVE m,stVec,dirVec,endVec |
---|
323 | Variable nTry,periodic |
---|
324 | |
---|
325 | Variable added,nx,ny,nz,nDimx,nDimy,nDimz |
---|
326 | |
---|
327 | nDimx=DimSize(m,0)-1 |
---|
328 | nDimy=DimSize(m,1)-1 |
---|
329 | nDimz=DimSize(m,2)-1 |
---|
330 | added=0 |
---|
331 | nx=stVec[0] |
---|
332 | ny=stVec[1] |
---|
333 | nz=stVec[2] |
---|
334 | |
---|
335 | do |
---|
336 | nx+=dirVec[0] //coordinates on next test point to try to add in d(xyz) direction |
---|
337 | ny+=dirVec[1] |
---|
338 | nz+=dirVec[2] |
---|
339 | //check for out-of bounds - if periodic ==1, reflect to opposite face of cube, else return |
---|
340 | if( (nx>nDimx) || (nx<0) ) |
---|
341 | if(periodic) |
---|
342 | nx = abs(nDimx+1 - abs(nx)) |
---|
343 | else |
---|
344 | break |
---|
345 | endif |
---|
346 | Endif |
---|
347 | if( (ny>nDimy) || (ny<0) ) |
---|
348 | if(periodic) |
---|
349 | ny = abs(nDimy+1 - abs(ny)) |
---|
350 | else |
---|
351 | break |
---|
352 | endif |
---|
353 | Endif |
---|
354 | if( (nz>nDimz) || (nz<0) ) |
---|
355 | if(periodic) |
---|
356 | nz = abs(nDimz+1 - abs(nz)) |
---|
357 | else |
---|
358 | break |
---|
359 | endif |
---|
360 | Endif |
---|
361 | //add the point or exit |
---|
362 | if(m[nx][ny][nz]) |
---|
363 | break //point is already occupied |
---|
364 | else |
---|
365 | m[nx][ny][nz] = 1 |
---|
366 | // nx +=1 //we'll back this off before returning |
---|
367 | // ny +=1 |
---|
368 | // nz +=1 |
---|
369 | added += 1 |
---|
370 | endif |
---|
371 | while(added<nTry) |
---|
372 | |
---|
373 | endVec[0] = nx - dirVec[0] //back up to the last point actually added |
---|
374 | endVec[1] = ny - dirVec[1] |
---|
375 | endVec[2] = nz - dirVec[2] |
---|
376 | return added |
---|
377 | End |
---|
378 | |
---|
379 | Proc PutRandomCylindersAtPoints(w,num,len,periodic) |
---|
380 | String w="mat" |
---|
381 | Variable num=1000,len=10,periodic=1 |
---|
382 | Prompt w,"matrix" |
---|
383 | Prompt num,"number of starting points" |
---|
384 | prompt len,"length of cylinders" |
---|
385 | prompt periodic"1=periodic, 0=non-periodic boundaries" |
---|
386 | |
---|
387 | $w=0 |
---|
388 | RandomFill3DMat($w,num) |
---|
389 | CylindersAtPoints($w,len,periodic) |
---|
390 | |
---|
391 | NumberOfPoints() |
---|
392 | end |
---|
393 | |
---|
394 | Proc PutXAxisCylindersAtPoints(w,num,rad,len,periodic,sobol) |
---|
395 | String w="mat" |
---|
396 | Variable num=100,rad=20,len=300,periodic=1,sobol=1 |
---|
397 | Prompt w,"matrix" |
---|
398 | Prompt num,"number of starting points" |
---|
399 | prompt rad,"radius of cylinders" |
---|
400 | prompt len,"length of cylinders" |
---|
401 | prompt periodic,"1=periodic, 0=non-periodic boundaries" |
---|
402 | Prompt sobol,"1=Sobol, 0=random" |
---|
403 | |
---|
404 | $w=0 |
---|
405 | X_CylindersAtPoints($w,num,rad,len,sobol,periodic) |
---|
406 | |
---|
407 | NumberOfPoints() |
---|
408 | end |
---|
409 | |
---|
410 | Proc PutXAxisCylindersSquare(w,rad,len,sep) |
---|
411 | String w="mat" |
---|
412 | Variable rad=20,len=300,sep=80 |
---|
413 | Prompt w,"matrix" |
---|
414 | prompt rad,"radius of cylinders" |
---|
415 | prompt len,"length of cylinders" |
---|
416 | prompt sep,"center-to-center separation of cylinders" |
---|
417 | |
---|
418 | $w=0 |
---|
419 | X_CylindersSquareGrid($w,rad,len,sep) |
---|
420 | |
---|
421 | NumberOfPoints() |
---|
422 | end |
---|
423 | |
---|
424 | Proc PutXAxisCylindersHexagonal(w,rad,len,sep) |
---|
425 | String w="mat" |
---|
426 | Variable rad=20,len=300,sep=80 |
---|
427 | Prompt w,"matrix" |
---|
428 | prompt rad,"radius of cylinders" |
---|
429 | prompt len,"length of cylinders" |
---|
430 | prompt sep,"center-to-center separation of cylinders" |
---|
431 | |
---|
432 | $w=0 |
---|
433 | X_CylindersHexagonalGrid($w,rad,len,sep) |
---|
434 | |
---|
435 | NumberOfPoints() |
---|
436 | end |
---|
437 | |
---|
438 | Function CylindersAtPoints(mat,len,periodic) |
---|
439 | Wave mat |
---|
440 | variable len //length in direction |
---|
441 | Variable periodic //==1 if periodic boundaries |
---|
442 | |
---|
443 | Make/O/B/N=3 stVec,dirVec,endVec |
---|
444 | |
---|
445 | ParseMatrix3D_rho(mat) |
---|
446 | Wave x3d=x3d |
---|
447 | Wave y3d=y3d |
---|
448 | Wave z3d=z3d |
---|
449 | |
---|
450 | Variable ii=0,num=numpnts(x3d) |
---|
451 | NVAR grid=root:FFT_T |
---|
452 | |
---|
453 | Variable nAdd |
---|
454 | |
---|
455 | for(ii=0;ii<num;ii+=1) |
---|
456 | dirVec[0] = PickDirection() |
---|
457 | dirVec[1] = PickDirection() |
---|
458 | dirVec[2] = PickDirection() |
---|
459 | stVec[0] = x3d[ii] |
---|
460 | stVec[1] = y3d[ii] |
---|
461 | stVec[2] = z3d[ii] |
---|
462 | nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic) |
---|
463 | //FillSphereRadius(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],1) |
---|
464 | endfor |
---|
465 | |
---|
466 | return(0) |
---|
467 | End |
---|
468 | |
---|
469 | /////////////// |
---|
470 | |
---|
471 | Function X_CylindersAtPoints(mat,num,rad,len,sobol,periodic) |
---|
472 | Wave mat |
---|
473 | variable num,rad,len //length in direction |
---|
474 | Variable periodic //==1 if periodic boundaries |
---|
475 | Variable sobol //==1 if sobol selection of points (2D) |
---|
476 | |
---|
477 | |
---|
478 | Variable np |
---|
479 | np = DimSize(mat,0) // assumes that all dimensions are the same |
---|
480 | |
---|
481 | // fill a 2D plane with points |
---|
482 | Make/O/B/N=(np,np) plane |
---|
483 | plane = 0 |
---|
484 | if(sobol) |
---|
485 | if(WaveExists($"SobolInit") == 0) |
---|
486 | Make/O/D/N=2 SobolInit |
---|
487 | SobolX(-1,SobolInit) |
---|
488 | else |
---|
489 | SobolPoints2D(plane,num) |
---|
490 | endif |
---|
491 | else |
---|
492 | RandomPoints2D(plane,num) |
---|
493 | endif |
---|
494 | // put it in the proper plane of the matrix |
---|
495 | mat[np/2][][] = plane[q][r] // in the YZ plane |
---|
496 | |
---|
497 | ParseMatrix3D_rho(mat) |
---|
498 | Wave x3d=x3d |
---|
499 | Wave y3d=y3d |
---|
500 | Wave z3d=z3d |
---|
501 | |
---|
502 | Variable ii=0 |
---|
503 | NVAR grid=root:FFT_T |
---|
504 | |
---|
505 | for(ii=0;ii<num;ii+=1) |
---|
506 | FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1) //cylinder 1 |
---|
507 | endfor |
---|
508 | |
---|
509 | return(0) |
---|
510 | End |
---|
511 | |
---|
512 | |
---|
513 | Function X_CylindersSquareGrid(mat,rad,len,sep) |
---|
514 | Wave mat |
---|
515 | variable rad,len //length of cylinders |
---|
516 | Variable sep // EDGE separation, in same units as cylinder |
---|
517 | |
---|
518 | |
---|
519 | NVAR grid=root:FFT_T |
---|
520 | Variable np,spacing |
---|
521 | np = DimSize(mat,0) // assumes that all dimensions are the same |
---|
522 | |
---|
523 | // fill a 2D plane with points |
---|
524 | Make/O/B/N=(np,np) plane |
---|
525 | plane = 0 |
---|
526 | |
---|
527 | |
---|
528 | spacing = round(sep/grid) // so it's an integer |
---|
529 | FillPlaneSquareGrid(plane,spacing) |
---|
530 | |
---|
531 | // put it in the proper plane of the matrix |
---|
532 | mat[np/2][][] = plane[q][r] // in the YZ plane |
---|
533 | |
---|
534 | ParseMatrix3D_rho(mat) |
---|
535 | Wave x3d=x3d |
---|
536 | Wave y3d=y3d |
---|
537 | Wave z3d=z3d |
---|
538 | |
---|
539 | Variable ii=0,num |
---|
540 | num = numpnts(x3d) |
---|
541 | |
---|
542 | for(ii=0;ii<num;ii+=1) |
---|
543 | FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1) //cylinder 1 |
---|
544 | endfor |
---|
545 | |
---|
546 | return(0) |
---|
547 | End |
---|
548 | |
---|
549 | Function X_CylindersHexagonalGrid(mat,rad,len,sep) |
---|
550 | Wave mat |
---|
551 | variable rad,len //length of cylinders |
---|
552 | Variable sep // EDGE separation, in same units as cylinder |
---|
553 | |
---|
554 | |
---|
555 | NVAR grid=root:FFT_T |
---|
556 | Variable np,spacing |
---|
557 | np = DimSize(mat,0) // assumes that all dimensions are the same |
---|
558 | |
---|
559 | // fill a 2D plane with points |
---|
560 | Make/O/B/N=(np,np) plane |
---|
561 | plane = 0 |
---|
562 | |
---|
563 | spacing = round(sep/grid) // so it's an integer |
---|
564 | FillPlaneHexagonal(plane,spacing) |
---|
565 | |
---|
566 | // put it in the proper plane of the matrix |
---|
567 | mat[np/2][][] = plane[q][r] // in the YZ plane |
---|
568 | |
---|
569 | ParseMatrix3D_rho(mat) |
---|
570 | Wave x3d=x3d |
---|
571 | Wave y3d=y3d |
---|
572 | Wave z3d=z3d |
---|
573 | |
---|
574 | Variable ii=0,num |
---|
575 | num = numpnts(x3d) |
---|
576 | |
---|
577 | for(ii=0;ii<num;ii+=1) |
---|
578 | FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1) //cylinder 1 |
---|
579 | endfor |
---|
580 | |
---|
581 | //// makes a crude core-shell cylinder |
---|
582 | // for(ii=0;ii<num;ii+=1) |
---|
583 | // FillXCylinder(mat,grid,rad-10,x3d[ii],y3d[ii],z3d[ii],len,3) //cylinder 1 |
---|
584 | // endfor |
---|
585 | // |
---|
586 | return(0) |
---|
587 | End |
---|
588 | |
---|
589 | |
---|
590 | |
---|
591 | Function FillPlaneSquareGrid(w,sp) |
---|
592 | Wave w |
---|
593 | Variable sp //spacing in pixels |
---|
594 | |
---|
595 | Variable np,ii,jj |
---|
596 | |
---|
597 | np = DimSize(w,0) |
---|
598 | |
---|
599 | for(ii=sp;ii<(np-sp/2);ii+=sp) |
---|
600 | for(jj=sp;jj<(np-sp/2);jj+=sp) |
---|
601 | w[ii][jj] = 1 |
---|
602 | endfor |
---|
603 | endfor |
---|
604 | |
---|
605 | |
---|
606 | return(0) |
---|
607 | End |
---|
608 | |
---|
609 | // sp is the |
---|
610 | Function FillPlaneHexagonal(w,sp) |
---|
611 | Wave w |
---|
612 | Variable sp //spacing in pixels |
---|
613 | |
---|
614 | Variable np,ii,jj,xx,yy,fill |
---|
615 | |
---|
616 | np = DimSize(w,0) |
---|
617 | |
---|
618 | fill = 1 |
---|
619 | // do even columns, then odds |
---|
620 | for(ii=sp;ii<(np-sp/2);ii+=(sp*sqrt(3)) ) |
---|
621 | for(jj=sp;jj<(np-sp/2);jj+=sp) |
---|
622 | w[ii][jj] = fill |
---|
623 | endfor |
---|
624 | endfor |
---|
625 | |
---|
626 | // now odds |
---|
627 | for(ii=(sp+sp*sqrt(3)/2);ii<(np-sp/2);ii+=(sp*sqrt(3))) |
---|
628 | for(jj=sp/2;jj<(np-sp/2);jj+=sp) |
---|
629 | w[ii][jj] = fill |
---|
630 | endfor |
---|
631 | endfor |
---|
632 | |
---|
633 | return(0) |
---|
634 | End |
---|
635 | |
---|
636 | Proc HexagonalCylinders() |
---|
637 | |
---|
638 | Variable xc,yc,zc,len,rad,aa,cosa,sina,meanRad,pd,meanLen |
---|
639 | |
---|
640 | Variable grid=root:FFT_T |
---|
641 | Variable npt=root:FFT_N |
---|
642 | |
---|
643 | xc=npt/2 |
---|
644 | yc=npt/2 |
---|
645 | zc=npt/2 |
---|
646 | |
---|
647 | meanRad=20 |
---|
648 | pd=0.01 |
---|
649 | meanLen=400 |
---|
650 | len=meanLen |
---|
651 | |
---|
652 | // aa=meanRad/grid*2 +40 //grid units, +2 gives 2 units of separation |
---|
653 | aa=meanRad/grid*2 +2 //grid units, +2 gives 2 units of separation |
---|
654 | cosa = round(0.5*aa) |
---|
655 | sina = round(0.866*aa) |
---|
656 | |
---|
657 | |
---|
658 | Variable ii,nPass,x1,x2,x3 |
---|
659 | ii=0 |
---|
660 | nPass=1 |
---|
661 | |
---|
662 | |
---|
663 | x1 = Gauss(meanRad-pd*meanRad,meanRad,pd*meanRad) |
---|
664 | x2 = Gauss(meanRad,meanRad,pd*meanRad) |
---|
665 | x3 = Gauss(meanRad+pd*meanRad,meanRad,pd*meanRad) |
---|
666 | |
---|
667 | Print x1,x2,x3 |
---|
668 | // |
---|
669 | /////////// average over 3 point distribution -- |
---|
670 | // //calculate and sum the three contributions |
---|
671 | // FastOp mat=0 |
---|
672 | // rad = meanRad-pd*meanRad |
---|
673 | //// len = meanLen + gnoise(pd*meanLen) |
---|
674 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
675 | // Execute "DoFFT()" |
---|
676 | // Duplicate/O iBin iBin_1g_128 |
---|
677 | // iBin_1g_128 *= x1 |
---|
678 | // |
---|
679 | // FastOp mat=0 |
---|
680 | // rad = meanRad |
---|
681 | //// len = meanLen + gnoise(pd*meanLen) |
---|
682 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
683 | // Execute "DoFFT()" |
---|
684 | // iBin_1g_128 += x2*iBin |
---|
685 | // |
---|
686 | // FastOp mat=0 |
---|
687 | // rad = meanRad+pd*meanRad |
---|
688 | //// len = meanLen + gnoise(pd*meanLen) |
---|
689 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
690 | // Execute "DoFFT()" |
---|
691 | // iBin_1g_128 += x3*iBin |
---|
692 | // |
---|
693 | // iBin_1g_128 /= (x1+x2+x3) //renormalize |
---|
694 | // |
---|
695 | ////// end 3 pt average |
---|
696 | |
---|
697 | ////////////////////// |
---|
698 | // paired cylinders |
---|
699 | ///////////////////////// |
---|
700 | // //calculate and sum the three contributions |
---|
701 | // FastOp mat=0 |
---|
702 | // rad = meanRad-pd*meanRad |
---|
703 | //// len = meanLen + gnoise(pd*meanLen) |
---|
704 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
705 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
706 | // Execute "DoFFT()" |
---|
707 | // Duplicate/O iBin iBin_2g_256_sep40 |
---|
708 | // iBin_2g_256_sep40 *= x1 |
---|
709 | // |
---|
710 | // FastOp mat=0 |
---|
711 | // rad = meanRad |
---|
712 | //// len = meanLen + gnoise(pd*meanLen) |
---|
713 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
714 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
715 | // Execute "DoFFT()" |
---|
716 | // iBin_2g_256_sep40 += x2*iBin |
---|
717 | // |
---|
718 | // FastOp mat=0 |
---|
719 | // rad = meanRad+pd*meanRad |
---|
720 | //// len = meanLen + gnoise(pd*meanLen) |
---|
721 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
722 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
723 | // Execute "DoFFT()" |
---|
724 | // iBin_2g_256_sep40 += x3*iBin |
---|
725 | // |
---|
726 | // iBin_2g_256_sep40 /= (x1+x2+x3) //renormalize |
---|
727 | |
---|
728 | ///////////////////////// |
---|
729 | // |
---|
730 | // |
---|
731 | //// |
---|
732 | // ii = 0 |
---|
733 | // do |
---|
734 | // FastOp mat=0 |
---|
735 | // rad = meanRad + gnoise(pd*meanRad) |
---|
736 | // print ii,rad |
---|
737 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) |
---|
738 | // Execute "DoFFT()" |
---|
739 | // if(ii==0) |
---|
740 | // Duplicate/O iBin iBin_1g_256 |
---|
741 | // else |
---|
742 | // iBin_1g_256 += iBin |
---|
743 | // endif |
---|
744 | // ii+=1 |
---|
745 | // while(ii<nPass) |
---|
746 | // iBin_1g_256 /= nPass |
---|
747 | // |
---|
748 | ///// |
---|
749 | // ii=0 |
---|
750 | // nPass=1 |
---|
751 | // do |
---|
752 | // FastOp mat=0 |
---|
753 | // rad = meanRad + gnoise(pd*meanRad) |
---|
754 | // print ii,rad |
---|
755 | // FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) //cylinder 1 |
---|
756 | // rad = meanRad + gnoise(pd*meanRad) |
---|
757 | // print ii,rad |
---|
758 | // FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
759 | // Execute "DoFFT()" |
---|
760 | // if(ii==0) |
---|
761 | // Duplicate/O iBin iBin_2g_256 |
---|
762 | // else |
---|
763 | // iBin_2g_256 += iBin |
---|
764 | // endif |
---|
765 | // ii+=1 |
---|
766 | // while(ii<nPass) |
---|
767 | // iBin_2g_256 /= nPass |
---|
768 | // |
---|
769 | ii=0 |
---|
770 | npass=1 |
---|
771 | do |
---|
772 | FastOp mat=0 |
---|
773 | rad = meanRad + gnoise(pd*meanRad) |
---|
774 | print ii,rad |
---|
775 | FillZCylinder(mat,grid,rad,xc,yc,zc,len,1) //cylinder 1 |
---|
776 | rad = meanRad + gnoise(pd*meanRad) |
---|
777 | print ii,rad |
---|
778 | FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) //cylinder 2 |
---|
779 | rad = meanRad + gnoise(pd*meanRad) |
---|
780 | print ii,rad |
---|
781 | FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1) //cylinder 3 |
---|
782 | Execute "DoFFT()" |
---|
783 | if(ii==0) |
---|
784 | Duplicate/O iBin iBin_3g_256 |
---|
785 | else |
---|
786 | iBin_3g_256 += iBin |
---|
787 | endif |
---|
788 | ii+=1 |
---|
789 | while(ii<nPass) |
---|
790 | iBin_3g_256 /= nPass |
---|
791 | // |
---|
792 | rad = meanRad + gnoise(pd*meanRad) |
---|
793 | print rad |
---|
794 | FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1) |
---|
795 | Execute "DoFFT()" |
---|
796 | Duplicate/O iBin iBin_2 |
---|
797 | /// |
---|
798 | rad = meanRad + gnoise(pd*meanRad) |
---|
799 | print rad |
---|
800 | FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1) |
---|
801 | Execute "DoFFT()" |
---|
802 | Duplicate/O iBin iBin_3 |
---|
803 | /// |
---|
804 | rad = meanRad + gnoise(pd*meanRad) |
---|
805 | print rad |
---|
806 | FillZCylinder(mat,grid,rad,xc+sina,yc-cosa,zc,len,1) |
---|
807 | Execute "DoFFT()" |
---|
808 | Duplicate/O iBin iBin_4 |
---|
809 | /// |
---|
810 | rad = meanRad + gnoise(pd*meanRad) |
---|
811 | print rad |
---|
812 | FillZCylinder(mat,grid,rad,xc,yc-aa,zc,len,1) |
---|
813 | Execute "DoFFT()" |
---|
814 | Duplicate/O iBin iBin_5 |
---|
815 | /// |
---|
816 | rad = meanRad + gnoise(pd*meanRad) |
---|
817 | print rad |
---|
818 | FillZCylinder(mat,grid,rad,xc-sina,yc-cosa,zc,len,1) |
---|
819 | Execute "DoFFT()" |
---|
820 | Duplicate/O iBin iBin_6 |
---|
821 | /// |
---|
822 | rad = meanRad + gnoise(pd*meanRad) |
---|
823 | print rad |
---|
824 | FillZCylinder(mat,grid,rad,xc-sina,yc+cosa,zc,len,1) |
---|
825 | Execute "DoFFT()" |
---|
826 | Duplicate/O iBin iBin_7 |
---|
827 | |
---|
828 | end |
---|
829 | |
---|
830 | |
---|
831 | Proc Vary_N_In_Direction() |
---|
832 | |
---|
833 | FFTEraseMatrixButtonProc("") |
---|
834 | ConnectedRodFill(mat,2,30000*4,1) |
---|
835 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
836 | Duplicate/O ival_XOP ival_XOP_2 |
---|
837 | Duplicate/O qval_XOP qval_XOP_2 |
---|
838 | |
---|
839 | FFTEraseMatrixButtonProc("") |
---|
840 | ConnectedRodFill(mat,3,30000*4,1) |
---|
841 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
842 | Duplicate/O ival_XOP ival_XOP_3 |
---|
843 | Duplicate/O qval_XOP qval_XOP_3 |
---|
844 | |
---|
845 | FFTEraseMatrixButtonProc("") |
---|
846 | ConnectedRodFill(mat,4,30000*4,1) |
---|
847 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
848 | Duplicate/O ival_XOP ival_XOP_4 |
---|
849 | Duplicate/O qval_XOP qval_XOP_4 |
---|
850 | |
---|
851 | FFTEraseMatrixButtonProc("") |
---|
852 | ConnectedRodFill(mat,5,30000*4,1) |
---|
853 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
854 | Duplicate/O ival_XOP ival_XOP_5 |
---|
855 | Duplicate/O qval_XOP qval_XOP_5 |
---|
856 | |
---|
857 | FFTEraseMatrixButtonProc("") |
---|
858 | ConnectedRodFill(mat,6,30000*4,1) |
---|
859 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
860 | Duplicate/O ival_XOP ival_XOP_6 |
---|
861 | Duplicate/O qval_XOP qval_XOP_6 |
---|
862 | |
---|
863 | FFTEraseMatrixButtonProc("") |
---|
864 | ConnectedRodFill(mat,7,30000*4,1) |
---|
865 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
866 | Duplicate/O ival_XOP ival_XOP_7 |
---|
867 | Duplicate/O qval_XOP qval_XOP_7 |
---|
868 | |
---|
869 | FFTEraseMatrixButtonProc("") |
---|
870 | ConnectedRodFill(mat,8,30000*4,1) |
---|
871 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
872 | Duplicate/O ival_XOP ival_XOP_8 |
---|
873 | Duplicate/O qval_XOP qval_XOP_8 |
---|
874 | |
---|
875 | FFTEraseMatrixButtonProc("") |
---|
876 | ConnectedRodFill(mat,9,30000*4,1) |
---|
877 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
878 | Duplicate/O ival_XOP ival_XOP_9 |
---|
879 | Duplicate/O qval_XOP qval_XOP_9 |
---|
880 | |
---|
881 | FFTEraseMatrixButtonProc("") |
---|
882 | ConnectedRodFill(mat,10,30000*4,1) |
---|
883 | DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5) |
---|
884 | Duplicate/O ival_XOP ival_XOP_10 |
---|
885 | Duplicate/O qval_XOP qval_XOP_10 |
---|
886 | |
---|
887 | End |
---|