算法:具有两个约束的三角形，每个角都在给定线上

Question

前段时间我在 math.stackexchange 上问了一个问题并得到了答案。我很难从该答案中得出算法，因为我的背景是设计，希望你们中的一些人能帮助我。

带视觉草图的原始问题和可能的答案在这里: https://math.stackexchange.com/questions/667432/triangle-with-two-constraints-each-corner-on-a-given-line

问题是:给定 3 条 3 维线(a、b 和 c)，它们在公共(public)点 S 和给定的点 B 上重合，我正在寻找 a 上的点 A 和点 C 上的点c 其中AB和BC的长度相同，角ABC为90度。

我将不得不用命令式语言实现这个算法，任何 C++、Java、命令式伪代码或类似代码都可以。

此外，同样欢迎使用不同的方法来解决这个问题。另外:感谢任何提示，如果完整的解决方案确实太耗时了!

Answer 1

两个关键公式是

(我已经在数学堆栈交换站点回复了公式的推导)

将第一个代入第二个最终得到一个四次方程，用封闭形式求解这个方程很烦人。因此，我在 Python 中使用了一个简单的数值求解器:

# function to solve (we look for t such that f(t)=0)
def f(t):
    s = (t*cB - B2) / (t*ac - aB)
    return s*s - 2*s*aB - t*t + 2*t*cB

# given f and an interval to search generates all solutions in the range
def solutions(f, x0, x1, n=100, eps=1E-10):
    X = [x0 + i*(x1 - x0)/(n - 1) for i in xrange(n)]
    Y = map(f, X)
    for i in xrange(n-1):
        if (Y[i]<0 and Y[i+1]>=0 or Y[i+1]<0 and Y[i]>=0):
            xa, xb = X[i], X[i+1]
            ya, yb = Y[i], Y[i+1]
            if (xb - xa) < eps:
                # Linear interpolation
                # 0 = ya + (x - xa)*(yb - ya)/(xb - xa)
                yield xa - ya * (xb - xa) / (yb - ya)
            else:
                for x in solutions(f, xa, xb, n, eps):
                    yield x

搜索算法在区间内对函数进行采样，当它发现两个相邻样本穿过 f=0 线时，会在这两个样本之间递归地重复搜索(除非区间大小低于指定限制，用以下函数逼近函数一条线并计算这种情况下的交叉点)。

我已经测试了生成随机问题并用 解决它们的算法

from random import random as rnd

for test in xrange(1000):
    a = normalize((rnd()-0.5, rnd()-0.5, rnd()-0.5))
    b = normalize((rnd()-0.5, rnd()-0.5, rnd()-0.5))
    c = normalize((rnd()-0.5, rnd()-0.5, rnd()-0.5))

    L = rnd() * 100
    B = tuple(x*L for x in b)

    aB = dot(a, B)
    cB = dot(c, B)
    B2 = dot(B, B)
    ac = dot(a, c)

sols = list(solutions(f, -1000., 1000.))

在某些情况下，解决方案是 0、1、2、3 或 4。例如问题

a = (-0.5900900304960981, 0.4717596600172049, 0.6551614908475357)
b = (-0.9831451620384042, -0.10306322574446096, 0.15100848274062748)
c = (-0.6250439408232388, 0.49902426033920616, -0.6002456660677057)
B = (-33.62793897729328, -3.5252208930692497, 5.165162011403056)

有四种不同的解决方案:

s = 57.3895941365 , t = -16.6969433689
    A = (-33.865027354189415, 27.07409541837935, 37.59945205363035)
    C = (10.436323283003153, -8.332179814593692, 10.022267893763457)
    |A - B| = 44.5910029061
    |C - B| = 44.5910029061
    (A - B)·(C - B) = 1.70530256582e-13

s = 43.619078237 , t = 32.9673082734
    A = (-25.739183207076163, 20.5777215193455, 28.577540327140607)
    C = (-20.606016281518986, 16.45148662649085, -19.78848391300571)
    |A - B| = 34.5155582156
    |C - B| = 34.5155582156
    (A - B)·(C - B) = 1.13686837722e-13

s = -47.5886624358 , t = 83.8222109697
    A = (28.08159526800866, -22.450411211385674, -31.17825902887765)
    C = (-52.39256507303229, 41.82931682916268, -50.313918854788845)
    |A - B| = 74.0747844969
    |C - B| = 74.0747844969
    (A - B)·(C - B) = 4.54747350886e-13

s = 142.883074325 , t = 136.634726869
    A = (-84.31387768560096, 67.4064705656035, 93.61148799140805)
    C = (-85.40270813540043, 68.1840435123674, -82.01440263735996)
    |A - B| = 124.189861967
    |C - B| = 124.189861967
    (A - B)·(C - B) = -9.09494701773e-13