这个解决方案本质上是Alex Martelli的解决方案，但我添加了一个Point和LineSegment类以使阅读更容易。我还调整了格式并添加了一些测试。
线段交点是错误的，但似乎对于计算线段距离并不重要。如果您对正确的线段交点测试感兴趣，请看这里：如何检测两条线段是否相交？
#!/usr/bin/env python
"""Calculate the distance between line segments."""
import math
class Point(object):
"""A two dimensional point."""
def __init__(self, x, y):
self.x = float(x)
self.y = float(y)
class LineSegment(object):
"""A line segment in a two dimensional space."""
def __init__(self, p1, p2):
assert isinstance(p1, Point), \
"p1 is not of type Point, but of %r" % type(p1)
assert isinstance(p2, Point), \
"p2 is not of type Point, but of %r" % type(p2)
self.p1 = p1
self.p2 = p2
def segments_distance(segment1, segment2):
"""Calculate the distance between two line segments in the plane.
>>> a = LineSegment(Point(1,0), Point(2,0))
>>> b = LineSegment(Point(0,1), Point(0,2))
>>> "%0.2f" % segments_distance(a, b)
'1.41'
>>> c = LineSegment(Point(0,0), Point(5,5))
>>> d = LineSegment(Point(2,2), Point(4,4))
>>> e = LineSegment(Point(2,2), Point(7,7))
>>> "%0.2f" % segments_distance(c, d)
'0.00'
>>> "%0.2f" % segments_distance(c, e)
'0.00'
"""
if segments_intersect(segment1, segment2):
return 0
# try each of the 4 vertices w/the other segment
distances = []
distances.append(point_segment_distance(segment1.p1, segment2))
distances.append(point_segment_distance(segment1.p2, segment2))
distances.append(point_segment_distance(segment2.p1, segment1))
distances.append(point_segment_distance(segment2.p2, segment1))
return min(distances)
def segments_intersect(segment1, segment2):
"""Check if two line segments in the plane intersect.
>>> segments_intersect(LineSegment(Point(0,0), Point(1,0)), \
LineSegment(Point(0,0), Point(1,0)))
True
"""
dx1 = segment1.p2.x - segment1.p1.x
dy1 = segment1.p2.y - segment1.p2.y
dx2 = segment2.p2.x - segment2.p1.x
dy2 = segment2.p2.y - segment2.p1.y
delta = dx2 * dy1 - dy2 * dx1
if delta == 0: # parallel segments
# TODO: Could be (partially) identical!
return False
s = (dx1 * (segment2.p1.y - segment1.p1.y) +
dy1 * (segment1.p1.x - segment2.p1.x)) / delta
t = (dx2 * (segment1.p1.y - segment2.p1.y) +
dy2 * (segment2.p1.x - segment1.p1.x)) / (-delta)
return (0 <= s <= 1) and (0 <= t <= 1)
def point_segment_distance(point, segment):
"""
>>> a = LineSegment(Point(1,0), Point(2,0))
>>> b = LineSegment(Point(2,0), Point(0,2))
>>> point_segment_distance(Point(0,0), a)
1.0
>>> "%0.2f" % point_segment_distance(Point(0,0), b)
'1.41'
"""
assert isinstance(point, Point), \
"point is not of type Point, but of %r" % type(point)
dx = segment.p2.x - segment.p1.x
dy = segment.p2.y - segment.p1.y
if dx == dy == 0: # the segment's just a point
return math.hypot(point.x - segment.p1.x, point.y - segment.p1.y)
if dx == 0:
if (point.y <= segment.p1.y or point.y <= segment.p2.y) and \
(point.y >= segment.p2.y or point.y >= segment.p2.y):
return abs(point.x - segment.p1.x)
if dy == 0:
if (point.x <= segment.p1.x or point.x <= segment.p2.x) and \
(point.x >= segment.p2.x or point.x >= segment.p2.x):
return abs(point.y - segment.p1.y)
# Calculate the t that minimizes the distance.
t = ((point.x - segment.p1.x) * dx + (point.y - segment.p1.y) * dy) / \
(dx * dx + dy * dy)
# See if this represents one of the segment's
# end points or a point in the middle.
if t < 0:
dx = point.x - segment.p1.x
dy = point.y - segment.p1.y
elif t > 1:
dx = point.x - segment.p2.x
dy = point.y - segment.p2.y
else:
near_x = segment.p1.x + t * dx
near_y = segment.p1.y + t * dy
dx = point.x - near_x
dy = point.y - near_y
return math.hypot(dx, dy)
if __name__ == '__main__':
import doctest
doctest.testmod()
这里是Java的解决方案（使用点检查轻松完成，因此可能不够高效）：
public static double getDistanceBetweenLineSegments(
PointDouble line1Start, PointDouble line1End,
PointDouble line2Start, PointDouble line2End) {
double result = 0;
// If they don't intersect, then work out the distance
if (!isLineIntersectingLine(line1Start, line1End, line2Start, line2End)) {
double p1 = getDistanceBetweenPointAndLine(line1Start, line2Start, line2End);
double p2 = getDistanceBetweenPointAndLine(line1End, line2Start, line2End);
double p3 = getDistanceBetweenPointAndLine(line2Start, line1Start, line1End);
double p4 = getDistanceBetweenPointAndLine(line2End, line1Start, line1End);
result = MathSafe.min(p1, MathSafe.min(p2, MathSafe.min(p3, p4)));
}
return result;
}
还有你需要的所有其他代码：
public class PointDouble {
private double x;
private double y;
public PointDouble(double x, double y) {
this.x = x;
this.y = y;
}
public double getX() {
return x;
}
public double getY() {
return y;
}
}
private static int relativeCCW(
double x1, double y1,
double x2, double y2,
double px, double py) {
x2 -= x1;
y2 -= y1;
px -= x1;
py -= y1;
double ccw = px * y2 - py * x2;
if (ccw == 0.0) {
ccw = px * x2 + py * y2;
if (ccw > 0.0) {
px -= x2;
py -= y2;
ccw = px * x2 + py * y2;
if (ccw < 0.0) {
ccw = 0.0;
}
}
}
return (ccw < 0.0) ? -1 : ((ccw > 0.0) ? 1 : 0);
}
public static boolean isLineIntersectingLine(
PointDouble line1Start, PointDouble line1End,
PointDouble line2Start, PointDouble line2End) {
return (
(relativeCCW(line1Start.getX(), line1Start.getY(), line1End.getX(), line1End.getY(), line2Start.getX(), line2Start.getY()) *
relativeCCW(line1Start.getX(), line1Start.getY(), line1End.getX(), line1End.getY(), line2End.getX(), line2End.getY()) <= 0)
&&
(relativeCCW(line2Start.getX(), line2Start.getY(), line2End.getX(), line2End.getY(), line1Start.getX(), line1Start.getY()) *
relativeCCW(line2Start.getX(), line2Start.getY(), line2End.getX(), line2End.getY(), line1End.getX(), line1End.getY()) <= 0));
}
public static double getDistanceBetweenPointAndLine(PointDouble pt, PointDouble linePt1, PointDouble linePt2) {
double lineX = linePt2.getX() - linePt1.getX();
double lineY = linePt2.getY() - linePt1.getY();
double dot = (pt.getX() - linePt1.getX()) * lineX + (pt.getY() - linePt1.getY()) * lineY;
double len_sq = lineX * lineX + lineY * lineY;
double param = -1;
double xx;
double yy;
if (len_sq != 0) {
param = dot / len_sq;
}
if (param < 0) {
xx = linePt1.getX();
yy = linePt1.getY();
}
else if (param > 1) {
xx = linePt2.getX();
yy = linePt2.getY();
}
else {
xx = linePt1.getX() + param * lineX;
yy = linePt1.getY() + param * lineY;
}
return MathSafe.hypot(pt.getX() - xx, pt.getY() - yy);
}
这里有一个使用Perl编写的代码，与Fnord的代码有一些不同：
1e-6值用于处理可能的浮点误差。如果您需要不同的公差，请调整或删除这些值：use strict;
use warnings;
use Math::Vector::Real;
use Math::Matrix;
# True if $a >= $b within $tolerance
sub approx_greater
{
my ($a, $b, $tolerance) = @_;
$tolerance //= 1e-6;
return ($a-$b) > -$tolerance;
}
# True if $a <= $b within $tolerance
sub approx_lesser
{
my ($a, $b, $tolerance) = @_;
$tolerance //= 1e-6;
return ($b-$a) > -$tolerance;
}
# True if $a == $b within $tolerance
sub approx_equal
{
my ($a, $b, $tolerance) = @_;
$tolerance //= 1e-6;
return abs($a-$b) < $tolerance;
}
# Returns shortest line: [ [x1,y1,z1], [x2,y2,z2], distance ].
# If askew lines cross (or nearly-intersect) then xyz1 and xyz2 are undefined
# and only distance is returned.
#
# Thank you to @Fnord: https://dev59.com/WXE85IYBdhLWcg3wViA4#18994296
sub line_intersect
{
# Map @_ as vectors:
my ($a0, $a1, $b0, $b1) = map { V(@{ $_ }) } @_;
my $A = ($a1-$a0);
my $B = ($b1-$b0);
my $magA = abs($A);
my $magB = abs($B);
# If length line segment:
if ($magA == 0 || $magB == 0)
{
return V(undef, undef, 0);
}
my $_A = $A / $magA;
my $_B = $B / $magB;
my $cross = $_A x $_B;
my $denom = $cross->norm2;
# If lines are parallel (denom=0) test if lines overlap.
# If they don't overlap then there is a closest point solution.
# If they do overlap, there are infinite closest positions, but there is a closest distance
#if ($denom == 0)
if (approx_equal($denom, 0))
{
my $d0 = $_A * ($b0-$a0);
my $d1 = $_A * ($b1-$a0);
# Is segment B before A?
#if ($d0 <= 0 && 0 >= $d1)
if (approx_lesser($d0, 0) && approx_greater(0, $d1))
{
if (abs($d0) < abs($d1))
{
return V($a0, $b0, abs($a0-$b0));
}
else
{
return V($a0, $b1, abs($a0-$b1));
}
}
# Is segment B after A?
#elsif ($d0 >= $magA && $magA <= $d1)
elsif (approx_greater($d0, $magA) && approx_lesser($magA, $d1))
{
if (abs($d0) < abs($d1))
{
return V($a1, $b0, abs($a1-$b0));
}
else
{
return V($a1, $b1, abs($a1-$b1));
}
}
else
{
# Segments overlap, return distance between parallel segments
return V(V(), V(), abs((($d0*$_A)+$a0)-$b0));
}
}
else
{
# Lines criss-cross: Calculate the projected closest points
my $t = ($b0 - $a0);
# Math::Matrix won't wirth with Math::Vector::Real
# even though they are blessed arrays,
# so convert them to arrays and back to refs:
my $detA = Math::Matrix->new([ [ @$t ], [ @$_B ], [ @$cross] ])->det;
my $detB = Math::Matrix->new([ [ @$t ], [ @$_A ], [ @$cross] ])->det;
my $t0 = $detA / $denom;
my $t1 = $detB / $denom;
my $pA = $a0 + ($_A * $t0); # Projected closest point on segment A
my $pB = $b0 + ($_B * $t1); # Projected closest point on segment A
if ($t0 < 0)
{
$pA = $a0;
}
elsif ($t0 > $magA)
{
$pA = $a1;
}
if ($t1 < 0)
{
$pB = $b0;
}
elsif ($t1 > $magB)
{
$pB = $b1;
}
# Clamp projection A
if ($t0 < 0 || $t0 > $magA)
{
my $dot = $_B * ($pA-$b0);
if ($dot < 0)
{
$dot = 0;
}
elsif ($dot > $magB)
{
$dot = $magB;
}
$pB = $b0 + ($_B * $dot)
}
# Clamp projection B
if ($t1 < 0 || $t1 > $magB)
{
my $dot = $_A * ($pB-$a0);
if ($dot < 0)
{
$dot = 0;
}
elsif ($dot > $magA)
{
$dot = $magA;
}
$pA = $a0 + ($_A * $dot)
}
return V($pA, $pB, abs($pA-$pB));
}
}
print "example: " . line_intersect(
[13.43, 21.77, 46.81 ], [27.83, 31.74, -26.60 ],
[77.54, 7.53, 6.22 ], [26.99, 12.39, 11.18 ]) . "\n" ;
print "contiguous: " . line_intersect(
[0, 0, 0 ], [ 0, 0, 1 ],
[0, 0, 1 ], [ 0, 0, 2 ],
) . "\n" ;
print "contiguous 90: " . line_intersect(
[0, 0, 0 ], [ 0, 0, 1 ],
[0, 0, 1 ], [ 0, 1, 1 ],
) . "\n" ;
print "colinear separate: " . line_intersect(
[0, 0, 0 ], [ 0, 0, 1 ],
[0, 0, 2 ], [ 0, 0, 3 ],
) . "\n" ;
print "cross: " . line_intersect(
[1, 1, 0 ], [ -1, -1, 0 ],
[-1, 1, 0 ], [ 1, -1, 0 ],
) . "\n" ;
print "cross+z: " . line_intersect(
[1, 1, 0 ], [ -1, -1, 0 ],
[-1, 1, 1 ], [ 1, -1, 1 ],
) . "\n" ;
print "full overlap1: " . line_intersect(
[2, 0, 0 ], [ 5, 0, 0 ],
[3, 0, 0 ], [ 4, 0, 0 ],
) . "\n" ;
print "full overlap2: " . line_intersect(
[3, 0, 0 ], [ 4, 0, 0 ],
[2, 0, 0 ], [ 5, 0, 0 ],
) . "\n" ;
print "partial overlap1: " . line_intersect(
[2, 0, 0 ], [ 5, 0, 0 ],
[3, 0, 0 ], [ 6, 0, 0 ],
) . "\n" ;
print "partial overlap2: " . line_intersect(
[3, 0, 0 ], [ 6, 0, 0 ],
[2, 0, 0 ], [ 5, 0, 0 ],
) . "\n" ;
print "parallel: " . line_intersect(
[3, 0, 0 ], [ 6, 0, 0 ],
[3, 0, 1 ], [ 6, 0, 1 ],
) . "\n" ;
example: {{20.29994361624, 26.5264817954106, 11.7875999397098}, {26.99, 12.39, 11.18}, 15.6513944955904}
contiguous: {{0, 0, 1}, {0, 0, 1}, 0}
contiguous 90: {{0, 0, 1}, {0, 0, 1}, 0}
colinear separate: {{0, 0, 1}, {0, 0, 2}, 1}
cross: {{-2.22044604925031e-16, -2.22044604925031e-16, 0}, {2.22044604925031e-16, -2.22044604925031e-16, 0}, 4.44089209850063e-16}
cross+z: {{-2.22044604925031e-16, -2.22044604925031e-16, 0}, {2.22044604925031e-16, -2.22044604925031e-16, 1}, 1}
full overlap1: {{}, {}, 0}
full overlap2: {{}, {}, 0}
partial overlap1: {{}, {}, 0}
partial overlap2: {{}, {}, 0}
parallel: {{}, {}, 1}
前面的答案看起来太复杂了，所以我决定在Python/Blender中实现一个简短的脚本来计算由两个2D端点标识的两个分段之间的最短距离。这是基于正交投影是最短距离的原理：
我正在寻找一种计算大量三维线条之间最短距离的方法。来自Fnord的代码是一个起点，但仅适用于单个批处理。我将其大部分移植到了多批处理中，如下所示。
import numpy as np
import time
def closest_distance_between_lines(a0, a1, b0, b1):
# equation for a line
number_of_lines = a0.shape[0]
A = a1 - a0
B = b1 - b0
# get the magnitude of each line
magA = np.broadcast_to(np.expand_dims(np.linalg.norm(A, axis=1), axis=1), (number_of_lines, 3))
magB = np.broadcast_to(np.expand_dims(np.linalg.norm(B, axis=1), axis=1), (number_of_lines, 3))
# get the unit-normalized lines
_A = A / magA
_B = B / magB
# get the perpendicular line
cross = np.cross(_A, _B, axis=1)
# normalize the perpendicular line
denom = np.linalg.norm(cross, axis=1)**2
# get the line between the start of the previous two lines
t = (b0 - a0)
stacked_matrices_A = np.stack([t, _B, cross], axis=1)
detA = np.linalg.det(stacked_matrices_A)
stacked_matrices_B = np.stack([t, _A, cross], axis=1)
detB = np.linalg.det(stacked_matrices_B)
t0 = detA/denom
t1 = detB/denom
t0 = np.broadcast_to(np.expand_dims(t0, axis=1), (number_of_lines, 3))
t1 = np.broadcast_to(np.expand_dims(t1, axis=1), (number_of_lines, 3))
# get the points that represent the closest projected distance between the two lines
pA = a0 + (_A * t0)
pB = b0 + (_B * t1)
return pA, pB, np.linalg.norm(pA-pB, axis=1)
number_of_samples = 10000000
a0 = np.random.rand(number_of_samples,3)
a1 = np.random.rand(number_of_samples,3)
b0 = np.random.rand(number_of_samples,3)
b1 = np.random.rand(number_of_samples,3)
start_time = time.time()
point_line_a, point_line_b, distance = closest_distance_between_lines(a0, a1, b0, b1)
end_time = time.time()
print("\nPoints pA on line A that are closest to line B with shape {}:\n{}".format(point_line_a.shape, point_line_a))
print("\nPoints pB on line B that are closest to line A with shape {}:\n{}".format(point_line_b.shape, point_line_b))
print("\nDistances between pA and pB with shape {}:\n{}".format(distance.shape, distance))
print("\n{:.2} seconds to compute closest distances between {:,} lines\n".format(end_time - start_time, number_of_samples))
这是上面的输出：
Points pA on line A that are closest to line B with shape (10000000, 3):
[[0.51075268 0.63261352 0.41417815]
[0.88225225 0.41163515 0.06090485]
[0.27712801 0.99045177 0.58932854]
...
[0.25982217 0.11225041 0.79015618]
[0.58653313 0.47864821 0.11680724]
[0.38297058 0.2251661 1.11088736]]
Points pB on line B that are closest to line A with shape (10000000, 3):
[[0.5084189 0.42417817 0.5847618 ]
[0.83041058 0.38914519 0.02384158]
[0.26068716 0.98567827 0.5010647 ]
...
[0.34356827 0.42162445 0.75820875]
[0.44523571 0.40278146 0.0014156 ]
[0.28498604 0.23670301 1.00712087]]
Distances between pA and pB with shape (10000000,):
[0.26935018 0.0675799 0.08990881 ... 0.32209679 0.19757518 0.14318364]
8.3 seconds to compute closest distances between 10,000,000 lines
import torch
import time
device = torch.device('cuda:0')
def closest_distance_between_lines(a0, a1, b0, b1):
number_of_lines = a0.shape[0]
# equation for a line
A = a1 - a0
B = b1 - b0
# get the magnitude of each line
magA = torch.broadcast_to(torch.unsqueeze(torch.linalg.norm(A, dim=1), dim=1), (number_of_lines, 3))
magB = torch.broadcast_to(torch.unsqueeze(torch.linalg.norm(B, dim=1), dim=1), (number_of_lines, 3))
# get the unit-normalized lines
_A = A / magA
_B = B / magB
# get the perpendicular line
cross = torch.cross(_A, _B, dim=1)
# normalize the perpendicular line
denom = torch.linalg.norm(cross, dim=1)**2
# get the line between the start of the previous two lines
t = (b0 - a0)
stacked_matrices_A = torch.stack([t, _B, cross], dim=1)
detA = torch.linalg.det(stacked_matrices_A)
stacked_matrices_B = torch.stack([t, _A, cross], dim=1)
detB = torch.linalg.det(stacked_matrices_B)
t0 = detA / denom
t1 = detB / denom
t0 = torch.broadcast_to(torch.unsqueeze(t0, dim=1), (number_of_lines, 3))
t1 = torch.broadcast_to(torch.unsqueeze(t1, dim=1), (number_of_lines, 3))
# get the points that represent the closest projected distance between the two lines
pA = a0 + (_A * t0)
pB = b0 + (_B * t1)
return pA, pB, torch.linalg.norm(pA-pB, dim=1)
number_of_samples = 10000000
a0 = torch.rand(number_of_samples,3).to(device=device)
a1 = torch.rand(number_of_samples,3).to(device=device)
b0 = torch.rand(number_of_samples,3).to(device=device)
b1 = torch.rand(number_of_samples,3).to(device=device)
start_time = time.time()
point_line_a, point_line_b, distance = closest_distance_between_lines(a0, a1, b0, b1)
end_time = time.time()
print("\nPoints pA on line A that are closest to line B with shape {}:\n{}".format(point_line_a.shape, point_line_a))
print("\nPoints pB on line B that are closest to line A with shape {}:\n{}".format(point_line_b.shape, point_line_b))
print("\nDistances between pA and pB with shape {}:\n{}".format(distance.shape, distance))
print("\n{:.2} seconds to compute closest distances between {:,} lines\n".format(end_time - start_time, number_of_samples))
在我的GPU上，从NumPy到Torch的运行时间从8.3秒加速到0.13秒。
享受吧。
这是一个基于@TreyReynolds的简洁优雅的Lua实现，并使用的C#实现。
private static Line ShortestBridge(Line a, Line b, bool clamp = true)
{
if (a.Length == 0 || b.Length == 0) throw new ArgumentException("The supplied lines must not have a length of zero!");
var eta = 1e-6;
var r = b.From - a.From;
var u = a.To - a.From;
var v = b.To - b.From;
var ru = Vector3.Dot(r, u);
var rv = Vector3.Dot(r, v);
var uu = Vector3.Dot(u, u);
var uv = Vector3.Dot(u, v);
var vv = Vector3.Dot(v, v);
var det = uu * vv - uv * uv;
float s, t;
if (det < eta * uu * vv) {
s = OptionalClamp01(ru / uu, clamp);
t = 0;
} else {
s = OptionalClamp01((ru * vv - rv * uv) / det, clamp);
t = OptionalClamp01((ru * uv - rv * uu) / det, clamp);
}
var S = OptionalClamp01((t * uv + ru) / uu, clamp);
var T = OptionalClamp01((s * uv - rv) / vv, clamp);
var A = a.From + S * u;
var B = b.From + T * v;
return new Line(A, B);
}
private static float OptionalClamp01(float value, bool clamp)
{
if (!clamp) return value;
if (value >1) return 1;
if (value <0) return 0;
return value;
}
如果你没有使用带有Line类的库，那么这就是最简单的实现方式，可以根据需要进行扩展...
using System.Numerics;
public class Line
{
public Line(Vector3 _From, Vector3 _To)
{
From=_From;
To=_To;
}
public Vector3 From {get;}
public Vector3 To {get;}
public float Length => (To-From).Length();
}
测试或使用如下：
public static void Main()
{
// define line 1
var a = new Vector3(-10,-10, 10);
var b = new Vector3(10,10,10);
var l1 = new Line(a,b);
// define line 2
var c = new Vector3(5,0,3);
var d = new Vector3(6,3,-1);
var l2 = new Line(c,d);
// get shortest bridging line between them
var result = ShortestBridge(l1,l2);
Console.WriteLine($"bridge is from {result.From}, to {result.To}.");
}