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

Resources

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 add(Vector offset) {
    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 rotate(float &rads) {
    Vector rotated;
    float s, c;
    s=sin(rads); c=cos(rads);
    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);
    s=sin(rads); c=cos(rads);
    y2= rotated.y*c + rotated.x*s;
    rotated.x = rotated.x*c - rotated.y*s;
    rotated.y = y2;
    rotated.set2_add(basept);
    return rotated;              }

/*  rotation in XZ axis */
Vector tilt(float &rads) {
    Vector rotated;
    float s, c;
    s=sin(rads); c=cos(rads);
    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);
    s=sin(rads); c=cos(rads);
    z2= rotated.z*c + rotated.x*s;
    rotated.x = rotated.x*c - rotated.z*s;
    rotated.z = z2;
    rotated.set2_add(basept);
    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;
    scaled.set2_add(basept);
    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;
    scaled.set2_add(basept);
    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

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…

Graph Threory

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];
    }