I am not a professional mathematics theorist. These are just my personal notes. There are probably errors here. You’ve been warned!

Annoying Symbols

• ∘ Function composition as in f(g(x)) or g(f(x)) which are f∘g and g∘f respectively. The composition symbol ∘  is encoded as U+2218 ∘ ring operator (HTML ∘) Also Hadamard product.

• (a,b] - interval from a to b excluding a and including b.

• ℝ Fancy R - real numbers, U+211D ℝ DOUBLE-STRUCK CAPITAL R (HTML ℝ)

• ℤ Fancy Z - integers

• ℕ Fancy N - natural numbers

• ℚ Fancy Q - rational numbers

• ℙ Fancy P - irrational numbers, or primes, or probability

• ∃ backwards E - "there exists"

• ∀ upsidedown A - "for every"

• U+27E8 ⟨ MATHEMATICAL LEFT ANGLE BRACKET \langle

• U+27E9 ⟩ MATHEMATICAL RIGHT ANGLE BRACKET \rangle Inner product of vectors A and B can be written as ⟨A,B⟩

• ⌊π⌋ = 3 - Floor brackets

• ⌈π⌉ = 4 - Ceiling brackets

• ¬ - Logical negation.

• AT means A, but with its rows swapped for columns.

• ⊗ - Tensor product {1 2 3 4} ⊗ {1 1 2} = {{1 1 2} {2 2 4} {3 3 6} {4 4 8}}

• s.t. - "such that"

Vector Math

I feel vaguely guilty building a transformation matrix by concatenating rotation / transform / scale instead of directly.

2014-05-06
— John Carmack

Frankly, I find it excessively difficult and pointless to follow the gobbledygook notations of people like Pascal, LaPlace, Euler, and their pals. I can do it, but it’s a lot of work. Why not cut the ambiguity and go right to code? For all practical purposes, this will suffice. I keep needing to do geometry things and I think about looking up how to do that and writing code. Well, I’ve already written this kind of code a dozen times so just having it in one convenient place is what I really need. Here is one of my vector libraries simplified a bit; this won’t compile, but it should be clear what it does and how. It is worth noting that this is optimized for clarity, not performance. I usually have two functions for most things, one to act on the object itself and one to return an object; I mostly left out the former functions (referred to with "set" by the other code). If you need performance or just intense fanciness, why not check out third party libraries like…

Vector(){ x = 0; y = 0; z = 0; }
Vector(float X, float Y, float Z){ x=X; y=Y; z=Z; }
~Vector(){ }

float valX(){ return x; }
float valY(){ return y; }
float valZ(){ return z; }

Vector inverse() { Vector I; I.x=-x; I.y=-y; I.z=-z; return I; }

float magnitude() { float m=sqrt(x*x + y*y + z*z); return m;}

float distance(Vector P2) { Vector P=(*this-P2); return P.magnitude(); }

Vector direction() {
float mag;
Vector unit;
mag = magnitude(); // this->magnitude() works too...
if (mag > 0) {
unit.x = valX() / mag;
unit.y = valY() / mag;
unit.z = valZ() / mag; }
else { unit.x = 0; unit.y = 0; unit.z =0;}
return unit; }

Vector sum;
sum.x = x + offset.x;
sum.y = y + offset.y;
sum.z = z + offset.z;
return sum; }

Vector subtract(Vector offset) {
Vector difference;
difference.x = x - offset.x;
difference.y = y - offset.y;
difference.z = z - offset.z;
return difference; }

Vector rotated;
float s, c;
rotated.x = x*c - y*s;
rotated.y= y*c + x*s;
rotated.z= z;
return rotated;        }

Vector rotate(float &rads, Vector &basept) {
Vector rotated;
float y2, s, c;
rotated = this->subtract(basept);
y2= rotated.y*c + rotated.x*s;
rotated.x = rotated.x*c - rotated.y*s;
rotated.y = y2;
return rotated;              }

/*  rotation in XZ axis */
Vector rotated;
float s, c;
rotated.x = x*c - y*s;
rotated.z= z*c + z*s;
rotated.y= y;
return rotated;        }

Vector tilt(float &rads, Vector &basept) {
Vector rotated;
float z2, s, c;
rotated = this->subtract(basept);
z2= rotated.z*c + rotated.x*s;
rotated.x = rotated.x*c - rotated.z*s;
rotated.z = z2;
return rotated;              }

