0
$\begingroup$

Related to this question: https://math.stackexchange.com/questions/407897/randomly-generate-point-on-shell-from-3-points-2-angles-with-uniform-angle-dis

I'm trying to reverse engineer the math-algorithms behind a complex vector geometry-generating code.

Can you please help me:

  1. Figure out the fundamental math operations used in this generation, based on the sample code below.
  2. Suggest why this (more complex) approach is used (better distribution? less trig ops.)?

The code generates points on the rim of an uniformally rotated (as in the two rotation angles phi and theta are drawn from flat distributions) of arbitrary length (vectLen) and with an arbitrary tip/attachment point (B) in 3D space.
(courtesy of Wolfram alpha)

It appears to be doing this via some sort of projection (see code regarding vector V) that uses the quadratic formula. This may be via the use of a conic:
Conic section

(Note, I translated this code from fortran, have not tested it yet.)

W.x = A.x - B.x;
W.y = A.y - B.y;
W.z = A.z - B.z;
W.unitize();

do
{
   V.x = rand(-1.0,1.0);
   V.y = rand(-1.0,1.0);
} while (V.len()>1.0);

a = W.x * W.x + W.y * W.y - 1.0;
b = -2.0 * ( V.x * W.x * W.z + V.y * W.y * W.z );
c = V.y * V.y * W.x * W.x + V.x * V.x * W.y * W.y +
   V.x * V.x * W.z * W.z + V.y * V.y * W.z * W.z -
   2.0 * V.x * V.y * W.x * W.y -  V.x * V.x - V.y * V.y;
d = b * b - 4.0 * a * c;

if ( d < 0.0 && d > -0.000001 )
   V.z = -b / ( 2.0*a );
else if ( rand(0,1.0) < 0.5 )
   V.z = ( -b + sqrt(d) ) / ( 2.0*a );
else
   V.z = ( -b - sqrt(d) ) / ( 2.0*a );

V.unitize();

for (unsigned int i = 0, i < numToDraw; i++)
{
   theta = rand(0, 2*M_PI);
   phi = rand(0, 2*M_PI);
   h = 1.0 / tan( phi );

   cs = cos(phi);
   s = sin(phi);

   U.x = cs * V.x + s * W.z * V.y - s * W.y * V.z;
   U.y = cs * V.y + s * W.x * V.z - s * W.z * V.x;
   U.z = cs * V.z + s * W.y * V.x - s * W.x * V.y;

   T.x = h * W.x + U.x;
   T.y = h * W.y + U.y;
   T.z = h * W.z + U.z;

   T.unitize();

   C[i].x = B.x + vectLen * T.x;
   C[i].y = B.y + vectLen * T.y;
   C[i].z = B.z + vectLen * T.z;
}

With respect to question #2, my guesses are as follows... It seems like it cuts down on the # of trig. function calls as there are only 3 used here, versus multiple ones for a standard rotation-matrix based solution.

It appears that the projection may also provide a more uniform distribution of coordinates while maintaining the angular distribution of the cones rotation about its two degrees of freedom (versus just drawing phi and theta and churning them through a series of matrix ops). Maybe I'm mistaken on that point, though... what do you think?

$\endgroup$
  • 1
    $\begingroup$ What does unitize () do? I would have thought that it normalizes a vector to unit length, but if so, a would always be zero? $\endgroup$ – joriki Jun 3 '13 at 17:57
  • $\begingroup$ Your guess regarding the function's functionality (as the name implies) is correct, however you assertion regarding a being zero is flawed. For a unit vector in 3D x^2+y^2+z^2=1.0, hence x^2+y^2-1.0=-(z^2). So it's the reflection across the xy plane of the z component of the unit form of W.... Perhaps you missed that it was using W and thought it was using V (in which case there was no z set yet, and your assertion would be correct assuming initialization to zero). $\endgroup$ – Jason R. Mick Jun 3 '13 at 19:02
  • 1
    $\begingroup$ Sorry, yes, I was distracted by V. $\endgroup$ – joriki Jun 3 '13 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.