If you have trouble viewing this, try the pdf of this post. You can download the documents in this post.

Computational geometry is an interesting and important topic for imaging in general and x-ray imaging in particular. In this post, I describe the basic formulas for perhaps the most fundamental geometric object—a straight line in three dimensions. The object can also be used in 2D to represent a line in an image.

I used the results described here in the code for my previous posts (here,here, and here) on a CT simulator. In particular, they used a line object. In this post I will derive the basic formulas for an object that encapsulates a three dimensional line. In following posts, I will describe C++ and Matlab implementations and show some applications.

Standard books[1] specify a line by any two points on it. I prefer a more structured definition, which has substantial advantages in implementation of computations as you will see in the examples. I specify a line by two vectors: a vector from the origin perpendicular to the line and a unit length vector direction vector **s**̂ along the line (see Fig. 1↓).

The most basic operation with a line is to construct it from two points. The details of the construction are shown in Fig. 1↓. The direction vector is
**s**̂ = (_{1} − _{2})/(|_{1} − _{2}|).
Note that the direction of **s**̂ is ambiguous depending on the ordering of the two points. But this is also the case with alternate definitions and the user needs to take this into account. The end of the normal vector is on the line so _{1} + *t*_{0}**s**̂. Since = is perpendicular to **s**̂
**s**̂ = 0 = _{1}⋅**s**̂ + *t*_{0}.
Solving, we see that ⋅*t*_{0} = − _{1}⋅**s**̂. With this value we can compute the line vector as _{1} + *t*_{0}**s**̂.

I specify a plane as two parameters: a unit vector **p**̂ perpendicular to the plane and the distance of the perpendicular to the plane *p*. I use two parameters because just using the perpendicular vector cannot represent a plane passing through the origin. I also have a convention that the *p* parameter is positive. If the computation results in a negative number, we can use its absolute value and negate**p**̂.

The problem is to compute the parameters of the line and **s**̂ from the two planes, **p**̂_{1}, *p*_{1} and **p**̂_{2}, *p*_{2} . First, since the line is in both planes, it is perpendicular to both **p**̂ vectors so it is parallel to their cross product
**s**̂ = (**p**̂_{1} × **p**̂_{2})/(|**p**̂_{1} × **p**̂_{2}|).

Also, is in the plane spanned by the **p**̂ vectors
*k*_{1}**p**̂_{1} + *k*_{2}**p**̂_{2}.
=

The end of the vector is in both planes so
**p**̂_{1} =
⋅
*p*_{1} =
*k*_{1} + *k*_{2}*p*_{12}
**p**̂_{2} =
⋅
*p*_{2} =
*k*_{1}*p*_{12} + *k*_{2}
where *p*_{12} = **p**̂_{1}⋅**p**̂_{2}. Solving these two simultaneous equations for *k*_{1} and *k*_{2}
*k*_{1} = (*p*_{1} − *p*_{2}*p*_{12})/(1 − *p*^{2}_{12})
*k*_{2} = ( − *p*_{1}*p*_{12} + *p*_{2})/(1 − *p*^{2}_{12})
Notice that if **p**̂_{1} = **p**̂_{2}, the planes are parallel, *p*_{12} = 1, and the equations break down.

Unlike the 2D case, lines in three dimensions usually do not intersect. Even if they do, the computation of the intersection point is subject to round off errors. I get around this by defining the intersection as the midpoint of the shortest line segment between the two lines as in Fig. 2↓. This is still undefined if the lines are parallel but that is to be expected. The general formula for points on the lines is _{i} = _{i} + *t*_{i}**s**̂_{i} *i* = 1, 2. The line segment between points on the lines is _{1} − ** **_{2} = + *t*_{1}**s**̂_{1} − *t*_{2}**s**̂_{2} where = _{1} − _{2}. The minimum length segment is perpendicular to both lines =
**s**̂_{1} =
⋅
0 =
**s**̂_{1} + *t*_{1} − *t*_{2}*s*_{12}
⋅
**s**̂_{2} =
⋅
0 =
**s**̂_{2} + *t*_{1}*s*_{12} − *t*_{2}
⋅
where *s*_{12} = **s**̂_{1}⋅**s**̂_{2}. Solving the equations
*t*_{1, min} = ( − *d*_{1} + *d*_{2}*s*_{12})/(1 − *s*^{2}_{12})
*t*_{2, min} = ( − *d*_{1}*s*_{12} + *d*_{2})/(1 − *s*^{2}_{12})
where *d*_{i} = ⋅**s**̂_{i} *i* = 1, 2. Substituting these in the general formula for points on a line gives the points on each line closest to the other, _{1, min}, _{2, min}. The “intersection” point is then
_{intersect} = 0.5(** **_{1, min} + _{2, min}).
The user can check the length of the line segment [_{1, min}, _{2, min}] to decide whether this is an intersection. The formula gives the correct point for the intersection in two dimensions.

A common operation is to compute the perpendicular distance of a point from a line. Suppose the point is at . A point on the line is _{L} = + *t***s**̂. The vector joining the original point and a point on the line is − _{L}. This vector is perpendicular to the line if =
**s**̂ = 0 = ⋅**s**̂ − *t*_{perp}

Solving, *t*_{perp} = **s**̂ and the distance is
*d* = | − − *t*_{perp}**s**̂|.

From elementary geometry a plane can be defined from a line and a point outside the line. The line segment from the point to the line vector, − , is in the plane as is the line **s**̂ vector. Therefore _{perp} = ( − ) × **s**̂ is perpendicular to the plane and the plane normal vector is
**p**̂ = (_{perp})/(|_{perp}|)
The line vector is in the plane so the plane distance parameter is *p* = ⋅**p**̂.

Another common operation is to find the intersection point of a line and a plane. As usual, a point on the line is _{L} = + *t***s**̂. At the intersection, _{L}⋅**p**̂ = *p* so

( + *t*_{intersect}**s**̂)⋅**p**̂ = *p*

and
*t*_{intersect} = (*p* − ⋅**p**̂)/(**s**̂⋅**p**̂).

Last edited Aug. 17, 2011

© 2011 by Aprend Technology and Robert E. Alvarez

Linking is allowed but reposting or mirroring is expressly forbidden.

[1] : *Geometric Tools for Computer Graphics*. Morgan Kaufmann, 2002.