content/developer-guide/annexes/design/JacobianUsage.html

<?xml version="1.0" encoding="UTF-8" ?> <!DOCTYPE html>  <html xmlns="http://www.w3.org/1999/xhtml" xmlns:xi="http://www.w3.org/2001/XInclude" xml:lang="en"> <head> <title>JacobianUsage</title> <meta charset="UTF-8"/> <link rel="stylesheet" type="text/css" href="../../book.css"/> </head> <body>  <section> <header> <h3 id="JacobianUsage">Usages of Jacobian matrix</h3> </header> <p> The Jacobian matrix (or transform derivative) was introduced in the <a href="#TransformDerivative">section about coordinate operations</a>. This section explains how derivatives are used by the Apache <abbr>SIS</abbr> implementation of some high-level operations. </p> <h4 id="DerivativeAndEnvelope">Transform derivatives applied to envelopes</h4> <p> Geographic information systems often need to transform an envelope from one <abbr>CRS</abbr> to another. A <a href="#EnvelopeTransform">previous section</a> has shown why a naive approach transforming the 4 corners is not sufficient. A simple way to reduce the problem could be to sample a large number of points on each envelope border. For example we could sample 40 points on each border, transform them to the target <abbr>CRS</abbr> and compute an envelope containing all the transformed points. But even with 40 points, the point which is closest to the minimum in above figure can still be above that minimum. Consequently we still miss a small portion of the projected shape. It may be tempting to declare this error is negligible, but a single missing pixel can be enough for causing artifacts such as thin black lines in a mosaic or missing “piece of pie” in a projection over a pole. A workaround could be to arbitrarily increase the envelope size by some percentage, but a percentage too high will cause a larger amount of unnecessary data processing (e.g. more data to load from a file) while a percentage too low will leave some artifacts. </p><p> Map projection derivatives offer a more efficient way to resolve this problem. The figure below reproduces the same projected shape than above figure but with exaggerated deformations. The Apache <abbr>SIS</abbr> algorithm projects the 4 corners as in the naive approach, but opportunistically computes also their derivatives (the Jacobian matrices). Between two corners and using derivative information, there is only one possible cubic curve. We can describe that curve by the <i>f(<var>x</var>)</i> = <var>C</var>₀ + <var>C</var>₁<var>x</var> + <var>C</var>₂<var>x</var>² + <var>C</var>₃<var>x</var>³ equation with <var>C</var> coefficients computed from corner positions and derivatives. This approximation (shown by the red curve in above figure) does not match exactly the desired curve (shown by the blue shape in above figure) but is close. We do not use directly the cubic curve minimum, but instead we take the longitude λ (for above example) where this minimum is located. This is shown by the green dashed line in above figure. This position is easy to compute when the <var>C</var> coefficients are known. Assuming that the λ coordinate of the cubic curve minimum is close to the λ coordinates of the projected shape minimum, we can project the point on the envelope border at that location instead of projecting blindly 40 points on that border. </p> <div class="row-of-boxes"> <div> <img src="../../../static/book/images/Approximation.png" alt="Cubic approximation of an envelope bounds"/> </div><div> Legend: <ul> <li><b>Blue:</b> the geometric shape of the envelope after projection. This is the shape from which to get a new envelope.</li> <li><b>Red</b> (with hash): The <var>y</var> = <var>C</var>₀ + <var>C</var>₁λ + <var>C</var>₂λ² + <var>C</var>₃λ³ approximation.</li> <li><b>Green</b> (dashed line): Position λ<sub>m</sub> of approximation minimum, obtained by resolving 0 = <var>C</var>₁ + 2<var>C</var>₂λ<sub>m</sub> + 3<var>C</var>₃λ<sub>m</sub>². The same cubic line can have two minimums.</li> </ul> </div> </div> <p> In practice Apache <abbr>SIS</abbr> uses 8 points: the 4 corners plus one point at the center of each border of the envelope to project. The center points are added as a safety for map projections having symmetric deformations above and below equator. According our tests, using only those 8 points together with derivatives as described above gives more accurate results than “brute force” approach with 160 points on the 4 envelope borders. Saving 150 points seems small compared to computer performances. But above discussion used a two-dimensional envelope as an example. The same discussion applies to <var>n</var>-dimensional envelopes as well. Apache <abbr>SIS</abbr> can apply this algorithm on envelopes with any number of dimensions up to 64. The performance gain offered by this algorithm compared to “brute force” approach increases exponentially with the number of dimensions. </p><p> Apache <abbr>SIS</abbr> implements the algorithm described in this section by the static <code>Envelopes.transform(CoordinateOperation, Envelope)</code> method. An alternative <code>Envelopes.transform(MathTransform, Envelope)</code> method is also available, but should be used only when the <code>CoordinateOperation</code> is unknown. The reason is because <code>MathTransform</code> objects does not contain any information about coordinate system axes, which prevent the <code>Envelopes.transform(…)</code> method to handle special cases such as envelopes containing a pole. </p> <h4 id="DerivativeAndRaster">Transform derivatives applied to rasters</h4> <p> A raster can be projected by preparing an initially empty raster which will contain the resampled pixel values. Then for each pixel in the <em>destination</em> raster, we use the <em>inverse</em> of the map projection (<var>T</var>⁻¹) for computing the coordinates of the corresponding pixel in source raster. The transformed coordinates may not fall at the center of a source raster pixel, so source pixel values usually need to be interpolated. </p> <div style="display:flex; flex-direction:column; align-items:center"> <div style="display:flex; justify-content:space-between; width:564px"> <div class="caption">Source image</div> <div class="caption">Destination image</div> </div> <img src="../../../static/book/images/RasterProjection.png" alt="Raster reprojection"/> </div> <p> However applying the inverse transform on all pixel coordinates can be relatively slow. We can accelerate raster reprojections by pre-computing an <dfn>interpolation grid</dfn>. This grid contains the coordinate values of inverse transform <var>T</var>⁻¹ for only a few points. The coordinate values of other points are obtained by bilinear interpolations between interpolation grid points. Historically, this algorithm was implemented in the <code>WarpGrid</code> object of <cite>Java Advanced Imaging</cite> library. The performance gain is yet better if the same interpolation grid is reused for many rasters having their pixels at the same coordinates. </p><p> A difficulty in using interpolation grids is to determine how many points we need for keeping errors (defined as the difference between interpolated positions and actual positions computed by <var>T</var>⁻¹) below some threshold value (for example ¼ of pixel size). A “brute force” approach is to begin with a grid of size 3×3, then increase the number of points iteratively: </p> <div class="row-of-boxes"> <div> <img src="../../../static/book/images/WarpGrid.png" alt="Warp grid"/> </div><div> Legend: <ul> <li><b>Blue dots:</b> first iteration (9 points).</li> <li><b>Green dots:</b> second iteration (25 points, including 16 news).</li> <li><b>Red dots:</b> third iteration (81 points, including 56 news).</li> </ul> Continuing… <ul> <li>Forth iteration: 289 points, including 208 news.</li> <li>Fifth iteration: 1089 points, including 800 news.</li> <li>Sixth iteration: 4225 points, including 3136 news.</li> <li>…</li> </ul> </div> </div> <p> We could stop dividing the grid when, after having computed new points, we verified that the differences between coordinates computed by <var>T</var>⁻¹ during last iteration and coordinates computed by bilinear interpolations for the same points are lower than the threshold value. Unfortunately this approach can only tell us <strong>after</strong> computing new points… that those new points were unnecessary. This is unfortunate considering that each iteration adds an amount of points approximately equal to 3 times the <em>sum</em> of the number of points of all previous iterations. </p><p> Map projection derivatives help to improve the speed of interpolation grid computations. Using the derivatives, we can estimate whether another iteration is necessary <strong>before</strong> to execute it. The basic idea is to check if the derivatives at two neighbor points define two tangent lines that are almost parallel. In such case, we assume that the transform between those two points is almost linear. In order to define numerically the meaning of “almost linear”, we compute the intersection of the two tangent lines as illustrated below (the blue arrows). Then we compute the distance between that intersection and the straight line connecting the two points (the dashed line in the figure below). </p> <p style="text-align:center"><img style="border: solid 1px" src="../../../static/book/images/IntersectionOfDerivatives.png" alt="Intersection of derivatives"/></p> <p> In “brute force” algorithm without derivatives, iteration stops when the distance between interpolated positions (blue dashed line) and actual projected positions (red curve) is less than the threshold value. This criteria requires to compute the projected positions before the decision can be made. But with an algorithm using derivatives, the projected positions (red curve) are replaced by the tangents intersection point. If the actual projected curve does not have too much deformations (which should be the case for pairs of neighbor points close enough), then the red curve stay between the blue dashed line and the tangents intersection point. By using that intersection point, we reduce greatly the number of points to transform (3× the sum of all previous iterations). </p> <h4 id="GetDerivative">Getting the derivative at a point</h4> <p> The algorithms discussed in previous section would be irrelevant if map projection derivatives were costly to compute. But analytic derivations of their formulas show that map projection and derivative formulas have many terms in common. Furthermore map projection formulas are simplified by Apache <abbr>SIS</abbr> implementation strategy, which takes out linear terms (delegated to affine transforms) so that <code>NormalizedProjection</code> subclasses can focus on the non-linear terms only. Map projection implementations in Apache <abbr>SIS</abbr> exploit those characteristics for computing Jacobian matrices (if required) together with projected points in a single operation, reusing many common terms between the two set of formulas. Example:</p> <pre><code>AbstractMathTransform projection = ...; // An Apache SIS map projection. double[] sourcePoint = {longitude, latitude}; // The geographic coordinate to project. double[] targetPoint = new double[2]; // Where to store the projection result. Matrix derivative = projection.<code class="SIS">transform</code>(sourcePoint, 0, targetPoint, 0, true);</code></pre> <p> If only the Jacobian matrix is desired (without the projected point), then the <code>MathTransform.derivative(DirectPosition)</code> method offers a more readable alternative. </p><p> Many operations on coordinate transforms can also be applied on transform derivatives: concatenation of a chain of transforms, inversion, taking a subset of the dimensions, <i>etc.</i> reverse projection formulas are usually more complicated than forward projections formulas, but thankfully their derivatives do not need to be computed: the Jacobian matrix of an inverse transform is the inverse of the Jacobian matrix of the forward transform. Calculation of reverse projections can be implemented as below: </p> <pre><code>@Override public Matrix derivative(DirectPosition p) throws TransformException { Matrix jac = inverse().derivative(transform(p)); return Matrices.inverse(jac); }</code></pre> </section> </body> </html>

content/developer-guide/annexes/design/JacobianUsage.html (197 lines of code) (raw):