Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java?rev=1707296&r1=1707295&r2=1707296&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java [UTF-8] Wed Oct 7 13:20:26 2015 @@ -57,56 +57,32 @@ import static org.apache.sis.math.MathFu * @author Martin Desruisseaux (Geomatys) * @author Rémi Maréchal (Geomatys) * @since 0.6 - * @version 0.6 + * @version 0.7 * @module * * @see Mercator * @see ObliqueMercator */ -public class TransverseMercator extends NormalizedProjection { +public class TransverseMercator extends ConformalProjection { /** * For cross-version compatibility. */ - private static final long serialVersionUID = -4717976245811852528L; + private static final long serialVersionUID = 8167193470189463354L; /** - * Whether to use the original formulas a published by EPSG, or their form modified using trigonometric identities. - * The modified form uses trigonometric identifies for reducing the amount of calls to the {@link Math#sin(double)} - * and similar method. The identities used are: + * Coefficients in the series expansion of the forward projection, + * depending only on {@linkplain #excentricity excentricity} value. + * The series expansion is of the following form: * - *
cf₂⋅f(2θ) + cf₄⋅f(4θ) + cf₆⋅f(6θ) + cf₈⋅f(8θ)* - * Hyperbolic formulas: - *
Consider those fields as final! + * They are not final only for the purpose of {@link #computeCoefficients()}.
*/ - private final double h1, h2, h3, h4, ih1, ih2, ih3, ih4; + private transient double cf2, cf4, cf6, cf8; /** * Creates a Transverse Mercator projection from the given parameters. @@ -159,7 +135,6 @@ public class TransverseMercator extends * Opportunistically use double-double arithmetic for computation of B since we will store * it in the denormalization matrix, and there is no sine/cosine functions involved here. */ - final double n; final DoubleDouble B; { // For keeping the 't' variable locale. /* @@ -168,35 +143,16 @@ public class TransverseMercator extends */ final DoubleDouble t = initializer.axisLengthRatio(); // t = b/a t.ratio_1m_1p(); // t = (1 - t) / (1 + t) - n = t.doubleValue(); + computeCoefficients(t.doubleValue()); /* * Compute B = (1 + n²/4 + n⁴/64) / (1 + n) */ B = new DoubleDouble(t); // B = n B.square(); B.series(1, 0.25, 1./64); // B = (1 + n²/4 + n⁴/64) - t.add(1,0); + t.add(1, 0); B.divide(t); // B = (1 + n²/4 + n⁴/64) / (1 + n) } - final double n2 = n * n; - final double n3 = n2 * n; - final double n4 = n2 * n2; - /* - * Coefficients for direct projection. - * Add the smallest values first in order to reduce rounding errors. - */ - h1 = ( 41. / 180)*n4 + ( 5. / 16)*n3 + (-2. / 3)*n2 + n/2; - h2 = ( 557. / 1440)*n4 + (-3. / 5)*n3 + (13. / 48)*n2; - h3 = ( -103. / 140)*n4 + (61. / 240)*n3; - h4 = (49561. / 161280)*n4; - /* - * Coefficients for inverse projection. - * Add the smallest values first in order to reduce rounding errors. - */ - ih1 = ( -1. / 360)*n4 + (37. / 96)*n3 + (-2. / 3)*n2 + n/2; - ih2 = (-437. / 1440)*n4 + ( 1. / 15)*n3 + ( 1. / 48)*n2; - ih3 = ( -37. / 840)*n4 + (17. / 480)*n3; - ih4 = (4397. / 161280)*n4; /* * Compute M₀ = B⋅(ξ₁ + ξ₂ + ξ₃ + ξ₄) and negate in anticipation for what will be needed * in the denormalization matrix. We opportunistically use double-double arithmetic but @@ -210,10 +166,10 @@ public class TransverseMercator extends final double Q = asinh(tan(φ0)) - excentricity * atanh(excentricity * sin(φ0)); final double β = atan(sinh(Q)); final DoubleDouble M0 = new DoubleDouble(); - M0.value = h4 * sin(8*β) - + h3 * sin(6*β) - + h2 * sin(4*β) - + h1 * sin(2*β) + M0.value = cf8 * sin(8*β) + + cf6 * sin(6*β) + + cf4 * sin(4*β) + + cf2 * sin(2*β) + β; M0.multiply(B); M0.negate(); @@ -234,6 +190,81 @@ public class TransverseMercator extends final MatrixSIS denormalize = context.getMatrix(false); denormalize.convertBefore(0, B, null); denormalize.convertBefore(1, B, M0); + /* + * When rewriting equations using trigonometric identities, some constants appear. + * For example sin(2θ) = 2⋅sinθ⋅cosθ, so we can factor out the 2 constant into the + * corresponding 'c' field. Note: this factorization can only be performed after + * the constructor finished to compute other constants. + */ + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { + // Multiplication by powers of 2 does not bring any additional rounding error. + cf4 *= 4; ci4 *= 4; + cf6 *= 16; ci6 *= 16; + cf8 *= 64; ci8 *= 64; + } + } + + /** + * Automatically invoked after deserialization for restoring transient fields. + */ + @Override + final void computeCoefficients() { + /* + * Double-double precision is necessary for computing 'n', otherwise we have rounding errors + * in the 3 last digits. Note that we still have sometime a 1 ULP difference compared to the + * 'n' value at serialization time. + */ + final DoubleDouble t = new DoubleDouble(1, 0); + t.subtract(excentricitySquared, 0); + t.sqrt(); + t.ratio_1m_1p(); + computeCoefficients(t.doubleValue()); + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { + // Same scaling than in the constructor. + cf4 *= 4; ci4 *= 4; + cf6 *= 16; ci6 *= 16; + cf8 *= 64; ci8 *= 64; + } + } + + /** + * Computes the transient fields after construction or deserialization. + * The n parameter is defined by the EPSG guide as: + * + *n = f / (2-f)+ * + * Where f is the flattening factor (1 - b/a). This equation can be rewritten as: + * + *
n = (1 - b/a) / (1 + b/a)+ * + * As much as possible, b/a should be computed from the map projection parameters. + * However if those parameters are not available anymore, then they can be computed + * from the excentricity as: + * + *
b/a = √(1 - ℯ²)+ * + * @param n The value of {@code f / (2-f)} where {@code f} is the flattening factor. + */ + private void computeCoefficients(final double n) { + final double n2 = n * n; + final double n3 = n2 * n; + final double n4 = n2 * n2; + /* + * Coefficients for the forward projections. + * Add the smallest values first in order to reduce rounding errors. + */ + cf2 = ( 41. / 180)*n4 + ( 5. / 16)*n3 + (-2. / 3)*n2 + n/2; + cf4 = ( 557. / 1440)*n4 + (-3. / 5)*n3 + (13. / 48)*n2; + cf6 = ( -103. / 140)*n4 + (61. / 240)*n3; + cf8 = (49561. / 161280)*n4; + /* + * Coefficients for the inverse projections. + * Add the smallest values first in order to reduce rounding errors. + */ + ci2 = ( -1. / 360)*n4 + (37. / 96)*n3 + (-2. / 3)*n2 + n/2; + ci4 = (-437. / 1440)*n4 + ( 1. / 15)*n3 + ( 1. / 48)*n2; + ci6 = ( -37. / 840)*n4 + (17. / 480)*n3; + ci8 = (4397. / 161280)*n4; } /** @@ -241,19 +272,15 @@ public class TransverseMercator extends */ TransverseMercator(final TransverseMercator other) { super(other); - h1 = other. h1; - h2 = other. h2; - h3 = other. h3; - h4 = other. h4; - ih1 = other.ih1; - ih2 = other.ih2; - ih3 = other.ih3; - ih4 = other.ih4; + cf2 = other.cf2; + cf4 = other.cf4; + cf6 = other.cf6; + cf8 = other.cf8; } /** * Returns the sequence of normalization → {@code this} → denormalization transforms - * as a whole. The transform returned by this method except (longitude, latitude) + * as a whole. The transform returned by this method expects (longitude, latitude) * coordinates in degrees and returns (x,y) coordinates in metres. * *
The non-linear part of the returned transform will be {@code this} transform, except if the ellipsoid
@@ -285,15 +312,13 @@ public class TransverseMercator extends
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
- final double λ = srcPts[srcOff];
- final double φ = srcPts[srcOff + 1];
-
- final double ℯsinφ = excentricity * sin(φ);
- final double Q = asinh(tan(φ)) - atanh(ℯsinφ) * excentricity;
+ final double λ = srcPts[srcOff ];
+ final double φ = srcPts[srcOff+1];
final double sinλ = sin(λ);
+ final double ℯsinφ = sin(φ) * excentricity;
+ final double Q = asinh(tan(φ)) - atanh(ℯsinφ) * excentricity;
final double coshQ = cosh(Q);
final double η0 = atanh(sinλ / coshQ);
-
/*
* Original formula: η0 = atanh(sin(λ) * cos(β)) where
* cos(β) = cos(atan(sinh(Q)))
@@ -306,16 +331,15 @@ public class TransverseMercator extends
*/
final double coshη0 = cosh(η0);
final double ξ0 = asin(tanh(Q) * coshη0);
-
/*
* Compute sin(2⋅ξ₀), sin(4⋅ξ₀), sin(6⋅ξ₀), sin(8⋅ξ₀) and same for cos, but using the following
* trigonometric identities in order to reduce the number of calls to Math.sin and cos methods.
*/
- final double sin_2ξ0, sin_4ξ0, sin_6ξ0, sin_8ξ0,
- cos_2ξ0, cos_4ξ0, cos_6ξ0, cos_8ξ0;
- if (ORIGINAL_FORMULA) {
- sin_2ξ0 = sin(2*ξ0);
- cos_2ξ0 = cos(2*ξ0);
+ final double sin_2ξ0 = sin(2*ξ0);
+ final double cos_2ξ0 = cos(2*ξ0);
+ final double sin_4ξ0, sin_6ξ0, sin_8ξ0,
+ cos_4ξ0, cos_6ξ0, cos_8ξ0;
+ if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
sin_4ξ0 = sin(4*ξ0);
cos_4ξ0 = cos(4*ξ0);
sin_6ξ0 = sin(6*ξ0);
@@ -323,28 +347,25 @@ public class TransverseMercator extends
sin_8ξ0 = sin(8*ξ0);
cos_8ξ0 = cos(8*ξ0);
} else {
- sin_2ξ0 = sin(2*ξ0); // sin(2⋅ξ₀);
- cos_2ξ0 = cos(2*ξ0); // cos(2⋅ξ₀)
final double sin2 = sin_2ξ0 * sin_2ξ0;
final double cos2 = cos_2ξ0 * cos_2ξ0;
- sin_4ξ0 = 2 * sin_2ξ0 * cos_2ξ0; // sin(4⋅ξ₀)
- cos_4ξ0 = cos2 - sin2; // cos(4⋅ξ₀)
- sin_6ξ0 = (3 - 4*sin2) * sin_2ξ0; // sin(6⋅ξ₀)
- cos_6ξ0 = (4*cos2 - 3) * cos_2ξ0; // cos(6⋅ξ₀)
- sin_8ξ0 = 4*cos_4ξ0 * (sin_2ξ0 * cos_2ξ0); // sin(8⋅ξ₀)
- cos_8ξ0 = 1 - 2*sin_4ξ0*sin_4ξ0; // cos(8⋅ξ₀)
+ sin_4ξ0 = sin_2ξ0 * cos_2ξ0; assert identityEquals(sin_4ξ0, sin(4*ξ0) / 2) : ξ0;
+ cos_4ξ0 = (cos2 - sin2) * 0.5; assert identityEquals(cos_4ξ0, cos(4*ξ0) / 2) : ξ0;
+ sin_6ξ0 = (0.75 - sin2) * sin_2ξ0; assert identityEquals(sin_6ξ0, sin(6*ξ0) / 4) : ξ0;
+ cos_6ξ0 = (cos2 - 0.75) * cos_2ξ0; assert identityEquals(cos_6ξ0, cos(6*ξ0) / 4) : ξ0;
+ sin_8ξ0 = sin_4ξ0 * cos_4ξ0; assert identityEquals(sin_8ξ0, sin(8*ξ0) / 8) : ξ0;
+ cos_8ξ0 = 0.125 - sin_4ξ0 * sin_4ξ0; assert identityEquals(cos_8ξ0, cos(8*ξ0) / 8) : ξ0;
}
-
/*
* Compute sinh(2⋅ξ₀), sinh(4⋅ξ₀), sinh(6⋅ξ₀), sinh(8⋅ξ₀) and same for cosh, but using the following
* hyperbolic identities in order to reduce the number of calls to Math.sinh and cosh methods.
* Note that the formulas are very similar to the above ones, with only some signs reversed.
*/
- final double sinh_2η0, sinh_4η0, sinh_6η0, sinh_8η0,
- cosh_2η0, cosh_4η0, cosh_6η0, cosh_8η0;
- if (ORIGINAL_FORMULA) {
- sinh_2η0 = sinh(2*η0);
- cosh_2η0 = cosh(2*η0);
+ final double sinh_2η0 = sinh(2*η0);
+ final double cosh_2η0 = cosh(2*η0);
+ final double sinh_4η0, sinh_6η0, sinh_8η0,
+ cosh_4η0, cosh_6η0, cosh_8η0;
+ if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
sinh_4η0 = sinh(4*η0);
cosh_4η0 = cosh(4*η0);
sinh_6η0 = sinh(6*η0);
@@ -352,18 +373,15 @@ public class TransverseMercator extends
sinh_8η0 = sinh(8*η0);
cosh_8η0 = cosh(8*η0);
} else {
- sinh_2η0 = sinh(2*η0); // sinh(2⋅η₀);
- cosh_2η0 = cosh(2*η0); // cosh(2⋅η₀)
final double sinh2 = sinh_2η0 * sinh_2η0;
final double cosh2 = cosh_2η0 * cosh_2η0;
- sinh_4η0 = 2 * sinh_2η0 * cosh_2η0; // sinh(4⋅η₀)
- cosh_4η0 = cosh2 + sinh2; // cosh(4⋅η₀)
- sinh_6η0 = (3 + 4*sinh2) * sinh_2η0; // sinh(6⋅η₀)
- cosh_6η0 = (4*cosh2 - 3) * cosh_2η0; // cosh(6⋅η₀)
- sinh_8η0 = 4*cosh_4η0 * (sinh_2η0 * cosh_2η0); // sinh(8⋅η₀)
- cosh_8η0 = 1 + 2*sinh_4η0*sinh_4η0; // cosh(8⋅η₀)
+ cosh_4η0 = (cosh2 + sinh2) * 0.5; assert identityEquals(cosh_4η0, cosh(4*η0) / 2) : η0;
+ sinh_4η0 = cosh_2η0 * sinh_2η0; assert identityEquals(sinh_4η0, sinh(4*η0) / 2) : η0;
+ cosh_6η0 = cosh_2η0 * (cosh2 - 0.75); assert identityEquals(cosh_6η0, cosh(6*η0) / 4) : η0;
+ sinh_6η0 = sinh_2η0 * (sinh2 + 0.75); assert identityEquals(sinh_6η0, sinh(6*η0) / 4) : η0;
+ cosh_8η0 = sinh_4η0 * sinh_4η0 + 0.125; assert identityEquals(cosh_8η0, cosh(8*η0) / 8) : η0;
+ sinh_8η0 = sinh_4η0 * cosh_4η0; assert identityEquals(sinh_8η0, sinh(8*η0) / 8) : η0;
}
-
/*
* Assuming that (λ, φ) ↦ Proj((λ, φ))
* where Proj is defined by: Proj((λ, φ)) : (η(λ, φ), ξ(λ, φ)).
@@ -371,36 +389,35 @@ public class TransverseMercator extends
* => (λ, φ) ↦ (η(λ, φ), ξ(λ, φ)).
*/
//-- ξ(λ, φ)
- final double ξ = h4 * sin_8ξ0 * cosh_8η0
- + h3 * sin_6ξ0 * cosh_6η0
- + h2 * sin_4ξ0 * cosh_4η0
- + h1 * sin_2ξ0 * cosh_2η0
+ final double ξ = cf8 * sin_8ξ0 * cosh_8η0
+ + cf6 * sin_6ξ0 * cosh_6η0
+ + cf4 * sin_4ξ0 * cosh_4η0
+ + cf2 * sin_2ξ0 * cosh_2η0
+ ξ0;
//-- η(λ, φ)
- final double η = h4 * cos_8ξ0 * sinh_8η0
- + h3 * cos_6ξ0 * sinh_6η0
- + h2 * cos_4ξ0 * sinh_4η0
- + h1 * cos_2ξ0 * sinh_2η0
+ final double η = cf8 * cos_8ξ0 * sinh_8η0
+ + cf6 * cos_6ξ0 * sinh_6η0
+ + cf4 * cos_4ξ0 * sinh_4η0
+ + cf2 * cos_2ξ0 * sinh_2η0
+ η0;
if (dstPts != null) {
- dstPts[dstOff ] = η;
- dstPts[dstOff + 1] = ξ;
+ dstPts[dstOff ] = η;
+ dstPts[dstOff+1] = ξ;
}
-
if (!derivate) {
return null;
}
- final double cosλ = cos(λ); //-- λ
- final double cosφ = cos(φ); //-- φ
- final double cosh2Q = coshQ * coshQ; //-- Q
+ final double cosλ = cos(λ); //-- λ
+ final double cosφ = cos(φ); //-- φ
+ final double cosh2Q = coshQ * coshQ; //-- Q
final double sinhQ = sinh(Q);
final double tanhQ = tanh(Q);
- final double cosh2Q_sin2λ = cosh2Q - sinλ * sinλ; //-- Qλ
- final double sinhη0 = sinh(η0); //-- η0
- final double sqrt1_thQchη0 = sqrt(1 - tanhQ * tanhQ * coshη0 * coshη0); //-- Qη0
+ final double cosh2Q_sin2λ = cosh2Q - (sinλ * sinλ); //-- Qλ
+ final double sinhη0 = sinh(η0); //-- η0
+ final double sqrt1_thQchη0 = sqrt(1 - (tanhQ * tanhQ) * (coshη0 * coshη0)); //-- Qη0
//-- dQ_dλ = 0;
final double dQ_dφ = 1 / cosφ - excentricitySquared * cosφ / (1 - ℯsinφ * ℯsinφ);
@@ -410,11 +427,10 @@ public class TransverseMercator extends
final double dξ0_dλ = sinhQ * sinhη0 * cosλ / (cosh2Q_sin2λ * sqrt1_thQchη0);
final double dξ0_dφ = (dQ_dφ * coshη0 / cosh2Q + dη0_dφ * sinhη0 * tanhQ) / sqrt1_thQchη0;
-
/*
* Assuming that Jac(Proj((λ, φ))) is the Jacobian matrix of Proj((λ, φ)) function.
*
- * So derivative Proj((λ, φ)) is defined by:
+ * So the derivative of Proj((λ, φ)) is defined by:
* ┌ ┐
* │ dη(λ, φ) / dλ, dη(λ, φ) / dφ │
* Jac = │ │
@@ -423,31 +439,31 @@ public class TransverseMercator extends
*/
//-- dξ(λ, φ) / dλ
final double dξ_dλ = dξ0_dλ
- + 2 * (h1 * (dξ0_dλ * cos_2ξ0 * cosh_2η0 + dη0_dλ * sinh_2η0 * sin_2ξ0)
- + 3 * h3 * (dξ0_dλ * cos_6ξ0 * cosh_6η0 + dη0_dλ * sinh_6η0 * sin_6ξ0)
- + 2 * (h2 * (dξ0_dλ * cos_4ξ0 * cosh_4η0 + dη0_dλ * sinh_4η0 * sin_4ξ0)
- + 2 * h4 * (dξ0_dλ * cos_8ξ0 * cosh_8η0 + dη0_dλ * sinh_8η0 * sin_8ξ0)));
+ + 2 * (cf2 * (dξ0_dλ * cos_2ξ0 * cosh_2η0 + dη0_dλ * sinh_2η0 * sin_2ξ0)
+ + 3 * cf6 * (dξ0_dλ * cos_6ξ0 * cosh_6η0 + dη0_dλ * sinh_6η0 * sin_6ξ0)
+ + 2 * (cf4 * (dξ0_dλ * cos_4ξ0 * cosh_4η0 + dη0_dλ * sinh_4η0 * sin_4ξ0)
+ + 2 * cf8 * (dξ0_dλ * cos_8ξ0 * cosh_8η0 + dη0_dλ * sinh_8η0 * sin_8ξ0)));
//-- dξ(λ, φ) / dφ
final double dξ_dφ = dξ0_dφ
- + 2 * (h1 * (dξ0_dφ * cos_2ξ0 * cosh_2η0 + dη0_dφ * sinh_2η0 * sin_2ξ0)
- + 3 * h3 * (dξ0_dφ * cos_6ξ0 * cosh_6η0 + dη0_dφ * sinh_6η0 * sin_6ξ0)
- + 2 * (h2 * (dξ0_dφ * cos_4ξ0 * cosh_4η0 + dη0_dφ * sinh_4η0 * sin_4ξ0)
- + 2 * h4 * (dξ0_dφ * cos_8ξ0 * cosh_8η0 + dη0_dφ * sinh_8η0 * sin_8ξ0)));
+ + 2 * (cf2 * (dξ0_dφ * cos_2ξ0 * cosh_2η0 + dη0_dφ * sinh_2η0 * sin_2ξ0)
+ + 3 * cf6 * (dξ0_dφ * cos_6ξ0 * cosh_6η0 + dη0_dφ * sinh_6η0 * sin_6ξ0)
+ + 2 * (cf4 * (dξ0_dφ * cos_4ξ0 * cosh_4η0 + dη0_dφ * sinh_4η0 * sin_4ξ0)
+ + 2 * cf8 * (dξ0_dφ * cos_8ξ0 * cosh_8η0 + dη0_dφ * sinh_8η0 * sin_8ξ0)));
//-- dη(λ, φ) / dλ
final double dη_dλ = dη0_dλ
- + 2 * (h1 * (dη0_dλ * cosh_2η0 * cos_2ξ0 - dξ0_dλ * sin_2ξ0 * sinh_2η0)
- + 3 * h3 * (dη0_dλ * cosh_6η0 * cos_6ξ0 - dξ0_dλ * sin_6ξ0 * sinh_6η0)
- + 2 * (h2 * (dη0_dλ * cosh_4η0 * cos_4ξ0 - dξ0_dλ * sin_4ξ0 * sinh_4η0)
- + 2 * h4 * (dη0_dλ * cosh_8η0 * cos_8ξ0 - dξ0_dλ * sin_8ξ0 * sinh_8η0)));
+ + 2 * (cf2 * (dη0_dλ * cosh_2η0 * cos_2ξ0 - dξ0_dλ * sin_2ξ0 * sinh_2η0)
+ + 3 * cf6 * (dη0_dλ * cosh_6η0 * cos_6ξ0 - dξ0_dλ * sin_6ξ0 * sinh_6η0)
+ + 2 * (cf4 * (dη0_dλ * cosh_4η0 * cos_4ξ0 - dξ0_dλ * sin_4ξ0 * sinh_4η0)
+ + 2 * cf8 * (dη0_dλ * cosh_8η0 * cos_8ξ0 - dξ0_dλ * sin_8ξ0 * sinh_8η0)));
//-- dη(λ, φ) / dφ
final double dη_dφ = dη0_dφ
- + 2 * (h1 * (dη0_dφ * cosh_2η0 * cos_2ξ0 - dξ0_dφ * sin_2ξ0 * sinh_2η0)
- + 3 * h3 * (dη0_dφ * cosh_6η0 * cos_6ξ0 - dξ0_dφ * sin_6ξ0 * sinh_6η0)
- + 2 * (h2 * (dη0_dφ * cosh_4η0 * cos_4ξ0 - dξ0_dφ * sin_4ξ0 * sinh_4η0)
- + 2 * h4 * (dη0_dφ * cosh_8η0 * cos_8ξ0 - dξ0_dφ * sin_8ξ0 * sinh_8η0)));
+ + 2 * (cf2 * (dη0_dφ * cosh_2η0 * cos_2ξ0 - dξ0_dφ * sin_2ξ0 * sinh_2η0)
+ + 3 * cf6 * (dη0_dφ * cosh_6η0 * cos_6ξ0 - dξ0_dφ * sin_6ξ0 * sinh_6η0)
+ + 2 * (cf4 * (dη0_dφ * cosh_4η0 * cos_4ξ0 - dξ0_dφ * sin_4ξ0 * sinh_4η0)
+ + 2 * cf8 * (dη0_dφ * cosh_8η0 * cos_8ξ0 - dξ0_dφ * sin_8ξ0 * sinh_8η0)));
return new Matrix2(dη_dλ, dη_dφ,
dξ_dλ, dξ_dφ);
@@ -463,20 +479,22 @@ public class TransverseMercator extends
final double[] dstPts, final int dstOff)
throws ProjectionException
{
- final double η = srcPts[srcOff ];
- final double ξ = srcPts[srcOff + 1];
+ final double η = srcPts[srcOff ];
+ final double ξ = srcPts[srcOff+1];
/*
* Following calculation of sin_2ξ, sin_4ξ, etc. is basically a copy-and-paste of the code in transform(…).
* Its purpose is the same than for transform(…): reduce the amount of calls to Math.sin(double) and other
* methods.
*/
- final double sin_2ξ, sin_4ξ, sin_6ξ, sin_8ξ,
- cos_2ξ, cos_4ξ, cos_6ξ, cos_8ξ,
- sinh_2η, sinh_4η, sinh_6η, sinh_8η,
- cosh_2η, cosh_4η, cosh_6η, cosh_8η;
- if (ORIGINAL_FORMULA) {
- sin_2ξ = sin(2*ξ);
- cos_2ξ = cos(2*ξ);
+ final double sin_2ξ = sin (2*ξ);
+ final double cos_2ξ = cos (2*ξ);
+ final double sinh_2η = sinh(2*η);
+ final double cosh_2η = cosh(2*η);
+ final double sin_4ξ, sin_6ξ, sin_8ξ,
+ cos_4ξ, cos_6ξ, cos_8ξ,
+ sinh_4η, sinh_6η, sinh_8η,
+ cosh_4η, cosh_6η, cosh_8η;
+ if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
sin_4ξ = sin(4*ξ);
cos_4ξ = cos(4*ξ);
sin_6ξ = sin(6*ξ);
@@ -484,8 +502,6 @@ public class TransverseMercator extends
sin_8ξ = sin(8*ξ);
cos_8ξ = cos(8*ξ);
- sinh_2η = sinh(2*η);
- cosh_2η = cosh(2*η);
sinh_4η = sinh(4*η);
cosh_4η = cosh(4*η);
sinh_6η = sinh(6*η);
@@ -493,40 +509,36 @@ public class TransverseMercator extends
sinh_8η = sinh(8*η);
cosh_8η = cosh(8*η);
} else {
- sin_2ξ = sin(2*ξ);
- cos_2ξ = cos(2*ξ);
final double sin2 = sin_2ξ * sin_2ξ;
final double cos2 = cos_2ξ * cos_2ξ;
- sin_4ξ = 2 * sin_2ξ * cos_2ξ;
- cos_4ξ = cos2 - sin2;
- sin_6ξ = (3 - 4*sin2) * sin_2ξ;
- cos_6ξ = (4*cos2 - 3) * cos_2ξ;
- sin_8ξ = 4*cos_4ξ * (sin_2ξ * cos_2ξ);
- cos_8ξ = 1 - 2*sin_4ξ*sin_4ξ;
+ sin_4ξ = sin_2ξ * cos_2ξ; assert identityEquals(sin_4ξ, sin(4*ξ) / 2) : ξ;
+ cos_4ξ = (cos2 - sin2) * 0.5; assert identityEquals(cos_4ξ, cos(4*ξ) / 2) : ξ;
+ sin_6ξ = (0.75 - sin2) * sin_2ξ; assert identityEquals(sin_6ξ, sin(6*ξ) / 4) : ξ;
+ cos_6ξ = (cos2 - 0.75) * cos_2ξ; assert identityEquals(cos_6ξ, cos(6*ξ) / 4) : ξ;
+ sin_8ξ = sin_4ξ * cos_4ξ; assert identityEquals(sin_8ξ, sin(8*ξ) / 8) : ξ;
+ cos_8ξ = 0.125 - sin_4ξ * sin_4ξ; assert identityEquals(cos_8ξ, cos(8*ξ) / 8) : ξ;
- sinh_2η = sinh(2*η);
- cosh_2η = cosh(2*η);
final double sinh2 = sinh_2η * sinh_2η;
final double cosh2 = cosh_2η * cosh_2η;
- sinh_4η = 2 * sinh_2η * cosh_2η;
- cosh_4η = cosh2 + sinh2;
- sinh_6η = (3 + 4*sinh2) * sinh_2η;
- cosh_6η = (4*cosh2 - 3) * cosh_2η;
- sinh_8η = 4*cosh_4η * (sinh_2η * cosh_2η);
- cosh_8η = 1 + 2*sinh_4η*sinh_4η;
+ cosh_4η = (cosh2 + sinh2) * 0.5; assert identityEquals(cosh_4η, cosh(4*η) / 2) : η;
+ sinh_4η = cosh_2η * sinh_2η; assert identityEquals(sinh_4η, sinh(4*η) / 2) : η;
+ cosh_6η = cosh_2η * (cosh2 - 0.75); assert identityEquals(cosh_6η, cosh(6*η) / 4) : η;
+ sinh_6η = sinh_2η * (sinh2 + 0.75); assert identityEquals(sinh_6η, sinh(6*η) / 4) : η;
+ cosh_8η = sinh_4η * sinh_4η + 0.125; assert identityEquals(cosh_8η, cosh(8*η) / 8) : η;
+ sinh_8η = sinh_4η * cosh_4η; assert identityEquals(sinh_8η, sinh(8*η) / 8) : η;
}
/*
* The actual inverse transform.
*/
- final double ξ0 = ξ - (ih4 * sin_8ξ * cosh_8η
- + ih3 * sin_6ξ * cosh_6η
- + ih2 * sin_4ξ * cosh_4η
- + ih1 * sin_2ξ * cosh_2η);
-
- final double η0 = η - (ih4 * cos_8ξ * sinh_8η
- + ih3 * cos_6ξ * sinh_6η
- + ih2 * cos_4ξ * sinh_4η
- + ih1 * cos_2ξ * sinh_2η);
+ final double ξ0 = ξ - (ci8 * sin_8ξ * cosh_8η
+ + ci6 * sin_6ξ * cosh_6η
+ + ci4 * sin_4ξ * cosh_4η
+ + ci2 * sin_2ξ * cosh_2η);
+
+ final double η0 = η - (ci8 * cos_8ξ * sinh_8η
+ + ci6 * cos_6ξ * sinh_6η
+ + ci4 * cos_4ξ * sinh_4η
+ + ci2 * cos_2ξ * sinh_2η);
final double β = asin(sin(ξ0) / cosh(η0));
final double Q = asinh(tan(β));
@@ -538,8 +550,8 @@ public class TransverseMercator extends
final double c = excentricity * atanh(excentricity * tanh(Qp));
Qp = Q + c;
if (abs(c - p) <= ITERATION_TOLERANCE) {
- dstPts[dstOff ] = asin(tanh(η0) / cos(β));
- dstPts[dstOff + 1] = atan(sinh(Qp));
+ dstPts[dstOff ] = asin(tanh(η0) / cos(β));
+ dstPts[dstOff+1] = atan(sinh(Qp));
return;
}
p = c;
@@ -581,8 +593,8 @@ public class TransverseMercator extends
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
- final double λ = srcPts[srcOff];
- final double φ = srcPts[srcOff + 1];
+ final double λ = srcPts[srcOff ];
+ final double φ = srcPts[srcOff+1];
final double sinλ = sin(λ);
final double cosλ = cos(λ);
final double sinφ = sin(φ);
Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/io/wkt/WKTParserTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/io/wkt/WKTParserTest.java?rev=1707296&r1=1707295&r2=1707296&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/io/wkt/WKTParserTest.java [UTF-8] (original)
+++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/io/wkt/WKTParserTest.java [UTF-8] Wed Oct 7 13:20:26 2015
@@ -533,7 +533,7 @@ public final strictfp class WKTParserTes
* AXIS[“longitude”,east,ORDER[2]],
* ANGLEUNIT[“degree”,0.0174532925199433]],
* VERTCRS[“NAVD88”,
- * VDATUM[“North American Vertical Datum 1983”],
+ * VDATUM[“North American Vertical Datum 1988”],
* CS[vertical,1],
* AXIS[“gravity-related height (H)”,up],
* LENGTHUNIT[“metre”,1]]]
Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java?rev=1707296&r1=1707295&r2=1707296&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java [UTF-8] (original)
+++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java [UTF-8] Wed Oct 7 13:20:26 2015
@@ -37,7 +37,7 @@ import org.junit.Test;
import static java.lang.StrictMath.*;
import static java.lang.Double.*;
-import static org.junit.Assert.*;
+import static org.apache.sis.test.Assert.*;
/**
@@ -49,7 +49,7 @@ import static org.junit.Assert.*;
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
* @since 0.6
- * @version 0.6
+ * @version 0.7
* @module
*/
@DependsOn(ConformalProjectionTest.class)
@@ -336,4 +336,23 @@ public final strictfp class LambertConic
tolerance = Formulas.LINEAR_TOLERANCE;
compareEllipticalWithSpherical(CoordinateDomain.GEOGRAPHIC_SAFE, 0);
}
+
+ /**
+ * Verifies that deserialized projections work as expected. This implies that deserialization
+ * recomputed the internal transient fields, especially the series expansion coefficients.
+ *
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a coordinate.
+ */
+ @Test
+ @DependsOnMethod("testLambertConicConformal1SP")
+ public void testSerialization() throws FactoryException, TransformException {
+ createNormalizedProjection(true, 40);
+ final double[] source = CoordinateDomain.GEOGRAPHIC_RADIANS_NORTH.generateRandomInput(TestUtilities.createRandomNumberGenerator(), 2, 10);
+ final double[] target = new double[source.length];
+ transform.transform(source, 0, target, 0, 10);
+ transform = assertSerializedEquals(transform);
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ verifyTransform(source, target);
+ }
}
Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java?rev=1707296&r1=1707295&r2=1707296&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java [UTF-8] (original)
+++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java [UTF-8] Wed Oct 7 13:20:26 2015
@@ -75,6 +75,7 @@ final strictfp class NoOp extends Confor
super(new Initializer(new DefaultOperationMethod(
Collections.singletonMap(DefaultOperationMethod.NAME_KEY, parameters.getDescriptor().getName()),
2, 2, parameters.getDescriptor()), parameters, Collections. This class works as a wrapper around a collection of identifiers. Because all operations
* are performed by an iteration over the collection elements, this implementation is suitable
* only for small maps (less than 10 elements). Given that objects typically have only one or
* two identifiers, this is considered acceptable. Because {@code null} elements are ignored, this method may return 0 even if {@link #isEmpty()}
+ * returns {@code false}. However this inconsistency should not happen in practice because
+ * {@link org.apache.sis.metadata.ModifiableMetadata} internal collection implementations
+ * do not allow null values. This is a helper method for {@link IdentifierMapWithSpecialCases#put(Citation, String)},
- * in order to be able to return the old value if that value was a {@link String} rather than
- * the specialized type. We do not return the string for the specialized case in order to avoid
- * the cost of invoking {@code toString()} on the specialized object (some may be costly). Such
- * call would be useless because {@code IdentifierMapWithSpecialCase} discard the value of this
- * method when it found a specialized type. If the backing identifier collection contains null entries, those entries will be ignored.
+ * If the backing collection contains many entries for the same authority, then only the first
+ * occurrence is included. This iterator supports the {@link #remove()} operation if the underlying collection
- * supports it. This iterator supports the {@link #remove()} operation if the underlying collection supports it. The map entries are used as a safety against duplicated authority values. The map values
* are non-null only after we iterated over an authority. Then the value is {@link Boolean#TRUE}
@@ -416,7 +394,7 @@ public class IdentifierMapAdapter extend
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.3
- * @version 0.3
+ * @version 0.7
* @module
*/
@SuppressWarnings("serial") // Not intended to be serialized.
@@ -439,11 +417,17 @@ public class IdentifierMapAdapter extend
private transient Citation authority;
/**
+ * {@code true} if the iterator should support the {@link #remove()} operation.
+ */
+ private final boolean isModifiable;
+
+ /**
* Creates a new iterator for the given collection of identifiers.
*/
- Iter(final Collection extends Identifier> identifiers) {
+ Iter(final Collection extends Identifier> identifiers, final boolean isModifiable) {
super(hashMapCapacity(identifiers.size()));
this.identifiers = identifiers.iterator();
+ this.isModifiable = isModifiable;
}
/**
@@ -518,6 +502,9 @@ public class IdentifierMapAdapter extend
*/
@Override
public void remove() throws IllegalStateException {
+ if (!isModifiable) {
+ throw new UnsupportedOperationException();
+ }
final Iterator extends Identifier> it = identifiers;
if (it == null || next != null) {
throw new IllegalStateException();
Modified: sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/jaxb/NonMarshalledAuthority.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/jaxb/NonMarshalledAuthority.java?rev=1707296&r1=1707295&r2=1707296&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/jaxb/NonMarshalledAuthority.java [UTF-8] (original)
+++ sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/jaxb/NonMarshalledAuthority.java [UTF-8] Wed Oct 7 13:20:26 2015
@@ -17,12 +17,18 @@
package org.apache.sis.internal.jaxb;
import java.util.Iterator;
+import java.util.List;
import java.util.ArrayList;
+import java.util.Map;
+import java.util.IdentityHashMap;
import java.util.Collection;
+import java.util.Collections;
import org.opengis.metadata.Identifier;
import org.opengis.metadata.citation.Citation;
import org.apache.sis.internal.simple.CitationConstant;
+import org.apache.sis.internal.util.CollectionsExt;
import org.apache.sis.internal.util.UnmodifiableArrayList;
+import org.apache.sis.util.collection.Containers;
import org.apache.sis.xml.IdentifierSpace;
@@ -57,7 +63,7 @@ import org.apache.sis.xml.IdentifierSpac
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.3
- * @version 0.6
+ * @version 0.7
* @module
*
* @see IdentifierSpace
@@ -73,6 +79,7 @@ public final class NonMarshalledAuthorit
* mirror the constants defined in the {@link IdentifierSpace} interface
* and {@link org.apache.sis.metadata.iso.citation.DefaultCitation} class.
*/
+ @SuppressWarnings("FieldNameHidesFieldInSuperclass")
public static final byte ID=0, UUID=1, HREF=2, XLINK=3, ISSN=4, ISBN=5;
// If more codes are added, please update readResolve() below.
@@ -103,6 +110,9 @@ public final class NonMarshalledAuthorit
* "special" identifiers (ISO 19139 attributes, ISBN codes...), which are recognized by
* the implementation class of their authority.
*
+ * This method is used for implementation of {@code getIdentifier()} methods (singular form)
+ * in public metadata objects. This method is used for implementation of {@code setIdentifier(Identifier)} methods
+ * in public metadata objects. This method is used for implementation of {@code getIdentifiers()} methods (plural form)
+ * in public metadata objects. Note that those methods override
+ * {@link org.apache.sis.xml.IdentifiedObject#getIdentifiers()}, which is expected to return
+ * all identifiers in normal (non-marshalling) usage. This method is invoked from {@code setIdentifiers(Collection This method is used for implementation of {@code setIdentifiers(Collection)} methods
+ * in public metadata objects. {@link org.apache.sis.internal.jaxb.SpecializedIdentifier} wraps {@link org.apache.sis.xml.XLink},
* {@link java.net.URI} and {@link java.util.UUID} as {@link org.opengis.metadata.Identifier} instances.
Modified: sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/collection/CodeListSet.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/collection/CodeListSet.java?rev=1707296&r1=1707295&r2=1707296&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/collection/CodeListSet.java [UTF-8] (original)
+++ sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/collection/CodeListSet.java [UTF-8] Wed Oct 7 13:20:26 2015
@@ -547,7 +547,6 @@ public class CodeListSet
+ *
+ *
*
- *
+ * This class performs opportunist null checks as an additional safety, but consistency is not
+ * guaranteed. See {@link #size()} for more information.
*
*
- *
- *
- * @param The type of elements in the storage (original) set.
* @param