CODEWITHSUNDEEP
javascriptmath3dtrigonometry
Michael Schade
Michael Schade

Reputation: 418

How to calculate a circle perpendicular to a vector?

Given two points (P1 and P2) in XYZ space, create a tube with a given radius. In order to do this, I need to calculate points for a circle around each of the two points, such that the circles are perpendicular to P1→P2 (and parallel to each other). The dx/dy/dz for one circle can be used to make other circles. The form of the code would look like:

function circle(radius, segments, P1, P2) {
  // 3D circle around the origin, perpendicular to P1P2
  var circle = [];  
  var Q = [P2[0] - P1[0], P2[1] - P1[1], P2[2] - P1[2]];
  for (var i = 0; i < segments; i++) {
    var theta = 2*Math.PI*segment/i;
    var dx = mysteryFunctionX(Q, theta, radius);
    var dy = mysteryFunctionY(Q, theta, radius);
    var dz = mysteryFunctionZ(Q, theta, radius);
    circle.push([dx, dy, dz]);
    }
  return circle;
  }

What is the calculation needed for each mystery function?

Upvotes: 3

Views: 2185

Answers (2)

Michael Schade
Michael Schade

Reputation: 418

Thank you Forward Ed and dmuir - that helped. Here is the code I made that seems to work:

function addTube(radius, segments, P1, P2) {
  // Q = P1→P2 moved to origin
  var Qx = P2[0] - P1[0];
  var Qy = P2[1] - P1[1];
  var Qz = P2[2] - P1[2];
  // Create vectors U and V that are (1) mutually perpendicular and (2) perpendicular to Q
  if (Qx != 0) {  // create a perpendicular vector on the XY plane
    // there are an infinite number of potential vectors; arbitrarily select y = 1
    var Ux = -Qy/Qx;
    var Uy = 1;
    var Uz = 0;
    // to prove U is perpendicular:
    // (Qx, Qy, Qz)·(Ux, Uy, Uz) = Qx·Ux + Qy·Uy + Qz·Uz = Qx·-Qy/Qx + Qy·1 + Qz·0 = -Qy + Qy + 0 = 0
    }
  else if (Qy != 0) {  // create a perpendicular vector on the YZ plane
    var Ux = 0;
    var Uy = -Qz/Qy;
    var Uz = 1;
    }
  else {  // assume Qz != 0; create a perpendicular vector on the XZ plane
    var Ux = 1;
    var Uy = 0;
    var Uz = -Qx/Qz;
    }
  // The cross product of two vectors is perpendicular to both, so to find V:
  // (Vx, Vy, Vz) = (Qx, Qy, Qz)×(Ux, Uy, Uz) = (Qy×Uz - Qz×Uy, Qz×Ux - Qx×Uz, Qx×Uy - Qy×Ux)
  var Vx = Qy*Uz - Qz*Uy;
  var Vy = Qz*Ux - Qx*Uz;
  var Vz = Qx*Uy - Qy*Ux;
  // normalize U and V:
  var Ulength = Math.sqrt(Math.pow(Ux, 2) + Math.pow(Uy, 2) + Math.pow(Uz, 2));
  var Vlength = Math.sqrt(Math.pow(Vx, 2) + Math.pow(Vy, 2) + Math.pow(Vz, 2));
  Ux /= Ulength;
  Uy /= Ulength;
  Uz /= Ulength;
  Vx /= Vlength;
  Vy /= Vlength;
  Vz /= Vlength;
  for (var i = 0; i < segments; i++) {
    var θ = 2*Math.PI*i/segments;  // theta
    var dx = radius*(Math.cos(θ)*Ux + Math.sin(θ)*Vx);
    var dy = radius*(Math.cos(θ)*Uy + Math.sin(θ)*Vy);
    var dz = radius*(Math.cos(θ)*Uz + Math.sin(θ)*Vz);
    drawLine(P1[0] + dx, P1[1] + dy, P1[2] + dz,  // point on circle around P1
             P2[0] + dx, P2[1] + dy, P2[2] + dz)  // point on circle around P2 
    }
  }

I'm sure there are many ways to shorten the code and make it more efficient. I created a short visual demo online using Three.JS, at http://mvjantzen.com/tools/webgl/cylinder.html

Upvotes: 1

dmuir
dmuir

Reputation: 4431

As pointed out in the link in Ed's post, if you have vectors u and v that are perpendicular to your axis Q, and to each other, and each of length 1 then the points 

 P + cos(theta)*u + sin(theta)*v

are, as theta goes between 0 and 2pi, the points on a circle with centre P on a plane perpendicular to Q.

It can be a bit tricky, given Q, to figure out what u and v should be. One way is to use Householder reflectors. It is straightforward to find a reflector that maps (1,0,0) say to a multiple of Q. If we apply this reflector to (0,1,0) and (0,0,1) we will get vectors u and v as required above. The algebra is a little tedious but the following C code does the job:

static  void    make_basis( const double* Q, double* u, double* v)
{
double  L = hypot( Q[0], hypot( Q[1], Q[2])); // length of Q
double  sigma = (Q[0]>0.0) ? L : -L;    // copysign( l, Q[0]) if you have it
double  h = Q[0] + sigma;   // first component of householder vector
double  beta = -1.0/(sigma*h);  // householder scale
    // apply to (0,1,0)'
double  f = beta*Q[1];
    u[0] = f*h;
    u[1] = 1.0+f*Q[1];
    u[2] = f*Q[2];
    // apply to (0,0,1)'
double  g = beta*Q[2];
    v[0] = g*h;
    v[1] = g*Q[1];
    v[2] = 1.0+g*Q[2];
}

Upvotes: 1

Related Questions