Vector viewvector(Vector &VV) {
if ((VV.x != 0) || (VV.y != 0)) {
Vector T, Xplus(1,0,0);
float r, t;
T.x=VV.x; T.y=VV.y; T.z=0;
r= ( T.y<0 ? T.angle(Xplus) : -T.angle(Xplus) );
t = (M_PI/2) -  atan(VV.z / sqrt(VV.x * VV.x + VV.y * VV.y) );
T = rotate(r);
T.set2_tilt(t);
r=-T.x; T.x=T.y; T.y=r; // turn 90deg- Z always points up
return T;}
else return *this; }

Vector perspective(float &lens, float &camdist) {
Vector T; float tz;
tz = camdist - z;
if (tz) {
T.x= (x * lens) / tz;
T.y= (y * lens) / tz;}
return T;  }

Vector scale(float &fact) {
Vector scaled;
scaled.x = x*fact;
scaled.y= y*fact;
scaled.z= z*fact;
return scaled;        }

Vector scale(float &fact, Vector &basept) {
Vector scaled;
scaled = this->subtract(basept);
scaled.x = scaled.x*fact;
scaled.y = scaled.y*fact;
scaled.z = scaled.z*fact;
return scaled;              }

Vector scale(Vector &fact3d) {
Vector scaled;
scaled.x = x * fact3d.x;
scaled.y = y * fact3d.y;
scaled.z = z * fact3d.z;
return scaled;        }

Vector scale(Vector &fact3d, Vector &basept) {
Vector scaled;
scaled = this->subtract(basept);
scaled.x = scaled.x * fact3d.x;
scaled.y = scaled.y * fact3d.y;
scaled.z = scaled.z * fact3d.z;
return scaled;              }

float dot(Vector &B) {
float dotprod;
dotprod = x * B.x + y * B.y + z * B.z;
return dotprod; }

