Tissot's indicatrix

Let a 3D point, ${\displaystyle {\hat {X}}}$, on an ellipsoid be parameterized as:

${\displaystyle {\hat {X}}(\lambda ,\phi )=\left[{\begin{matrix}N\cos {\lambda }\cos {\phi }\\-N(1-e^{2})\sin {\phi }\\N\sin {\lambda }\cos {\phi }\end{matrix}}\right]}$

where ${\displaystyle (\lambda ,\phi )}$ are longitude and latitude, respectively, and ${\displaystyle N}$ is a function of the equatorial radius, ${\displaystyle R}$, and eccentricity, ${\displaystyle e}$:

${\displaystyle N={\frac {R}{\sqrt {1-e^{2}\sin ^{2}(\phi )}}}}$

The element of distance on the sphere, ${\displaystyle ds}$ is defined by the first fundamental form:

${\displaystyle ds^{2}={\begin{bmatrix}d\lambda &d\phi \end{bmatrix}}{\begin{bmatrix}E&F\\F&G\end{bmatrix}}{\begin{bmatrix}d\lambda \\d\phi \end{bmatrix}}}$

whose coefficients are defined as:

${\displaystyle {\begin{aligned}&E={\frac {\partial {\hat {X}}}{\partial \lambda }}{\boldsymbol {\cdot }}{\frac {\partial {\hat {X}}}{\partial \lambda }}\\&F={\frac {\partial {\hat {X}}}{\partial \lambda }}{\boldsymbol {\cdot }}{\frac {\partial {\hat {X}}}{\partial \phi }}\\&G={\frac {\partial {\hat {X}}}{\partial \phi }}{\boldsymbol {\cdot }}{\frac {\partial {\hat {X}}}{\partial \phi }}\\\end{aligned}}}$

Computing the necessary derivatives gives:

${\displaystyle {\frac {\partial {\hat {X}}}{\partial \lambda }}=\left[{\begin{matrix}-N\sin {\lambda }\cos {\phi }\\0\\N\cos {\lambda }\cos {\phi }\end{matrix}}\right]\qquad \qquad {\frac {\partial {\hat {X}}}{\partial \phi }}=\left[{\begin{matrix}-M\cos {\lambda }\sin {\phi }\\-M\cos {\phi }\\M\sin {\lambda }\sin {\phi }\end{matrix}}\right]}$

where ${\displaystyle M}$ is a function of the equatorial radius, ${\displaystyle R}$, and the ellipsoid eccentricity, ${\displaystyle e}$:

${\displaystyle M={\frac {R(1-e^{2})}{(1-e^{2}\sin ^{2}(\phi ))^{\frac {3}{2}}}}}$

Substituting these values into the first fundamental form gives the formula for elemental distance on the ellipsoid:

${\displaystyle ds^{2}=\left(N\cos {\phi }\right)^{2}d\lambda ^{2}+M^{2}d\phi ^{2}}$

This result relates the measure of distance on the ellipsoid surface as a function of the spherical coordinate system.

Recall that the purpose of Tissot's indicatrix is to relate how distances on the sphere change when mapped to a planar surface. Specifically, the desired relation is the transform ${\displaystyle {\mathcal {T}}}$ that relates differential distance along the bases of the spherical coordinate system to differential distance along the bases of the Cartesian coordinate system on the planar map. This can be expressed by the relation:

${\displaystyle {\begin{bmatrix}dx\\dy\end{bmatrix}}={\mathcal {T}}{\begin{bmatrix}ds(\lambda ,0)\\ds(0,\phi )\end{bmatrix}}}$

where ${\displaystyle ds(\lambda ,0)}$ and ${\displaystyle ds(0,\phi )}$ represent the computation of ${\displaystyle ds}$ along the longitudinal and latitudinal axes, respectively. Computation of ${\displaystyle ds(\lambda ,0)}$ and ${\displaystyle ds(0,\phi )}$ can be performed directly from the equation above, yielding:

${\displaystyle {\begin{aligned}&ds(\lambda ,0)=N\cos(\phi )d\lambda \\&ds(0,\phi )=Md\phi \end{aligned}}}$

For the purposes of this computation, it is useful to express this relationship as a matrix operation:

${\displaystyle {\begin{bmatrix}d\lambda \\d\phi \end{bmatrix}}=K{\begin{bmatrix}ds(\lambda ,0)\\ds(0,\phi )\end{bmatrix}},\qquad K={\begin{bmatrix}{\frac {1}{N\cos {\phi }}}&0\\0&{\frac {1}{M}}\end{bmatrix}}}$

Now, in order to relate the distances on the ellipsoid surface to those on the plane, we need to relate the coordinate systems. From the chain rule, we can write:

${\displaystyle {\begin{bmatrix}dx\\dy\end{bmatrix}}=J{\begin{bmatrix}d\lambda \\d\phi \end{bmatrix}}}$

where J is the Jacobian matrix:

${\displaystyle J={\begin{bmatrix}{\frac {\partial x}{\partial \lambda }}&{\frac {\partial x}{\partial \phi }}\\{\frac {\partial y}{\partial \lambda }}&{\frac {\partial y}{\partial \phi }}\end{bmatrix}}}$

Plugging in the matrix expression for ${\displaystyle d\lambda }$ and ${\displaystyle d\phi }$ yields the definition of the transform ${\displaystyle {\mathcal {T}}}$ represented by the indicatrix:

${\displaystyle {\begin{bmatrix}dx\\dy\end{bmatrix}}=JK{\begin{bmatrix}ds(\lambda ,0)\\ds(0,\phi )\end{bmatrix}}}$

${\displaystyle {\mathcal {T}}=JK}$

This transform ${\displaystyle {\mathcal {T}}}$ encapsulates the mapping from the ellipsoid surface to the plane. Expressed in this form, SVD can be used to parcel out the important components of the local transformation.

In order to extract the desired distortion information, at any given location in the spherical coordinate system, the values of ${\displaystyle K}$ can be computed directly. The Jacobian, ${\displaystyle J}$, can be computed analytically from the mapping function itself, but it is often simpler to numerically approximate the values at any location on the map using central differences. Once these values are computed, SVD can be applied to each transformation matrix to extract the local distortion information. Remember that, because distortion is local, every location on the map will have its own transformation.

Recall the definition of SVD:

${\displaystyle \mathrm {SVD} ({\mathcal {T}})=U\Lambda V^{T}}$

It is the decomposition of the transformation, ${\displaystyle {\mathcal {T}}}$, into a rotation in the source domain (i.e. the ellipsoid surface), ${\displaystyle V^{T}}$, a scaling along the basis, ${\displaystyle \Lambda }$, and a subsequent second rotation, ${\displaystyle U}$. For understanding distortion, the first rotation is irrelevant, as it rotates the axes of the circle but has no bearing on the final orientation of the ellipse. The next operation, represented by the diagonal singular value matrix, scales the circle along its axes, deforming it to an ellipse. Thus, the singular values represent the scale factors along axes of the ellipse. The first singular value provides the semi-major axis, ${\displaystyle a}$, and the second provides the semi-minor axis, ${\displaystyle b}$, which are the directional scaling factors of distortion. Scale distortion can be computed as the area of the ellipse, ${\displaystyle ab}$, or equivalently by the determinant of ${\displaystyle {\mathcal {T}}}$. Finally, the orientation of the ellipse, ${\displaystyle \theta }$, can be extracted from the first column of ${\displaystyle U}$ as:

${\displaystyle \theta =\arctan \left({\frac {u_{1,0}}{u_{0,0}}}\right)}$