Geometrical Primitives

In this section, geometrical primitives including

point,
vector,
line,
plane,
spheres, and
ellipsoids

and their mathematical operations, are introduced.

Point¶

In geometry, a point is an abstract idealization of an exact position, with no size or mass, in space. For this ource, this space is the [Cartesian space]Cartesian coordinate system, thus the position of a point is its coordinates.

Visualization¶

Visualization using computer graphics is a very helpful way to

understand the geometrical relationships, and
validation the results of a numerical solution of an algorithm.

Here, we introduce Visualization Toolkit, or VTK.

The example below depicts a single point as a sphere in a 3D Cartesian coordinate system. The full source code is here.

Step 1: Import External Libraries¶

Working in Python, we need to import all the neessary functions and external libraries first. This includes numpy and vtk filters:

import math
import numpy as np

# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonTransforms import vtkTransform
from vtkmodules.vtkFiltersSources import vtkSphereSource
from vtkmodules.vtkRenderingAnnotation import vtkAxesActor
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkPolyDataMapper,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
    vtkRenderer
)

Step 2: Rendering Pipeline¶

We need to set up a rendering pipeline. In VTK, this includes:

A vtkRenderWindow is a window in a graphical user interface (GUI) where a renderer draw their images,
A vtkRenderer is an object that controls the rendering process of objects, and
A vtkRenderWindowInteractor provides an interaction mechanism better users (i.e. mouse inputs) and the render window

# Visualization Pipeline

# a renderer and render window
ren = vtkRenderer()
renWin = vtkRenderWindow()
renWin.SetWindowName('AISE 4025: point and axes')
renWin.SetSize( 640, 480 )
renWin.AddRenderer( ren )

# an interactor
iren = vtkRenderWindowInteractor()
iren.SetRenderWindow( renWin )

Step 3: Geometrical Primitives¶

A geometrical primitives such as a point and a line is created by

a source, i.e. its mathematical description
a mapper, an object that translates the source to,
an actor, the final object that is to be rendered in the renderer. For geometrical objects, they are typically a set of connected triangles

colors = vtkNamedColors()

# Create a sphere to represent a point
point = np.array([0, 5, 0])         # position of a point

# use a sphere with a finite size to visualiza a point, since a point, by definition, has no size
sphereSource = vtkSphereSource()
sphereSource.SetCenter( point.tolist() )
sphereSource.SetRadius( 0.5 )       # an arbitrary size

# create a mapper
sphereMapper = vtkPolyDataMapper()
sphereMapper.SetInputConnection( sphereSource.GetOutputPort() )

# create an actor
sphereActor = vtkActor()
sphereActor.SetMapper( sphereMapper )
sphereActor.GetProperty().SetColor( colors.GetColor3d('Red') )

We can add a 3D axes to help us orient

# add an 3D Axes
axes = vtkAxesActor()
axes.SetTotalLength(10,10,10)
axes.AxisLabelsOff()

# properties of the axes labels can be set as follows
# this sets the x axis label to red
axes.GetXAxisCaptionActor2D().GetCaptionTextProperty().SetColor(colors.GetColor3d('Red'));
axes.GetYAxisCaptionActor2D().GetCaptionTextProperty().SetColor(colors.GetColor3d('Green'));
axes.GetZAxisCaptionActor2D().GetCaptionTextProperty().SetColor(colors.GetColor3d('Blue'));

Step 4: Rendering¶

Add all actors to the renderer so the renderer knows what to draw.

# Add the actors to the scene
ren.AddActor( sphereActor )
ren.AddActor( axes )
ren.SetBackground( colors.GetColor3d( 'MidnightBlue') ) # the color can be specified as a name
ren.SetBackground( .1, .2, .4 )                         # or as RGB
ren.SetBackground( 1, 1, 1 )                            # this is white
ren.ResetCamera()

renWin.Render()
iren.Start()

In this particular example, the point is loccated at $(0, 5, 0)$ , i.e. on the $y$ -axis. We draw a red sphere at that location, where each of the $x$ -, $y$ -, and $z$ -axis is shown as red, green, blue lines with arrow.

A point on a Cartesian Coordinate System — Figure 1:A point visualized as a red sphere.

Draw onto the render window and start the user interaction.

Vector and Unit Vector¶

While a point specifies the position of an object in the {wiki:Cartesian_coordinate_system}, a vector is used to specify the orientation of an object. As such, a vector is not fixed on the {wiki:Cartesian_coordinate_system} at a given location.

It is often useful to think of the position of a point as a vector, originated from the origin to the point, instead.

For mathmatical convenience, a unit vector, a vector of length 1, is often used to represent the orientation.

A unit vector is a vector with length one:

That is,

length(P)=\sqrt{x^2+y^2+z^2} = 1

(1)

"""
@author: Prof. Elvis C. S. Chen

January 8, 2026
"""

import math
import numpy as np 
import numpy.matlib

p = np.array([ 1, 0, 0 ])      # This is a ROW vector
p = p.reshape(-1, 1)           # this is a COLUMN vector
p = np.array([ [0],[1],[0] ])  # alternative
print("The 3x1 colume vector:\n", p)
print(p.shape)
print("The length, or norm, of this vector is: ", np.linalg.norm(p))       # Length of the vector can be calculated using numpy

The 3x1 colume vector:
 [[0]
 [1]
 [0]]
(3, 1)
The length, or norm, of this vector is:  1.0

Scaling of a Vector¶

A vector can be scaled up (i.e. lengthening) or down (i.e. shortening) by a scalar factor $s$

A scaled vector — Figure 3:A vector $P$ beins scaled as $q$ .

If

$s>1$ , the vector is being lengthened,
$s<1$ , the vector is being shortened,
$s=1$ , nothing happens,
$s=0$ , something is probably wrong,
$s<0$ , the orientation of the vector is changed, in addition to being lengthened/shortened.

# demonstration of scaling of a vector
s = 0.75 # Scaling the vector by a factor of 0.75
q = s * p
print("The original vector p is\n", p, "\nand the scaling factor is:", s, "\n")
print("q is the scaled vector of p: \n", q, "\nwith the norm of the vector being:", np.linalg.norm(q))

The original vector p is
 [[0]
 [1]
 [0]] 
and the scaling factor is: 0.75 

q is the scaled vector of p: 
 [[0.  ]
 [0.75]
 [0.  ]] 
with the norm of the vector being: 0.75

Obviously, scaling of a vector changes its length accordingly.

Normalizing a Vector¶

Normalizing a vector meaning to change the vector length to 1:

A normalized vector — Figure 4:A vector $P$ beins normalized as as $q$ .

In Figure 4, the vector $q$ is the normalized version of $P$ , as the length of $q$ is 1; both $q$ and $P$ are originated from the origin, pointing in the same direction.

p = np.array([ 3, 4, 5 ])      # This is a ROW vector
p = p.reshape(-1, 1)           # this is a COLUMN vector

q = p/np.linalg.norm(p)
print("The original vector p\n", p, "\nhas a length of", np.linalg.norm(p), "\n")
print("The normalizes vesion q\n", q, "\nhas a length of", np.linalg.norm(q), "\n")

The original vector p
 [[3]
 [4]
 [5]] 
