Description
I am working on creating a wireframe sphere that showcases latitude and longitude lines. A visual representation of what I'm attempting to achieve can be found here. To be clear, I'm looking for an orthographic projection, not a perspective projection.
To define a couple things, I'm saying $\theta$ is the rotation from the $x$-axis to the $y$-axis, $\phi$ is the rotation from the $z$-axis to the $y$-axis, and finally $\psi$ is the rotation from the $x$-axis to the $z$-axis.
My question
TL;DR my question is this: What is the $\mathbb R^2$ equation of a circle "in" $\mathbb R^3$ that is rotated $(\theta,\phi,\psi)$ about the origin and has offset $d$ that is the distance from the center of the circle to the origin that is projected onto the $xz$-axis?
Constraints
I am making this using a computer-generated shader. Meaning, in order to draw said object, each pixel will have to be calculated to "check" if it is apart of the circle once projected. For this reason, I am constrained to non-paraterized functions. (I.E. my check would be, as a simple case, does pix_location_x * pix_location_x + pix_location_y * pix_location_y == 1?)
Initial Attempts
Idea #1, using ellipses
I've began by breaking down this problem into separate smaller parts. I first tried drawing rotations on $(\theta,\phi,\psi)$ against the $(x,y,z)$ axis using ellipse equations; since ellipses appear to be rotated circles.
The reason I want to do this is because I can then combine multiple "circles" with offset $\theta$ and offset $center$ to mimic latitude and longitude lines.
My first attempt on the first part (the latitude lines) has fallen short; (I have come up with the following desmos project (The $x$-axis is the 3D $x$-axis, while the $y$-axis is the 3D $z$-axis).)
I attempted to use the following equation:
$$\frac{(x\cos(\psi) – y\sin(\psi))^2}{(\frac{\mod(2\theta, 360)}{360}-0.5)^2}+\frac{(x\sin(\psi) – y\cos(\psi))^2}{(\frac{\mod(2\phi, 360)}{360}-0.5)^2}=1$$
This equation seems to successfully rotate about $\theta$, $\phi$, and $\psi$ only when they are the only rotations taking place. It seems that each is trying to rotate about a fixed $(x,y,z)$ axis system where as $(\theta,\phi,\psi)$ actually changes the direction of that axis system, and so certain combinations seem odd.
Idea #2, using projections and intersections
I havent't actually gone ahead and done this yet, but my next idea was to use the standard sphere:
$$x^2+y^2+z^2=1$$ and intersect a plane rotated about $(\theta,\phi,\psi)$ on the origin. Then take the intersection of the two surfaces, and project that curve onto the $xz$-axis for my drawing.
This however seems overboard, unless I can simplify the expression before hand and get some closed form solution, but I think at that point I am solving my ellipse problem from above.
Any help or pointers would be much appreciated!
Best Answer
Working in homogeneous coordinates can be convenient for problems like these. Planes in $\mathbb R^3$ (really the projective space $\mathbb{RP}^3$) can be represented by vectors, while quadric surfaces and conics can be represented by symmetric matrices. All of the various transformations involved are linear, so the necessary computations can usually be constructed as a sequence of matrix multiplications and the occasional inversion (or pseudo-inverse).
Going with your second approach, we intersect a plane with the unit sphere and then project this intersection orthogonally onto the $x$-$z$ plane. In homogeneous coordinates, the unit sphere has the matrix equation $\mathbf X^T\mathtt Q\mathbf X=0$, where $\mathtt Q=\operatorname{diag}(1,1,1,-1)$. The final projection matrix is quite simple, too. It’s just $$\mathtt P = \begin{bmatrix}1&0&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix},$$ which simply discards the $y$-coordinate of a point. I would probably negate the first column so that we’re looking down the positive $y$-axis toward the image plane instead of viewing it “from behind,” but that’s not a critical detail. The outline of the sphere’s projection has the matrix $\mathtt P\mathtt Q\mathtt P^T = \operatorname{diag}(1,1,-1)$, as one would expect, but that’s not quite what we’re after.†
It looks like you know how to compute a normal $\tilde{\mathbf N} = (n_x,n_y,n_z)^T$ and offset $d$ for the intersecting plane to produce its normal equation $\tilde{\mathbf N}\cdot(x,y,z)^T + d = 0$, so I’ll just take them as given. The plane is then represented by the homogeneous vector $\mathbf\pi=\left(\tilde{\mathbf N}^T,d\right)^T$. A point $\mathbf x$ on the image plane back-projects to the ray $\mathtt P^T\mathbf x+\lambda\mathbf E_2$, which intersects the slicing plane at $$\mathbf X = (\mathbf E_2^T\mathbf\pi)(\mathtt P^T\mathbf x) - (\mathtt P^T\mathbf x)^T\mathbf\pi\mathbf E_2.\tag1$$ (This expression is obtained from the Plücker matrix of the back-projected ray; $\mathbf E_2$ is the basis vector $(0,1,0,0)^T$.) Substituting this into the equation of the sphere and fiddling with it a bit eventually produces the matrix $$\mathtt C = \begin{bmatrix} n_x^2+n_y^2 & n_xn_z & n_xd \\ n_xn_z & n_z^2+n_y^2 & n_zd \\ n_xd & n_yd & d^2-n_y^2 \end{bmatrix} \tag2$$ for the projected image of the intersection curve, which corresponds to the implicit Cartesian equation $$(n_x^2+n_y^2)x^2 + 2n_xn_z xy + (n_z^2+n_y^2)y^2 +2d(n_xx+n_zy) + (d^2-n_y^2) = 0. \tag3$$ In practice, you should probably check for $|\mathbf x^T\mathtt C\mathbf x| \lt \varepsilon$ for some fixed small positive $\varepsilon$, both to account for truncation errors and because it’s unlikely that many pixel centers will fall exactly on this curve.
With this approach, you’ll need to deal separately with the degenerate case $n_y=0$. The intersecting plane is parallel to the projection direction, so the image of the intersection should be a line segment, but the equation that results from the above is that of a (double) line.
Here’s an example plot from Mathematica using $(1,1,1)^T$ as the globe’s axis. This is the normal for the stack of planes used to generate the lines of latitude and the rotation axis for the the planes that generate the longitude.
† Actually, that conjugation to get the outline should be applied to the dual quadric and produces the dual of the outline, i.e., $\mathtt C^*=\mathtt P\mathtt Q^*\mathtt P^T$, but the unit sphere and resulting unit circle are self-dual, so I glossed over that complication.
Addendum: The above development works, with a couple of small modifications, for any quadric surface and any central projection onto a plane. With some other type of surface, you can work with the back-projection of the image point and directly test to see if that lies on the surface. We can do this here, too: expanding (1) in coordinates the image point $\mathbf x = (x,y,1)^T$ back-projects to $(n_y x,n_x x+n_z y+d,n_y y,n_y)$, and by using the homogeneous equation $X^2+Y^2+Z^2-W^2=0$ of the sphere directly you can avoid dehomogenizing. Computing this looks like it might be more efficient than using equation (3). You can verify that the equation of the ellipse that you get by expanding the homogeneous sphere equation is equivalent to (3).
For this simple case of a sphere and orthogonal projection, an equation of the ellipse can be developed directly, too. Switching now to a unit normal, $d$ is the signed distance from the plane to the origin, so the plane-sphere intersection is a circle with center $\tilde{\mathbf C}=-d\tilde{\mathbf N}$ and radius $r=\sqrt{1-d^2}$. The projection foreshortens along the projection of $\tilde{\mathbf N}$ by a factor of $\cos\phi$, where $\phi$ is the dihedral angle between the cutting and image planes. This is simply $n_y$. We thus get an ellipse centered at the image of $\tilde{\mathbf C}$ with half-axis lengths $r$ and $r\cos\phi$ and minor axis parallel to the image of $\tilde{\mathbf N}$. The equation of the minor axis is therefore $n_z x-n_x y = 0$ and that of the major axis is $n_x x+n_z z+d(n_x^2+n_z^2)=0$. Dividing each of these equation by $\sqrt{n_x^2+n_z^2}$ to normalize, we get the equation $${(n_z x-n_x y)^2 \over (1-d^2)(n_x^2+n_z^2)} + {(n_x x+n_z z+d(n_x^2+n_z^2))^2 \over n_y^2(1-d^2)(n_x^2+n_z^2)} = 1$$ for the ellipse. There are three degenerate cases: $n_y=0$ was already covered, for $d=\pm1$ the sphere-plane intersection is a point, and for $\tilde{\mathbf N}=(0,1,0)^T$, the cutting plane is parallel to the image plane, so the projected intersection is the circle $x^2+y^2=1-d^2$.