// Calculates the Vertex of two straight lines define by the vectors base and dir
//
// g: x1 = base1 + l * dir1
// h: x2 = base2 + m * dir2 , where l,m are real numbers
//
// 1. are g and h
// parallel / identical, i.e. are dir1 and dir2 linear dependent?
//
// /-
// |
// | = 0 linear dependent, no unique solution, returning dummy
// => cross product : dir1 x dir2 = -|
// | != 0 linear independent
// |
// \-
//
// 2. are g and h
// skew or do they have a crossing point, i.e are dir1, dir2 and (base1 - base2) linear dependent ?
//
// /-
// |
// | = 0 linear dependent
// | g and h are intersecting
// | calculating vertex as point of intersection
// |
// => determinant: det[ dir1, dir2, base1-base2]= -|
// | != 0 linear independent
// | g and h are skew
// | calulating vertex as point of closest approach
// |
// \-
//
// 3.
// (a) calculating intersection point
// NOTE: [a,b,c] is a MATRIX with a, b, c as column vectors
//
// g ^ h:
// x1 = x2
// --> base1 + l * dir2 = base2 + m * dir2
// --> base1 - base2 = m * dir2 - l * dir1
//
// --> base1 - base2 = [dir2, -dir1] * (m)
// (l)
//
// overdetermined set of linear equation: but we know that a cross point must exist.
// so just taking the first two (x,y) components for solving this set.
//
// solving with CRAMER's RULE
//
// m = det[(base1- base2), -dir1] / det[dir2, -dir1]
// l = det[ dir2, (base1- base2)] / det[dir2, -dir1] (for completeness, but not used)
//
// vertex = base2 + m * dir2
//
// (b) calculating point of closest approach
//
// from the equations of the straight lines of g and h
// g: x1 = base1 + l * dir1
// h: x2 = base2 + m * dir2
//
// you can construct the following planes:
//
// E1: e1 = base1 + a * dir1 + b * (dir1 x dir2)
// E2: e2 = base2 + s * dir2 + t * (dir1 x dir2)
//
// now the intersection point of E1 with g2 = P1
// and the intersection point of E2 with g1 = P2
//
// form the base points of the perpendicular to both straight lines.
//
// The point of closes approach is the middle point between P1 and P2:
//
// vertex = (p2 - p1)/2
//
// E1 ^ g2:
//
// e1 = x2
// --> base1 + a * dir1 + b * (dir1 x dir2) = base2 + m * dir2
// --> base1 - base2 = m * dir2 - a * dir1 - b * (dir1 x dir2)
// (m)
// --> [ dir2, -dir1, -(dir1 x dir2)] (a) = base1 - base2
// (b)
//
// again using CRAMER's RULE you can find the solution for m (a,b, not used)
//
// using the rules for converting determinants:
//
// D12 = det [dir2, -dir1, -(dir1 x dir2)]
// = det [dir2, dir1, (dir1 x dir2)]
//
// Dm = det [base1 - base2, -dir1, -(dir1 x dir2)]
// = det [base1 - base2, dir1, (dir1 x dir2)]
//
// m = Dm/D12
//
// P1: p1 = x2(m)
// = base2 + Dm/D12 * dir2
//
// E2 ^ g1:
//
// e2 = x1
// --> base2 + s * dir2 + t * (dir1 x dir2) = base1 + l * dir1
// --> base2 - base1 = l * dir1 - s * dir2 - t * (dir1 x dir2)
// (m)
// --> [ dir1, -dir2, -(dir1 x dir2)] (a) = base2 - base1
// (b)
//
// again using CRAMER's RULE you can find the solution for m (a,b, not used)
//
// using the rules for converting determinants:
//
// D21 = det [dir1, -dir2, -(dir1 x dir2)]
// = det [dir1, dir2, (dir1 x dir2)]
// = -det [dir2, dir1, (dir1 x dir2)]
// = -D12
//
// Dl = det [base2 - base1, -dir2, -(dir1 x dir2)]
// = det [base2 - base1, dir1, (dir1 x dir2)]
// = -det [base1 - base2, dir1, (dir1 x dir2)]
//
// l = Dl/D21
// = - Dl/D12
//
// P2: p2 = x1(m)
// = base1 - Dl/D12 * dir1
//
//
// vertex = 1/2 * (p2 - p1)
// = 1/2 * (base1 - base2) + 1/2 * 1/D12 ( Dm * dir2 - Dl * dir1)
//
//
...Back to
Vertex Finding
--
PeterZumbruch - 23 Feb 2004