has a length of 7.0710678118654755 

The normalizes vesion q
 [[0.42426407]
 [0.56568542]
 [0.70710678]] 
has a length of 0.9999999999999999

Sum of Vectors¶

Just like scalars, vectors can be added and substracted. Adding two vectors means that you contatenate the end of one vector to the origin of another:

In Figure 5, the vector $p_3$ is the sum of $p_1$ and $p_2$ :

p1 = np.array([ 1, 2, 3 ])
p2 = np.array([ 4, 5, 6 ])
p1 = p1.reshape(-1,1)
p2 = p2.reshape(-1,1)

p3 = p1 + p2

print("the sum of vector 1\n", p1, "\nand vector 2\n", p2, "\nis\n", p3)

the sum of vector 1
 [[1]
 [2]
 [3]] 
and vector 2
 [[4]
 [5]
 [6]] 
is
 [[5]
 [7]
 [9]]

For visualization, consult the following tutorial.

Distance between two Points¶

The distance between two points $A$ and $B$ is the length of the vector $\overrightarrow{AB}$ (or $\overrightarrow{BA}$ ):

print(np.linalg.norm(p2-p1))
print(np.linalg.norm(p1-p2))

5.196152422706632
5.196152422706632

Cross Product between two Vectors¶

Cross product between vectors $a$ and $b$ results in a third vector $c$ that is perpendicular to both $u$ and $v$ , with the direction of vector $q$ that obeys the RHR. That is, if the right-hand index finger points in the direction of the vector $a$ , and the right-hand middle finger points in the direction of the vector $b$ , then the thumb points to the direction of $c=a \times b$ . Cross product is denoted by the symbol $\times$ .

Right hand rule — Figure 7:Right-hand rule, image courtesy of Wikimedia Commons here, accessed January 9, 2026.

This directionality implies $a \times b \ne b \times a$ : cross product is not communitive. In fact, the cross product is anticommutative:

a \times b = - b \times a

(2)

Cross product is distributive, however:

a \times (b + c) = a \times b + a \times c

(3)

The cross product between $a=(a_1, a_2, a_3)$ and $b=(b_1,b_2,b_3)$ can be expressed as the formal determinant:

\begin{align} a \times b & = & \begin{vmatrix} i & j & k \\ a_1 & a_2 & a_3\\ b_1 & b_2 & b_3\\ \end{vmatrix} \\ & = & (a_2 b_3 - a_3 b_2)i + (a_3 b_1 - a_1 b_3)j + (a_1 b_2 - a_3 b_1)k \end{align}

(4)

As such, the magnitude of the resulting vector $c = a \times b$ is:

\left| c\right| = \left| a \right| \left| b \right| \sin{\alpha}

(5)

where $\alpha$ is the angle between vectors $a$ and $b$ .

Using this property, the cross product can be used to check if two vectors are parallel or antiparallel: in these scenarios, the cross product results in a zero vector with a length of 0:

a = np.array([ 1, 0, 0 ])     # this is a row vector
b = np.array([ 1, 0, 0 ])
a = a.reshape(-1,1)           # this is now a column vector
b = b.reshape(-1,1)

c = np.cross(a, b, axis=0)
print("The cross prouct is:\n",c)
print("with the length being:", np.linalg.norm(c))

# what does this line do?
print(np.asin( np.linalg.norm(c) / ((np.linalg.norm(a))*(np.linalg.norm(b))) ) * 180.0 / np.pi )

The cross prouct is:
 [[0]
 [0]
 [0]]
with the length being: 0.0
0.0

You should experiment with the above code, by using different values for $a$ and $b$ , and convince yourself that:

Indeed the cross product vector $c$ is perpendicular with both $a$ and $b$ , and
Indeed you can infer the angle between $a$ and $b$

Cross Product of two vectors — Figure 8:The cross product between the red vector with the green vector produces a blue vector that is perpendicular to both the red and green vectors.

For visualization, consult the following tutorial.

Orthonormal Basis from Three Points¶

One utility of cross product often encountered in CAS is to define a unique Cartesian coordinate system based on 3 points.

Orthonormal basis from 3 points — Figure 9:A orthonormal basis of $i,j,k$ can be created from 3 points $A$ , $B$ , and $C$ using vector subtraction and cross product, image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

The axes of the Orthonormal Basis¶

Given three points $A$ , $B$ , and $C$ , the process of creating a orthonormal basis while keeping the RHR are:

Create the axis $i$ by $B-A$ ,
Create the axis $k$ by $i \times (C-A)$ , and
Create the axis $j$ by $k \times i$
Normalize each axis

Note the order of which axis is computed, and try to understand why by visualizing. Normalization of the axes should be performed before taking the cross product to ensure numerical stability.

A = np.array([ 0, 0, 0 ])     # this is a row vector
B = np.array([ 2, 0, 0 ])
C = np.array([ 3, 4, 0 ])
A = A.reshape(-1,1)           # this is now a column vector
B = B.reshape(-1,1)
C = C.reshape(-1,1)

i = (B-A)/np.linalg.norm(B-A)

k = np.cross(i, (C-A)/np.linalg.norm(C-A), axis=0)
k = k/np.linalg.norm(k)

j = np.cross(k, i, axis=0)

print("Basis axis-i, i.e., x:\n", i)
print("Basis axis-j, i.e., y:\n", j)
print("Basis axis-k, i.e., z:\n", k)

Basis axis-i, i.e., x:
 [[1.]
 [0.]
 [0.]]
Basis axis-j, i.e., y:
 [[0.]
 [1.]
 [0.]]
Basis axis-k, i.e., z:
 [[0.]
 [0.]
 [1.]]

The Origin of the Orthonormal Basis¶

As a Cartesian coordinate system has an origin, i.e. a common point where all directional axes intersect at, we must define one.

There are two practical choices:

Desginate one of the three point as the origin, and associate all axes with it, or
use the Center of Gravity (COG), i.e. the centroid of the 3 points.

Mathematically, the origin of the newly define Cartesian coordinate system can be anywhere. The COG is a reasonable choice, calculated as the unweighted average of the points:

O = \frac{A+B+C}{3}

(6)

Use case for reference frame — Figure 10:Use case for defining a reference frame including tracking objects, image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

Defining a reference Cartesian coordinate system from three points is used in many techniques in CAS; mainly for tracking objects.

Reference Frame from points — Figure 11:A reference Cartesian frame $(i,j,k)$ defined from points $A$ , $B, and$ C $, where the origin$ O$ is located at the Centre of Gravity, image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

For visualization, consult the following tutorial.

Area of a Triangle¶

Visualization in CAS uses techniques from computer graphics including surface wire-frame or Volume rendering. The most common representation of an object with complex geometry is to represent it as a set of connected triangles via a process called polygon triangulation. In this regard, area of a triangle, and the summation of connected triangles, yields additional information such as the surface area of an organ.

Geometrical shapes represented via a series of connected triangles — Figure 12:Use case for triangulation: computing the surface dose of triangulated representation of anatomical structures (e.g. rectum wall) in radiotherapy, image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

The area of a triangle is half of the area of a parallelogram with the same base length and height.