float angle(Vector &B) {
float ang;
ang = magnitude() * B.magnitude();
if ( ang == 0 ) ang = 0;
else ang = acos( dot(B) / ang ) ;
if (!(ang == ang)) return 0; /* nan check!!!
else return ang; }

Vector midpoint(Vector &B) {
return (((*this - B) * .5) + B); }

float angle(Vector B, Vector basept) {
set2_subtract(basept); B.set2_subtract(basept);
return angle(B); }

Vector cross(Vector &B) {
Vector C;
C.x = y * B.z - z * B.y;
C.y = z * B.x - x * B.z;
C.z = x * B.y - y * B.x;
return C;  }

Geometry

Miscellaneous geometry notes. Possibly erroneous! Double check!

Point On A Line Perpendicular To Another Point

Points A and B are line endpoints, Point P is the point to get a perpendicular from off the line, Point N is the point on the line that is where P lines up at 90deg.

NA dot NP = 0
NB dot NP = 0
AN dot BN = 1

(Ax-Nx)(Px-Nx)+(Ay-Ny)(Py-Ny)+(Az-Nz)(Pz-Nz)=0
(Bx-Nx)(Px-Nx)+(By-Ny)(Py-Ny)+(Bz-Nz)(Pz-Nz)=0
(Ax-Nx)(Bx-Nx)+(Ay-Ny)(BY-Ny)+(Az-Nz)(Bz-Nz)=1

(Ax-Nx)(Px-Nx)+(Ay-Ny)(Py-Ny)+(Az-Nz)(Pz-Nz)=0
AxPx-NxPx-NxAx+Nx^2 + AyPy-NyPy-NyAy+Ny^2 + AzPz-NzPz-NzAz+Nz^2 = 0
BxPx-NxPx-NxBx+Nx^2 + ByPy-NyPy-NyBy+Ny^2 + BzPz-NzPz-NzBz+Nz^2 = 0
BxAx-NxAx-NxBx+Nx^2 + ByAy-NyAy-NyBy+Ny^2 + BzAz-NzAz-NzBz+Nz^2 = 1

-AxPx+NxPx+NxAx-Nx^2 = AyPy-NyPy-NyAy+Ny^2 + AzPz-NzPz-NzAz+Nz^2
-Nx^2 = AyPy-NyPy-NyAy+Ny^2 + AzPz-NzPz-NzAz+Nz^2 + AxPx+ NxPx +NxAx
Nx^2 = -AyPy+NyPy+NyAy-Ny^2 - AzPz+NzPz+NzAz-Nz^2 - AxPx- NxPx -NxAx

+AxPx -NxPx -NxAx -AyPy +NyPy +NyAy -Ny^2 - AzPz +NzPz +NzAz -Nz^2 -AxPx -NxPx -NxAx +AyPy -NyPy -NyAy +Ny^2 +AzPz -NzPz -NzAz +Nz^2 = 0

+AxPx -AyPy +NyPy +NyAy -Ny^2 - AzPz +NzPz +NzAz -Nz^2 -AxPx +AyPy -NyPy -NyAy +Ny^2 +AzPz -NzPz -NzAz +Nz^2 = NxPx +NxAx +NxPx +NxAx

2(Px +Ax)Nx = +AxPx -AyPy +NyPy +NyAy -Ny^2 - AzPz +NzPz +NzAz -Nz^2 -AxPx +AyPy -NyPy -NyAy +Ny^2 +AzPz -NzPz -NzAz +Nz^2

Nx = [(+AxPx -AyPy +NyPy +NyAy -Ny^2 - AzPz +NzPz +NzAz -Nz^2 -AxPx +AyPy -NyPy -NyAy +Ny^2 +AzPz -NzPz -NzAz +Nz^2)/2(Px+Ax)]

BxPx-[(+AxPx -AyPy +NyPy +NyAy -Ny^2 - AzPz +NzPz +NzAz -Nz^2 -AxPx +AyPy -NyPy -NyAy +Ny^2 +AzPz -NzPz -NzAz +Nz^2)/2(Px+Ax)]Px-[(+AxPx -AyPy +NyPy +NyAy -Ny^2 - AzPz +NzPz +NzAz -Nz^2 -AxPx +AyPy -NyPy -NyAy +Ny^2 +AzPz -NzPz -NzAz +Nz^2)/2(Px+Ax)]Bx+[(+AxPx -AyPy +NyPy +NyAy -Ny^2 - AzPz +NzPz +NzAz -Nz^2 -AxPx +AyPy -NyPy -NyAy +Ny^2 +AzPz -NzPz -NzAz +Nz^2)/2(Px+Ax)]^2 + ByPy-NyPy-NyBy+Ny^2 + BzPz-NzPz-NzBz+Nz^2 = 0

Midpoint Of A Line

Line with endpoint vectors A and B:

midpoint= ((A - B) * .5) + B

Intersection Of Line And Plane

• N=(Nx,Ny,Nz) normal vector of plane

• A=(Ax,Ay,Az) one endpoint of the line in question

• B=(Bx,By,Bz) other endpoint of the line in question

• P=(Px,Py,Pz) the point of intersection on plane along line

• Px=(Bx-Ax)alpha+Ax alpha=ratio of P along line AB

• alpha=(Px-Ax)/(Bx-Ax)=(Py-Ay)/(By-Ay)=(Pz-Az)/(Bz-Az)

Three Equations For Three Variables Px,Py,Pz

1. NxPx + NyPy + NzPz = 0   or N dot P=0
2. Py=(Px-Ax)(By-Ay)/(Bx-Ax) + Ay
3. Pz=(Px-Ax)(Bz-Az)/(Bx-Ax) + Az

0=NxPx+Ny((Px-Ax)(By-Ay)/(Bx-Ax)+Ay)+Nz((Px-Ax)(Bz-Az)/(Bx-Ax)+Az)

<after much fussy and punctilious rearranging>

Ny(AyBx-AxBy) + Nz(AzBx-AxBz)
Px = ---------------------------------
Nx(Ax-Bx) + Ny(Ay-By) + Nz(Az-Bz)

Intersection Of Two 2-d Line Segments

A nice website that calculates intersectoins.

Line 1 is (Ax,Ay) to (Bx,By)
Line 2 is (Cx,Cy) to (Dx,Dy)
y1 = m1*x1+b1
y2 = m2*x2+b2

Where the lines intersect, x1=x2 and y1=y2, so:

m1*x+b1 = m2*x+b2
(m1-m2)*x = b2-b1
x = (b2-b1)/(m1-m2)

Ah, but what is intercept b and slope m?

m1 = (By-Ay)/(Bx-Ax)
m2 = (Dy-Cy)/(Dx-Cx)

And now the intercept b?

Ay = m1*Ax + b1
Cy = m2*Cx + b2
b1 = Ay - (m1*Ax)
b2 = Cy - (m2*Cx)
b1 = Ay - ((By-Ay)/(Bx-Ax)*Ax)
b2 = Cy - ((Dy-Cy)/(Dx-Cx)*Cx)

Recall:

x = (b2-b1)/(m1-m2)
y = m1*x + b1   (or (y = m2*x + b2) at intersection)

Fill in b, m and the intersection point is x,y:

x = {[Cy-((Dy-Cy)/(Dx-Cx)*Cx)]-[Ay-((By-Ay)/(Bx-Ax)*Ax)]}  /
{[(By-Ay)/(Bx-Ax)]-[(Dy-Cy)/(Dx-Cx)]}
y = [(By-Ay)/(Bx-Ax)]*x + [Ay - ((By-Ay)/(Bx-Ax)*Ax)]

Of course this should be set up to catch special cases like parallel lines.

Intersection Of A Plane And A Z Vector At Arbitrary X,Y

 Note UNTESTED

This is a special case of the line/plane intersection. It is a simplified case with a particular use in Hidden Line Removal (HLR). When the view is lined up along the Z axis, and a point’s X,Y is within an opaque triangle’s X,Y borders, then it is important to know whether the point is nearer to the viewer than the opaque triange (visable) or farther (hidden). The Z, or nearness to viewer, of the point is known. How near the point on the triangle is at X,Y is not automatically known; the triangle could be quite oblique to the view.

Setup:

Coordinates of the triangular face:

(T1x,T1y,T1z) (T2x,T2yT2z) (T3x,T3y,T3z)

Coordinates of the intersecting point:

(Ax,Ay,Az)  Ax and Ay are known constants, Az is the variable

First move the whole mess so that a triangle vertex is 0,0,0:

T2x -= T1x; T2y -= T1y; T2z -= T1z (i.e. T2x=T2x-T1x, etc.)
T3x -= T1x; T3y -= T1y; T3z -= T1z
Ax -= T1x; Ay -= T1y

T1 now is effectively (0,0,0)

The triangle needs a useful handle in the form of a normal vector at a vertex (the new origin vertex T1 is a good choice):

N = T2 crossproduct T3
Nx= T2y*T3z - T2z*T3y
Ny= T2x*T3z - T2z*T3x
Nz= T2x*T3y - T2y*T3x

If N is normal to the tri and A is on the tri, then they must be perpendicular and satisfy this equation:

A dotproduct N = 0
And that's it. One unknown, one equation. Onward:
Nx*Ax + Ny*Ay + Nz*Az = 0
Nz*Az = -Nx*Ax - Ny*Ay
Az = (-Nx*Ax - Ny*Ay)/Nz

Now put the thing back where it was (only Z is necessary):

Az += T1z;

So, the complete deal with everything included should be:

Az = {[-[(T2y-T1y)*(T3z-T1z) - (T2z-T1z)*(T3y-T1y)]*(Ax-T1x) -
[(T2x-T1x)*(T3z-T1z) - (T2z-T1z)*(T3x-T1x)]*(Ay-T1y)]  /
[(T2x-T1x)*(T3y-T1y) - (T2y-T1y)*(T3x-T1x)]} + T1z

Intersection Of A Line And An Intersecting Z Vector

Calculates Z for a point of known X,Y on a 3-d line with known endpoints. The goal of this is to determine which of 2 lines is in front of the other. If two lines are found to intersect in x-y dimensions (i.e. a Z aligned view) at X,Y, then what are the points for Z for each of the lines. This approach looks at each line one at a time. So there is the x-y screen intersection point, (Px,Py), and a line to find the depth of contact of, (Ax,Ay)→(Bx,By). The object is to find Pz for this line. Then find the Z at this point for the other line and compare.

There seem to be two simple approaches to this - one slightly simpler than the other. Since AB is a line, things are linear. This means that the Z value will change with respect to X or Y. So one approach is to drop X or Y since only one is needed to relate to the Z. Here’s that approach:

Ignoring Y              Ignoring X
Ax-Bx/Az-Bz=Ax-Px/Az-Pz     Ay-By/Az-Bz=Ay-Py/Az-Pz
Az-Pz=(Ax-Px)*(Az-Bz)/Ax-Bx Az-Pz=(Ay-Py)*(Az-Bz)/Ay-By
Pz=Az-(Ax-Px)*(Az-Bz)/Ax-Bx Pz=Az-(Ay-Py)*(Az-Bz)/Ay-By

The problem with this approach is that if the Y is ignored, and the line is mostly along the Y axis, then the chance of precision errors gets big. Of course if the line is completely on the Y axis, there is a division by zero problem. One approach is to wisely choose which of the above formulas to use by this heuristic:

if (Ax-Bx) > (Ay-By) then (Ignore Y) else (Ignore X)

The grand complete solution is to not use progress along the line in some specific axis, but rather the 2-d distance from the starting point as the independent variable. This means that when the point has travelled halfway along the line in XY, it has travelled half of the way in Z. This is opposed to saying that when it has travelled half of the X displacement, then it has also moved half of the Z. This may be true, but the X displacement might be 0. Here’s the full formula:

sqrt[(Ax-Bx)^2+(Ay-By)^2]/(Az-Bz) = sqrt[(Ax-Px)^2+(Ay-Py)^2]/(Az-Pz)
Pz=Az-{sqrt[(Ax-Px)^2+(Ay-Py)^2]*(Az-Bz)/sqrt[(Ax-Bx)^2+(Ay-By)^2]}

Yuck…

Linear Algebra and Matrix Operations

Matrix Multiplication
/* Input: */
/* A= matrix of dimension x by y. */
/* B= matrix of dimension y by z. */
/* Output: */
/* C= matrix of dimension x by z
/* where C[i][j] is dot product of ith row of A and jth column of B */

for (i=1;i<=x;i++)
for (j=1;j<=z;j++) {
C[i][j]= 0;
for (k=1;k<=y;k++)
C[i][j]+= A[i][k] * B[k][j];
}

Here is a small comprehensible 2d matrix library from Sebastian Thrun (and Ernesto P. Adorio for the Cholesky functions). I find this far easier to understand than the obfuscatory typography gloop found in math books. Not only does this clearly explain how matrix operations work in unambiguous terms, you can actually use this library directly to do useful things.

#!/usr/bin/python3
from math import *

class matrix:

def __init__(self, value):
self.value = value
self.dimx = len(value)
self.dimy = len(value)
if value == [[]]:
self.dimx = 0

def zero(self, dimx, dimy):
assert dimx >= 1 or dimy >= 1, "Invalid size of matrix"
self.dimx = dimx
self.dimy = dimy
self.value = [[0 for row in range(dimy)] for col in range(dimx)]

def identity(self, dim):
assert dim >= 1, "Invalid size of matrix"
self.dimx = dim
self.dimy = dim
self.value = [[0 for row in range(dim)] for col in range(dim)]
for i in range(dim):
self.value[i][i] = 1

def show(self):
for i in range(self.dimx):
print(self.value[i])
print(' ')

assert self.dimx == other.dimx or self.dimy == other.dimy, \
"Matrices must be of equal dimensions to add"
res = matrix([[]])
res.zero(self.dimx, self.dimy)
for i in range(self.dimx):
for j in range(self.dimy):
res.value[i][j] = self.value[i][j] + other.value[i][j]
return res

def __sub__(self, other):
assert self.dimx == other.dimx or self.dimy == other.dimy, \
"Matrices must be of equal dimensions to subtract"
res = matrix([[]])
res.zero(self.dimx, self.dimy)
for i in range(self.dimx):
for j in range(self.dimy):
res.value[i][j] = self.value[i][j] - other.value[i][j]
return res

def __mul__(self, other):
assert self.dimy == other.dimx, \
"Matrices must be m*n and n*p to multiply"
res = matrix([[]])
res.zero(self.dimx, other.dimy)
for i in range(self.dimx):
for j in range(other.dimy):
for k in range(self.dimy):
res.value[i][j] += self.value[i][k] * other.value[k][j]
return res

def transpose(self):
res = matrix([[]])
res.zero(self.dimy, self.dimx)
for i in range(self.dimx):
for j in range(self.dimy):
res.value[j][i] = self.value[i][j]
return res

def Cholesky(self, ztol=1.0e-5):
# Computes the upper triangular Cholesky factorization of
# a positive definite matrix.
res = matrix([[]])
res.zero(self.dimx, self.dimx)
for i in range(self.dimx):
S = sum([(res.value[k][i])**2 for k in range(i)])
d = self.value[i][i] - S
if abs(d) < ztol:
res.value[i][i] = 0.0
else:
assert d >= 0.0, "Matrix not positive-definite"
res.value[i][i] = sqrt(d)
for j in range(i+1, self.dimx):
S = sum([res.value[k][i] * res.value[k][j] for k in range(self.dimx)])
if abs(S) < ztol:
S = 0.0
res.value[i][j] = (self.value[i][j] - S)/res.value[i][i]
return res

def CholeskyInverse(self):
# Computes inverse of matrix given its Cholesky upper Triangular
# decomposition of matrix.
res = matrix([[]])
res.zero(self.dimx, self.dimx)
# Backward step for inverse.
for j in reversed(range(self.dimx)):
tjj = self.value[j][j]
S = sum([self.value[j][k]*res.value[j][k] for k in range(j+1, self.dimx)])
res.value[j][j] = 1.0/tjj**2 - S/tjj
for i in reversed(range(j)):
res.value[j][i] = res.value[i][j] = -sum([self.value[i][k]*res.value[k][j] for k in range(i+1, self.dimx)])/self.value[i][i]
return res

def inverse(self):
aux = self.Cholesky()
res = aux.CholeskyInverse()
return res

def __repr__(self):
return repr(self.value)