[RESOLVED] Closest distance between two line segments

[Simon666]

Does anyone know a good algorithm for calculating the closest points/distance between two line segments? I use some pretty general code:

http://geomalgorithms.com/a07-_distance.html

which seems widely used and advertized on the web and works in most cases but seems to often fail horribly when line segments are nearly parallel. I've been messing with the SMALL_NUM value for division overflow to no avail. The calculated distance can still vary widely when nearly parallel.

I managed to isolate a specific incident where this happens in my code. The distance between segments P1P2 and Q1Q2 changed abruptly in one timestep from 1.05 mm to 0.90 mm (yarn radius = 1 mm), causing abrupt compression spikes. In reality the distance in the original timestep is definitely also around 0.90 mm but is not calculated as such. I find that the values of s and t (s=0 for P1, 1 for P2, t=0 for Q1, 1 for Q2) for the closest points are originally 0 and 0 (as well as in the previous time steps) and then change abruptly to 0 and around 0.29 in the new time step. What it should be, I still need to check out.

P1 = (0.012711 ,-0.000688 ,-0.001097);
P2 = (0.012895 ,-0.000686 ,-0.001133);
Q1 = (0.012676 ,-0.000689 ,-0.000999);
Q2 = (0.012859 ,-0.000687 ,-0.001034);

P1new = (0.012712, -0.000689, -0.001095);
P2new = (0.012895, -0.000687, -0.001131);
Q1new = (0.012676, -0.000690, -0.000996);
Q2new = (0.012859, -0.000688, -0.001032);

Results when calculating self contact of 15000 line segments in a few tens of fibers, one big mess:



Anyone know of a better algorithm/corrections?

[Simon666]

I have decided to contact Dr. David H. Eberly, which is referenced as the source upon which the author in the first link I provided based himself upon. Dr. Eberly was kind enough to write an extensive reply. If I may quote part of it:

The algorithm Dan Sunday uses is different from mine (that he makes mention about at his web page). They are equivalent theoretically, and both are correct algorithms when using real-valued arithmetic. The problem shows up in an implementation of the algorithm when using floating-point arithmetic. Both algorithms rely on classifying whether two segments are parallel or not parallel. When floating-point arithmetic is used, rounding errors can cause two nonparallel segments to be processed in the parallel clause. Or two parallel segments can be processed in the nonparallel clause. The question is how the misclassification affects the results.

The parallel/nonparallel classification reduces to testing whether a 2x2 determinant is zero/not-zero. I just looked at Dan's code. He uses a floating-point tolerance on the determinant. The tolerance is 1e-8, apparently an arbitrary choice. He also mixes 'float' and 'double' in his dist3D_Segment_to_Segment function. I had tolerances in my Wild Magic source code, but removed them in my GTEngine code.

(...)

I have support for arbitrary precision arithmetic. For SegP/SegQ, the determinant is nonzero. It uses more bits of precision than 'double', (...) If you want, I can send to you a sample source code program that uses the arbitrary precision arithmetic of my library. The program will compute the squared distance between segments using exact arithmetic, switching back to 'float' or 'double' only at the very end to compute the distance (sqrt cannot be performed with arbitrary precision arithmetic).

I already asked for that sample source code and am yet to hear back from him. Meanwhile, I have downloaded the Geometric Tools engine source code from:

http://www.geometrictools.com/Downloads/Downloads.html

And looked at the code for segment to segment distance calculation, adjusted it to fit into my code and now it works fine using double's everywhere. I am still going to look into how to implement this "arbitrary precision arithmetic", I want to ask whether and where this arithmetic is necessary and still have some questions about some issues. I have however currently a temporary solution that seems to be working a lot better already.