:label: fig_area_triangle :alt: Area of a triangle :align: center :width: 80%

The area of a triangle as half of area of a parallelogram that has the same base length and height, image courtesy of Wikimedia Commons here, accessed January 9, 2026.

If a triangle is presented by three points $b$, $c$, with the angle $\alpha$ between them, The base of the triangle can be taken as $c$ and the height of the triangle, $h$, is equal to the length of $b$ multiplied by $\sin{\alpha}$:

```{figure} ./../images/002/fig_triangle_vectors.png
:label: fig_area_triangle_vectors
:alt: Area of a triangle
:align: center
:width: 80%

A triangle represented by 2 vectors $b$ and $c$, separated by an angle $\alpha$. The **base** is $c$ and the **height** is $h$.

That is, the area of the triangle is calculated as: $$

\begin{align} area & = & \frac{1}{2} \cdot base \cdot height \\ & = & \frac{1}{2} \left| c \right| \cdot \left| h \right| \\ & = & \frac{1}{2} \left| c \right| \cdot \left| b \right| \cdot \sin{\alpha} \\ \end{align}

(7)

$$

But $\left| c \right| \cdot \left| b \right| \cdot \sin{\alpha}$ is the length of the cross product $b \times c$ ! Thus, the area of a triangle can be computed as:

AREA = \frac{1}{2} \left| b \times c \right|

(8)

A = np.array([ 3, 0, 0 ])     # this is a row vector
B = np.array([ 0, 4, 0 ])
A = A.reshape(-1,1)           # this is now a column vector
B = B.reshape(-1,1)

Area = 0.5 * np.linalg.norm( np.cross(A, B, axis=0) )

print( "Area of the triangle is:", Area)

Area of the triangle is: 6.0

For visualization and other implementations, consult the following tutorial.

Dot product of Two Vectors¶

Dot product, or the inner product, of two vectors, is the sum of the element-wise product. The final value is a scalar and not a vector:

dot(u,v) = u \cdot v = (x_1 x_2 + y_1 y_2 + z_1 z_2) = \left| u \right| \left| v \right| \cos{\alpha}

(9)

where $\alpha$ is the angle between these two vectors.

Dot product between two vectors — Figure 13:Illustration showing how to find the angle between vectors using the dot product, image courtesy of Wikimedia Commons here, accessed January 9, 2026.

Note the properties of dot product:

Commutative, $a b = b a$
Not Associative (but why?): $(a \cdot b) c \ne a (b \cdot c)$ , but
Distributive with addition: $(a+b) \cdot c = a \cdot c + b \cdot c$

It can be trivially shown that the length of a vector is just the square root of the dot product with itself:

length(v) = \sqrt{x^2 + y^2 + z^2} = \sqrt{dot(v,v)}

(10)

Because the dot product is related to the cosine of the angle between two vectors, we can use it to evaluate if two vectors are orthogonal with eather other: if their dot product is 0, implying that $\cos{\alpha}$ is 0, these two vectors must be perpendicular:

A = np.array([ 3, 0, 0 ])     # this is a row vector
B = np.array([ 0, 4, 0 ])
A = A.reshape(-1,1)           # this is now a column vector
B = B.reshape(-1,1)

# np.dot(A,B) ## won't work as numpy assume row vector. Instead, use this:
print("The dot product is:", A.T@B)

The dot product is: [[0]]

Angle between Two Vectors¶

Given two vectors $u$ and $v$ with the angle $\alpha$ in between, we now have 2 ways to calculate the angles:

Via dot product

\begin{align} \frac{u \cdot v}{\left| u \right| \left| v \right|} & = & \cos{\alpha} \\ \alpha & = & \arccos{\alpha}\\ \end{align}

(11)

or via cross product

\begin{align} \frac{\left| u \times v \right|}{\left| u \right| \left| v \right|} & = & \sin{\alpha} \\ \alpha & = & \arcsin{\alpha} \end{align}

(12)

Given a situation, which one should you use? One recommendation is to consider the slope of the $\sin$ and $\cos$ function and avoid the flat area of the curve:

Figure 14:Sine and cosine as function of the angle, image courtesy of Wikimedia Commons here, accessed January 9, 2026.

Around $0\degree$ use sine, i.e. cross product,
Around $90\degree$ use cosine, i.e. dot product

A more robust method is to use both the cross and dot products:

atan2( norm( cross(u,v)), dot(u,v))

Here is an example:

a = 1e-10                               # a small number
u = 4.0 * np.array([[1],                # a non-unit vector along the x-axis
                    [0],
                    [0]])             
v = 5.0 * np.array([[ math.cos(a) ],    # a non-unit vector different from u
                    [ math.sin(a) ],    # by a very small angle
                      [0]]) 

dotproduct = ((u.T@v)[0])[0]
crossproduct = np.linalg.cross( u, v, axis=0 )
print( "The angle in radian between 2 vectors is:\n", a )
print( "The angle computed using dot product (cosine):\n",
      math.acos( dotproduct / ((np.linalg.norm(u))*(np.linalg.norm(v)))) )
print( "The angle computed using cross product (sine):\n",
      math.asin( np.linalg.norm(crossproduct) / ((np.linalg.norm(u))*(np.linalg.norm(v)))) )
print( "The angle computed using both cross and dot products: \n",
      math.atan2( np.linalg.norm( crossproduct ), dotproduct ))

The angle in radian between 2 vectors is:
 1e-10
The angle computed using dot product (cosine):
 0.0
The angle computed using cross product (sine):
 1e-10
The angle computed using both cross and dot products: 
 1e-10

An use case is to determine the Cobb angle, a measurement of the bending disorders of the vertebral column, for patients with scoliosis based on an x-ray image:

Cobb angle to measure scoliosis — Figure 15:Cobb angle is determined from an frontal view of surface-rendered CT volume based on the angle between two vectors, image courtesy of Prof. Gabor Fichtinger of Queen’s University, Canada.

Vector Representation of a Line¶

There are many ways to represent a line. In fact, the appropriate representation often depends on mathematical convenience and the definition of a line. In CAS, a line could mean

an oriented line: a line with a point on the line (i.e. the origin) and a direction (often represented by a unit vector)
a line segment: a line with length
- can be represented by the two end points of the line segment, or
- can be represented by one of the end point plus a vector of the known length

Line representation — Figure 16:A line can be presented by a fixed point $P$ , a unit dictional vector $v$ , and a parameter $t$ , where another point $L$ on the lie is $L=P+t \cdot v$ .

Parameterized Line¶

Let $P$ be a point on a line, and a unit vector $v$ be the direction (orientation) of the line. Another point $L$ on the line can be parameterized as:

L = P + t \cdot v

(13)

where $t$ is a scalar, effectively specifying how far $L$ is from $P$ . Specifically, $L$ is calculated as

\begin{align} L_x & = & P_x + t * v_x\\ L_y & = & P_y + t * v_y\\ L_z & = & P_z + t * v_z\\ \end{align}

(14)

The orientation vector $v$ does not need to be a unit vector, but mathematically convenient if it is. The range of $t$ specify the type of the line:

Infinite line: $t \in \{ -\infty, \infty \}$
Ray: $t \in \{ 0, \infty \}$
Line Segment: $t \in \{ t_{min}, t_{max} \}$

To understand what the parameterizing parameter $t$ is, let’s rearrange the equation of the line:

\begin{align} L & = & P + t \cdot v\\ L-P & = & t \cdot v \\ \left| L - P \right| & = & t \left| v \right|\\ \end{align}

(15)

Because $\left| v \right| = 1$ , i.e. $v$ was normalized, then

\left| L - P \right| = t

(16)

That is, if $v$ is a unit vector, the parameter $t$ measures the distance between the fixed point $P$ and the running point $L$ .

P = np.array([[0],          # a fixed point on the line
              [0],
              [0]])  
v = np.array([[1],          # orientation of the line as a unit vector
              [0],
              [0]])
t = 3.0                     # line parameterization
L = P + t * v               # another point on the line
print( L )

[[3.]
 [0.]
 [0.]]

Line Defined by Two Points¶

Give two points $P$ and $L$ on a line, the orientation of the line can be determined by

dir = L-P

(17)

The orientation unit vector $v$ is

v = \frac{dir}{\left| dir \right|}

(18)

That is, the orientaiton of the line is determined by substracting one point from another; the order of which is important as it determines the direction of the line. The orientation vector is normalized for mathematical convenience.

P = np.array([[0],          # a fixed point on the line
              [0],
              [0]])  
L = np.array([[3],          # another point on the same line
              [0],
              [0]])
dir = L-P                   # orientation of the line
v = dir/np.linalg.norm(dir) # normalized orientation, i.e. a unit vector
print(v)

[[1.]
 [0.]
 [0.]]

Distance of a Point from a Line¶

The distance between a point $A$ and a line $L$ is defined as the shortest distance beween $A$ to any points on $L$ . Suppose a line is represented by the fixed point on the line $P$ and a unit orientation vector $v$ . The length $c=\left| \overrightarrow{AP} \right|$ is the distance between point $A$ and $P$ and is the length of the hypotenuse of a right-angle triangle. The projection of $\overrightarrow{AP}$ onto $L$ is the adjacent of the right-angle triangle, with the angle denoted as $\alpha$ . The distance between the point $A$ and line $L$ is the length of the opposite of the right-angle triangle.

Distance between a point and a line — Figure 18:The distance between a point $A$ and line $L$ can be calculated based a triangle being formed.

Using this right-angle triangle, one can calculate the distance $d$ using:

\begin{align} c & = & \left| A - P \right| \\ a & = & \left| v \cdot (A-P) \right| \\ d^2 & = & c^2 - a^2\\ d^2 & = & (c+a)(c-a)\\ d & = & \sqrt{(c+a)(c-a)}\\ \end{align}

(19)

P = np.array([[0],              # a fixed point on the line
              [0],
              [0]])
v = np.array([[1],              # a unit vector representing the orientation of the line
              [0],
              [0]])
A = np.array([[3],              # a point NOT on the line
              [4],
              [0]])

c = np.linalg.norm( A-P )
a = np.linalg.norm( v.T@(A-P) )
d = math.sqrt( (c+a)*(c-a) )    # just a more efficient way of calculating c*c-a*a
print(d)

4.0

This formulation calculates the distance between a point $A$ and an line $L$ without knowing which point on the $L$ is the closest to $A$ . However, since we know $a$ , the length of the adjacent of the right-angle triangle, we can determine which point on $L$ is closest to A:

Q = P+a*v

(20)

Q = P + a*v                     # Q is a point on line L that is closest to point A
print( Q )
print( np.linalg.norm(A-Q))     # The length of the vector AQ

[[3.]
 [0.]
 [0.]]
4.0

A use case in CAS: suppose we want to insert a biopsy needle (a line) into the center of the prostate (a point) to take a sample; one may wish to know how far away the biopsy needle is to the intented target.

Use case for the distance between a point and a line — Figure 19:Screen capture in ultrasound-guided prostate biopsy navigation. Computing the needle placement error as the distance between the planned target (a point) and the actual needle path observed in ultrasound (a line), image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

Plane¶

A plane is represented by a point $A$ on the plane and the surface normal $n$ of the plane. For mathematical convenience, the normal of the plane $n$ is normalized with a length of 1.

Given 3 points $A$ , $B$ , and $C$ on a plane, the plane normal vector $n$ can be computed by their cross product:

n = normalize( (B-A) \times (C-A) )

(21)

A = np.array([[0],              # a point on the plane
              [0],
              [0]])
B = np.array([[1],              # a point on the plane
              [0],
              [0]])
C = np.array([[0],              # a point on the plane
              [1],
              [0]])

dir = np.linalg.cross( (B-A), (C-A), axis=0 )
n = dir/np.linalg.norm(dir)
print( "plane normal vector is:\n", n)

plane normal vector is:
 [[0.]
 [0.]
 [1.]]

Distance from a Plane¶

The distance between a point $P$ from a plane $(n, A)$ , where $n$ is a unit vector specifying the normal of the plane, and $A$ is a point on the plane, is the length of the projection of $\overrightarrow{AP}$ to the surface normal $n$ :

D = (L-A) \cdot n

(22)

Distance between a plane and a point — Figure 22:The distance between a point $P$ and plane $(n, A)$ can be calculated as a dot product.

Notice the plane normal $n$ has a direction:

if $d>0$ , then $P$ is above the plane,
if $d<0$ , then $P$ is below the plane,
if $d=0$ , then $P$ is on the plane.

A = np.array([[0],                 # a point on the plane
              [0],
              [0]])
n = np.array([[0],                 # plane normal
              [0],
              [1]])
P = np.array([[3],                 # a point 
              [4],
              [5]])

d = (n.T@(P-A))[0][0]
print("The distance is:", d)

The distance is: 5

Intersection of Line and Plane¶

The intersection between a line, $L=P+vt$ , and a plane, $(n,A)$ can be calculated as the following.

Line and plane intersection — Figure 23:The intersection between a parameterized line $P+vt$ to a plane $(n,A)$ .

Assuming the line is parameterized with $t$ :

L=P+v t

(23)

When $L$ intersects with the plane, i.e. $L$ is on the plane, then

(L-A) \cdot n = 0

(24)

Substitute Eqn.(23) into Eqn.(24):

(P+v t -A) \cdot n = 0

(25)

Using the distributive property of dot product, one obtains:

(P-A)\cdot n + (v \cdot n) t = 0

(26)

Rearrange for $t$ :

t = \frac{(A-P) \cdot n}{ v \cdot n}

(27)

Substitute $t$ back to Eqn.(23) to obtain the intersection point $L=P+v t$ .

For visualization and other implementations, consult the following tutorial.

One of the use case involves the use of X-ray or Fluoroscopy where X-ray, originated from a point source, is treated as a line, and the image detector is the plane.

Use case for line and plane intersection — Figure 24:The intersection between an X-ray to the detector as a use case for calculating the intersection between a line and a plane, image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

Intersection of Lines in 3D¶

Intersection of lines in 3D is common in motion analysis, often involved the usa of stereo camera of bi-plane fluoroscopy. In Figure 25, the X-ray source is treated as a Pinhole camera, and the principle of triangulation is used to localize an object.

Bi-plane X-ray — Figure 25:Biplane X-ray target reconstruction, image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

Let two lines be represented parametrically:

\begin{align} L_1 & = & P_1 + t_1 v_1\\ L_2 & = & P_2 + t_2 v_2\\ \end{align}

(28)

Expanding for each $(x,y,z)$ components, one obtains:

\begin{align} L_{1x} & = & P_{1x} + t_{1} v_{1x}\\ L_{1y} & = & P_{1y} + t_{1} v_{1y}\\ L_{1z} & = & P_{1z} + t_{1} v_{1z}\\ \end{align}

(29)

and

\begin{align} L_{2x} & = & P_{2x} + t_{2} v_{2x}\\ L_{2y} & = & P_{2y} + t_{2} v_{2y}\\ L_{2z} & = & P_{2z} + t_{2} v_{2z}\\ \end{align}

(30)

where

\begin{align} t_1 & \in & (-\infty, \infty) \ t_2 & \in & (-\infty, \infty) \\ \end{align}

(31)

If lines intersect with each other, they share a common point in space. That is,

\begin{align} L_1 & = & L_2 \\ L_1 - L_2 & = & 0 \end{align}

(32)

for some $L_1$ and $L_2$ , and

\begin{align} P_{1x} - P_{2x} + t_1 v_{1x} - t_2 v_{2x} & = & 0\\ P_{1y} - P_{2y} + t_1 v_{1y} - t_2 v_{2y} & = & 0\\ P_{1z} - P_{2z} + t_1 v_{1z} - t_2 v_{2z} & = & 0\\ \end{align}

(33)

There $t_1$ and $t_2$ are unknown. You should note that these are 3 equations with 2 unknowns, thus, in general, two lines do not intersect in 3D.

Figure 26:Intersection of 2 lines in 3D.

In other words, 2 lines intersect with each other if one of these three equations cancels out. That is, if one of the equations can be expressed as a linear combination of the other two: they lie on a common plane. However, this does not happen in general, either because two lines in fact do not intersect with each other, or they do but numerically the equations above cannot be satisfied due to digital representation.

But we can in fact calculate a “symbolic intersection”: find the shortest distance between two lines and the midway point.

The symbolic intersection point $M$ is a point that is equal distance from both lines. It can be computed in 3 steps:

Find a line $L_3$ that is perpendicular to both lines,
Find the intersection of line $L_3$ with both lines $L_1$ and $L_2$ , and
Find the middle point $M$ on $L_3$

Symbolic intersection of 2 lines in 3D. — Figure 27:Symbolic intersection, $M$ , of 2 lines in 3D. $v_3$ is perpendicular to both $v_1$ and $v_2$ and $M$ is mid-way from $L_1$ and $L_2$ with a distance $\frac{d}{2}$ .

The distance between the mid-way point $M$ is denoted as the Residual Eerror Metric (REM), which is $\frac{d}{2}$ , and $d$ is the distance between lines $L_1$ and $L_2$ .

Line 3 $L_3=P_3 + t_3 v_3$ is perpendicular to the other two lines, thus, its directional vector $v_3$ can be computed as the cross product of the other two directional vectors:

v_3 = normalize( v_1 \times v_2 )

(34)

The length of $v_3$ is normalized to be 1, again, for mathematical convenience.

As such, the fixed point $P_3$ and the parameterized point $L_3$ on this line can be designate as

P_3 = L_2

(35)

L_3 = L_1

(36)

Because the equations for all three lines are:

L_1 = P_1 + t_1 v_1

(37)

L_2 = P_2 + t_2 v_2

(38)

L_3 = P_3 + t_3 v_3

(39)

Substituting Eqn.(35) and Eqn.(36) into Eqn.(39) yields:

\begin{align} L_3 & = & P_3 + t_3 v_3 \\ P_1 + t_1 v_1 & = & P_2 + t_2 v_2 + t_3 v_3 \end{align}

(40)

After rearrangement:

P_1 - P_2 = - t_1 v_1 + t_2 v_2 + t_3 v_3

(41)

and expanding it to each of the 3 dimensions:

\begin{align} P_{1x} - P_{2x} & = & - t_1 v_{1x} + t_2 v_{2x} + t_3 v_{3x} \\ P_{1y} - P_{2z} & = & - t_1 v_{1y} + t_2 v_{2y} + t_3 v_{3y} \\ P_{1z} - P_{2z} & = & - t_1 v_{1z} + t_2 v_{2z} + t_3 v_{3z} \\ \end{align}

(42)

where $t_1$ , $t_2$ and $t_3$ are unknown. Now there are 3 equations and 3 unknown, implying that a unique solution does exist.

These equations can be represented using matrices:

\begin{bmatrix} P_{1x} - P_{2x} \\ P_{1y} - P_{2y} \\ P_{1z} - P_{2z} \end{bmatrix} = \begin{bmatrix} - v_{1x} & v_{2x} & v_{3x}\\ - v_{1y} & v_{2y} & v_{3y}\\ - v_{1z} & v_{2z} & v_{3z}\\ \end{bmatrix} \begin{bmatrix} t_{1} \\ t_{2} \\ t_{3} \\ \end{bmatrix}

(43)

The unknown $t_1$ , $t_2$ and $t_3$ can be solved using

Gaussian elimination
Algebraic substitution, or
Matrix inversion (shown below)

$$

\begin{align} \begin{bmatrix} P_{1x} - P_{2x} \\ P_{1y} - P_{2y} \\ P_{1z} - P_{2z} \end{bmatrix} & = & \begin{bmatrix} - v_{1x} & v_{2x} & v_{3x}\\ - v_{1y} & v_{2y} & v_{3y}\\ - v_{1z} & v_{2z} & v_{3z}\\ \end{bmatrix} \begin{bmatrix} t_{1} \\ t_{2} \\ t_{3} \\ \end{bmatrix} \\ \begin{bmatrix} - v_{1x} & v_{2x} & v_{3x}\\ - v_{1y} & v_{2y} & v_{3y}\\ - v_{1z} & v_{2z} & v_{3z}\\ \end{bmatrix}^{-1} \begin{bmatrix} P_{1x} - P_{2x} \\ P_{1y} - P_{2y} \\ P_{1z} - P_{2z} \end{bmatrix} & = & \begin{bmatrix} - v_{1x} & v_{2x} & v_{3x}\\ - v_{1y} & v_{2y} & v_{3y}\\ - v_{1z} & v_{2z} & v_{3z}\\ \end{bmatrix}^{-1} \begin{bmatrix} - v_{1x} & v_{2x} & v_{3x}\\ - v_{1y} & v_{2y} & v_{3y}\\ - v_{1z} & v_{2z} & v_{3z}\\ \end{bmatrix} \begin{bmatrix} t_{1} \\ t_{2} \\ t_{3} \\ \end{bmatrix} \\ \begin{bmatrix} - v_{1x} & v_{2x} & v_{3x}\\ - v_{1y} & v_{2y} & v_{3y}\\ - v_{1z} & v_{2z} & v_{3z}\\ \end{bmatrix}^{-1} \begin{bmatrix} P_{1x} - P_{2x} \\ P_{1y} - P_{2y} \\ P_{1z} - P_{2z} \end{bmatrix} & = & I \begin{bmatrix} t_{1} \\ t_{2} \\ t_{3} \\ \end{bmatrix}\\ \begin{bmatrix} - v_{1x} & v_{2x} & v_{3x}\\ - v_{1y} & v_{2y} & v_{3y}\\ - v_{1z} & v_{2z} & v_{3z}\\ \end{bmatrix}^{-1} \begin{bmatrix} P_{1x} - P_{2x} \\ P_{1y} - P_{2y} \\ P_{1z} - P_{2z} \end{bmatrix} & = & \begin{bmatrix} t_{1} \\ t_{2} \\ t_{3} \\ \end{bmatrix} \end{align}

(44)

$$

Once $t_1$ , $t_2$ , and $t_3$ are calculated, substitute them back to

\begin{align} L_1 & = & P_1 + t_1 v_1 \\ L_2 & = & P_2 + t_2 v_2 \end{align}

(45)

to solve for $L_1$ and $L_2$ . Note that

t_3 = \left| L_1 - L_2 \right|

(46)

and

M = \frac{L_1+L_2}{2}

(47)

The REM, which is the distance of $M$ from either $L_1$ or $L_2$ is:

REM = \left| L_1 - M \right| = \left| L_2 - M \right| = \frac{t_3}{2}

(48)

P1 = np.array([[0],                 # a point on z=0 plane
              [0],
              [0]])
v1 = np.array([[1],                 # line orientation +x
              [0],
              [0]])
P2 = np.array([[0],                 # a point on z=5 plane
              [0],
              [5]])
v2 = np.array([[0],                 # line orientation +y
              [1],
              [0]])

v3 = np.linalg.cross( v1, v2, axis=0 )


V = np.zeros((3,3))                 # construct the 3x3 matrix
V[:,0] = -1*v1.reshape(1,3)
V[:,1] =    v2.reshape(1,3)
V[:,2] =    v3.reshape(1,3)

P = P1-P2
t = np.matmul( np.linalg.inv(V), P)

# solve for L1 and L2
L1 = P1 + t[0] * v1
L2 = P2 + t[1] * v2
M = 0.5 * (L1+L2)
REM = .5 * t[2][0]
print("The symbolic intersection is: \n", M,
      "\nwith REM of ", REM)

The symbolic intersection is: 
 [[0. ]
 [0. ]
 [2.5]] 
with REM of  -2.5

Symbolic Intersection of $n$ Lines in 3D¶

One utility of computing the intersection of $n$ lines in 3D is to find the centre of the volume in a reconstructed CT volume.

Rotational X-ray Reconstruction. — Figure 29:Rotational X-ray Reconstruction, where center of the imaged volume is the intersection of multiple lines in 3D.

First, write the equation for all ines

\begin{align} L_1 & = & P_1 + t_1 v_1\\ L_2 & = & P_2 + t_2 v_2\\ & \cdots & \\ L_n & = & P_n + t_n v_n\\ \end{align}

(49)

The the problem of computing the symbolic intersection can be posed as a least-square optimization with the cost function of:

S = {d_1}^2 + {d_2}^2 + \cdots + {d_n}^2

(50)

and compute $M$ where $S$ is the minimum.

Intersection of N lines — Figure 30:Intersection of $n$ lines in 3D where the point $M$ is the symbolic intersection, having the shortest sum of distance from all $n$ lines.

Throughout the optimization process, one should always be mindful of deriving and monitoring error metrics including Root Mean Square (RMS), average (Avg), and standard deviation (STD), and use them a means to detect and remove Outlier. For example, an outlier can be defined as:

\left| Avg - d_i \right| > (k \cdot STD )

(51)

or

\left| RMS - d_i \right| > (k \cdot STD )

(52)

(typically $k=2$ or $k=1$ ). Once outlier is detected, they are moved from the population and the optimization is recomputed, so are the error metrics. This process is repeated until no outlier is left, or the input data is completely removed (a bad situation).

One easy way to approximate the symbolic intersection point $M$ is to compute pairwise intersections for all pairs $(M_{12}, M_{13}, \cdots, M_{1n}, M_{23}, \cdots)$ and approximate $M$ as the average, i.e. COG, of these pairwise intersections:

Approximation of the intersection of N lines — Figure 31:Approximation of the symbolic intersection of $n$ lines in 3D as the centroid of the pair-wise line intersections.

If the input is decent, i.e. no outliers and all lines come very close, this simple average is almost as good as the proper least square minimum.

In fact, there is indeed no outlier, a closed-form solutoin exist. Let $P_i$ be the fixed point on $i^{th}$ line and $v_i$ be the unit vector along the $i^{th}$ line, then $M$ can be calculated as:

M = \left( \sum_i I - v_i {v_i}^T \right)^{-1} \left( \sum_i \left(I - v_i {v_i}^T \right) P_i \right)

(53)

P1 = np.array([[5],                 # 
              [0],
              [0]])
v1 = np.array([[1],                 # 
              [0],
              [0]])
P2 = np.array([[0],                 # 
              [0],
              [9]])
v2 = np.array([[0],                 # 
              [1],
              [0]])
P3 = np.array([[1],                 # 
              [1],
              [0]])
v3 = np.array([[0],                 # 
              [0],
              [1]])

# initialization
I = np.identity(3)
Sum1 = np.zeros((3,3))
Sum2 = np.zeros((3,1))
# 1st point
Sum1 += I-v1@v1.T
Sum2 += np.matmul(I-(v1@v1.T), P1)

# second point
Sum1 += I-v2@v2.T
Sum2 += np.matmul(I-(v2@v2.T), P2)

# third point
Sum1 += I-v3@v3.T
Sum2 += np.matmul(I-(v3@v3.T), P3)

# solution
M = np.matmul(np.linalg.inv(Sum1), Sum2)
print(M)

[[0.5]
 [0.5]
 [4.5]]

This is a good time to learn about iterative solution vs a closed-form solution. In mathematics, a closed-form solution is often optimal with regarding to certain error metric, in this case, the Euclidean distance. However, closed-form solution are sensitive to outlier.

On the other hand, iterative solution (i.e. via optimization) may not reach a solution that is optimal with respect to the given error metric. However, it allows rooms to integrate outlier detection and removal naturally during the iterative process.

Cartesian Equations of Sphere and Ellipsoid¶

Spheres and ellipsoid are often used to approximate the size of anatomical structures, such as a tumour, in 3D from medical imaging.

Sphere and ellipsoid to determine sizing. — Figure 32:Sphere and ellipsoid are used to determine the size of anatomical structures in 3D, image courtesy of Prof. Gabor Fichtinger at Queen’s University, Canada.

A sphere is represented by its center $C(x_0, y_0, z_0)$ and the radius $R$ , with the corresponding equation of:

(x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 = R^2

(54)

which can be used to test if a point $P=(x,y,z)$ is on the surface of the sphere or not.

An ellipsoidis represented by its center $C(x_0, y_0, z_0)$ and the lengths of principal axes $(a,b,c)$ :

mathematically represented as

\frac{(x-x_0)^2}{a^2}+\frac{(y-y_0)^2}{b^2}+\frac{(z-z_0)^2}{c^2}=1

(55)

Intersection of Line and Sphere/Ellipsoid¶

Gamma rays enter head in radiosurgery — Figure 35:An ellipsoid is defined by its center $C$ and the lengths of principal axes.

Let the line be presented by $L=P+t \cdot v$ the $P$ being the fixed point on the line, $v$ a unit vector representing the direction of the line, and $t$ being the scalar.

L = P + t \cdot v

(56)

And the equation of the sphere is:

\left| Q - C \right| = R^2

(57)

where $Q$ is a point on the surface of the sphere, $C$ is the center of the sphere and $R$ is the radius. If the point $L$ on the line intersects with the sphere, $L$ must also satisfy the equation of the sphere:

\begin{align} \left| Q - C \right|^2 & = & R^2 \\ \left| P + t \cdot v - C \right|^2 & = & R^2 \\ \end{align}

(58)

which can be evaluated using dot product:

(P + t \cdot v - C) \cdot (P + t \cdot v - C)= R^2

(59)

Using the distributive property of dot product and rearrange the terms:

t^2 (v \cdot v) + 2 t \left[ v \cdot (P-C) \right] + (P-C)\cdot(P-C) - R^2 = 0

(60)

which has the form of a quadratic formula:

a t^2 + b t + c = 0

(61)

where $$

\begin{align} a & = & v \cdot v = \left| v \right|^2\\ b & = & 2 \left[ v \cdot (P-C) \right] \\ c & = & (P-C)\cdot(P-C) - R^2\\ \end{align}

(62)

$$

and $t$ can be solved as:

d=\frac{-2 \left[ v \cdot (P-C) \right] \pm \sqrt{ (2 \left[ v \cdot (P-C) \right])^2 - 4 (\left| v \right|^2)((P-C)\cdot(P-C) - R^2)}}{2 \left| v \right|^2}

(63)

Because $v$ is a unit vector, the above solution can be simplified as:

d=-\left[ v \cdot (P-C) \right] \pm ( \left[ v \cdot (P-C) \right])^2 - ((P-C)\cdot(P-C) - R^2)

(64)

Eqn.(63) is worth investigating: let $\Delta$ be the term inside the square root:

If $\Delta < 0$ , there is no solution: the line does not intersect the sphere,
If $\Delta = 0$ , there is exactly one solution: the line is tangent to the sphere,
If $\Delta > 0$ , there are 2 solutions: the line intersect the sphere in 2 locations (i.e. going through)

C = np.array([[0],                  # centre of a sphere
              [0],
              [0]])
R = 5                               # radius of the sphere
v = np.array([[-1],                  # orientation of a line (unit vector) 
              [0],
              [0]])
P = np.array([[10],                  # a fixed point on the line
              [0],
              [0]])

PC=P-C
print(PC)
determinant = ((v.T@PC)*(v.T@PC) - PC.T@PC + R*R)[0][0]
print(determinant)
a = (-1 * v.T@PC)[0][0]

if determinant < 0:
    print("no solution")
elif determinant == 0:
    t = a
    L = P + t * v
    print("one intersection at: \n", L)
else:
    t = a + math.sqrt(determinant)
    L = P + t * v
    print("one intersection at: \n", L)
    t = a - math.sqrt(determinant)
    L = P + t * v
    print("one intersection at: \n", L)

[[10]
 [ 0]
 [ 0]]
25
one intersection at: 
 [[-5.]
 [ 0.]
 [ 0.]]
one intersection at: 
 [[5.]
 [0.]
 [0.]]

For visualization and other implementations, consult the following tutorial.

Ellipsoid¶

Similiarly, the intersection between a line and an ellipsoid can be expressed as an equation with a quadratic form. For completeness, let the line be presented as:

\begin{align} L & = &P+t \cdot v\\ L_x & = &P_x+t \cdot v_x\\ L_y & = &P_y+t \cdot v_y\\ L_z & = &P_z+t \cdot v_z\\ \end{align}

(65)

where $P$ is a fixed point on the line, $v$ be the unit orientation vector, $t$ a scalar, and $L$ is the parameterized point on the same line.

An ellipsoid is defined by its center $C$ and the lengths of principal axes, mathematically represented as

\frac{(x-x_0)^2}{a^2}+\frac{(y-y_0)^2}{b^2}+\frac{(z-z_0)^2}{c^2}=1

(66)

The point $L$ on the line is on the surface of the ellipsoid when

\begin{align} \frac{(L_x-x_0)^2}{a^2}+\frac{(L_y-y_0)^2}{b^2}+\frac{(L_z-z_0)^2}{c^2}&=&1\\ \frac{((P_x+t \cdot v_x)-x_0)^2}{a^2}+\frac{((P_y+t \cdot v_y)-y_0)^2}{b^2}+\frac{((P_z+t \cdot v_z)-z_0)^2}{c^2}&=&1 \end{align}

(67)

Expand and simplifying:

\begin{align} \frac{(P_x-x_0)^2+(t v_x)^2 + 2(P_x-x_0)(t v_x)}{a^2}+\frac{(P_y-y_0)^2+ (t v_y)^2 + 2(P_y-y_0)(t v_y)}{b^2}+\frac{(P_z-z_0)^2+(t v_z)^2+2(P_z-z_0)(t v_z)}{c^2}&=&1\\ \left[ \frac{{v_x}^2}{a^2} + \frac{{v_y}^2}{b^2} +\frac{{v_z}^2}{c^2} \right] t^2 + \left[ 2\frac{(P_x-x_0)v_x}{a^2} + 2\frac{(P_y-y_0)v_y}{b^2} + 2\frac{(P_z-z_0)v_z}{c^2} \right] t + \left[ \frac{(P_x - x_0)^2}{a^2} + \frac{(P_y - y_0)^2}{b^2} + \frac{(P_z - z_0)^2}{c^2} - 1 \right]&=&0 \end{align}

(68)

which is a quadratic equation of $t$ . Similar to the case of intersection between a line and a sphere, there will be three cases:

The line does not intersect with the ellipsoid,
The line intersects with the ellipsoid at one location, i.e. the line is tangent to the ellipsoid, or
The line intersects with the ellipsoid at two locations.

Below, you should try to implement this equation and come up with a evaluation case to verify your implementation.

For visualization and other implementations, consult the following tutorial.

Fitting Points to Sphere in 3D¶

One if the critical use case in CAS for fitting points to a sphere in 3D is tool-tip calibration where, when a tracking sensor is attached to a (as an example) to a surgical drill, the location of the drill tip can be determined with respect to the tracking sensor via pivot calibration.

In pivot calibration, where the tip of the surgical drill remains stationary but the body of the surgical drill pivot about this stationary point, the path of any point on the body of the surgical drill lies on the surface of a sphere (with an error margin $d_i$ ). Thus, the problem of fitting a set of measured points $P_i$ to a sphere centered at $c(x_0,y_0,z_0)$ with a radius $r$ can be treated as a least square optimization, with the error to minimize being the sum of square distance of all measured points from the surface of the sphere:

S={d_1}^2 +{d_2}^2 +{d_3}^2 + \cdots {d_n}^2

(69)

and it to find an optimal $(R, x_0, y_0, z_0)$ . Akin with other optimization approach, you should always compute an error metric, for examples RMS, Avg, and STD, and set rules for outliers to maintain data integrity.

Many closed for solutions also exist and are suitable if there is no data outliers. Al-Sharadqah & Chernov (2009) provides a good background materials from a pure mathematical point of view (but for fitting circle in 2D). A mathod utilizing the Singular Value Decomposition (SVD) with hyperaccuracy and no essential bias is coded in Matlab below.

function Par = HyperSVD( XYZ )
%
% Algebraic sphere fit with "hyperaccuracy" 
% i.e. with zero essential bias
%
%
% INPUTS:  XYZ - 3xn matrix of 3D points that are (noisy) measurements of
%                points laying on a sphere
%
% OUTPUT:  Par = [ a b c R ] is the fitting sphere, 
%              centered at ( a, b, c ) with radius R
%
% Based on the paper by
% 
%     Error analysis for circle fitting algorithms
%     A. Al-Sharadqah and N. Chernov
%     Electronic Journal of Statistics
%     Vol. 3, 2009
%     page 886--911
%
% Elvis C.S. Chen
% chene@robarts.ca
%
% Jan. 30, 2013
%
% heavily based on the original 2D circle fit by Chernov
%

% assume row vectors
centroid = mean( XYZ, 2 );

% de-mean and calculate the squares
X = XYZ(1,:) - centroid(1);
Y = XYZ(2,:) - centroid(2);
Z = XYZ(3,:) - centroid(3);
ZZ = X.*X + Y.*Y + Z.*Z;

% construct the matrix Z based on equation (21)
ZZZ = [ ZZ' X' Y' Z' ones( length(Z),1) ];
[U,S,V] = svd( ZZZ, 0 );

if ( S(5,5)/S(1,1) < 1e-12 ) % singular case
    A = V(:,5);
else
   % consult page 904 for an algorithm describing the following code
   R = mean( ZZZ ); 
   
   % construct the constrain matrix N based on equation (87)
   %
   % this is stupid, R(2:4) are zero since it is de-mean already
   %N = [ 8*R(1) 4*R(2) 4*R(3) 4*R(4) 2; ...
   %      4*R(2) 1      0      0      0; ...
   %      4*R(3) 0      1      0      0; ...
   %      4*R(4) 0      0      1      0; ...
   %      2      0      0      0      0 ];
   %N = eye(5);
   %N(1,1) = 8*R(1);
   %N(1,5) = 2;
   %N(5,1) = 2;
   %N(5,5) = 0;
   
   % we can build the inverse of N directly
   invN = eye(5);
   invN(1,1) = 0;
   invN(1,5) = 0.5;
   invN(5,1) = 0.5;
   invN(5,5) = -2*R(1);
   W = V * S * V';
   %[E,D] = eig( W * inv(N) * W );
   [E,D] = eig( W * invN * W );
   [Dsort,ID] = sort(diag(D));
   Astar=E(:,ID(2)); % find the smallest positive eigen value
   A = W\Astar;
end

% recover the natural parameter based on equation (12)
Par = -.5 * A(2:4)'/A(1) + centroid'; % center
Par = [ Par, sqrt( A(2)^2 + A(3)^2 + A(4)^2 - 4*A(1)*A(5))/abs(A(1))/2 ];
end

Fitting Points to line in 3D¶

Fitting a set of measured points in 3D to a line effectively finds the orientation of the line. Using the tracked surgical drills as an example,

The pivot calibration determines the drill tip location respecto to the tracking sensor, and
fitting points measured along the drill bit and fitting them to aline allows the orientation of the drill bit to be determined.

Figure 37:Fitting points to line in 3D, where $P_i$ are observed points, $P$ is the **average**, i.e. centre of gravity or centroid, of all $P_i$ , $c_i$ are the unit orientation vector from $P$ to each $P_i$ , and $v$ is the average of all $c_i$ .

The core concept for an easy solution is that, assuming the input data is decent in the sense that

there are sufficiently large number of observed points, and
there are little outliers, then the average is just as good as the least square minimum. This can be achieved by computing both an average point $P$ and an average orientation vector $v$ :

\begin{align} P & = & Avg(P_1, P_2, \cdots, P_n)\\ v & = & Avg(c_1, c_2, \cdots, c_n)\\ \end{align}

(70)

One should normalize all $c_1, c_2, \cdots, c_n$ and flip them so they have the same genrral direction before computing $v$ .

As an optimization problem, the same precaution should be applied. One should always compute a set of error metric including but not limited to:

RMS $(d_1,d_2,\cdots,d_n)$
Avg $(d_1,d_2,\cdots,d_n)$
STD $(d_1,d_2,\cdots,d_n)$

and set rules for outliers, such as

$\|Avg-d_i \| > k*STD$ or
$\|RMS-d_i \| > k*STD$ or and remove these outliers and recompute the optimization and error metrics. Repeat until no outlier is left or the input has been exausted.

Closed form solution via PCA¶

A closed form solution can be achieved using the Principal component analysis (PCA).

Thoughts on Averaging vs. Optimization¶

Optimization tend to yield better results than simple averaging,
Optimization are computationally more expensive than simple averaging,
Optimization (especially multi-variate non-linear ones) can be finicky, and they can be stuck in local minima and fail to accept acceptable results,
Most optimizations require an initialization (a guess serving as a starting point): most of the time an estimate such as averaging can serve as one,
If the input is decent, i.e. the cardinality of the dataset is sufficiently large and the cardinality of the outlier is relatively small (or you found a way to remove these outlier), the optimization may not improve much on the initial guess -- so you get by with simple averaging,
The good news is: for many least-square optimizations, such as fitting lines, planes, and quadratic surfaces, the optimum is guaranteed to exist and can be found without an initial guess,
Many of these fitting solutions exists on the internet,
Many closed-form solutions also exist, but they assume the absence of outliers!

Point¶

Visualization¶

Step 1: Import External Libraries¶

Step 2: Rendering Pipeline¶

Step 3: Geometrical Primitives¶

Step 4: Rendering¶

Vector and Unit Vector¶

Scaling of a Vector¶

Normalizing a Vector¶

Sum of Vectors¶

Distance between two Points¶

Cross Product between two Vectors¶

Orthonormal Basis from Three Points¶

The axes of the Orthonormal Basis¶

The Origin of the Orthonormal Basis¶

Area of a Triangle¶

Dot product of Two Vectors¶

Angle between Two Vectors¶

Vector Representation of a Line¶

Parameterized Line¶

Line Defined by Two Points¶

Distance of a Point from a Line¶

Plane¶

Distance from a Plane¶

Intersection of Line and Plane¶

Intersection of Lines in 3D¶

Symbolic Intersection of nnn Lines in 3D¶

Cartesian Equations of Sphere and Ellipsoid¶

Intersection of Line and Sphere/Ellipsoid¶

Ellipsoid¶

Fitting Points to Sphere in 3D¶

Fitting Points to line in 3D¶

Closed form solution via PCA¶

Thoughts on Averaging vs. Optimization¶

Symbolic Intersection of $n$ Lines in 3D¶