How to compute the 2D equations of 3D circular arcs?

Question

I'd like to obtain these equations for the ellipses produced by the perspective projections of (3-dimensionally transformed) circles.

This is useful for rendering in 2D contexts which provide curve primitives. I'm using HTML5's canvas, so I get Beziers, arcs, and quadratic curves.

See here:

The projection of a sphere outside of the plane of projection is an ellipse because the view is a cone (silhouette of a sphere is a circle).

However if I want to draw my sphere using circular wireframes, that projection-cone is no longer a circular cone. So it's not your traditional conic section anymore.

How to deal with this?

Answer 1

I wouldn't bother with heavy maths to get that ellipse's equation. Here's what I recommend:

• Discretize your wireframe sphere (the usual way, using lattitude/longitude, or with a geodesic grid, or fancier). This gives you a list of 3D vertices;
• Transform them to screen space with the usual model-view-projection matrix;
• Then just render line segments between those transformed vertices.

Of course you've got line segments instead of arcs, but drawing arcs in an HTML5 canvas kills performance, so I wouldn't recommend it. Just subdivide your sphere until the result looks OK.

There's a bunch of 3D engines for Javascript that implement software rendering. For instance, you could have a look at Three.js for inspiration.

That's if you want to do it all yourself, on the CPU. You could also use WebGL, there's already a tutorial for that. Browser support is not very broad yet, but it's supposed to be the future.

Answer 2

However if I want to draw my sphere using circular wireframes, that projection-cone is no longer a circular cone. So it's not your traditional conic section anymore.

The perspective projection of a 3d circular disk is also an ellipse, as long as the disk is entirely in front of the camera. (If the disk crosses the camera plane, it projects into a hyperbola).

It is only slightly more complicated to compute the ellipse parameters of a projected disk than to compute the ellipse parameters of a projected sphere. And as you've pointed out, there can be greater utility, efficiency, and accuracy in using curve primitives for rendering when curve primitives are available, especially when using a vector based renderer like SVG rather than raster based renderer.

There are several methods for computing the ellipse parameters of a perspective-projected circle, using geometric constructions, or using analytic / algebraic methods. I've found the algebraic methods easier to implement.

Perhaps the simplest algebraic method is to put your disk into matrix form [1], run it through the perspective transform [2], and then extract the center and axes from the result [3]. A simple alternative to the first two steps is to compute the projected ellipse equation by taking 5 sample points from your disk, projecting them to the image plane, and converting the 5 points directly into the ellipse equation coefficients [4], then proceed with extracting the center and axes from the coefficients [3].

There is a nice complete code example of computing the ellipse using a transform on ShaderToy [5] which has a blog post attached describing the math and an interesting application of it, undistorting text. [6]

One constructive geometry method is to perspective-project a square that encloses your 3d disk, and find the ellipse inscribed by the projected quad [7]. It uses known geometric relationships between the quad and the ellipse (such as the tangents of the contact points) to find the ellipse parameters. The geometric relationships of parts of the ellipse have been explored in depth [8]

And just in case you want to go the simpler route of drawing the sphere's silhouette ellipse (rather than wireframe disks), Inigo Quilez wrote a nice article deriving the math [9] and he provides a complete implementation [10].

[1] https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections

[2] https://math.stackexchange.com/a/2899066/616990

[3] https://en.wikipedia.org/wiki/Ellipse#General_ellipse

[4] https://math.stackexchange.com/a/163931/616990

[5] https://www.shadertoy.com/view/MlGGz3

[6] https://mzucker.github.io/2016/10/11/unprojecting-text-with-ellipses.html

[7] https://web.archive.org/web/20190924033522/http://www.math.nus.edu.sg/~graeme/Geometry/Circle_in_Perspective/construct_axes.html

[8] https://archive.org/details/constructivegeom00eagluoft/page/96/mode/2up (pages 96-150)

[9] https://www.iquilezles.org/www/articles/sphereproj/sphereproj.htm

[10] https://www.shadertoy.com/view/XdBGzd

Answer 3

I already answered the question How to project a 3D circle/ellipse to 2D? by providing the JavaScript implementation below.

Christopher Brierley Jones describes the process on his website how an ellipse can be fit perfectly into any convex quadrilateral.

So one way to compute the 2D equations of 3D circular arcs, without discretizing into line segments first, is to project the bounding rectangle and the extract the ellipse parameters from the resulting convex quadrilateral as described by Jones.

You can see a live demo here or a video here

Now if you want to project a sphere instead of a plain circle from 3d to 2d the steps would be:

1. create a circle with the same center as the sphere and the same radius but oriented with it's normal towards the camera plane, i.e. that the circle's plane is parallel to the projection plane As DMGregory pointed out the correct orientation for the circle face is not parallel to the projection plane but instead facing the camera center, as illustrated in this diagram.
2. so instead: create a circle with the same center as the sphere and the same radius but oriented with it's normal towards the camera position
3. create a bounding square around this circle
4. project the 4 vertices of this square with the camera projection matrix to obtain a quadrilateral
5. reconstruct the ellipse parameters by applying the procedure below to the resulting quadrilateral

