Calculating an AABB for a transformed ellipse

Question

I am looking to compute the axis-aligned bounding box (AABB) of a 2D ellipse on which a tranformation matrix was applied (rotation, scale, translation, etc.)

Something similar to this solution : Calculating an AABB for a transformed sphere

So far, it doesn't seem to work for 2D ellipses.

This is what I got (in pseudo-code) :

Matrix M; // Transformation matrix (already existing)
Matrix C = new Matrix( // Conic matrix
    radiusX, 0,         0,
    0,       radiusY,   0,
    0,       0,         -1
);

Matrix MT = M.transpose();
Matrix CI = C.inverse();

Matrix R = M*CI*MT;

int minX = (R13 + sqrt(R13^2 - (R11 * R33))) / R33;
int minY = (R23 + sqrt(R23^2 - (R22 * R33))) / R33;

// maxX etc...
// Build AABB Rectangle out of min & max...

Simple demo of the current behavior

radiusX = 2    
radiusY = 2                              // To keep it simple, M is identity
                                         // (no transformation on the ellipse)
M = /1 0 0\                              // /M11 M21 M31\ 
    |0 1 0|                              // |M12 M22 M32| Transform matrix format
    \0 0 1/                              // \0   0   1  /

C = /2 0  0\                             // C as conic
    |0 2  0|
    \0 0 -1/

CI =/0.5 0   0\                          // CI as dual conic
    |0   0.5 0|
    \0   0  -1/

R = /1 0 0\ * /0.5 0   0\ * /1 0 0\      // R = M*CI*MT
    |0 1 0|   |0   0.5 0|   |0 1 0|
    \0 0 1/   \0   0  -1/   \0 0 1/

  = /0.5 0   0\                          // /R11 R12 R13\
    |0   0.5 0|                          // |R12 R22 R23| (R is symmetric)
    \0   0  -1/                          // \R13 R23 R33/

minX = (0 + sqrt(0^2 - (0.5 * -1))) / -1

     = -0.7071                           // Should be -2

                                         // Also, using R = MIT*C*MI
                                         // leads to -1.4142

Solution (using dual conic matrix)

Matrix M;
Matrix C = new Matrix(
    1/radiusX^2, 0,           0,
    0,           1/radiusY^2, 0,
    0,           0,           -1
);
Matrix MT = M.transpose();
Matrix CI = C.inverse();
Matrix R = M*CI*MT;

int minX = (R13 + sqrt(R13^2 - (R11 * R33))) / R33;
int minY = (R23 + sqrt(R23^2 - (R22 * R33))) / R33;

Final solution (no direct use of conic matrix)

Here's a simplified version.

Matrix M;

int xOffset = sqrt((M11^2 * radiusX^2) + (M21^2 * radiusY^2));
int yOffset = sqrt((M12^2 * radiusX^2) + (M22^2 * radiusY^2));

int centerX = (M11 * ellipse.x + M21 * ellipse.y) + M31; // Transform center of 
int centerY = (M12 * ellipse.x + M22 * ellipse.y) + M32; // ellipse using M
                                                         // Most probably, ellipse.x = 0 for you, but my implementation has an actual (x,y) AND a translation
int xMin = centerX - xOffset;
int xMax = centerX + xOffset;

int yMin = centerY - yOffset;
int yMax = centerY + yOffset;

Answer 1

From dual conic

So you state that M is a transformation matrix. But what does it transform, is it points or lines? I assume points. How do you represent points, as a row vector so that the point is on the left and the matrix on the right, or as a column vector so that the matrix is on the left and the point on the right of a multiplication? I'll assume column vectors. So a transformation would be p' = M*p for some point p.

Next is C. The way you write it, that's an ellipse but not with the radii you are using. A point lies on the ellipse if it satisfies (x/radiusX)^2 + (y/radiusY)^2 = 1 so the values on the main diagonal have to be (1/radiusX^2, 1/radiusY^2, -1). I repeatedly missed this mistake in pervious revisions of my answer.

Next you combine these things. Suppose CP were the primal conic, i.e. the conic as a set of points. Then you'd obtain the transformed version by doing MT.inverse()*CP*M.inverse(). The reason is because you apply M.inverse() to every point and then check whether it lies on the original conic. But you are not using M.inverse(), you are using M. This indicates that you try to transform a dual conic. If M transforms points, then MT.inverse() transforms lines, so M*CD*MT is the correct transformation if CD is a dual conic.

And if R is a dual conic, then your formulas are correct. So perhaps the main problem with your code is the fact that you forgot to use inverse radii in the matrix C.

From primal conic

When I read your post for the first time, I assumed R would describe a set of points, i.e. that a point (x,y) lies on that ellipse if (x,y,1)*R*(x,y,1).transpose()=0. Based on this, I did come up with formulas for the AABB without using the dual conic. I'm not saying that this is simpler, particularly not if you have matrix inversion available as a building block. But I'll still leave it here for reference. Keep in mind that the R in this paragraph is a different one from the one used in your code example.

For my approach, consider that R*(1,0,0) (which is simply the first column of R) is some vector (a,b,c) which you can interpret as a definition of a line ax+by+c=0. Intersect that line with the conic and you get the points where the tangents are horizontal, which are the extrema in y direction. Do the same for R*(0,1,0) (i.e. the seond column) to find extrema in the x direction.

The key idea here is that R*p computes the polar line for some point p, so we are constructing the polar line for the point at infinity in x resp. y direction. That polar line will intersect the conic in those points where the tangents through p touch the conic, which in this case would be horizontal resp. vertical tangents since parallel lines intersect at infinity.

If I do the above computation symbolically, I get the following formulas:

xmin, xmax = (R13*R22^2 - R12*R22*R23 ± sqrt(R13^2*R22^4 - 2*R12*R13*R22^3*R23 + R11*R22^3*R23^2 + (R12^2*R22^3 - R11*R22^4)*R33))/(R12^2*R22 - R11*R22^2)
ymin, ymax = (R11*R12*R13 - R11^2*R23 ± sqrt(R11^3*R13^2*R22 - 2*R11^3*R12*R13*R23 + R11^4*R23^2 + (R11^3*R12^2 - R11^4*R22)*R33))/(R11^2*R22 - R11*R12^2)

These expressions can certainly be simplified, but it should get you started. Feel free to edit this post if you reformulate this to something simpler, or easier to read, or whatever.