When computing inverse projection, temporarily disable the setting of out-of-domain coordinates to NaN.
We use the Inverter class as a flag for identifying those situations.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
index 8504238..ff93401 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
@@ -30,7 +30,7 @@
* └ ┘</pre></blockquote>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.4
+ * @version 1.1
*
* @see Matrix2
* @see Matrix3
@@ -39,8 +39,7 @@
* @since 0.4
* @module
*/
-@SuppressWarnings("CloneableClassWithoutClone") // No field in this class needs clone.
-public final class Matrix1 extends MatrixSIS {
+public class Matrix1 extends MatrixSIS {
/**
* Serial number for inter-operability with different versions.
*/
@@ -260,6 +259,14 @@
}
/**
+ * {@inheritDoc}
+ */
+ @Override
+ public Matrix1 clone() {
+ return (Matrix1) super.clone();
+ }
+
+ /**
* Returns {@code true} if the specified object is of type {@code Matrix1} and
* all of the data members are equal to the corresponding data members in this matrix.
*
@@ -268,7 +275,7 @@
*/
@Override
public boolean equals(final Object object) {
- if (object instanceof Matrix1) {
+ if (object != null && object.getClass() == getClass()) {
final Matrix1 that = (Matrix1) object;
return Numerics.equals(m00, that.m00);
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
index 04e0f51..c3f455a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
@@ -31,7 +31,7 @@
* └ ┘</pre></blockquote>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.4
+ * @version 1.1
*
* @see Matrix1
* @see Matrix3
@@ -40,8 +40,7 @@
* @since 0.4
* @module
*/
-@SuppressWarnings("CloneableClassWithoutClone") // No field in this class needs clone.
-public final class Matrix2 extends MatrixSIS {
+public class Matrix2 extends MatrixSIS {
/**
* Serial number for inter-operability with different versions.
*/
@@ -284,6 +283,14 @@
}
/**
+ * {@inheritDoc}
+ */
+ @Override
+ public Matrix2 clone() {
+ return (Matrix2) super.clone();
+ }
+
+ /**
* Returns {@code true} if the specified object is of type {@code Matrix2} and
* all of the data members are equal to the corresponding data members in this matrix.
*
@@ -292,7 +299,7 @@
*/
@Override
public boolean equals(final Object object) {
- if (object instanceof Matrix2) {
+ if (object != null && object.getClass() == getClass()) {
final Matrix2 that = (Matrix2) object;
return Numerics.equals(this.m00, that.m00) &&
Numerics.equals(this.m01, that.m01) &&
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java
index 6edfb8f..d246bc0 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix3.java
@@ -30,7 +30,7 @@
* └ ┘</pre></blockquote>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.6
+ * @version 1.1
*
* @see Matrix1
* @see Matrix2
@@ -39,8 +39,7 @@
* @since 0.4
* @module
*/
-@SuppressWarnings("CloneableClassWithoutClone") // No field in this class needs clone.
-public final class Matrix3 extends MatrixSIS {
+public class Matrix3 extends MatrixSIS {
/**
* Serial number for inter-operability with different versions.
*/
@@ -298,4 +297,12 @@
swap = m02; m02 = m20; m20 = swap;
swap = m12; m12 = m21; m21 = swap;
}
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public Matrix3 clone() {
+ return (Matrix3) super.clone();
+ }
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java
index 8093b87..2f14f39 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix4.java
@@ -31,7 +31,7 @@
* └ ┘</pre></blockquote>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.6
+ * @version 1.1
*
* @see Matrix1
* @see Matrix2
@@ -40,8 +40,7 @@
* @since 0.4
* @module
*/
-@SuppressWarnings("CloneableClassWithoutClone") // No field in this class needs clone.
-public final class Matrix4 extends MatrixSIS {
+public class Matrix4 extends MatrixSIS {
/**
* Serial number for inter-operability with different versions.
*/
@@ -336,4 +335,12 @@
swap = m13; m13 = m31; m31 = swap;
swap = m23; m23 = m32; m32 = swap;
}
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public Matrix4 clone() {
+ return (Matrix4) super.clone();
+ }
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java
index 0005bb1..db87112 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java
@@ -69,7 +69,7 @@
* Like SIS, Vecmath is optimized for small matrices of interest for 2D and 3D graphics.</p>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.8
+ * @version 1.1
* @since 0.4
* @module
*/
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Inverter.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Inverter.java
new file mode 100644
index 0000000..d928119
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Inverter.java
@@ -0,0 +1,105 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
+
+import static java.lang.Math.abs;
+
+
+/**
+ * A temporary Jacobian matrix where to write the derivative of a forward projection.
+ * This Jacobian matrix is used for calculation of inverse projection when no inverse
+ * formulas is available, or when the inverse formula is too approximate (for example
+ * because eccentricity is too high). This class processes as below:
+ *
+ * <ol>
+ * <li>Given a first estimation of (λ,φ), compute the forward projection (x′,y′)
+ * for that estimation together with the Jacobian matrix at that position.</li>
+ * <li>Compute the errors compared to the specified (x,y) values.</li>
+ * <li>Convert that (Δx,Δy) error into a (Δλ,Δφ) error using the inverse of the Jacobian matrix.</li>
+ * <li>Correct (λ,φ) and continue iteratively until the error is small enough.</li>
+ * </ol>
+ *
+ * This algorithm described in EPSG guidance note for {@link Orthographic} projection
+ * but can be applied to any map projection, not only orthographic.
+ *
+ * <p>This algorithm is defined in a {@link Matrix2} subclass for allowing map projection
+ * implementations to use {@code if (derivative instanceof Inverter)} check for detecting
+ * when a {@code transform} method is invoked for the purpose of an inverse projection.</p>
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ *
+ * @ee <a href="https://issues.apache.org/jira/browse/SIS-478">SIS-478</a>
+ *
+ * @since 1.1
+ * @module
+ */
+@SuppressWarnings({"CloneableClassWithoutClone", "serial"})
+final class Inverter extends Matrix2 {
+ /**
+ * Creates a new matrix initialized to identity.
+ */
+ Inverter() {
+ }
+
+ /**
+ * Computes the inverse of the given projection. Before this method is invoked,
+ * the {@code dstPts[dstOff]} and {@code dstPts[dstOff+1]} array elements must
+ * contain the initial λ and φ estimations respectively. After method returns,
+ * the same array elements contain the refined λ and φ values.
+ *
+ * <p>Note: restricted to {@link Orthographic} projection for now,
+ * but may be generalized to any projection in a future version.</p>
+ *
+ * @param projection the forward projection for which to compute an inverse projection.
+ * @param x the easting value from {@code srcPts[srcOff]}.
+ * @param y the northing value from {@code srcPts[srcOff+1]}.
+ * @param dstPts the array where to refine the (λ,φ) values.
+ * @param dstOff index of the λ element in the {@code dstPts} array.
+ * @throws ProjectionException if the iterative algorithm does not converge.
+ */
+ final void inverseTransform(final Orthographic projection, final double x, final double y,
+ final double[] dstPts, final int dstOff) throws ProjectionException
+ {
+ double λ = dstPts[dstOff ];
+ double φ = dstPts[dstOff+1];
+ for (int it=NormalizedProjection.MAXIMUM_ITERATIONS; --it >= 0;) {
+ final double cosφ = projection.transform(dstPts, dstOff, dstPts, dstOff, this);
+ final double ΔE = x - dstPts[dstOff ];
+ final double ΔN = y - dstPts[dstOff+1];
+ final double D = (m01 * m10) - (m00 * m11); // Determinant.
+ final double Δλ = (m01*ΔN - m11*ΔE) / D;
+ final double Δφ = (m10*ΔE - m00*ΔN) / D;
+ dstPts[dstOff ] = λ += Δλ;
+ dstPts[dstOff+1] = φ += Δφ;
+ /*
+ * Following condition uses ! for stopping iteration on NaN values.
+ * We do not use Math.max(…) because we want to check NaN separately
+ * for φ and λ.
+ */
+ if (!(abs(Δφ) > NormalizedProjection.ANGULAR_TOLERANCE ||
+ abs(Δλ*cosφ) > NormalizedProjection.ANGULAR_TOLERANCE))
+ {
+ return;
+ }
+ }
+ throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
+ }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java
index 9e3d024..87b9705 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java
@@ -25,7 +25,6 @@
import org.apache.sis.referencing.operation.matrix.Matrix2;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.transform.ContextualParameters;
-import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.Workaround;
@@ -150,13 +149,16 @@
* <div class="note"><b>Implementation note:</b>
* in other map projections, we use a different class for ellipsoidal formulas.
* But the orthographic projection is a bit different; for this one it is more
- * convenient to use a flag.</div>
+ * convenient to use {@code if} statements.</div>
*
+ * @param derivative where to store the Jacobian matrix, or {@code null} if none.
+ * If this matrix is an {@link Inverter} instance, we take that as a flag
+ * meaning to not set out-of-range results to NaN.
* @return {@code cos(φ)}, useful for error tolerance check.
*/
- private double transform(final double[] srcPts, final int srcOff,
- final double[] dstPts, final int dstOff,
- final Matrix2 derivative)
+ final double transform(final double[] srcPts, final int srcOff,
+ final double[] dstPts, final int dstOff,
+ final Matrix2 derivative)
{
final double λ = srcPts[srcOff ];
final double φ = srcPts[srcOff+1];
@@ -174,7 +176,7 @@
* then the point is on the opposite hemisphere and should be clipped.
*/
double rν;
- if (cosc >= 0) {
+ if (cosc >= 0 || derivative instanceof Inverter) {
x = cosφ * sinλ; // Snyder (20-3)
y = mℯ2_cosφ0*sinφ - sinφ0*cosφ_cosλ; // Snyder (20-4) × (1 – ℯ²)
rν = 1;
@@ -203,7 +205,7 @@
derivative.m00 = cosφ_cosλ / rν; // ∂E/∂λ
derivative.m01 = -sinφ*sinλ * ρ; // ∂E/∂φ
derivative.m10 = sinφ0 * x; // ∂N/∂λ
- derivative.m11 = (cosφ0 * cosφ + sinφ0_sinφ*cosλ) * ρ; // ∂N/∂φ
+ derivative.m11 = (cosφ0*cosφ + sinφ0_sinφ*cosλ) * ρ; // ∂N/∂φ
}
return cosφ;
}
@@ -238,29 +240,9 @@
*
* See https://issues.apache.org/jira/browse/SIS-478
*/
- final Matrix2 J = new Matrix2(); // Jacobian matrix.
- double λ = dstPts[dstOff ];
- double φ = dstPts[dstOff+1];
- for (int it=0; it<MAXIMUM_ITERATIONS; it++) {
- final double cosφ = transform(dstPts, dstOff, dstPts, dstOff, J);
- final double ΔE = x - dstPts[dstOff ];
- final double ΔN = y - dstPts[dstOff+1];
- final double D = (J.m01 * J.m10) - (J.m00 * J.m11); // Determinant.
- final double Δλ = (J.m01*ΔN - J.m11*ΔE) / D;
- final double Δφ = (J.m10*ΔE - J.m00*ΔN) / D;
- dstPts[dstOff ] = λ += Δλ;
- dstPts[dstOff+1] = φ += Δφ;
- /*
- * Following condition uses ! for stopping iteration on NaN values.
- * We do not use Math.max(…) because we want to check NaN separately
- * for φ and λ.
- */
- if (!(abs(Δφ) > ANGULAR_TOLERANCE || abs(cosφ*abs(Δλ)) > ANGULAR_TOLERANCE)) {
- dstPts[dstOff] = IEEEremainder(dstPts[dstOff], PI);
- return;
- }
- }
- throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
+ final Inverter j = new Inverter(); // Jacobian matrix.
+ j.inverseTransform(this, x, y, dstPts, dstOff);
+ dstPts[dstOff] = IEEEremainder(dstPts[dstOff], PI);
}
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/OrthographicTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/OrthographicTest.java
index 42fc7c5..434d55b 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/OrthographicTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/OrthographicTest.java
@@ -16,6 +16,7 @@
*/
package org.apache.sis.referencing.operation.projection;
+import org.opengis.util.FactoryException;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.referencing.operation.transform.CoordinateDomain;
import org.apache.sis.internal.referencing.provider.MapProjection;
@@ -144,4 +145,18 @@
verifyInDomain(CoordinateDomain.GEOGRAPHIC_RADIANS_SOUTH, 753524735);
verifyDerivative(toRadians(5), toRadians(-85));
}
+
+ /**
+ * Tests the <cite>"Orthographic"</cite> (EPSG:9840) projection method.
+ * This test is defined in GeoAPI conformance test suite.
+ *
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a coordinate.
+ *
+ * @see org.opengis.test.referencing.ParameterizedTransformTest#testOrthographic()
+ */
+ @Test
+ public void runGeoapiTest() throws FactoryException, TransformException {
+ createGeoApiTest(new org.apache.sis.internal.referencing.provider.Orthographic()).testOrthographic();
+ }
}