// W,X,Y,Z are {x:number,y:number} vertices of the convex qudriliteral.
function quadToEllipse(W,X,Y,Z) {
    // Reconstruct matrix that transforms the unit square ((-1,-1), (1,1)) into quad (W,X,Y,Z)
    const m00 =  X.x * Y.x * Z.y - W.x * Y.x * Z.y - X.x * Y.y * Z.x + W.x * Y.y * Z.x - 
                W.x * X.y * Z.x + W.y * X.x * Z.x + W.x * X.y * Y.x - W.y * X.x * Y.x;
    const m01 =  W.x * Y.x * Z.y - W.x * X.x * Z.y - X.x * Y.y * Z.x + X.y * Y.x * Z.x - 
                W.y * Y.x * Z.x + W.y * X.x * Z.x + W.x * X.x * Y.y - W.x * X.y * Y.x;
    const m02 =  X.x * Y.x * Z.y - W.x * X.x * Z.y - W.x * Y.y * Z.x - X.y * Y.x * Z.x + 
                W.y * Y.x * Z.x + W.x * X.y * Z.x + W.x * X.x * Y.y - W.y * X.x * Y.x;
    const m10 =  X.y * Y.x * Z.y - W.y * Y.x * Z.y - W.x * X.y * Z.y + W.y * X.x * Z.y - 
                X.y * Y.y * Z.x + W.y * Y.y * Z.x + W.x * X.y * Y.y - W.y * X.x * Y.y;
    const m11 = -X.x * Y.y * Z.y + W.x * Y.y * Z.y + X.y * Y.x * Z.y - W.x * X.y * Z.y - 
                W.y * Y.y * Z.x + W.y * X.y * Z.x + W.y * X.x * Y.y - W.y * X.y * Y.x;
    const m12 =  X.x * Y.y * Z.y - W.x * Y.y * Z.y + W.y * Y.x * Z.y - W.y * X.x * Z.y - 
                X.y * Y.y * Z.x + W.y * X.y * Z.x + W.x * X.y * Y.y - W.y * X.y * Y.x;
    const m20 =  X.x * Z.y - W.x * Z.y - X.y * Z.x + W.y * Z.x - X.x * Y.y + W.x * Y.y + X.y * Y.x - W.y * Y.x;
    const m21 =  Y.x * Z.y - X.x * Z.y - Y.y * Z.x + X.y * Z.x + W.x * Y.y - W.y * Y.x - W.x * X.y + W.y * X.x;
    const m22 =  Y.x * Z.y - W.x * Z.y - Y.y * Z.x + W.y * Z.x + X.x * Y.y - X.y * Y.x + W.x * X.y - W.y * X.x;

    // invert matrix
    const determinant = +m00*(m11*m22-m21*m12) -m01*(m10*m22-m12*m20) +m02*(m10*m21-m11*m20);
    if(determinant == 0) return null;
    const invdet = 1/determinant; 
    const J =  (m11*m22-m21*m12)*invdet;
    const K = -(m01*m22-m02*m21)*invdet;
    const L =  (m01*m12-m02*m11)*invdet;
    const M = -(m10*m22-m12*m20)*invdet;
    const N =  (m00*m22-m02*m20)*invdet;
    const O = -(m00*m12-m10*m02)*invdet;
    const P =  (m10*m21-m20*m11)*invdet;
    const Q = -(m00*m21-m20*m01)*invdet;
    const R =  (m00*m11-m10*m01)*invdet;

    // extract ellipse coefficients from matrix
    const a = J*J + M*M - P*P;
    const b = J*K + M*N - P*Q;
    const c = K*K + N*N - Q*Q;
    const d = J*L + M*O - P*R;
    const f = K*L + N*O - Q*R;
    const g = L*L + O*O - R*R;

    // deduce ellipse center from coefficients
    const centerX = (c*d - b*f) / (b*b - a*c);
    const centerY = (a*f - b*d) / (b*b - a*c);

    // deduce ellipse radius from coefficients
    const radiusA = Math.sqrt(2*(a*f*f + c*d*d + g*b*b - 2*b*d*f - a*c*g)/((b*b - a*c) * (Math.sqrt((a-c)*(a-c) + 4*b*b) - (a+c))));
    const radiusB = Math.sqrt(2*(a*f*f + c*d*d + g*b*b - 2*b*d*f - a*c*g)/((b*b - a*c) * (-Math.sqrt((a-c)*(a-c) + 4*b*b) - (a+c))));

    // deduce ellipse rotation from coefficients
    let angle = 0;
    if(b==0 && a <= c) {
        angle = 0;
    } else if(b == 0 && a >= c) {
        angle = Math.PI / 2;
    } else if(b != 0 && a > c) {
        angle = Math.PI / 2 + 0.5 * (Math.PI / 2 - Math.atan2((a-c), (2*b)));
    } else if(b != 0 && a <= c) {
        angle = Math.PI / 2 + 0.5 * (Math.PI / 2 - Math.atan2((a-c), (2*b)));
    }
            
    return {
        centerX, centerY, radiusA, radiusB, angle
    }