Calculating the closest point to a set of lines

Questions
,
1

Hi,

as titled, I am trying to calculate the closest point to a set of lines in 3D

I saw here that this guy pointed out some Mathematica code.

I put it into the Mathematica Kernel and I got a very long result.. so I decided to write an algorithm by myself.

Now I would like, instead, given four known lines (each of them defined by two points), to get the closest point..

The lines are the following:

Line[] lines = new Line[]{
new Line(new Vec3(77.007f, 32.46f, 177.79f), new Vec3(139.398f, -309.133f, -116.29f)),
//
new Line(new Vec3(-79.504f, 31.969f, -179.158f), new Vec3(31.886f, -367.752f, 99.796f)),
//
new Line(new Vec3(-77.401f, 32.855f, 174.843f), new Vec3(25.944f, -350.134f, -129.523f)),
//
new Line(new Vec3(80.781f, 32.154f, -176.942f), new Vec3(-149.572f, -326.383f, 84.57f))
};

But I have no knowledge of Mathematica, could someone help me?

2

Hi Giuseppe,

… so I decided to write an algorithm by myself.

Excellent idea!

As a start you could have written at least the definition of the lines in WL … Well, anyway, here my solution, it is pretty straightforward:

ClearAll["Global`*"]
p = {p1, p2, p3};
g = {g1, g2, g3}/Sqrt[g1^2 + g2^2 + g3^2];
(* function to calculate distance between line and point: *)    
lineDist[line : {g1_, g2_, g3_}, point : {p1_, p2_, p3_}] = 
  FullSimplify[Norm[p\[Cross]g], Thread[{p1, p2, p3, g1, g2, g3} \[Element] Reals]];

(* your four lines *)    
v = {{{77.007, 32.46, 177.79}, {139.398, -309.133, -116.29}},
   {{-79.504, 31.969, -179.158}, {31.886, -367.752, 99.796}},
   {{-77.401, 32.855, 174.843}, {25.944, -350.134, -129.523}},
   {{80.781, 32.154, -176.942}, {-149.572, -326.383, 84.57}}};

(* search for the point, where the max-distance is minimized *)
sol = Minimize[Max[lineDist[#1, p - #2] & @@@ ({#2 - #1, #1} & @@@ v)], p]

Graphics3D[{{Red, Line[v]}, Sphere[p /. sol[[2]], sol[[1]]]}, SphericalRegion -> True]

The result is:

pDistOf4Lines.gif

Regards Henrik

3

Hi Henrik,

thanks for your code, one guy gave me this function:

void MathUtil::CalcBeamsClosestPoint(pfSeg* beams, int numBeams, pfVec3& p)
{
    float m00 = 0.0f;
    float m01 = 0.0f;
    float m02 = 0.0f;
    float m11 = 0.0f;
    float m12 = 0.0f;
    float m22 = 0.0f;
    float v0 = 0.0f;
    float v1 = 0.0f;
    float v2 = 0.0f;

    for (int i = 0; i < numBeams; i++) {
	    pfSeg& beam = beams[i];
	    pfVec3& dir = beam.dir;
	    pfVec3& pos = beam.pos;
	    const float ax = dir[0];
	    const float ay = dir[1];
	    const float az = dir[2];
	    const float bx = pos[0];
	    const float by = pos[1];
	    const float bz = pos[2];

	    m00 += ay * ay + az * az;
	    m11 += az * az + ax * ax;
	    m22 += ax * ax + ay * ay;
	    m01 += -ax * ay;
	    m02 += -az * ax;
	    m12 += -ay * az;
	    v0 += az * (az * bx - ax * bz) - ay * (ax * by - ay * bx);
	    v1 += ax * (ax * by - ay * bx) - az * (ay * bz - az * by);
	    v2 += ay * (ay * bz - az * by) - ax * (az * bx - ax * bz);
    }
    matrix33f m(vector3f(m00, m01, m02), vector3f(m01, m11, m12), vector3f(m02, m12, m22));
    float det = m.det();
    assert(det != 0.0f);
    vector3f pt = inverse(m) * vector3f(v0, v1, v2);
    p.set(pt.x, pt.y, pt.z);
}

it works with simple lines such as

       new Line(new Vec3(0, 0, 0), new Vec3(1, 1, 1)),

       new Line(new Vec3(2, 0, 0), new Vec3(1, 1, 1)),

but not with the ones I gave you previously (I obtain a different result).. I don’t know why, I can’t get anyway the theory behind his code

By the way, the formula is this one at the end http://math.stackexchange.com/a/61721/64943

I tried to wrote the algorithm, but I didn’t make it

I should come to a form such as:

A x = b

but I don’t have idea how to deal with sums and dot products..