For most people, anisotropic filtering is something you just flick on for instant gratification and then proceed to forget about it because the hardware does all the work automagically. Gamers crank it up to the maximum magic number of 16 in settings because it seemingly has no performance impact and programmers read that setting and just write it into the sampler's description. They know it makes things less blurry at oblique angles, but seldom how it does it.
I've been asked to make the problem more approachable numerous times; to provide a guided derivation which doesn't make people want to jump off a cliff. Specifically, the request was to skirt any mention of signal processing which is a tall order. I've been reluctant for many moons because postulating the problem clearly for people without a strong background is a non-trivial affair... But let's give it a shot.
Basics of texture mapping
Gotta be easy, right?
Early on, rasterizers only dealt with generating procedural detail across the surface of rasterized triangles. To add more realism, texture mapping was introduced. The idea was to encode surface details into a 2D array, an image, a texture that was stored in memory and accessed at runtime as needed. The contents of it can be anything you want, from material information to a drawing of Hello Kitty. One "wraps" it around the geometry by creating a mapping between it and the vertices. This is a non-trivial, challenging and even artistic process which involves creatively flattening a 3D mesh ("unwrapping") into a normalized 2D space with minimal seams, then storing those 2D coordinates of flattened vertices on the vertices as an attribute which we colloquially refer to as "UVs". Naturally, when accessing a texture, the normalized texture coordinates are first scaled by the dimensions of the texture being accessed. Depending on texture addressing mode, texture coordinates leaving the normalized range can be well-defined by allowing wrapping around (enabling tiling), edge clamping etc.
Vertices, and their attributes, get projected to the screen during rasterization. But before we can use these coordinates to actually access any texture, we need to interpolate them from the corners of the triangle being rasterized to the pixel we're currently processing. Correctly interpolating all the vertex attributes is of paramount importance, we have to account for perspective. Simply put, under the perspective transform, uniform steps in the image plane (iterating pixel by pixel of uniform width and height) do not necessarily map to uniform steps in the plane containing the triangle being interpolated. Triangle vertices in view space can have high differences in depth and, thus, we have to effectively interpolate linearly in the plane of the triangle. You can find my derivation here, but you can easily see the gist of it in the diagram below, the higher the tilt with respect to the image plane, the smaller the projection of the triangle. Naively interpolating on the projected triangle is thus wrong as $s \neq t$ .
By finding the appropriate transform for said interpolation and fully simplifying it, we discover that the inverse of depth interpolates linearly in screen space. Knowing how to correctly interpolate $z$ allows us to solve for interpolated vertex attributes because the interpolation factor $t = (z - z_0) / (z_1 - z_0)$ is shared by any attribute $a$ being interpolated: $t = (a - a_0) / (a_1 - a_0)$. By equating the two, you can derive $a$. Some rearranging ahead, turns out vertex attributes divided by the associated vertex's depth interpolate linearly in screen space as $a/z = (1-t)a_0/z_0 + ta_1/z_1$, recovering the final value comes down to multiplying by the interpolated $z$. As texture coordinates are vertex attributes, we have our means for interpolating texture coordinates to pixel centers.
So, what now? Do you just point sample the texture with the correctly interpolated texture coordinates? Hahah, no. That's your one-way ticket to Aliasville. To understand why not, consider your screen. A two-dimensional grid of square pixels (assume they're square for simplicity). Each pixel has a non-zero square area. Each, thus, must faithfully represent the area of a triangle projected to it. The projected triangle area may be covered by thousands of texels all competing for the attention of the poor single pixel, depending on the geometric relationship (distance, orientation) between the image plane and the plane containing the triangle. Huge area far away or tilted in view space projects to a tiny area. In other words, the partial derivatives of the texture coordinates wrt screen coordinates at the particular $(x,y)$ are huge.
Theoretically (and ideally), we can project the screen space pixel square into texture space where it morphs into a quadrilateral which covers the texels that desire to be considered. The goal is to iterate over all the texels covered by this quadrilateral (and not one more because ✨blur✨) and filter them into a single value to be used for the pixel being lit up. As you might imagine, this is computationally expensive to do at runtime. Any kind of direct convolution like this is, by definition.
In the 80s [notably Heckbert '86], a number of solutions performed the convolution offline, but with a twist. The unwieldy quadrilaterals are hard to handle mathematically, determining whether a texel is inside or out is non-trivial. Hence, Heckbert and those before him used to replace the projection with a local affine approximation, which approximates the unwieldy quadrilateral with a comparatively-tame parallelogram or an ellipse (if you modeled pixels as circles). This was still extremely expensive for realtime use, though, but it was doable. They used a weight LUT to actually figure out how each of the texels was weighted for the final filtered value to be computed. We can't afford the convolution at runtime, at all. We need to figure out a way to precompute this. At texture creation time, preferably.
Inventing "multum-in-parvo" maps
The multi-isotropic compromise
We can approximate a projection of the screen space pixel square (or circle) to a parallelogram (or ellipse) in texture space with a linear affine approximation, the details of which are the subject of next section. While that gives us the boundaries on the texels that require filtering, it's not particularly useful to us unless we can somehow prefilter the texture. It's not feasible to cover every possible "projected" parallelogram as you can see in the sample texture below:
Those are just a few examples and it is immediately obvious that no lookup structure can be devised here. What we could do is impose some restrictions on the filter shape, for example for it to be isotropic (same in every direction) and axis-aligned:
This kind of assumption makes creating the lookup structure a lot simpler, but it will not handle anisotropy well as a square naively approximating the parallelogram-shaped filter is going to capture texels outside of the filter shape, introducing overblurring in that direction. But if you think about it, we're not limited to picking just a single square. We can approximate a parallelogram-shaped filter by doing multiple square ones placed along the line of anisotropy (the "major axis" of the parallelogram):
To facilitate this isotropic prefilter structure, known as a mipmap, we create a hierarchy of textures where each level is quarter (or half if you're looking at dimensions individually) the resolution of the previous level and the data stored in it are the filtered values of all the $2 \times 2$ blocks of the previous level. This goes all the way to a level which is effectively a $1\times 1$ texture representing the filtered value of a $\textsf{w} \times \textsf{h}$ block of the texture. Since it's square and axis-aligned, you only need the center position (basically, the texture coordinates) and a single value to access the correct levels, the width of a square of texels we're filtering.
Reducing by $4$ on each level means it reduces by $2$ in each dimension. Given a width $w_f$ of square texels we're filtering (determined in the next section), we can take the $\log_2$ of $w_f$ to get a zero-based "index" into the levels. In quotes because it can have fractional parts, for example $\log(17) = 4.087$ meaning that we have to perform two bilinear samples on levels $4$ and $5$ and then use the fractional part $0.087$ to get the final filtered value, process called trilinear filtering. Final example, say $w_f$ turns out to be $2$, it means we want to access filtered $2 \times 2$ blocks. $\log_2(2)$ is $1$, meaning we want the first level where each texel represents filtered $2\times 2$ blocks of texels of the source texture (zeroth is the source texture) and so on. As we seldom want the actual blocks that were prefiltered and closer to those in-between, we perform bilinear filtering on the $4$ nearest texels each representing a filtered block of the source texture (size of block depends on level computed).
Rounding up texels
We're not going to hurt them, we just want to talk to them.
In the previous section, we found a way to prefilter our texture data and, from the diagrams, know that we need the line of anisotropy to determine the direction along which multiple isotropic mipmap samples are to be taken and the filter width $w_f$ from which we can compute the level(s) at which the mipmap will be sampled. According to the diagram, the line of anisotropy lies on the major axis of the ellipse drawn inside of our parallelogram of interest while the width $w_f$ can be derived from the distance from the extreme points on the ellipse along the minor axis. As you can clearly see, the ellipse can replace our parallelogram in the calculations and greatly simplify them as we only care about the extreme points on the ellipse (which are shared by the parallelogram). Thus, from now on, we'll use a circle instead of a square in all our calculations because ellipses intrinsically admit the notion of a major and minor axis.
As for the projection itself, it is reasonable to assume that we can approximate it with a locally-correct linear (affine) transform. To formulate such a transform, we need to determine how the texture coordinates change in the pixel's vicinity (the relationship between the two spaces). In other words, we're interested in computing the partial derivatives $\partial u/ \partial x$, $\partial u/ \partial y$, $\partial v/ \partial x$ and $\partial v/ \partial y$. Grouping these gives us two vectors indicating how $u$ and $v$ change in the XY plane around the point $(x,y)$, $\mathbf{dX}$ and $\mathbf{dY}$.
We don't actually compute partial derivatives proper on the GPU, pixel shaders instead execute in $2\times 2$ blocks which enables us to recover approximate screenspace partial derivatives via finite differencing (just subtracting), sufficient for practical applications dealing with discrete signals. By arranging our first-order partial derivative vectors as columns of a $2 \times 2$ matrix, we get a linear transform which approximately takes us to into texture space. This construct, when used more accurately, is known as the Jacobian (matrix):
$$\begin{split}
\Large { \mathbf{J} = \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial v}{\partial x} \\ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} \end{bmatrix}} \\
\end{split}$$
This linear transform, again, essentially approximates the reverse projection transform locally, we know how $(u,v)$ changes locally with respect to $(x,y)$. It is an approximation because we're assuming it'll be valid for more than just the point it was calculated on, ie the assumption of linearity is what is approximative.
So how do circles and ellipses work, again?
We've all dealt with circles at some point in our life, which are just a rotationally invariant kind of ellipse. The simplest formulation is usually $x^2 + y^2 = r^2$ which is trivial to understand. If you divide by $r^2$, you arrive at $(x/r)^2 + (y/r)^2 = 1$. You can see how the $r$ and $1$ on the right side constrict $x$ and $y$ isotropically from going beyond it. If you replace $r$ with $a$ and $b$, you can control the horizontal and vertical extent, effectively creating an ellipse centered on $(0,0)$ whose axes are aligned with the axes of the coordinate system:
$$\begin{split}
\Large{\left(\frac{x}{a}\right)^2 + \left(\frac{y}{b}\right)^2 = 1} \\
\end{split}$$
This form is a bit bare-bones, but it gets the point across. The full, implicit equation for an arbitrary ellipse is:
$$\begin{split}
ax^2 + bxy + cy^2 + dx + ey - f = 0 \\
\end{split}$$
$d$ and $e$ terms make the ellipse "leave" the origin and we're not particularly interested in them, so we drop them:
$$\begin{split}
ax^2 + bxy + cy^2 - f = 0 \\
\end{split}$$
It can be converted into matrix form if you transform it a bit:
$$\begin{split}
ax^2 + bxy + cy^2 &= f \\
ax^2 + \frac{b}{2}xy + \frac{b}{2}xy + cy^2 &= \\
x(ax + \frac{b}{2}y) + y(\frac{b}{2}x + cy) &= \\
\begin{bmatrix} x & y \end{bmatrix} \begin{bmatrix} a & b/2 \\ b/2 & c \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} &= \\
\end{split}$$
We have now shown that the following holds:
$$\begin{split}
ax^2 + bxy + cy^2 = \mathbf{p}\mathbf{q}\mathbf{p}^\mathsf{t} = f \\
\end{split}$$
This is neat because we can apply transformations to our implicit ellipse parametrization as if it were any other primitive. However, there is a twist, we apply the inverse of whatever transform we have in mind. Intuition for this can be found in the way you translate something like a line to the right by actually moving the entire thing left on the $x$ term. Recall that $x$ and $y$ effectively go from $-\infty$ to $\infty$ and there's no transforming them. If we want to apply the transform $\mathbf{M}$ to facilitate $\mathbf{p'} = \mathbf{p}\mathbf{M}$, we plug in the inverse into a hypothetical ellipse $\mathbf{p'}\mathbf{Q'}\mathbf{p'}^\mathsf{T} = f'$ to yield:
$$\begin{split}
\mathbf{p'}\mathbf{Q'}\mathbf{p'}^\mathsf{T} &= f'\\
\mathbf{pM^{-1}} \mathbf{Q'} \left(\mathbf{pM^{-1}}\right)^\mathsf{T} &= f'\\
\mathbf{pM^{-1}} \mathbf{Q'} \mathbf{M^{-1}}^\mathsf{T}\mathbf{p}^\mathsf{T} &= f' = f\\
\end{split}$$
We then have the transformed coefficient matrix specified as $\mathbf{Q} = \mathbf{M^{-1}} \mathbf{Q'} \mathbf{M^{-1}}^\mathsf{T}$. Also note that $f = f'$. Surely now you can see that one can start with a circle coefficient matrix (the identity matrix) and transform it into any ellipse. Further still, we can take a more generic, parametric formulation of the circle such as $\mathbf{p} = (\cos t, \sin t)$ and transform it arbitrarily then use $\mathbf{p}\mathbf{p}^{\textsf{T}} = \cos^2 t + \sin^2 t = 1$ to recover $\mathbf{Q}$ and $f$ of the implicit formulation.
Okay, now what?
What we're most interested in, though, is recovering the major and minor axis of the ellipse once we transform the parametric unit circle (representing our pixel) to the implicit ellipse using the finite difference relation clueing us in as to how texture space's $(u,v)$ changes locally with respect to screen space's $(x,y)$, our peasant-constructed Jacobian computed from $2\times 2$ blocks of pixels via finite differencing.
First things first, from the parametric unit circle to the implicit target ellipse. As mentioned, we have the Jacobian peasantly computed, so it comes down to merely transforming the points with it and recovering $\mathbf{Q}$ to get our $a, b, c, f$ coefficients:
$$\mathbf{p}\mathbf{p}^{\textsf{T}} = \cos^2 t + \sin^2 t = 1$$
$$\mathbf{p'M^{-1}} \mathbf{M^{-1}}^\mathsf{T}\mathbf{p'}^\mathsf{T} = 1 \iff \mathbf{p} = \mathbf{p'M^{-1}}$$
$$\mathbf{M^{-1}} \mathbf{M^{-1}}^\mathsf{T} = \frac{\begin{bmatrix} v_y & -v_x \\ -u_y & u_x \end{bmatrix} \begin{bmatrix} v_y & -u_y \\ -v_x & u_x \end{bmatrix}}{(u_xv_y - u_yv_x)^2} \iff M = \begin{bmatrix} u_x & v_x \\ u_y & v_y \end{bmatrix}$$
$f$ is in plain sight, it's just the denominator. The rest of them are a matrix multiply away, just map them to the $\mathbf{Q}$ entries:
$$\begin{split}
a &= v_x^2 + v_y^2 \\
b &= -2(u_xv_x + u_yv_y) \\
c &= u_x^2 + u_y^2 \\
f &= (u_xv_y - u_yv_x)^2 \\
\end{split}$$
While we now have our circle "projected" to an ellipse in texture space, we need to find its major and minor axis that we'll be using for filtering (specifically, farthest points along those axes). The minor will be used to select the LOD of the mipmap and the major is our line of anisotropy we'll be sampling along at the LOD selected by the minor axis. Think back to our parametric circle prior to transform, it's a locus of points, yes? As we applied our linear transformation to it, the Jacobian, it elongated into an ellipse. There are four points on that ellipse which remained on the same two lines from the origin under the transform we performed, as if they were scaled by some value along a normalized direction which could replace the transform in that case. Yep, it's our major and minor direction "end points". This is a trivial (albeit tedious) eigenvalue and eigenvector determination problem.
Any ellipse centered around the origin can be viewed as the result of non-uniformly scaling and rotating a circle, meaning we can equate our Jacobian transform with the following:
$$\mathbf{M} = \mathbf{SR} = \begin{bmatrix} u_x & v_x \\ u_y & v_y \end{bmatrix} = \begin{bmatrix} a & 0 \\ 0 & b \end{bmatrix}\begin{bmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix}$$
The scaling is orthogonal while the rotation is orthonormal. As you know, the latter means rows are both mutually orthogonal and of unit length which makes its inverse equal to its transpose. Given that we define the transformed $\mathbf{Q}$ here as $\mathbf{M^{-1}}\mathbf{I}\mathbf{M^{-1}}^\mathsf{T}$, we have:
$$\mathbf{Q} = \mathbf{(SR)^{-1}}\mathbf{(SR)^{-1}}^\mathsf{T} = \mathbf{R}^{-1} \mathbf{S}^{-1} \mathbf{S^{-1}}^{\textsf{T}} \mathbf{R} = \mathbf{R}^{-1} \mathbf{S^{-2}} \mathbf{R}$$
Any well-behaved $n\times n$ matrix $\mathbf{A}$ encoding a linear transform will have $n$ eigenvalues and we should be able to find $n$ linearly independent eigenvectors. If you were to place all the eigenvectors as rows of a square matrix $\mathbf{E}$ and multiplied it with the associated matrix $\mathbf{A}$ the result would be a matrix where all the rows were multiplied with respective eigenvalues $\lambda_i$... As if you had built a diagonal square matrix containing eigenvalues $\mathbf{D}$ and multiplied it with $\mathbf{E}$ yielding $\mathbf{EA} = \mathbf{DE} \iff \mathbf{EAE^{-1}} = \mathbf{D}$. Going the "other way", we have $\mathbf{E^{-1}DE} = \mathbf{A}$. This leaves us with the following equivalence:
$$\mathbf{Q} = \mathbf{R}^{-1} \mathbf{S^{-2}} \mathbf{R} = \mathbf{E}^{-1} \mathbf{D} \mathbf{E}$$
Hence, $\mathbf{R} = \mathbf{E}$ and $\mathbf{S} = \mathbf{D}^{-1/2}$. Computing the eigenvalues for a generic $\mathbf{Q}$ amounts to:
$$\begin{bmatrix} x & y \end{bmatrix}\begin{bmatrix} (a - \lambda) & b/2 \\ b/2 & (c - \lambda) \end{bmatrix} = 0$$
$$(a - \lambda)(c - \lambda) - \frac{b^2}{4} = 0$$
$$\lambda^2 - (a+c)\lambda + \frac{4ac - b^2}{4} = 0$$
$$\lambda_{1/2} = \frac{a + c \pm \sqrt{(a+c)^2 - 4ac + b^2}}{2}$$
$$\lambda_{1/2} = \frac{a + c \pm \sqrt{(a-c)^2 + b^2}}{2}$$
It seems the square root and the $a \pm c$ terms will be following us around, so we might as well give them names. In the interest of keeping the result succinct, we can introduce the following substitutions $q = a + c, p = a - c, t = \sqrt{ p^2 + b^2}$. With $\mathbf{D}$, we easily resolve $\mathbf S$:
$$\mathbf{S} = \begin{bmatrix} \sqrt{\frac{2}{q + t}} & 0 \\ 0 & \sqrt{\frac{2}{q - t}} \end{bmatrix} \iff \mathbf{D} = \begin{bmatrix} \frac{q + t}{2} & 0 \\ 0 & \frac{q - t}{2} \end{bmatrix}$$
Going forward, remember that direction is our main concern, and as it aliases $\mathbf{R}$, we want them to be unit length and mutually orthogonal (a higher requirement than linear independence). Again, direction remains pristine. Thus, all shared magnitude-varying scalar factors will be dropped unceremoniously.
$$\begin{bmatrix} x & y \end{bmatrix}\begin{bmatrix} a & b/2 \\ b/2 & c \end{bmatrix} = \lambda_{1/2}\begin{bmatrix} x & y \end{bmatrix}$$
$$ax + \frac{b}{2}y = \lambda_{1/2}x \space \Rightarrow \space y = \frac{\pm t - p}{b}x$$
$$\frac{b}{2}x + cy = \lambda_{1/2}y \space \Rightarrow \space y = \frac{b}{\pm t + p}x$$
This nets us a pair of eigenvectors for each eigenvalue (remember the $\pm$):
$$\mathbf{e_{\lambda_{1/2}}} = \begin{bmatrix} b & \pm t - p \end{bmatrix}, \begin{bmatrix} \pm t + p & b \end{bmatrix}$$
Our goal is to normalize these vectors before moving forward, so let's do the first one of the first pair first... See what shakes out:
$$\begin{split}
n_{\lambda_1,1} &= \sqrt{ (\pm t - p)^2 + b^2 } \\
&= \sqrt{ t^2 - 2tp + p^2 + b^2 } \\
&= \sqrt{ 2t^2 - 2tp } \\
&= \sqrt{2t}\sqrt{ t - p } \\
\end{split}$$
Similarly, you will find that $n_{\lambda_1,2} = \sqrt{2t}\sqrt{ t + p }$. Further, $n_{\lambda_2,i} = n_{\lambda_1,i}$ as the squaring action terminates the $\pm$ associated with $t$. But the biggest realization here is that $b$ doesn't contribute as a free term meaning that it has no bearing on the direction. We're then free to reuse the two distinct length calculations and recombine them into a single vector. Then, we can express the other vector as an implicit 90° rotation (as you know, swap the coordinates of the first one, negate one of them). This simplification will cost us some implicit sign information, but we can recover that from the sign of $b$. We have:
$$\mathbf{e_{\lambda_{1}}} = \begin{bmatrix} \sqrt{2t}\sqrt{ t + p } & \textsf{sgn}(b)\sqrt{2t}\sqrt{ t - p } \end{bmatrix}$$
$$\mathbf{e_{\lambda_{2}}} = \begin{bmatrix} -\textsf{sgn}(b)\sqrt{2t}\sqrt{ t - p } & \sqrt{2t}\sqrt{ t + p } \end{bmatrix}$$
As we isolated a scalar factor from all of them, we can drop it because it doesn't vary the direction:
$$\mathbf{e_{\lambda_{1}}} = \begin{bmatrix} \sqrt{ t + p } & \textsf{sgn}(b)\sqrt{ t - p } \end{bmatrix}$$
$$\mathbf{e_{\lambda_{2}}} = \begin{bmatrix} -\textsf{sgn}(b)\sqrt{ t - p } & \sqrt{ t + p } \end{bmatrix}$$
In this state, it becomes clear that the normalization factor of both is equivalent and it's $\sqrt{2t}$ which we will just prepend below. We set the two vectors as rows of $\mathbf{R}$ as it's equal to $\mathbf{E}$. Then we multiply it with $\mathbf{S}$. The rows of the resulting matrix represent the final directions for the major and minor axis.
$$\begin{split}
\mathbf{M} = \mathbf{SR} &= \sqrt{\frac{1}{2t}} \begin{bmatrix} \sqrt{\frac{2}{q +t}} & 0 \\ 0 & \sqrt{\frac{2}{q - t}} \end{bmatrix} \begin{bmatrix} \sqrt{t + p} & \textsf{sgn}(b)\sqrt{t - p} \\ -\textsf{sgn}(b)\sqrt{t - p} & \sqrt{t + p} \end{bmatrix} \\
\end{split}$$
Therefore, our final expression for $\mathbf{M}$ is:
$$\begin{split}
\large{ M = \begin{bmatrix} \sqrt{\frac{t+p}{t(q +t)}} & \textsf{sgn}(b)\sqrt{\frac{t-p}{t(q +t)}} \\ -\textsf{sgn}(b)\sqrt{\frac{t-p}{t(q - t)}} & \sqrt{\frac{t+p}{t(q - t)}} \end{bmatrix} }\\
\end{split}$$
From this we can calculate the farthest points of the ellipse along the normalized major and minor axis directions by factoring in $f$, to finalizing our boundary computation (if you can't figure out how $f$ works here, you're reading this article prematurely):
$$\mathbf{dX_{TS}} = \begin{bmatrix} \sqrt{f\frac{t+p}{t(q +t)}} & \textsf{sgn}(b)\sqrt{f\frac{t-p}{t(q +t)}} \end{bmatrix}$$
$$\mathbf{dY_{TS}} = \begin{bmatrix} -\textsf{sgn}(b)\sqrt{f\frac{t-p}{t(q - t)}} & \sqrt{f\frac{t+p}{t(q - t)}} \end{bmatrix}$$
These points can be interpreted as the transformed $\mathbf{dX}$ and $\mathbf{dY}$ vectors, hence the name. Doesn't really matter what you name them, though. Whichever of the two is longer can be thought of as defining the line of anisotropy. The shorter one, its length is more important. Specifically, you need to double it to get the width of the square filter we'll be using, the $\log_2$ of which nets us the mipmap level selection value. From there, $n$ trilinearly filtered samples are taken along the line of anisotropy (defined by the other vector) and filtered into a single value. $n$ is conservative, it is dynamically determined by the ratio of the anisotropy with a maximum defined by the user (and hardware). It tries to do it with least samples, unless the count of "least samples" is greater than the maximum anisotropy set by the user / programmer. At which point, it is clamped. The spacing of the samples is determined by the hardware, how much the axis-aligned squares end up overlapping affects weighing.
And in case you were wondering, yup, this is pretty much how Direct3D, OpenGL and Vulkan perform LOD and Line Of Anisotropy determination per specification. To properly get a feel for ellipses, I suggest you use these equations we derived on a graphing tool which are usually available online, my favorite is Desmos.
