The Power of the Vector Dot Product and Cross Product in (ℝ^3)

John McGee
Wolfram Technical

Join the Conversation #WolframTechConf

Initialization

eps=10^(-10);  (* the zero thereshold *)

setZeroThreshold[neweps_:10^(-10)]:=If[neweps>0&&neweps<1,eps=neweps]

(* determine the length of the projection of u onto v *)
scalarProject[u_,v_]:=With[{nv=Norm[v]},If[nv<eps,0,u.v/nv]]

(* determine the length of the projection of u onto the unit vector v *)
scalarProjectUnit[u_,vh_]:=u.vh

(* distance of p from the line's anchor point q line q + k v *)
distAlongLine[p_,{q_,v_}]:=(p-q).v

(* order the points pnts with respect to thier paramter k on the line q + k v *)
orderPointsOnLine[pnts_, {q_, v_}] := Ordering[((# - q).v) & /@ pnts]

(* give the vector projection of u onto v, the component of u in the direction of v *)
vectorProject[u_,v_]:=((u.v)/(v.v)) v

(* determine the vector projection of u onto the unit vector vh *)
vectorProjectUnit[u_,vh_]:=(u.vh) vh

(* project the point p onto the parametric line q + k v*)
projectPointLine[p_,{q_,vh_}]:=q+vectorProjectUnit[p-q,vh]

(* compute the perpendicular distance between the point p and the parametric line p + q v *)
distancePointLine[p_,{q_,vh_}]:=Sqrt[(p-q).(p-q)-((p-q).vh)^2]

distancePointPlane[p_,{q_,vh_}]:=(p-q).vh

projectPointPlane[p_,{q_,vh_}]:=p-distancePointPlane[p,{q,vh}] vh

intersectLinePlane[{p_,uh_},{q_,vh_}]:=With[{k=((q-p).vh)/(vh.uh)},p+k uh]

intersectPlanePlane[{q1_,vh1_},{q2_,vh2_}]:=Module[{w,nw,uh,p},
  (* the line of intersection is perpendicular to both plane normals *)
  w=Cross[vh1,vh2];nw=Norm[w];
  (* if cross product vanishes the planes are parallel, no line of intersection *)
  If[nw<eps,Return[{}]];
  (* create a direction from q2 directly to the intersection point *)
  uh=Normalize[Cross[vh2,w]];
  p=intersectLinePlane[{q2,uh},{q1,vh1}];
 {p,w}
]

intersectPlaneCylinder[{q_,vh_},{p_,uh_,r_}]:=Module[{cp,wh,zh,a,b},
  cp=intersectLinePlane[{p,uh},{q,vh}];
  wh=Normalize[Cross[vh,uh]]; (* the minor axis direction *)
  zh=Cross[vh,wh]; (* the major axis direction *)
  b=r;a=r/(vh.uh); (* the semi-major and semi-minor axis lengths *)
  {cp,zh,wh,a,b}
]

Join the Conversation #WolframTechConf

Outline

● Definitions and Properties of the dot product (OverVector[u])·(OverVector[v]) in (ℝ^n)

● Special behavior of (ℝ^3), including the cross product (OverVector[u])×(OverVector[v])

● Basic geometric operations based on the dot project

	vector and scalar projection, distance computations
	
● Geometric constructions using the dot and cross product

● Applications of these basic methods

Vectors in (ℝ^n)

A vector in (ℝ^n) is an ordered n-tuple v=<(Subscript[v, 1]),(Subscript[v, 2]),...,(Subscript[v, n])>, with addition and scalar multiplication defined component-wise as:

	u+v=<(Subscript[u, 1])+(Subscript[v, 1]), (Subscript[u, 2])+(Subscript[v, 2]),...,(Subscript[u, n])+(Subscript[v, n])>
	
	s u=<s (Subscript[u, 1]),s (Subscript[u, 2]),...,s (Subscript[u, n])>

With these definitions (ℝ^n) becomes a vector space, e.g. an Abelian group under addition, with scalar multiplication as a field acting on the group.  In particular we get two distributive properties

	s (u+v)=s u+s v	(s+r) v = s v + r v	
	
With these operations we get a algebra of vectors.

The Geometric Vector Dot Product

The dot product of two vectors u,v in (ℝ^n) is proportional to the product of the length (norm) of the vectors with the cosine of the angle between them.

	u·v=||u||×||v||×cos(θ)
	
in particular for unit vectors (OverHat[u]),(OverHat[v]), e.g.||(OverHat[u])||=||(OverHat[v])||=1,

	(OverHat[u])·(OverHat[v])=cos(θ)
	
This is a coordinate-free representation of the dot product

	note that u·v=0 when cos(θ)=0, e.g. when θ=π/2,-π/2
	
	also (OverHat[u])·(OverHat[v])=1 when θ=0 and (OverHat[u])·(OverHat[v])=-1 when θ=π

The Algebraic Vector Dot Product

The dot product of two vectors u,v in (ℝ^n) is given by

	u·v=(∑_(i=1)^n)(Subscript[u, i]) (Subscript[v, i])
	
in particular for vectors u=<(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])>,v=<(Subscript[x, 2]),(Subscript[y, 2]),(Subscript[v, 2])>

	u·v=(Subscript[x, 1 ])(Subscript[x, 2])+(Subscript[y, 1])(Subscript[y, 2])+(Subscript[z, 1]) (Subscript[z, 2])
	
With the definition ||u((||)^2)=u·u, and d(u,v)=||u-v||=(Sqrt[(u-v)·(u-v)]),
	  (ℝ^n) becomes a metric space

This means that we can determine the angle between two vectors as

	cos(θ)=((u·v)/(||u|| ×||v||));0<θ<π
	
where we can use the coordinates of u,v to perform the computations.

Scalar and Vector Projection

The scalar projection of the vector u onto the vector v gives the length of the component of u that lies in the direction of v.

	a=u·v/||v|| = u·(OverHat[v])
	
The vector projection of the vector u onto the vector v gives the component of u that lies in the direction of v.

	w=a (OverHat[v])= (u·(v/(||v|| )))(OverHat[v])= ((u·v)/(v·v))v
	
Note: in this form we do not need to compute the norm of v

The Vector Cross Product in (ℝ^3)

For vectors u,v∈(ℝ^3), the vector w=u×v is a vector perpendicular to both u and v with length

	||w||=||u||×||v||×sin(θ) ,
	
		 where θ is the angle between u and v
	 
	 note that w is the zero vector when sin(θ)=0
	 
for vectors u=<(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])>,v=<(Subscript[x, 2]),(Subscript[y, 2]),(Subscript[v, 2])>

	w=<(Subscript[y, 1]) (Subscript[z, 2])-(Subscript[y, 2]) (Subscript[z, 1]),(Subscript[x, 2]) (Subscript[z, 1])-(Subscript[x, 1]) (Subscript[z, 2]),(Subscript[x, 1]) (Subscript[y, 2])-(Subscript[x, 2]) (Subscript[y, 1])>

“Proof” that w⊥u,v

u={(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[z, 1])};v={(Subscript[x, 2]),(Subscript[y, 2]),(Subscript[z, 2])};
w=u\[Cross]v

{w.u,w.v}//Simplify (* u\[Cross]v is ⊥ to both u and v *)

Basic Applications

Project a point p onto a line with parametric form q+k (OverHat[v])

	r= q+a (OverHat[v])   where a is the scalar projection of (p-q) onto (OverHat[v])
				
	
	
Distance between a point p and a line with parametric form q+k (OverHat[v])

	d=(Sqrt[||p-q((||)^2)-(a^2)])  where a is the scalar projection of (p-q) onto (OverHat[v])

Plane Applications

Distance between point p, and a plane passing through point q with normal vector (OverHat[v])

	the plane equation is (OverHat[v])·x=(OverHat[v])·q=d
	
	note: if (OverHat[v])=<a,b,c>, x=<x,y,z> this becomes a x+b y+c z =d

	h=(p-q)·(OverHat[v]), i.e. the scalar projection of (p-q) onto (OverHat[v])
	
	distancePointPlane[p_, {q_, vh_}] := (p - q).vh
		
Project a point p onto a plane passing through point q with normal vector (OverHat[v])

	w=p-h (OverHat[v])
	
	projectPointPlane[p_,{q_,vh_}] := p - ((p - q).vh) vh
	
Intersect a line p+k (OverHat[u]) with plane (OverHat[v])·x=(OverHat[v])·q

	intersectLinePlane[{p_,uh_},{q_,vh_}] :=
	   With[{k = ((q - p).vh)/(vh.uh)}, p + k uh]

Intersect Plane with Cylinder

The center of the ellipse is at the intersection of the cylinder axis with the plane
The minor axis direction is perpendicular to both the cylinder axis and the plane normal
The major axis direction is perpendicular to both the plane axis and the minor axis
The semi-minor axis length b is the same as the cylinder radius
The semi-major axis length a=a/cos(θ)=a/((OverHat[v])·(OverHat[u]))

intersectPlaneCylinder[{q_,vh_},{p_,uh_,r_}]:=Module[{cp,wh,zh,a,b},
   cp = intersectLinePlane[{p,uh},{q,vh}];   (* the center point *)
   wh = Normalize[Cross[vh,uh]];   (* the minor axis direction *)
   zh = Cross[vh ,wh];     (* the major axis direction *)
   b = r;  a = r/(vh.uh);    (* the semi-major and semi-minor axis lengths *)
   {cp,zh,wh,a,b}
]
      Ellipse {a Cos(θ) (Subscript[v, x])+b Sin(θ) (Subscript[v, y])|0≤θ<2 π}

Block Alignment Problem

Connect Block R, Cube 0, Face 4 to Block Y cube 0 face 3 with W rotation

(1) Label the blocks with face vectors

(2) Select from the cube group (S3) the four matrices M that satisfy
	(M·(Subscript[v, 2])),·(Subscript[v, 1])=-1

rotmatset = Select[allmats, (#.face2vec).face1vec == -1 &];

	choose the matrix that aligns the two faces North to North
	
		rotmat = Select[rotmatset, ((#.v12).v11==1)& ];
	
	rotate the second block about its mating face vector by the given compass point

		rotangle=orient/.{“N”->0,”W”->π/2,”S”->π,”E”->-π/2};
		
		orientmat=RotationMatrix[rotangle, fv22];
	
	Connect Block R, Cube 0, Face 4 to Block Y cube 0 face 3 with W rotation

(3) Translate the faces into alignment

Application - Intersect Plane Mesh

Efficiently find the loops of intersection between the plane and the mesh

(1) Select only the polygons that cross the plane using DistancePointPlane on the vertices

d=(p-q)·v;   d>0 when p is above the plane d<0 when p is below

	Select[vertices, SomeTrue[Sign[(q-p).v==1]]

(2) Intersect this (small) set of polygons with the plane using IntersectLinePlane
	and order the resulting line segments into loops

(3) Show the result

(4) Similar approach gives intersect cylinder - mesh

Simulation of an Interferometer

Simulate a Fabry-Pérot Interferometer by modeling reflection and refraction

 
 
(*  vi is the incident direction
   vn is the direction normal to the interface surface
   ni is the refractive index of the incident material
   nr is the refractive index of the refractive material
   returns the refracted vector *)
 refract[vi_ ,vn_, ni_, nr_]:=Module[{th1, th2},
	th1 = angleBetween[vi, vn];
	th2 = ArcSin[ni Sin[th1]/nr];
	RotationMatrix[th2].(-vn)
]

(* reflect the vector vi at the surface with normal vector vn *)
reflect[vi_ ,vn_] := ReflectionMatrix[vn].vi

A Development Method Note

To optimize a computational function that calls one or more lower level functions:  e.g. k(x,y)=f(g(x.y)+h(x×y))

(1) call the outer function with symbolic vectors and scalars

(2) apply FullSimplify to the symbolic result

(3) re-define the outer function using the simplified result

(4) Compile and/or parallelize

Conclusions and Future Work

The Wolfram Language has a rich set of 3D geometry functions
	
Simple functions based on Dot and Cross can have better performance for specialized applications than the more general

	RegionIntersection, LinearSolve 
	
Continue the development of the vector geometry package
 * optional use of ParallelMap, etc. for large sets of points/vectors
 * integrate more functions from classical geometry