This is an automated email from the ASF dual-hosted git repository. desruisseaux pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit eeb214ce1ad8bbedb30a10689298823da13cfd65 Author: Martin Desruisseaux <martin.desruisse...@geomatys.com> AuthorDate: Tue Apr 30 12:01:22 2024 +0200 Rewrite the Equidistant Cylindrical equations using Clenshaw summation. --- .../projection/EquidistantCylindrical.java | 151 +++++++++++++++------ .../apache/sis/referencing/ClenshawSummation.java | 44 +++++- 2 files changed, 146 insertions(+), 49 deletions(-) diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/EquidistantCylindrical.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/EquidistantCylindrical.java index de2edd02df..00d7883f4e 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/EquidistantCylindrical.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/EquidistantCylindrical.java @@ -25,6 +25,7 @@ import org.opengis.parameter.ParameterDescriptor; import org.opengis.referencing.operation.Matrix; import org.opengis.referencing.operation.MathTransform; import org.opengis.referencing.operation.OperationMethod; +import org.apache.sis.util.Debug; import org.apache.sis.util.Workaround; import org.apache.sis.util.privy.DoubleDouble; import org.apache.sis.parameter.Parameters; @@ -52,17 +53,22 @@ public class EquidistantCylindrical extends NormalizedProjection { */ private static final long serialVersionUID = 6598912212024442670L; + /** + * Whether to allow the use of Clenshaw summation. This flag should be {@code true} for performance reason, + * as it reduces the number of calls to trigonometric functions. A value of {@code false} causes this class + * to use the formulas as published by EPSG guidance notes instead, which may be useful if we suspect a bug + * in our application of Clenshaw summation. + */ + @Debug + private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true; + /** * Coefficients for the forward (f) and inverse (i) projection. * This is used in a trigonometric series with multiple angles. - * - * <h4>Missed optimization</h4> - * We should replace the multiple angles by power polynomials using the Clenshaw summation algorithm. - * However the {@code org.apache.sis.referencing.ClenshawSummation} class (in test packages) that we - * used for this purpose in other map projections is limited to 6 terms, and we have 7 terms here. */ - private final double cf0, cf2, cf4, cf6, cf8, cf10, cf12, cf14, - ci2, ci4, ci6, ci8, ci10, ci12, ci14; + private final double c0, + cf1, cf2, cf3, cf4, cf5, cf6, cf7, + ci1, ci2, ci3, ci4, ci5, ci6, ci7; /** * Creates an Equidistant Cylindrical projection from the given parameters. @@ -99,33 +105,71 @@ public class EquidistantCylindrical extends NormalizedProjection { final DoubleDouble sx = initializer.rν2(sin(φ1)).sqrt().multiply(cos(φ1), false); final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION); denormalize.convertBefore(0, sx, null); - - // Coefficients for the forward projection. - final double e2 = eccentricitySquared; - final double e4 = e2*e2; - final double e6 = e2*e4; - final double e8 = e4*e4; - cf0 = 1 - e2*(1./4 + e2*( 3./64 + e2*( 5./256 + e2*(175./16384 + e2*( 441./65536 + e2*( 4851./1048576 + e2*( 14157./4194304 ))))))); - cf2 = - e2*(3./8 + e2*( 3./32 + e2*(45./1024 + e2*(105./4096 + e2*(2205./131072 + e2*( 6237./524288 + e2*(297297./33554432))))))); - cf4 = e4*(15./256 + e2*(45./1024 + e2*(525./16384 + e2*(1575./65536 + e2*(155925./8388608 + e2*(495495./33554432)))))); - cf6 = - e6*(35./3072 + e2*(175./12288 + e2*(3675./262144 + e2*( 13475./1048576 + e2*(385385./33554432))))); - cf8 = + e8*(315./131072 + e2*(2205./524288 + e2*( 43659./8388608 + e2*(189189./33554432)))); - cf10 = - (e4*e6)*( 693./1310720 + e2*( 6237./5242880 + e2*(297297./167772160))); - cf12 = (e6*e6)*( 1001./8388608 + e2*( 11011./33554432)); - cf14 = - (e6*e8)*( 6435./234881024); - - // Coefficients for the inverse projection. - final double n = initializer.axisLengthRatio().ratio_1m_1p().doubleValue(); - final double n2 = n*n; - final double n3 = n*n2; + /* + * Coefficients for the forward projection. The formulas are available in two variants: + * as published by EPSG guidance notes, or modified with Clenshaw summation for reducing + * the number of calls to trigonometric functions. Coefficients for the modified form were + * computed by `org.apache.sis.referencing.ClenshawSummation.equidistantCylindrical(boolean)` + * in the test packages. + */ + final double e2 = eccentricitySquared; + final double e4 = e2*e2; + final double e6 = e2*e4; + final double e8 = e4*e4; + final double e10 = e4*e6; + final double e12 = e6*e6; + final double e14 = e6*e8; + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { + // Formulas modified with Clenshaw summation. + c0 = ( -14157./ 4194304)*e14 + (-4851./1048576)*e12 + ( -441./ 65536)*e10 + (-175./16384)*e8 + ( -5./256)*e6 + (-3./ 64)*e4 + (-1./4)*e2 + 1; + cf1 = ( 16159./18350080)*e14 + ( -77./ 327680)*e12 + ( -273./ 81920)*e10 + ( -35./ 3072)*e8 + (-25./768)*e6 + (-3./ 32)*e4 + (-3./8)*e2; + cf2 = ( 75075./ 8388608)*e14 + (35805./2097152)*e12 + ( 4095./131072)*e10 + (1785./32768)*e8 + ( 45./512)*e6 + (15./128)*e4; + cf3 = (-464893./18350080)*e14 + (-6083./ 163840)*e12 + (-2037./ 40960)*e10 + (-175./ 3072)*e8 + (-35./768)*e6; + cf4 = ( 145145./ 4194304)*e14 + (39655./1048576)*e12 + ( 2205./ 65536)*e10 + ( 315./16384)*e8; + cf5 = (-480051./18350080)*e14 + (-6237./ 327680)*e12 + ( -693./ 81920)*e10; + cf6 = ( 11011./ 1048576)*e14 + ( 1001./ 262144)*e12; + cf7 = ( -6435./ 3670016)*e14; + } else { + // Formulas as published by EPSG guidance notes. + c0 = 1 + (-1./4)*e2 + (-3./ 64)*e4 + ( -5./ 256)*e6 + (-175./ 16384)*e8 + ( -441./ 65536)*e10 + ( -4851./1048576)*e12 + ( -14157./ 4194304)*e14; + cf1 = + (-3./8)*e2 + (-3./ 32)*e4 + (-45./1024)*e6 + (-105./ 4096)*e8 + (-2205./ 131072)*e10 + ( -6237./ 524288)*e12 + (-297297./ 33554432)*e14; + cf2 = (15./256)*e4 + ( 45./1024)*e6 + ( 525./ 16384)*e8 + ( 1575./ 65536)*e10 + (155925./8388608)*e12 + ( 495495./ 33554432)*e14; + cf3 = (-35./3072)*e6 + (-175./ 12288)*e8 + (-3675./ 262144)*e10 + (-13475./1048576)*e12 + (-385385./ 33554432)*e14; + cf4 = + ( 315./131072)*e8 + ( 2205./ 524288)*e10 + ( 43659./8388608)*e12 + ( 189189./ 33554432)*e14; + cf5 = ( -693./1310720)*e10 + ( -6237./5242880)*e12 + (-297297./167772160)*e14; + cf6 = ( 1001./8388608)*e12 + ( 11011./ 33554432)*e14; + cf7 = ( -6435./234881024)*e14; + } + /* + * Coefficients for the inverse projection, available in two variants. + * See the forward case above. + */ + final double n1 = initializer.axisLengthRatio().ratio_1m_1p().doubleValue(); + final double n2 = n1*n1; + final double n3 = n1*n2; final double n4 = n2*n2; - ci2 = n*(3./2 + n2*(-27./32 + n2*(269./512 + n2*( -6607./24576)))); - ci4 = n2*( 21./16 + n2*( -55./32 + n2*( 6759./4096))); - ci6 = n3*(151./96 + n2*(-417./128 + n2*( 87963./20480))); - ci8 = n4*(1097./512 + n2*( -15543./2560)); - ci10 = (n3*n2)*(8011./2560 + n2*( -69119./6144)); - ci12 = (n3*n3)*( 293393./61440); - ci14 = (n3*n4)*(6845701./860160); + final double n5 = n2*n3; + final double n6 = n3*n3; + final double n7 = n3*n4; + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { + // Formulas modified with Clenshaw summation. + ci1 = (-5112013./215040)*n7 + ( 553./ 80)*n5 + (-29./12)*n3 + (3./2)*n1; + ci2 = ( 143969./ 2560)*n6 + ( -1537./128)*n4 + ( 21./ 8)*n2; + ci3 = ( 3074943./ 8960)*n7 + (-32373./640)*n5 + (151./24)*n3; + ci4 = ( -386651./ 1920)*n6 + ( 1097./ 64)*n4; + ci5 = (-2927011./ 3584)*n7 + ( 8011./160)*n5; + ci6 = ( 293393./ 1920)*n6; + ci7 = ( 6845701./ 13440)*n7; + } else { + // Formulas as published by EPSG guidance notes. + ci1 = (3./2)*n1 + (-27./32)*n3 + ( 269./ 512)*n5 + ( -6607./ 24576)*n7; + ci2 = ( 21./16)*n2 + ( -55./ 32)*n4 + ( 6759./ 4096)*n6; + ci3 = (151./96)*n3 + (-417./ 128)*n5 + ( 87963./ 20480)*n7; + ci4 = (1097./ 512)*n4 + ( -15543./ 2560)*n6; + ci5 = (8011./2560)*n5 + ( -69119./ 6144)*n7; + ci6 = ( 293393./ 61440)*n6; + ci7 = (6845701./860160)*n7; + } } /** @@ -161,17 +205,32 @@ public class EquidistantCylindrical extends NormalizedProjection { final boolean derivate) throws ProjectionException { final double φ = srcPts[srcOff+1]; + final double sinθ = sin(2*φ); + final double cosθ = cos(2*φ); if (dstPts != null) { - dstPts[dstOff] = srcPts[srcOff]; - dstPts[dstOff+1] = cf14*sin(14*φ) + cf12*sin(12*φ) + cf10*sin(10*φ) + cf8*sin(8*φ) - + cf6*sin( 6*φ) + cf4*sin( 4*φ) + cf2*sin( 2*φ) + cf0*φ; + final double y; + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { + y = ((((((cf7*cosθ + cf6)*cosθ + cf5)*cosθ + cf4)*cosθ + cf3)*cosθ + cf2)*cosθ + cf1)*sinθ + c0*φ; + } else { + y = cf7*sin(14*φ) + cf6*sin(12*φ) + cf5*sin(10*φ) + cf4*sin(8*φ) + + cf3*sin( 6*φ) + cf2*sin( 4*φ) + cf1*sinθ + c0*φ; + } + dstPts[dstOff ] = srcPts[srcOff]; + dstPts[dstOff+1] = y; } if (!derivate) { return null; } final var derivative = new Matrix2(); - derivative.m11 = cf14*cos(14*φ)*14 + cf12*cos(12*φ)*12 + cf10*cos(10*φ)*10 + cf8*cos(8*φ)*8 - + cf6*cos( 6*φ)*6 + cf4*cos( 4*φ)*4 + cf2*cos( 2*φ)*2 + cf0; + final double d; + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { + d = (((((( cf7*cosθ + cf6)*cosθ + cf5)*cosθ + cf4)*cosθ + cf3)*cosθ + cf2)*cosθ + cf1)*cosθ + - (((((6*cf7*cosθ + 5*cf6)*cosθ + 4*cf5)*cosθ + 3*cf4)*cosθ + 2*cf3)*cosθ + cf2)*(sinθ*sinθ); + } else { + d = 7*cf7*cos(14*φ) + 6*cf6*cos(12*φ) + 5*cf5*cos(10*φ) + 4*cf4*cos(8*φ) + + 3*cf3*cos( 6*φ) + 2*cf2*cos( 4*φ) + cf1*cosθ; + } + derivative.m11 = 2*d + c0; return derivative; } @@ -186,9 +245,17 @@ public class EquidistantCylindrical extends NormalizedProjection { final double[] dstPts, final int dstOff) throws ProjectionException { - final double μ = srcPts[srcOff+1] / cf0; - dstPts[dstOff] = srcPts[srcOff]; - dstPts[dstOff+1] = ci14*sin(14*μ) + ci12*sin(12*μ) + ci10*sin(10*μ) + ci8*sin(8*μ) - + ci6*sin( 6*μ) + ci4*sin( 4*μ) + ci2*sin( 2*μ) + μ; + final double μ = srcPts[srcOff+1] / c0; + final double sinθ = sin(2*μ); + final double cosθ = cos(2*μ); + final double y; + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { + y = ((((((ci7*cosθ + ci6)*cosθ + ci5)*cosθ + ci4)*cosθ + ci3)*cosθ + ci2)*cosθ + ci1)*sinθ; + } else { + y = ci7*sin(14*μ) + ci6*sin(12*μ) + ci5*sin(10*μ) + ci4*sin(8*μ) + + ci3*sin( 6*μ) + ci2*sin( 4*μ) + ci1*sinθ; + } + dstPts[dstOff ] = srcPts[srcOff]; + dstPts[dstOff+1] = y + μ; } } diff --git a/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/ClenshawSummation.java b/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/ClenshawSummation.java index 2831ca0cd3..b8f57cd585 100644 --- a/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/ClenshawSummation.java +++ b/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/ClenshawSummation.java @@ -159,16 +159,16 @@ public final class ClenshawSummation { */ final void appendTo(final StringBuilder b) { boolean more = false; - for (int i=terms.length; --i >= 0;) { - final Term t = terms[i]; + for (int p = terms.length; --p >= 0;) { // `p` is power minus one. + final Term t = terms[p]; if (t != null) { if (more) { b.append(" + "); } t.appendTo(b); - b.append(" * η"); - if (i != 0) { - b.append(i+1); + b.append("*η"); + if (p != 0) { + b.append(p+1); } more = true; } @@ -299,9 +299,9 @@ public final class ClenshawSummation { } /** - * Returns a string representation for debugging purpose. + * Returns a string representation of this object, for debugging purposes or for showing the result. * - * @return a string representation of this object before summation. + * @return a string representation of this object before or after summation. */ @Override public String toString() { @@ -574,6 +574,36 @@ public final class ClenshawSummation { ); } + /** + * Coefficients for Equidistant Cylindrical (EPSG:1028) map projection. + * + * @param inverse {@code false} for the forward projection, or {@code true} for the inverse projection. + * + * @return Equidistant Cylindrical equation. + */ + public static ClenshawSummation equidistantCylindrical(final boolean inverse) { + if (inverse) { + return new ClenshawSummation( + new Coefficient(t(3, 2), null, t(-27, 32), null, t(269, 512), null, t(-6607, 24576)), + new Coefficient(null, t(21, 16), null, t( -55, 32), null, t(6759, 4096)), + new Coefficient(null, null, t(151, 96), null, t(-417, 128), null, t(87963, 20480)), + new Coefficient(null, null, null, t(1097, 512), null, t(-15543, 2560)), + new Coefficient(null, null, null, null, t(8011, 2560), null, t(-69119, 6144)), + new Coefficient(null, null, null, null, null, t( 293393, 61440)), + new Coefficient(null, null, null, null, null, null, t(6845701, 860160)) + ); + } + return new ClenshawSummation( + new Coefficient(t(-3, 8), t(-3, 32), t(-45, 1024), t(-105, 4096), t(-2205, 131072), t( -6237, 524288), t(-297297, 33554432)), + new Coefficient(null, t(15, 256), t( 45, 1024), t( 525, 16384), t( 1575, 65536), t( 155925, 8388608), t(495495, 33554432)), + new Coefficient(null, null, t(-35, 3072), t(-175, 12288), t(-3675, 262144), t( -13475, 1048576), t(-385385, 33554432)), + new Coefficient(null, null, null, t( 315, 131072), t( 2205, 524288), t( 43659, 8388608), t(189189, 33554432)), + new Coefficient(null, null, null, null, t(-693, 1310720), t(-6237, 5242880), t(-297297, 167772160)), + new Coefficient(null, null, null, null, null, t( 1001, 8388608), t( 11011, 33554432)), + new Coefficient(null, null, null, null, null, null, t( -6435, 234881024)) + ); + } + /** * Computes the Clenshaw summation of a series and display the code to standard output. * This method has been used for generating the Java code implemented in methods such as