source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Cylinder_2D_Fill.ipf @ 798

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

I'm not sure this is a great idea, but I'm putting the FFT / Debye sphere work that I have completed into SVN. It's all really rough - the math, I believe is correct, but the interface if really, really rough. But it's not going to develop without help.

File size: 14.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4Proc ConnectDots2D(maxNumConn,thick)
5        Variable maxNumConn=1,thick=1
6       
7        Variable num=numpnts(x2d)
8        Make/O/N=(num) numConnection2D=0
9        Make/O/T/N=(num) connectedTo2D=""
10        fConnectDots2D(maxNumConn,thick)
11end
12
13// connect the dots, no periodic boundary conditions
14Function fConnectDots2D(maxNumConn, thick)
15        Variable maxNumConn, thick
16       
17       
18        WAVE x2d=x2d
19        WAVE y2d=y2d
20        WAVE numConnection2D = numConnection2D
21        WAVE/T connectedTo2D = connectedTo2D
22        WAVE plane=plane
23        Variable num=numpnts(x2d),ii=0
24        Variable matSize=DimSize(plane,0),jj
25
26        Variable nnX,nnY,nnInd,nnDist           //return values of the neareset neighbor
27
28        do
29                for(ii=0;ii<num;ii+=1)
30                        if(numConnection2D[ii] < maxNumConn)
31                                FindNearestNeighbor2D(x2d,y2d,ii,maxNumConn,nnX,nnY,nnInd,nnDist)
32                               
33                                numConnection2D[nnInd] += 1
34                                numConnection2D[ii] += 1
35                                connectedTo2D[nnInd] += num2str(ii)+","
36                                connectedTo2D[ii] += num2str(nnInd)+","
37                               
38                                ConnectPoints2D(plane, x2d[ii],y2d[ii],x2d[nnInd],y2d[nnInd],thick)
39                        endif
40                endfor
41                //without rotation, a very closed network structure tends to form (bias?)
42                // rotate (10?) points so that we're not starting from the same place every time
43                // -- but the "connectedTo" information is LOST...
44                //Rotate 10, x2d,y2d, numConnection,connectedTo
45                //Print sum(numConnection,-inf,inf), num*maxNumConn
46        while(sum(numConnection2D,-inf,inf) < num*maxNumConn)
47
48End
49
50
51///////////////////
52// find the nearest point to the input point
53//
54// The nearest point that is returned must satisfy the input criteria
55// (1) not already connected to input point
56// (2) not already connected to the maximum number of interconnects
57//
58// Input:
59//      xw,yw: waves of identical length that are the XY pairs for the nodes to connect
60//      pt: index of the trial point
61//      nConn: wave of # of connections of the XY pair
62//      connTo: text wave of (comma-delimited) indexes of the connected points
63//      matSize: dimension of the plane (assumed square)
64//      maxN: maximum number of connections
65//
66// Output:
67//      (change) nConn, connTo (at each end)
68//      Index of the XY pair to connect (as pass by reference)
69//      (?) Distance
70//
71//
72Function FindNearestNeighbor2D(xw,yw,startPt,maxConn,nnX,nnY,nnInd,nnDist)
73        Wave xw,yw
74        Variable startPt,maxConn,&nnX,&nnY,&nnInd,&nnDist
75       
76        Variable num,matSize,numToFind
77       
78        //make waves to hold the answers
79        num=numpnts(xw)
80        Make/O/N=(num) distWave2D,distIndex2D
81       
82        //now calculate the distances
83        Variable ii,dist,testPt
84       
85        //this is probably the slowest step by far in the process...(tiny XOP?)
86        for(ii=0;ii<num;ii+=1)
87                distWave2D[ii] = (xw[ii]-xw[startPt])^2 +  (yw[ii]-yw[startPt])^2               //d^2
88        endfor
89
90        MakeIndex distWave2D distIndex2D                //then distWave[distIndex[ii]] will loop through in order of distance
91        WAVE numConnection2D = numConnection2D
92        WAVE/T connectedTo2D = connectedTo2D
93        WAVE plane = plane
94       
95        for(ii=1;ii<num;ii+=1)          //[0] point is the test point, dist == 0
96                testPt =  distIndex2D[ii]               //index
97                if(numConnection2D[testPt] < maxConn )          //can test pt accept another connection?
98                        if(WhichListItem(num2str(testPt),connectedTo2D[startPt],",",0) == -1) // not already connected to the starting point?
99        //                      Print ii,testPt
100        //                      Printf "nearest point is (%d,%d), dist^2 = %g\r",xw[testPt],yw[testPt],distWave[testPt]
101                                nnX = xw[testPt]                //return values as pass-by-reference
102                                nnY = yw[testPt]
103                                nnInd = testPt
104                                nnDist = distWave2D[testPt]
105
106                                return(0)               //found a point, return
107                        endif
108                endif
109        endfor
110
111        return (1)              //did not find a point, return error
112End
113
114Proc Test2D()
115        Make/O/N=(100,100) plane
116        DoWindow/F Plane_View
117        if(V_flag==0)
118                Display/K=1 /W=(40,44,330,334)
119                DoWindow/C Plane_View
120                AppendImage  root:plane
121                ModifyGraph width={Plan,1,bottom,left},height={Plan,1,left,bottom}
122        endif
123        plane=0
124        RandomPoints2D(plane,10)
125End
126
127Proc Erase2D()
128        plane=0
129End
130
131
132//for the 3D equivalent, see RandomFill3DMat(mat,num)
133//
134Function RandomPoints2D(w,num)
135        Wave w
136        variable num            //number of spheres to add
137       
138        Variable row,col,ii,xt,yt,zt,fail=0
139       
140        row=DimSize(w,0)               
141        col=DimSize(w,1)
142       
143        ii=0
144        do
145                xt=trunc(abs(enoise(row)))              //distr betw (0,npt)
146                yt=trunc(abs(enoise(col)))
147                if( w[xt][yt] == 0 )
148                        w[xt][yt] = 1
149                        ii+=1           //increment number of spheres actually added
150                        //Print "point ",ii
151                else
152                        fail +=1
153                endif
154        while(ii<num)
155        Print "failures = ",fail
156       
157        ParseMatrix2D(w)                //convert to XY pairs
158       
159        return(0)
160End
161
162// the Sobol sequence MUST be initlalized before passing to thie routine
163Function SobolPoints2D(w,num)
164        Wave w
165        variable num            //number of spheres to add
166       
167        Variable row,col,ii,xt,yt,zt,fail=0
168       
169        Make/O/D/N=2 Sobol2D
170       
171        // initialize Sobol
172//      SobolX(-1,Sobol2D)
173        row=DimSize(w,0)               
174        col=DimSize(w,1)
175       
176        for(ii=0;ii<num;ii+=1)
177                SobolX(2,Sobol2D)
178                xt = Sobol2D[0] *row
179                yt = Sobol2D[1] *col
180                w[xt][yt] = 1
181        endfor
182       
183        ParseMatrix2D(w)                //convert to XY pairs
184       
185        return(0)
186End
187
188// parses the 2d matrix to get xy coordinates
189// CHANGE - ? Make as byte waves to save space
190//
191// XY waves are hard-wired to be "x2D,y2D"
192// overwrites previous waves
193//
194//3D equivalent is ParseMatrix3D_rho(mat)
195//
196Function ParseMatrix2D(mat)
197        Wave mat
198       
199        Variable nptx,npty,ii,jj,num
200       
201        nptx=DimSize(mat,0)
202        npty=DimSize(mat,1)     
203        Make/O/N=0 x2D,y2D,ptNum
204        num=0
205       
206        for(ii=0;ii<nptx;ii+=1)
207                for(jj=0;jj<npty;jj+=1)
208                        if(mat[ii][jj])
209                                num+=1
210                                InsertPoints num,1,x2D,y2D,ptNum
211                                x2D[num-1]=ii
212                                y2D[num-1]=jj
213                                ptNum[num-1]=num-1
214                        endif
215                endfor
216        endfor
217       
218        return(0)
219End
220
221Function ConnectPoints2D(w, x1,y1,x2,y2,thick)
222        Wave w
223        Variable  x1,y1,x2,y2,thick
224       
225        //find all of the points within dist=thick of the line defined by p1 and p2
226        // and give these points (in w) the value 1
227       
228        //find the equation of the line Ax+By+C=0 (or mx - y + b = 0)
229        Variable A,B,C
230        A = (y2-y1)/(x2-x1)
231        B = -1
232        C = y1 - A*x1
233        if (x1 == x2)   //special case of a vertical line
234                A=1
235                B=0
236                C=-x1
237        endif
238        // search through the box defined by p1 and p2
239        Variable ii,jj,dist
240       
241        for(ii=min(x1,x2);ii<=max(x1,x2);ii+=1)
242                for(jj=min(y1,y2);jj<=max(y1,y2);jj+=1)
243                        //find distance between current point and P1P2
244                        dist = Pt_Line_Distance2D(ii,jj,A,B,C)
245                        //print dist
246                        if(dist<=thick)
247                                w[ii][jj] = 1           //add the point
248                        endif
249                endfor
250        endfor
251       
252End
253
254// connect two points, maybe needing to apply periodic conditions
255// wave w is the 2d matrix (the plane)
256//
257// P1 is in, P2 MAY be a ghost point
258//
259//
260// INCOMPLETE
261Function ConnectPoints2DPeriodic(w, x1,y1,x2,y2,thick,list,jj)
262        Wave w
263        Variable  x1,y1,x2,y2,thick
264        String list             //full list of points
265        Variable jj             //index of list we're using
266       
267        //find the intersection with the edge of the box
268        //
269        // then do two fills - from p1 to edge, then from reflection of edge to p2
270        Variable xProj,yProj,xIn,yIn,wDim
271        String keyStr,xyStr
272       
273        keyStr = "XY"+num2str(jj+1)
274        xyStr = StringByKey(keyStr,list,"=",";")
275        xProj = str2num(StringFromList(0, xyStr ,","))
276        yProj = str2num(StringFromList(1, xyStr ,","))
277        wDim = DimSize(w,0)
278         xIn = (xProj >= wDim) ? 0 : 1          //== 0 if x out of bounds, 1 if in- bounds
279         xIn = (xProj < 0) ? 0 : 1
280         yIn = (yProj >= wDim) ? 0 : 1          //== 0 if y out of bounds, 1 if in- bounds
281         yIn = (yProj < 0) ? 0 : 1
282        //find the equation of the line Ax+By+C=0 (or mx - y + b = 0)
283        Variable A,B,C
284       
285        if(xIn && yIn)          //normal connection
286                A = (y2-y1)/(x2-x1)
287                B = -1
288                C = y1 - A*x1
289                if (x1 == x2)   //special case of a vertical line
290                        A=1
291                        B=0
292                        C=-x1
293                endif
294                ConnectPoints2D(w, x1,y1,x2,y2,thick)           //do the normal connection
295                return(1)
296        endif
297        //x,y or both are out-of bounds, find the intersecting point
298        Variable xInt,yInt
299        xInt = (x1+xProj)/2
300        yInt = (y1+yProj)/2
301       
302        //*****still need to find the intersecting point and do the two fills
303       
304       
305End
306
307//distance between a point (x3,y3) and a line Ax + By + C = 0
308Function Pt_Line_Distance2D(x3,y3,A,B,C)
309        Variable x3,y3,A,B,C
310       
311        Variable dist
312       
313        dist = A*x3+B*y3+C
314        dist /= sqrt(A*A+B*B)
315       
316        return(abs(dist))
317End
318
319
320Proc ConnectDots2DPeriodic()
321       
322        Variable num=numpnts(x2d),ii=0,thick=1,numToFind=1,nei=0
323        Variable matSize=DimSize(plane,0),jj
324        String list="",keyStr="",ghostStr=""
325        ii=0
326        do
327                list = FindNearestNeighborPeriodic(x2d,y2d,ii,2,matSize,numToFind)
328                print list
329                jj=0
330                do
331                        //get point from list
332                        keyStr = "N"+num2str(jj+1)
333                        nei = NumberByKey(keyStr, list , "=",";" )
334//                      keyStr = "G"+num2str(jj+1)
335//                      ghostStr = StringByKey(keyStr,list,"=",";")
336                        ConnectPoints2D(plane, x2D[ii],y2d[ii],x2D[nei],y2D[nei],thick)
337//                      ConnectPoints2DPeriodic(plane, x2D[ii],y2d[ii],x2D[nei],y2D[nei],thick,list,jj)
338                        jj+=1
339                while(jj<numToFind)
340                ii+=1
341        while(ii<num)
342//      while(ii<13)            //do only point zero
343End
344
345//
346// ??? is this 3D or 2D  ?????
347//
348// find the numToFind nearest points to the input pt
349//NEED to pass in maximum dimension of the matrix (assume it's square!)
350// a KW=val string is returned
351// returns: N1=...point number of 1st neighbor
352//                      G1 = ...ghost type of neighbor n
353//                      D1.... distance
354//                      XY1.... X,Y pair (comma delimited of the actual point, may be a ghost)
355Function/S FindNearestNeighborPeriodic(xw,yw,pt,num,matSize,numToFind)
356        Wave xw,yw
357        Variable pt,num,matSize,numToFind
358               
359        // look over a box +/-num around pt
360        // using periodic conditions
361        String ghostStr="",retStr="",pointList="",xyStr=""
362        Variable boxDim
363       
364        // get a list of local points
365        boxDim = num
366        do
367                pointList=ListPointsInBox(xw,yw,pt,boxDim,matSize)
368                print pointList
369                if ( ItemsInList(pointList ,";") >= numToFind )
370                        break
371                endif
372                boxDim += 5
373        while(boxDim < matSize)
374        Print "boxDim = ",boxDim
375       
376        //make waves to hold the answers
377        Make/O/N=(ItemsInList(pointList ,";")) ptWave,distWave
378        Make/O/T/N=(ItemsInList(pointList ,";")) ghostWave,xyWave
379        //now calculate the distances
380        Variable ii,dist,testPt
381
382        for(ii=0;ii<ItemsInList(pointList ,";");ii+=1)
383                testPt = str2num(StringFromList(ii,pointList,";"))
384                if(testPt != pt)
385                        dist = Point2PointPeriodic(xw[pt],yw[pt],xw[testPt],yw[testPt],matSize,ghostStr,xyStr)
386                        ptWave[ii] = testPt
387                        distWave[ii] = dist
388                        ghostWave[ii] = ghostStr
389                        xyWave[ii] = xyStr
390                endif
391        endfor
392        //sort by distance
393        Sort distWave distWave,ptWave,ghostWave,xyWave
394        // build the return string
395        for(ii=0;ii<numToFind;ii+=1)
396                retStr += "N"+num2str(ii+1)+"="+num2str(ptWave[ii])+";"         //point number
397                retStr += "D"+num2str(ii+1)+"="+num2str(distWave[ii])+";"       //distance
398                retStr += "G"+num2str(ii+1)+"="+ghostWave[ii]+";"
399                retStr += "XY"+num2str(ii+1)+"="+xyWave[ii]+";"                         //xy point, may be ghost
400        endfor
401       
402        return (retStr)
403End
404
405Function/S ListPointsInBox(xw,yw,pt,num,matSize)
406        Wave xw,yw
407        Variable pt,num,matSize
408       
409        String list=""
410        Variable left,right,top,bottom
411        Variable ghostLeft,ghostRight,ghostTop,ghostBottom
412        Variable ii
413       
414        left = xw[pt] - num
415        if (left<0)
416                ghostLeft=1
417                left=0
418        endif
419        right = xw[pt]+num
420        if(right>matSize-1)
421                ghostRight=1
422                right = matSize-1
423        Endif
424        bottom = yw[pt] - num
425        if(bottom<0)
426                ghostBottom=1
427                bottom=0
428        endif
429        top = yw[pt] + num
430        if(top>matSize-1)
431                ghostTop=1
432                top= matSize-1
433        Endif
434       
435        //do the bit in the normal matrix using FindPointsInPoly
436        Make/O/N=4 tempx={left,right,right,left}
437        Make/O/N=4 tempy={bottom,bottom,top,top}
438        list = GetPointsInBox(xw,yw,tempx,tempy,list)
439       
440        //do the "ghost regions" - there will be either 0,1, or 3
441        //do each individually
442        if(ghostLeft)
443                Make/O/N=4 tempx={xw[pt]-num+matSize,matSize,matSize,xw[pt]-num+matSize}
444                Make/O/N=4 tempy={bottom,bottom,top,top}
445                list = GetPointsInBox(xw,yw,tempx,tempy,list)
446        endif
447        if(ghostRight)
448                Make/O/N=4 tempx={0,xw[pt] +num - matSize, xw[pt] +num - matSize,0}
449                Make/O/N=4 tempy={bottom,bottom,top,top}
450                list = GetPointsInBox(xw,yw,tempx,tempy,list)
451        endif
452        if(ghostBottom)
453                Make/O/N=4 tempx={left,right,right,left}
454                Make/O/N=4 tempy={yw[pt]-num+matSize,yw[pt]-num+matSize,matSize,matSize}
455                list = GetPointsInBox(xw,yw,tempx,tempy,list)
456        endif
457        if(ghostTop)
458                Make/O/N=4 tempx={left,right,right,left}
459                Make/O/N=4 tempy={0,0,yw[pt] +num - matSize,yw[pt] +num - matSize}
460                list = GetPointsInBox(xw,yw,tempx,tempy,list)
461        endif
462        // then do the diagonal, if necessary
463        if(ghostLeft && ghostBottom)   
464                Make/O/N=4 tempx={xw[pt]-num+matSize,matSize,matSize,xw[pt]-num+matSize}
465                Make/O/N=4 tempy={yw[pt]-num+matSize,yw[pt]-num+matSize,matSize,matSize}
466                list = GetPointsInBox(xw,yw,tempx,tempy,list)
467        endif
468        if(ghostLeft && ghostTop)       
469                Make/O/N=4 tempx={xw[pt]-num+matSize,matSize,matSize,xw[pt]-num+matSize}
470                Make/O/N=4 tempy={0,0,yw[pt] +num - matSize,yw[pt] +num - matSize}
471                list = GetPointsInBox(xw,yw,tempx,tempy,list)
472        endif
473        if(ghostRight && ghostBottom)   
474                Make/O/N=4 tempx={0,xw[pt] +num - matSize, xw[pt] +num - matSize,0}
475                Make/O/N=4 tempy={yw[pt]-num+matSize,yw[pt]-num+matSize,matSize,matSize}
476                list = GetPointsInBox(xw,yw,tempx,tempy,list)
477        endif
478        if(ghostRight && ghostTop)
479                Make/O/N=4 tempx={0,xw[pt] +num - matSize, xw[pt] +num - matSize,0}
480                Make/O/N=4 tempy={0,0,yw[pt] +num - matSize,yw[pt] +num - matSize}
481                list = GetPointsInBox(xw,yw,tempx,tempy,list)
482        endif
483        //remove the original point
484        list=RemoveFromList(num2str(pt), list ,";")
485        return(list)
486End
487
488Function/S GetPointsInBox(xw,yw,tempx,tempy,list)
489        Wave xw,yw,tempx,tempy
490        String list
491       
492        FindPointsInPoly xw,yw,tempx,tempy
493        Wave W_inPoly=W_inPoly
494       
495        Variable ii
496        //scan through all the points
497        for(ii=0;ii<numpnts(xw);ii+=1)
498                if(W_inPoly[ii] == 1)
499                        list += num2str(ii) + ";"
500                endif
501        endfor
502       
503        return list
504End
505
506//find the distance between two points, allowing for periodic boundaries
507// on a square grid of maxSize x maxSize
508//
509// consider x1 y1 as the fixed point
510//
511Function Point2PointPeriodic(x1,y1,x2,y2,maxSize,ghostStr,xyStr)
512        Variable x1,y1,x2,y2,maxSize
513        String &ghostStr,&xyStr
514       
515        Variable dist,dx,dy,tx,ty
516        String tempStr=""
517       
518        dx = abs(x1-x2)
519        dy = abs(y1-y2)
520       
521        // if dx is more than half of maxsize, pass through x (left) (right)
522        tx = x2
523        if (dx >= maxSize/2)
524                if(x1 <= maxSize/2)
525                        // ghost x2 to the left
526                        tx = x2 - maxSize
527                        tempStr += "left"
528                else
529                        //ghost x2 to the right
530                        tx = x2 + maxSize
531                        tempStr += "right"
532                endif
533        endif
534        // if dy is more than half of maxSize, pass through y (bottom) (top)
535        ty = y2
536        if (dy >= maxSize/2)
537                if(y1 <= maxSize/2)
538                        // ghost y2 to the bottom
539                        ty = y2 - maxSize
540                        tempStr += "bottom"
541                else
542                        //ghost y2 to the top
543                        ty = y2 + maxSize
544                        tempStr += "top"
545                endif
546        endif
547        //calculate the distance
548        dist = sqrt( (tx-x1)^2 + (ty-y1)^2 )
549//      Print "point2 = ",tx,ty,dist
550//      Print tempStr
551        ghostStr = tempStr
552        xyStr = num2str(tx)+","+num2str(ty)
553        return (dist)
554end
Note: See TracBrowser for help on using the repository browser.