Functions | Variables
test_distance Namespace Reference

Functions

def test_dAlgo ()
 

Variables

int _epsilon = 1
 
string OK = '\033[92m'
 
string NO = '\033[91m'
 
string BLUE = '\033[94m'
 
string ENDC = '\033[0m'
 

Function Documentation

def test_distance.test_dAlgo ( )

Definition at line 17 of file test_distance.py.

17 def test_dAlgo():
18 
19  debug()
20  debug(BLUE + "Precision Being Required to Consider Two numbers Equal: {0:.2e}".format(_epsilon) + ENDC)
21  debug()
22 
23  # number of times to test each function
24  tests = 10000
25 
26  # import Distance Algo
27  dAlgo = geoalgo.GeoAlgo()
28 
29  try:
30 
31  # test distance from point to infinite line
32  info('Testing Point & Infinite Line Distance')
33  totSuccess = 0
34  sqdistT = 0
35  closestT = 0
36  for y in range(tests):
37  success = 1
38  # generate a random point (will be point on line closest to external point
39  p1 = geoalgo.Vector(random(),random(),random())
40  # second point on line
41  p2 = geoalgo.Vector(random(),random(),random())
42  # generate random line
43  l = geoalgo.Line(p1,p2)
44  # generate a point in a direction transverse to the line
45  transX = random()
46  transY = random()
47  # now ensure perpendicularity by having dot product be == 0
48  dirct = l.Pt2()-l.Pt1()
49  transZ = (-transX*dirct[0]-transY*dirct[1])/dirct[2]
50  vectTranslate = geoalgo.Vector(transX,transY,transZ)
51  # generate point away starting from p1
52  pt = geoalgo.Vector(p1+vectTranslate)
53  # answer will be vectTranslate sq. length
54  answer = vectTranslate.SqLength()
55  # closest point will be p1
56  tim = time()
57  a1 = dAlgo.SqDist(pt,l)
58  sqdistT += (time()-tim)
59  a2 = dAlgo.SqDist(l,pt)
60  tim = time()
61  pAnswer1 = dAlgo.ClosestPt(pt,l)
62  closestT += (time()-tim)
63  pAnswer2 = dAlgo.ClosestPt(l,pt)
64  if not ( abs(answer-a1) < _epsilon ) : success = 0
65  if not ( abs(answer-a2) < _epsilon ) : success = 0
66  for x in xrange(3):
67  if not ( abs(p1[x]-pAnswer1[x]) < _epsilon) : success = 0
68  if not ( abs(p1[x]-pAnswer2[x]) < _epsilon) : success = 0
69  totSuccess += success
70  if ( float(totSuccess)/tests < 1):
71  info(NO + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
72  else:
73  info(OK + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
74  info("Time for SqDist : {0:.3f} us".format(1E6*sqdistT/tests))
75  info("Time for ClosestPt : {0:.3f} us".format(1E6*closestT/tests))
76 
77  # test distance from point to segment
78  info('Testing Point & LineSegment Distance')
79  totSuccess = 0
80  sqdistT_out = 0.
81  closestT_out = 0.
82  sqdistT_in = 0.
83  closestT_in = 0.
84  for y in range(tests):
85  success = 1
86  # generate a random segment
87  l = geoalgo.LineSegment(random(),random(),random(),random(),random(),random())
88  # get the segment length
89  lLen = l.Dir().Length()
90  # get the segment direction
91  d = l.Dir()
92  # generate a segment parallel to it
93  # to get a line parallel to the first,
94  # translate the Start and End points
95  # by some amount in a direction perpendicular
96  # to that of the original line
97  transX = random()
98  transY = random()
99  # now ensure perpendicularity by having dot product be == 0
100  transZ = (-transX*d[0]-transY*d[1])/d[2]
101  vectTranslate = geoalgo.Vector(transX,transY,transZ)
102  p1 = l.Start()+vectTranslate
103  p2 = l.End()+vectTranslate
104  # parallel segment:
105  lPar = geoalgo.LineSegment(p1,p2)
106 
107  # first, test a point that is "before start point"
108  # distance to this point should be distance between lines
109  # in quadrature with how much further from start point we go
110  dist = random()
111  # direction vector outwards from line
112  dirOut = lPar.Dir()*(-1*dist/lPar.Dir().Length())
113  pTest = p1+dirOut
114  answer = dirOut.SqLength()+vectTranslate.SqLength()
115  tim = time()
116  a1 = dAlgo.SqDist(pTest,l)
117  sqdistT_out += (time() - tim)
118  a2 = dAlgo.SqDist(l,pTest)
119  if not ( abs(answer-a1) < _epsilon): success = 0
120  if not ( abs(answer-a2) < _epsilon): success = 0
121  tim = time()
122  point1 = dAlgo.ClosestPt(pTest,l)
123  closestT_out += (time() - tim)
124  point2 = dAlgo.ClosestPt(l,pTest)
125  for x in xrange(3):
126  if not ( (l.Start()[x]-point1[x]) < _epsilon ) : success = 0
127  if not ( (l.Start()[x]-point2[x]) < _epsilon ) : success = 0
128 
129  # now, test a point inside the segment.
130  # distance between point & segment should be
131  # pre-determined distance between the two segments
132  dist = random()
133  dirIn = lPar.Dir()*dist #dist ensures < distance of full segment
134  pTest = p1+dirIn
135  answer = vectTranslate.SqLength()
136  tim = time()
137  a1 = dAlgo.SqDist(pTest,l)
138  sqdistT_in += (time() - tim)
139  a2 = dAlgo.SqDist(l,pTest)
140  if not (abs(answer-a1) < _epsilon): success = 0
141  if not (abs(answer-a2) < _epsilon): success = 0
142  pAns = l.Start()+dirIn
143  tim = time()
144  point1 = dAlgo.ClosestPt(pTest,l)
145  closestT_in += (time() - tim)
146  point2 = dAlgo.ClosestPt(l,pTest)
147  for x in xrange(3):
148  if not ( (pAns[x]-point1[x]) < _epsilon ) : success = 0
149  if not ( (pAns[x]-point2[x]) < _epsilon ) : success = 0
150 
151  # now test a point beyond the segment
152  dist = random()
153  dirOut = lPar.Dir()*(dist/lPar.Dir().Length())
154  pTest = p2+dirOut
155  answer = dirOut.SqLength()+vectTranslate.SqLength()
156  if not ( abs(answer-dAlgo.SqDist(pTest,l)) < _epsilon): success = 0
157  if not ( abs(answer-dAlgo.SqDist(l,pTest)) < _epsilon): success = 0
158  point1 = dAlgo.ClosestPt(pTest,l)
159  point2 = dAlgo.ClosestPt(pTest,l)
160  for x in xrange(3):
161  if not ( (l.End()[x]-point1[x]) < _epsilon ) : success = 0
162  if not ( (l.End()[x]-point2[x]) < _epsilon ) : success = 0
163 
164  if (success == 1) : totSuccess += 1
165  if ( float(totSuccess)/tests < 1):
166  info(NO + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
167  else:
168  info(OK + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
169  info("Time for SqDist (Pt Out of Segment) : {0:.3f} us".format(1E6*sqdistT_out/tests))
170  info("Time for ClosestPt (Pt Out of Segment): {0:.3f} us".format(1E6*closestT_out/tests))
171  info("Time for SqDist (Pt In Segment) : {0:.3f} us".format(1E6*sqdistT_in/tests))
172  info("Time for ClosestPt (Pt In Segment) : {0:.3f} us".format(1E6*closestT_in/tests))
173 
174  # test Point to HalfLine distance
175  debug('Testing Point & HalfLine Distance')
176  success = 1
177  totSuccess = 0
178  sqdistT_out = 0.
179  closestT_out = 0.
180  sqdistT_in = 0.
181  closestT_in = 0.
182  for x in range(tests):
183  # generate a random segment
184  l = geoalgo.HalfLine(random(),random(),random(),random(),random(),random())
185  # generate a segment parallel to it
186  # to get a line parallel to the first,
187  # translate the Start and End points
188  # by some amount in a direction perpendicular
189  # to that of the original line
190  transX = random()
191  transY = random()
192  # now ensure perpendicularity by having dot product be == 0
193  transZ = (-transX*l.Dir()[0]-transY*l.Dir()[1])/l.Dir()[2]
194  vectTranslate = geoalgo.Vector(transX,transY,transZ)
195  # create the new translated & parallel hal-fline
196  lPar = geoalgo.HalfLine(l.Start()+vectTranslate, l.Dir())
197  # first, test a point that is "before start point"
198  # distance to this point should be distance between lines
199  # in quadrature with how much further from start point we go
200  dist = random()
201  # direction vector outwards from line
202  dirOut = lPar.Dir()*(-1*dist)
203  pTest = lPar.Start()+dirOut
204  answer = dirOut.SqLength()+vectTranslate.SqLength()
205  tim = time()
206  a1 = dAlgo.SqDist(pTest,l)
207  tim = time() - tim
208  sqdistT_out += tim
209  a2 = dAlgo.SqDist(l,pTest)
210  if not ( abs(answer-a1) < _epsilon): success = 0
211  if not ( abs(answer-a2) < _epsilon): success = 0
212  tim = time()
213  point1 = dAlgo.ClosestPt(pTest,l)
214  tim = time() - tim
215  closestT_out += tim
216  point2 = dAlgo.ClosestPt(l,pTest)
217  for x in xrange(3):
218  if not ( (l.Start()[x]-point1[x]) < _epsilon ) : success = 0
219  if not ( (l.Start()[x]-point2[x]) < _epsilon ) : success = 0
220 
221  # now, test a point inside the segment.
222  # distance between point & segment should be
223  # pre-determined distance between the two segments
224  dist = random()
225  dirIn = lPar.Dir()*dist #dist ensures < distance of full segment
226  pTest = lPar.Start()+dirIn
227  answer = vectTranslate.SqLength()
228  pAns = l.Start()+dirIn
229  tim = time()
230  a1 = dAlgo.SqDist(pTest,l)
231  tim = time() - tim
232  sqdistT_in += tim
233  a2 = dAlgo.SqDist(l,pTest)
234  if not ( abs(answer-a1) < _epsilon): success = 0
235  if not ( abs(answer-a2) < _epsilon): success = 0
236  tim = time()
237  point1 = dAlgo.ClosestPt(pTest,l)
238  tim = time() - tim
239  closestT_in += tim
240  point2 = dAlgo.ClosestPt(l,pTest)
241  for x in xrange(3):
242  if not ( (pAns[x]-point1[x]) < _epsilon ) : success = 0
243  if not ( (pAns[x]-point2[x]) < _epsilon ) : success = 0
244 
245  if (success == 1) : totSuccess += 1
246  if ( float(totSuccess)/tests < 1):
247  info(NO + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
248  else:
249  info(OK + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
250  info("Time for SqDist (Pt Out of Segment) : {0:.3f} us".format(1E6*sqdistT_out/tests))
251  info("Time for ClosestPt (Pt Out of Segment): {0:.3f} us".format(1E6*closestT_out/tests))
252  info("Time for SqDist (Pt In Segment) : {0:.3f} us".format(1E6*sqdistT_in/tests))
253  info("Time for ClosestPt (Pt In Segment) : {0:.3f} us".format(1E6*closestT_in/tests))
254 
255  # test Distance between two Infinite Lines
256  debug('Testing Inf Line & Inf Line Distance')
257  totSuccess = 0
258  sqdistT = 0.
259  closestT = 0.
260  for y in range(tests):
261  success = 1
262  l1 = geoalgo.Line(random(),random(),random(),random(),random(),random())
263  # take a point a fixed distance away
264  # generate a random direction in the plane
265  # perpendicular to the plane connecting
266  # the point and the line
267  # the distance between the two lines
268  # should be the fixed amount selected previously
269  # use half-way point to do calculations
270  p1 = (l1.Pt2()+l1.Pt1())/2
271  # get line direction
272  d1 = (l1.Pt2()-l1.Pt1())
273  # move in a random direction perpendicular to line
274  dirx = random()
275  diry = random()
276  dirz = (-dirx*d1[0]-diry*d1[1])/d1[2]
277  vectTranslate = geoalgo.Vector(dirx,diry,dirz)
278  # need to re-orient in some random direction on this plane
279  # this direction has to be perpendicular to both
280  # the line's direction as well as to the direction
281  # of the segment uniting the two lines
282  # use cross-product
283  vectRotate = vectTranslate.Cross(d1)
284  l2 = geoalgo.Line(p1+vectTranslate,p1+vectTranslate+vectRotate)
285  # now calculate distance. Should be vectTranslate.Length()
286  answer = vectTranslate.SqLength()
287  tim = time()
288  a1 = dAlgo.SqDist(l1,l2)
289  tim = time() - tim
290  sqdistT += tim
291  # expect the closest points on both lines to be p1 & l2.Pt1()
292  ptL1 = geoalgo.Vector()
293  ptL2 = geoalgo.Vector()
294  a2 = dAlgo.SqDist(l1,l2,ptL1,ptL2)
295  if not (abs(answer-a1) < _epsilon): success = 0
296  if not (abs(answer-a2) < _epsilon) : success = 0
297  for x in xrange(3):
298  if not ( abs(ptL1[x]-p1[x]) < _epsilon ) : success = 0
299  if not ( abs(ptL2[x]-l2.Pt1()[x]) < _epsilon ) : success = 0
300 
301  if (success == 1) : totSuccess += 1
302  if ( float(totSuccess)/tests < 1):
303  info(NO + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
304  else:
305  info(OK + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
306  info("Time for SqDist : {0:.3f} us".format(1E6*sqdistT/tests))
307 
308  # test Distance between two Half-Infinite Lines
309  debug('Testing Half-Inf Line & Half-Inf Line Distance')
310  totSuccess = 0
311  sqdistT_in = 0
312  timesIN = 0
313  sqdistT_out = 0
314  timesOUT = 0
315  for y in range(tests):
316  success = 1
317  l1 = geoalgo.HalfLine(random(),random(),random(),random(),random(),random())
318  # take a point a fixed distance away
319  # then generate a new half-line starting
320  # at the same point (translated) and with
321  # the same (or opposite) direction.
322  # But rotated outward a bit
323  p1 = l1.Start()
324  # get line direction
325  d1 = l1.Dir()
326  # move in a random direction perpendicular to line
327  dirx = random()
328  diry = random()
329  dirz = (-dirx*d1[0]-diry*d1[1])/d1[2]
330  vectTranslate = geoalgo.Vector(dirx,diry,dirz)
331  # pick some random direction (aligned with 1st or not)
332  aligned = -1
333  if (random() < 0.5) : aligned = 1
334  l2 = geoalgo.HalfLine(p1+vectTranslate,d1+vectTranslate*aligned)
335  # now calculate distance.
336  # if aligned == -1 distance should be == 0 (lines intersect)
337  # if aligned == 1 then the two start points are the closest points
338  if (aligned == -1):
339  timesIN += 1
340  answer = 0.
341  tim = time()
342  a1 = dAlgo.SqDist(l1,l2)
343  tim = time() - tim
344  L1 = geoalgo.Vector()
345  L2 = geoalgo.Vector()
346  a2 = dAlgo.SqDist(l1,l2,L1,L2)
347  sqdistT_in += tim
348  if not (abs(answer-a1) < _epsilon): success = 0
349  if not (abs(answer-a2) < _epsilon): success = 0
350  for x in xrange(3):
351  if not (L1[x]-L2[x] < _epsilon) : success = 0
352  if (aligned == 1):
353  timesOUT += 1
354  answer = vectTranslate.SqLength()
355  tim = time()
356  a1 = dAlgo.SqDist(l1,l2)
357  tim = time() - tim
358  L1 = geoalgo.Vector()
359  L2 = geoalgo.Vector()
360  a2 = dAlgo.SqDist(l1,l2,L1,L2)
361  sqdistT_out += tim
362  if not (abs(answer-a1) < _epsilon): success = 0
363  if not (abs(answer-a2) < _epsilon): success = 0
364  for x in xrange(3):
365  if not (L1[x]-l1.Start()[x] < _epsilon) : success = 0
366  if not (L2[x]-l2.Start()[x] < _epsilon) : success = 0
367 
368  if (success == 1) : totSuccess += 1
369  if ( float(totSuccess)/tests < 1):
370  info(NO + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
371  else:
372  info(OK + "Success: {0}%".format(100*float(totSuccess)/tests) + ENDC)
373  info("Time for SqDist (OUT) : {0:.3f} us".format(1E6*sqdistT_out/timesOUT))
374  info("Time for SqDist (IN) : {0:.3f} us".format(1E6*sqdistT_in/timesIN))
375 
376 
377  # test Distance between Half-Line and Line Segment
378  debug('Testing Half Line & Line Segment Distance')
379  totSuccess1 = 0
380  sqdistT1 = 0
381  totSuccess2 = 0
382  sqdistT2 = 0
383  totSuccess3 = 0
384  sqdistT3 = 0
385  totSuccess4 = 0
386  sqdistT4 = 0
387  for y in range(tests):
388  success1 = 1
389  success2 = 1
390  success3 = 1
391  success4 = 1
392  l1 = geoalgo.HalfLine(random(),random(),random(),random(),random(),random())
393  # take a point a fixed distance away
394  # test multiple scenarios:
395  # 1) segment "away" from half-line but same direction
396  # -> closest points are half-line start and segment end
397  # 2) segment "away" from half-line and rotated
398  # -> closest points are half-line start and segment rotation pivot
399  # 3) segment "close to" half-line and tilted outwards
400  # -> closest points are point on half-line and segment start
401  # 4) segment "close to" half-line and rotated
402  # -> closest points are point on half-line and segment rotation pivot
403  p1 = l1.Start()
404  # get line direction
405  d1 = l1.Dir()
406  # move in a random direction perpendicular to line
407  dirx = random()
408  diry = random()
409  dirz = (-dirx*d1[0]-diry*d1[1])/d1[2]
410  vectTranslate = geoalgo.Vector(dirx,diry,dirz)
411  # 1) generate line-segment "behind" half-line
412  # use unit-vector nature of Dir() to go back slightly more
413  dist = random()
414  seg = geoalgo.LineSegment(p1+vectTranslate-d1*dist,p1+vectTranslate-d1*dist-d1)
415  # distance should be dist
416  # closest points should be seg.Start() and l1.Start()
417  answer = vectTranslate*vectTranslate+dist*dist
418  tim = time()
419  a1 = dAlgo.SqDist(l1,seg)
420  tim = time() - tim
421  sqdistT1 += tim
422  L1 = geoalgo.Vector()
423  L2 = geoalgo.Vector()
424  a2 = dAlgo.SqDist(l1,seg,L1,L2)
425  if not (abs(answer-a1) < _epsilon): success1 = 0
426  if not (abs(answer-a2) < _epsilon): success1 = 0
427  for x in xrange(3):
428  if not (L1[x]-l1.Start()[x] < _epsilon) : success1 = 0
429  if not (L2[x]-seg.Start()[x] < _epsilon) : success1 = 0
430  if (success1 == 1) : totSuccess1 += 1
431  # 2) generate line segment away but rotated
432  # rotate in a direction perpendicular to both line direction and translation direction
433  vectRotate = vectTranslate.Cross(d1)
434  # choose pivot point
435  pivot = p1+vectTranslate-d1*dist
436  seg = geoalgo.LineSegment(pivot-vectRotate*random(),pivot+vectRotate*random())
437  # distance should be pivot distance to half-line start
438  # closest points should be pivot and line start point
439  answer = vectTranslate*vectTranslate+dist*dist
440  tim = time()
441  a1 = dAlgo.SqDist(l1,seg,L1,L2)
442  sqdistT2 += (time()-tim)
443  a2 = dAlgo.SqDist(seg,l1)
444  if not (abs(answer-a1) < _epsilon): success2 = 0
445  if not (abs(answer-a2) < _epsilon): success2 = 0
446  for x in xrange(3):
447  if not (L1[x]-l1.Start()[x] < _epsilon) : success2 = 0
448  if not (L2[x]-pivot[x] < _epsilon) : success2 = 0
449  if (success2 == 1) : totSuccess2 += 1
450  # 3) generate line segment next to but tilted outwards
451  seg = geoalgo.LineSegment(p1+vectTranslate+d1*dist,p1+vectTranslate+d1*dist+vectTranslate*random())
452  # distance should be vectTranslate
453  # closest point should be segment start and point on line "d1*dist" away from line-start
454  answer = vectTranslate*vectTranslate
455  tim = time()
456  a1 = dAlgo.SqDist(l1,seg,L1,L2)
457  sqdistT3 += (time()-tim)
458  a2 = dAlgo.SqDist(seg,l1)
459  if not (abs(answer-a1) < _epsilon): success3 = 0
460  if not (abs(answer-a2) < _epsilon): success3 = 0
461  ptLine = l1.Start()+d1*dist
462  for x in xrange(3):
463  if not (L1[x]-ptLine[x] < _epsilon) : success3 = 0
464  if not (L2[x]-seg.Start()[x] < _epsilon) : success3 = 0
465  if (success3 == 1) : totSuccess3 += 1
466  # 4) generate line segment next to line but rotated
467  pivot = p1+vectTranslate+d1*dist
468  seg = geoalgo.LineSegment(pivot-vectRotate*random(),pivot+vectRotate*random())
469  # distance should be vectTranslate
470  # closest point should be pivot on segment and d1*dist away from start for line
471  answer = vectTranslate*vectTranslate
472  tim = time()
473  a1 = dAlgo.SqDist(l1,seg,L1,L2)
474  sqdistT4 += (time()-tim)
475  a2 = dAlgo.SqDist(seg,l1)
476  if not (abs(answer-a1) < _epsilon): success4 = 0
477  if not (abs(answer-a2) < _epsilon): success4 = 0
478  ptLine = l1.Start()+d1*dist
479  for x in xrange(3):
480  if not (L1[x]-ptLine[x] < _epsilon) : success4 = 0
481  if not (L2[x]-pivot[x] < _epsilon) : success4 = 0
482  if (success4 == 1) : totSuccess4 += 1
483 
484  if ( float(totSuccess1)/tests < 1):
485  info(NO + "Success: {0}%".format(100*float(totSuccess1)/tests) + ENDC)
486  else:
487  info(OK + "Success: {0}%".format(100*float(totSuccess1)/tests) + ENDC)
488  info("Time for SqDist (Case 1) : {0:.3f} us".format(1E6*sqdistT1/tests))
489  if ( float(totSuccess2)/tests < 1):
490  info(NO + "Success: {0}%".format(100*float(totSuccess2)/tests) + ENDC)
491  else:
492  info(OK + "Success: {0}%".format(100*float(totSuccess2)/tests) + ENDC)
493  info("Time for SqDist (Case 2) : {0:.3f} us".format(1E6*sqdistT2/tests))
494  if ( float(totSuccess3)/tests < 1):
495  info(NO + "Success: {0}%".format(100*float(totSuccess3)/tests) + ENDC)
496  else:
497  info(OK + "Success: {0}%".format(100*float(totSuccess3)/tests) + ENDC)
498  info("Time for SqDist (Case 3) : {0:.3f} us".format(1E6*sqdistT3/tests))
499  if ( float(totSuccess4)/tests < 1):
500  info(NO + "Success: {0}%".format(100*float(totSuccess4)/tests) + ENDC)
501  else:
502  info(OK + "Success: {0}%".format(100*float(totSuccess4)/tests) + ENDC)
503  info("Time for SqDist (Case 4) : {0:.3f} us".format(1E6*sqdistT4/tests))
504 
505  # test Distance between Half-Line and Line Segment
506  debug('Testing Line Segment & Line Segment Distance')
507  totSuccess1 = 0.
508  sqdistT1 = 0
509  totSuccess2 = 0.
510  sqdistT2 = 0
511  totSuccess3 = 0.
512  sqdistT3 = 0
513  totSuccess4 = 0.
514  sqdistT4 = 0
515  for y in range(tests):
516  success1 = 1
517  success2 = 1
518  success3 = 1
519  success4 = 1
520  seg1 = geoalgo.LineSegment(random(),random(),random(),random(),random(),random())
521  # take a point a fixed distance away
522  # test multiple scenarios:
523  # 1) segment "away" from half-line but same direction
524  # -> closest points are half-line start and segment end
525  # 2) segment "away" from half-line and rotated
526  # -> closest points are half-line start and segment rotation pivot
527  # 3) segment "close to" half-line and tilted outwards
528  # -> closest points are point on half-line and segment start
529  # 4) segment "close to" half-line and rotated
530  # -> closest points are point on half-line and segment rotation pivot
531  p1 = seg1.Start()
532  # get line direction
533  p2 = seg1.End()
534  d = p2-p1
535  length = d.Length()
536  # move in a random direction perpendicular to line
537  dirx = random()
538  diry = random()
539  dirz = (-dirx*d[0]-diry*d[1])/d[2]
540  vectTranslate = geoalgo.Vector(dirx,diry,dirz)
541  # 1) generate line-segment "behind" half-line
542  # use unit-vector nature of Dir() to go back slightly more
543  dist = random()
544  seg = geoalgo.LineSegment(p1+vectTranslate-d*dist,p1+vectTranslate-d*dist-d)
545  #seg = geoalgo.LineSegment(p1+vectTranslate,p1+vectTranslate-d)
546  # distance should be dist
547  # closest points should be seg.Start() and seg1.Start()
548  answer = vectTranslate*vectTranslate+dist*dist*d.SqLength()
549  tim = time()
550  a1 = dAlgo.SqDist(seg1,seg)
551  tim = time() - tim
552  sqdistT1 += tim
553  L1 = geoalgo.Vector()
554  L2 = geoalgo.Vector()
555  a2 = dAlgo.SqDist(seg1,seg,L1,L2)
556  if not (abs(answer-a1) < _epsilon): success1 = 0
557  if not (abs(answer-a2) < _epsilon): success1 = 0
558  for x in xrange(3):
559  if not (L1[x]-seg1.Start()[x] < _epsilon) : success1 = 0
560  if not (L2[x]-seg.Start()[x] < _epsilon) : success1 = 0
561  if (success1 == 1) : totSuccess1 += 1
562  # 2) generate line segment away but rotated
563  # rotate in a direction perpendicular to both line direction and translation direction
564  vectRotate = vectTranslate.Cross(d)
565  # choose pivot point
566  pivot = p1+vectTranslate-d*dist
567  seg = geoalgo.LineSegment(pivot-vectRotate*random(),pivot+vectRotate*random())
568  # distance should be pivot distance to half-line start
569  # closest points should be pivot and line start point
570  answer = vectTranslate*vectTranslate+dist*dist*d.SqLength()
571  tim = time()
572  a1 = dAlgo.SqDist(seg1,seg,L1,L2)
573  sqdistT2 += (time()-tim)
574  a2 = dAlgo.SqDist(seg,seg1)
575  if not (abs(answer-a1) < _epsilon): success2 = 0
576  if not (abs(answer-a2) < _epsilon): success2 = 0
577  for x in xrange(3):
578  if not (L1[x]-seg1.Start()[x] < _epsilon) : success2 = 0
579  if not (L2[x]-pivot[x] < _epsilon) : success2 = 0
580  if (success2 == 1) : totSuccess2 += 1
581  # 3) generate line segment next to but tilted outwards
582  distin = length*random() # ensures that we do not pass the segment
583  seg = geoalgo.LineSegment(p1+vectTranslate+d*distin,p1+vectTranslate+d*distin+vectTranslate*random())
584  # distance should be vectTranslate
585  # closest point should be segment start and point on line "d*distin" away from line-start
586  answer = vectTranslate*vectTranslate
587  tim = time()
588  a1 = dAlgo.SqDist(seg1,seg,L1,L2)
589  sqdistT3 += (time()-tim)
590  a2 = dAlgo.SqDist(seg,seg1)
591  if not (abs(answer-a1) < _epsilon): success3 = 0
592  if not (abs(answer-a2) < _epsilon): success3 = 0
593  ptLine = seg1.Start()+d*distin
594  for x in xrange(3):
595  if not (L1[x]-ptLine[x] < _epsilon) : success3 = 0
596  if not (L2[x]-seg.Start()[x] < _epsilon) : success3 = 0
597  if (success3 == 1) : totSuccess3 += 1
598  # 4) generate line segment next to line but rotated
599  pivot = p1+vectTranslate+d*distin
600  seg = geoalgo.LineSegment(pivot-vectRotate*random(),pivot+vectRotate*random())
601  # distance should be vectTranslate
602  # closest point should be pivot on segment and d*distin away from start for line
603  answer = vectTranslate*vectTranslate
604  tim = time()
605  a1 = dAlgo.SqDist(seg1,seg,L1,L2)
606  sqdistT4 += (time()-tim)
607  a2 = dAlgo.SqDist(seg,seg1)
608  if not (abs(answer-a1) < _epsilon): success4 = 0
609  if not (abs(answer-a2) < _epsilon): success4 = 0
610  ptLine = seg1.Start()+d*distin
611  for x in xrange(3):
612  if not (L1[x]-ptLine[x] < _epsilon) : success4 = 0
613  if not (L2[x]-pivot[x] < _epsilon) : success4 = 0
614  if (success4 == 1) : totSuccess4 += 1
615 
616  if ( float(totSuccess1)/tests < 1):
617  info(NO + "Success: {0}%".format(100*float(totSuccess1)/tests) + ENDC)
618  else:
619  info(OK + "Success: {0}%".format(100*float(totSuccess1)/tests) + ENDC)
620  info("Time for SqDist (Case 1) : {0:.3f} us".format(1E6*sqdistT1/tests))
621  if ( float(totSuccess2)/tests < 1):
622  info(NO + "Success: {0}%".format(100*float(totSuccess2)/tests) + ENDC)
623  else:
624  info(OK + "Success: {0}%".format(100*float(totSuccess2)/tests) + ENDC)
625  info("Time for SqDist (Case 2) : {0:.3f} us".format(1E6*sqdistT2/tests))
626  if ( float(totSuccess3)/tests < 1):
627  info(NO + "Success: {0}%".format(100*float(totSuccess3)/tests) + ENDC)
628  else:
629  info(OK + "Success: {0}%".format(100*float(totSuccess3)/tests) + ENDC)
630  info("Time for SqDist (Case 3) : {0:.3f} us".format(1E6*sqdistT3/tests))
631  if ( float(totSuccess4)/tests < 1):
632  info(NO + "Success: {0}%".format(100*float(totSuccess4)/float(tests)) + ENDC)
633  else:
634  info(OK + "Success: {0}%".format(100*float(totSuccess4)/float(tests)) + ENDC)
635  info("Time for SqDist (Case 4) : {0:.3f} us".format(1E6*sqdistT4/tests))
636 
637 
638  # test Distance between Point and AABox
639  debug('Testing Point and AABox Distance/Closest Point')
640  success1 = 1
641  totSuccess1 = 0.
642  sqdistT1 = 0
643  closestT1 = 0
644  success2 = 1
645  totSuccess2 = 0.
646  sqdistT2 = 0
647  closestT2 = 0
648  success3 = 1
649  totSuccess3 = 0.
650  sqdistT3 = 0
651  closestT3 = 0
652  success4 = 1
653  totSuccess4 = 0.
654  sqdistT4 = 0
655  closestT4 = 0
656  for y in xrange(tests):
657  success1 = 1
658  success2 = 1
659  success3 = 1
660  success4 = 1
661  # create a simple cubic box from 0,0,0 to 1,1,1
662  b = geoalgo.AABox(0,0,0,1,1,1)
663  # various cases to test:
664  # case 1) Point inside box
665  # case 2) Point outside box
666  # 1) POINT INSIDE BOX
667  p = geoalgo.Vector(random(),random(),random())
668  # if point not recognized as inside -> error!
669  if not ( b.Contain(p) ) : success1 = 0
670  dTop = 1-p[2]
671  dBot = p[2]
672  dLeft = p[0]
673  dRight = 1-p[0]
674  dFront = p[1]
675  dBack = 1-p[1]
676  # find smallest of these
677  dMin = dTop
678  if ( dBot < dMin ) : dMin = dBot
679  if ( dLeft < dMin ) : dMin = dLeft
680  if ( dRight < dMin ) : dMin = dRight
681  if ( dFront < dMin ) : dMin = dFront
682  if ( dBack < dMin ) : dMin = dBack
683  answer = dMin*dMin
684  tim = time()
685  a1 = dAlgo.SqDist(p,b)
686  sqdistT1 += (time()-tim)
687  a2 = dAlgo.SqDist(b,p)
688  # closest point
689  tim = time()
690  pt1 = dAlgo.ClosestPt(b,p)
691  closestT1 += (time()-tim)
692  pt2 = dAlgo.ClosestPt(p,b)
693  if not (abs(answer-a1) < _epsilon): success1 = 0
694  if not (abs(answer-a2) < _epsilon): success1 = 0
695  for x in xrange(3):
696  if not (pt1[x]-p[x] < _epsilon) : success1 = 0
697  if not (pt2[x]-p[x] < _epsilon) : success1 = 0
698  if (success1 == 1) : totSuccess1 += 1
699  # 2) POINT OUT OF BOX
700  # pick a side that is should exit on at random
701  pick = random()
702  side = 0
703  if ( (pick > 0.33) and (pick < 0.67) ) : side = 1
704  if ( pick > 0.67 ) : side = 2
705  # pick a direction to overflow the box by (-1 = back or +1 = front)
706  direction = 1
707  if ( random() < 0.5 ) : direction = -1
708  if ( side == 0 ) :
709  p = geoalgo.Vector(random()+direction,random(),random())
710  elif ( side == 1 ) :
711  p = geoalgo.Vector(random(),random()+direction,random())
712  else :
713  p = geoalgo.Vector(random(),random(),random()+direction)
714  # if point not recognized as outside -> error!
715  if ( b.Contain(p) ) : success1 = 0
716  # go through cases and find min distance
717  if ( (side == 0) and (direction == 1) ):
718  #overflow on +x direction
719  dMin = p[0]-1
720  pMin = geoalgo.Vector(1,p[1],p[2])
721  if ( (side == 0) and (direction == -1) ):
722  #overflow on -x direction
723  dMin = -p[0]
724  pMin = geoalgo.Vector(0,p[1],p[2])
725  if ( (side == 1) and (direction == 1) ):
726  #overflow on +y direction
727  dMin = p[1]-1
728  pMin = geoalgo.Vector(p[0],1,p[2])
729  if ( (side == 1) and (direction == -1) ):
730  #overflow on -y direction
731  dMin = -p[1]
732  pMin = geoalgo.Vector(p[0],0,p[2])
733  if ( (side == 2) and (direction == 1) ):
734  #overflow on +z direction
735  dMin = p[2]-1
736  pMin = geoalgo.Vector(p[0],p[1],1)
737  if ( (side == 2) and (direction == -1) ):
738  #overflow on -z direction
739  dMin = -p[2]
740  pMin = geoalgo.Vector(p[0],p[1],0)
741  answer = dMin*dMin
742  tim = time()
743  a1 = dAlgo.SqDist(p,b)
744  sqdistT2 += (time()-tim)
745  a2 = dAlgo.SqDist(b,p)
746  # closest point
747  tim = time()
748  pt1 = dAlgo.ClosestPt(b,p)
749  closestT2 += (time()-tim)
750  pt2 = dAlgo.ClosestPt(p,b)
751  #info("Point: [{0}, {1}, {2}]".format(p[0],p[1],p[2]))
752  #info("Expected: {0}. Got: {1}".format(answer,a1))
753  #info("Expected: {0}. Got: {1}".format(answer,a2))
754  if not (abs(answer-a1) < _epsilon): success2 = 0
755  if not (abs(answer-a2) < _epsilon): success2 = 0
756  for x in xrange(3):
757  #info("\tExpected: {0}. Got: {1}".format(p[x],pt1[x]))
758  #info("\tExpected: {0}. Got: {1}".format(p[x],pt2[x]))
759  if not (pt1[x]-pMin[x] < _epsilon) : success2 = 0
760  if not (pt2[x]-pMin[x] < _epsilon) : success2 = 0
761  #info("Success: {0}".format(success2))
762  if (success2 == 1) : totSuccess2 += 1
763 
764  if ( float(totSuccess1)/tests < 1):
765  info(NO + "Success: {0}%".format(100*float(totSuccess1)/tests) + ENDC)
766  else:
767  info(OK + "Success: {0}%".format(100*float(totSuccess1)/tests) + ENDC)
768  info("Time for SqDist (Case 1) : {0:.3f} us".format(1E6*sqdistT1/tests))
769  info("Time for ClosestPt (Case 1) : {0:.3f} us".format(1E6*closestT1/tests))
770  if ( float(totSuccess2)/tests < 1):
771  info(NO + "Success: {0}%".format(100*float(totSuccess2)/tests) + ENDC)
772  else:
773  info(OK + "Success: {0}%".format(100*float(totSuccess2)/tests) + ENDC)
774  info("Time for SqDist (Case 2) : {0:.3f} us".format(1E6*sqdistT2/tests))
775  info("Time for ClosestPt (Case 2) : {0:.3f} us".format(1E6*closestT2/tests))
776 
777  except Exception:
778  error('geoalgo::DistanceAlgo unit test failed.')
779  print traceback.format_exception(*sys.exc_info())[2]
780  return 1
781 
782  info('geoalgo::DistanceAlgo unit test complete.')
783  return 0
784 
tca::Length
float Length(const PFPStruct &pfp)
Definition: PFPUtils.cxx:3303
geoalgo::GeoAlgo
Algorithm to compute various geometrical relation among geometrical objects. In particular functions ...
Definition: GeoAlgo.h:43
reco_momentum_tuples.float
float
Definition: reco_momentum_tuples.py:12
test_distance.test_dAlgo
def test_dAlgo()
Definition: test_distance.py:17
geoalgo::LineSegment
Representation of a simple 3D line segment Defines a finite 3D straight line by having the start and ...
Definition: GeoLineSegment.h:25
format
static bool format(QChar::Decomposition tag, QString &str, int index, int len)
Definition: qstring.cpp:11496
geoalgo::AABox
Representation of a 3D rectangular box which sides are aligned w/ coordinate axis. A representation of an Axis-Aligned-Boundary-Box, a simple & popular representation of 3D boundary box for collision detection. The concept was taken from the reference, Real-Time-Collision-Detection (RTCD), and in particular Ch. 4.2 (page 77): .
Definition: GeoAABox.h:34
error
error
Definition: include.cc:26
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
generate_datataset.info
info
Definition: generate_datataset.py:142
geoalgo::Line
Representation of a 3D infinite line. Defines an infinite 3D line by having 2 points which completely...
Definition: GeoLine.h:27
python.larbatch_posix.debug
string debug
Definition: larbatch_posix.py:142
geoalgo::HalfLine
Representation of a 3D semi-infinite line. Defines a semi-infinite 3D line by having a start point (P...
Definition: GeoHalfLine.h:30

Variable Documentation

int test_distance._epsilon = 1
private

Definition at line 9 of file test_distance.py.

string test_distance.BLUE = '\033[94m'

Definition at line 14 of file test_distance.py.

string test_distance.ENDC = '\033[0m'

Definition at line 15 of file test_distance.py.

string test_distance.NO = '\033[91m'

Definition at line 13 of file test_distance.py.

string test_distance.OK = '\033[92m'

Definition at line 12 of file test_distance.py.

Generated by   doxygen 1.8.11