001 //$HeadURL: svn+ssh://rbezema@svn.wald.intevation.org/deegree/base/tags/2.1/src/org/deegree/model/csct/ct/GeocentricTransform.java $
002 /*---------------- FILE HEADER ------------------------------------------
003
004 This file is part of deegree.
005 Copyright (C) 2001 by:
006 EXSE, Department of Geography, University of Bonn
007 http://www.giub.uni-bonn.de/exse/
008 lat/lon GmbH
009 http://www.lat-lon.de
010
011 It has been implemented within SEAGIS - An OpenSource implementation of OpenGIS specification
012 (C) 2001, Institut de Recherche pour le D�veloppement (http://sourceforge.net/projects/seagis/)
013 SEAGIS Contacts: Surveillance de l'Environnement Assist�e par Satellite
014 Institut de Recherche pour le D�veloppement / US-Espace
015 mailto:seasnet@teledetection.fr
016
017
018 This library is free software; you can redistribute it and/or
019 modify it under the terms of the GNU Lesser General Public
020 License as published by the Free Software Foundation; either
021 version 2.1 of the License, or (at your option) any later version.
022
023 This library is distributed in the hope that it will be useful,
024 but WITHOUT ANY WARRANTY; without even the implied warranty of
025 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
026 Lesser General Public License for more details.
027
028 You should have received a copy of the GNU Lesser General Public
029 License along with this library; if not, write to the Free Software
030 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
031
032 Contact:
033
034 Andreas Poth
035 lat/lon GmbH
036 Aennchenstr. 19
037 53115 Bonn
038 Germany
039 E-Mail: poth@lat-lon.de
040
041 Klaus Greve
042 Department of Geography
043 University of Bonn
044 Meckenheimer Allee 166
045 53115 Bonn
046 Germany
047 E-Mail: klaus.greve@uni-bonn.de
048
049
050 ---------------------------------------------------------------------------*/
051 package org.deegree.model.csct.ct;
052
053 // OpenGIS dependencies (SEAGIS)
054 import java.io.Serializable;
055
056 import javax.media.jai.ParameterList;
057 import javax.media.jai.util.Range;
058
059 import org.deegree.model.csct.cs.Ellipsoid;
060 import org.deegree.model.csct.resources.css.ResourceKeys;
061 import org.deegree.model.csct.resources.css.Resources;
062 import org.deegree.model.csct.units.Unit;
063
064 /**
065 * Transforms three dimensional geographic points to geocentric
066 * coordinate points. Input points must be longitudes, latitudes
067 * and heights above the ellipsoid.
068 *
069 * @version 1.00
070 * @author Frank Warmerdam
071 * @author Martin Desruisseaux
072 */
073 class GeocentricTransform extends AbstractMathTransform implements Serializable {
074 /**
075 * Serial number for interoperability with different versions.
076 */
077 private static final long serialVersionUID = -3352045463953828140L;
078
079 /**
080 * Cosine of 67.5 degrees.
081 */
082 private static final double COS_67P5 = 0.38268343236508977;
083
084 /**
085 * Toms region 1 constant.
086 */
087 private static final double AD_C = 1.0026000;
088
089 /**
090 * Semi-major axis of ellipsoid in meters.
091 */
092 private final double a;
093
094 /**
095 * Semi-minor axis of ellipsoid in meters.
096 */
097 private final double b;
098
099 /**
100 * Square of semi-major axis (@link #a}�).
101 */
102 private final double a2;
103
104 /**
105 * Square of semi-minor axis ({@link #b}�).
106 */
107 private final double b2;
108
109 /**
110 * Eccentricity squared.
111 */
112 private final double e2;
113
114 /**
115 * 2nd eccentricity squared.
116 */
117 private final double ep2;
118
119 /**
120 * <code>true</code> if geographic coordinates
121 * include an ellipsoidal height (i.e. are 3-D),
122 * or <code>false</code> if they are strictly 2-D.
123 */
124 private final boolean hasHeight;
125
126 /**
127 * The inverse of this transform.
128 * Will be created only when needed.
129 */
130 private transient MathTransform inverse;
131
132 /**
133 * Construct a transform.
134 *
135 * @param ellipsoid The ellipsoid.
136 * @param hasHeight <code>true</code> if geographic coordinates
137 * include an ellipsoidal height (i.e. are 3-D),
138 * or <code>false</code> if they are strictly 2-D.
139 */
140 protected GeocentricTransform( final Ellipsoid ellipsoid, final boolean hasHeight ) {
141 this( ellipsoid.getSemiMajorAxis(), ellipsoid.getSemiMinorAxis(), ellipsoid.getAxisUnit(),
142 hasHeight );
143 }
144
145 /**
146 * Construct a transform.
147 *
148 * @param semiMajor The semi-major axis length.
149 * @param semiMinor The semi-minor axis length.
150 * @param units The axis units.
151 * @param hasHeight <code>true</code> if geographic coordinates
152 * include an ellipsoidal height (i.e. are 3-D),
153 * or <code>false</code> if they are strictly 2-D.
154 */
155 protected GeocentricTransform( final double semiMajor, final double semiMinor,
156 final Unit units, final boolean hasHeight ) {
157 this.hasHeight = hasHeight;
158 a = Unit.METRE.convert( semiMajor, units );
159 b = Unit.METRE.convert( semiMinor, units );
160 a2 = a * a;
161 b2 = b * b;
162 e2 = ( a2 - b2 ) / a2;
163 ep2 = ( a2 - b2 ) / b2;
164 checkArgument( "a", a, Double.MAX_VALUE );
165 checkArgument( "b", b, a );
166 }
167
168 /**
169 * Check an argument value. The argument must be greater
170 * than 0 and finite, otherwise an exception is thrown.
171 *
172 * @param name The argument name.
173 * @param value The argument value.
174 * @param max The maximal legal argument value.
175 */
176 private static void checkArgument( final String name, final double value, final double max )
177 throws IllegalArgumentException {
178 if ( !( value >= 0 && value <= max ) ) // Use '!' in order to trap NaN
179 {
180 throw new IllegalArgumentException(
181 Resources.format(
182 ResourceKeys.ERROR_ILLEGAL_ARGUMENT_$2,
183 name, new Double( value ) ) );
184 }
185 }
186
187 /**
188 * Converts geodetic coordinates (longitude, latitude, height) to
189 * geocentric coordinates (x, y, z) according to the current ellipsoid
190 * parameters.
191 */
192 public void transform( double[] srcPts, int srcOff, double[] dstPts, int dstOff, int numPts ) {
193 transform( srcPts, srcOff, dstPts, dstOff, numPts, false );
194 }
195
196 /**
197 * Implementation of geodetic to geocentric conversion. This
198 * implementation allows the caller to use height in computation.
199 */
200 private void transform( final double[] srcPts, int srcOff, final double[] dstPts, int dstOff,
201 int numPts, boolean hasHeight ) {
202 int step = 0;
203 final int dimSource = getDimSource();
204 hasHeight |= ( dimSource >= 3 );
205 if ( srcPts == dstPts && srcOff < dstOff && srcOff + numPts * dimSource > dstOff ) {
206 step = -dimSource;
207 srcOff -= ( numPts - 1 ) * step;
208 dstOff -= ( numPts - 1 ) * step;
209 }
210 while ( --numPts >= 0 ) {
211 final double L = Math.toRadians( srcPts[srcOff++] ); // Longitude
212 final double P = Math.toRadians( srcPts[srcOff++] ); // Latitude
213 final double h = hasHeight ? srcPts[srcOff++] : 0; // Height above the ellipsoid (metres).
214
215 final double cosLat = Math.cos( P );
216 final double sinLat = Math.sin( P );
217 final double rn = a / Math.sqrt( 1 - e2 * ( sinLat * sinLat ) );
218
219 dstPts[dstOff++] = ( rn + h ) * cosLat * Math.cos( L ); // X
220 dstPts[dstOff++] = ( rn + h ) * cosLat * Math.sin( L ); // Y
221 dstPts[dstOff++] = ( rn * ( 1 - e2 ) + h ) * sinLat; // Z
222 srcOff += step;
223 dstOff += step;
224 }
225 }
226
227 /**
228 * Converts geodetic coordinates (longitude, latitude, height) to
229 * geocentric coordinates (x, y, z) according to the current ellipsoid
230 * parameters.
231 */
232 public void transform( final float[] srcPts, int srcOff, final float[] dstPts, int dstOff,
233 int numPts ) {
234 int step = 0;
235 final int dimSource = getDimSource();
236 final boolean hasHeight = ( dimSource >= 3 );
237 if ( srcPts == dstPts && srcOff < dstOff && srcOff + numPts * dimSource > dstOff ) {
238 step = -dimSource;
239 srcOff -= ( numPts - 1 ) * step;
240 dstOff -= ( numPts - 1 ) * step;
241 }
242 while ( --numPts >= 0 ) {
243 final double L = Math.toRadians( srcPts[srcOff++] ); // Longitude
244 final double P = Math.toRadians( srcPts[srcOff++] ); // Latitude
245 final double h = hasHeight ? srcPts[srcOff++] : 0; // Height above the ellipsoid (metres).
246
247 final double cosLat = Math.cos( P );
248 final double sinLat = Math.sin( P );
249 final double rn = a / Math.sqrt( 1 - e2 * ( sinLat * sinLat ) );
250
251 dstPts[dstOff++] = (float) ( ( rn + h ) * cosLat * Math.cos( L ) ); // X
252 dstPts[dstOff++] = (float) ( ( rn + h ) * cosLat * Math.sin( L ) ); // Y
253 dstPts[dstOff++] = (float) ( ( rn * ( 1 - e2 ) + h ) * sinLat ); // Z
254 srcOff += step;
255 dstOff += step;
256 }
257 }
258
259 /**
260 * Converts geocentric coordinates (x, y, z) to geodetic coordinates
261 * (longitude, latitude, height), according to the current ellipsoid
262 * parameters. The method used here is derived from "An Improved
263 * Algorithm for Geocentric to Geodetic Coordinate Conversion", by
264 * Ralph Toms, Feb 1996.
265 */
266 protected final void inverseTransform( final double[] srcPts, int srcOff,
267 final double[] dstPts, int dstOff, int numPts ) {
268 int step = 0;
269 final int dimSource = getDimSource();
270 final boolean hasHeight = ( dimSource >= 3 );
271 boolean computeHeight = hasHeight;
272
273 if ( srcPts == dstPts && srcOff < dstOff && srcOff + numPts * dimSource > dstOff ) {
274 step = -dimSource;
275 srcOff -= ( numPts - 1 ) * step;
276 dstOff -= ( numPts - 1 ) * step;
277 }
278 while ( --numPts >= 0 ) {
279 final double x = srcPts[srcOff++];
280 final double y = srcPts[srcOff++];
281 final double z = srcPts[srcOff++];
282
283 // Note: The Java version of 'atan2' work correctly for x==0.
284 // No need for special handling like in the C version.
285 // No special handling neither for latitude. Formulas
286 // below are generic enough, considering that 'atan'
287 // work correctly with infinities (1/0).
288
289 // Note: Variable names follow the notation used in Toms, Feb 1996
290 final double W2 = x * x + y * y; // square of distance from Z axis
291 final double W = Math.sqrt( W2 ); // distance from Z axis
292 final double T0 = z * AD_C; // initial estimate of vertical component
293 final double S0 = Math.sqrt( T0 * T0 + W2 ); // initial estimate of horizontal component
294 final double sin_B0 = T0 / S0; // sin(B0), B0 is estimate of Bowring aux variable
295 final double cos_B0 = W / S0; // cos(B0)
296 final double sin3_B0 = sin_B0 * sin_B0 * sin_B0; // cube of sin(B0)
297 final double T1 = z + b * ep2 * sin3_B0; // corrected estimate of vertical component
298 final double sum = W - a * e2 * ( cos_B0 * cos_B0 * cos_B0 ); // numerator of cos(phi1)
299 final double S1 = Math.sqrt( T1 * T1 + sum * sum ); // corrected estimate of horizontal component
300 final double sin_p1 = T1 / S1; // sin(phi1), phi1 is estimated latitude
301 final double cos_p1 = sum / S1; // cos(phi1)
302
303 final double longitude = Math.toDegrees( Math.atan2( y, x ) );
304 final double latitude = Math.toDegrees( Math.atan( sin_p1 / cos_p1 ) );
305 final double height;
306
307 dstPts[dstOff++] = longitude;
308 dstPts[dstOff++] = latitude;
309 if ( computeHeight ) {
310 final double rn = a / Math.sqrt( 1 - e2 * ( sin_p1 * sin_p1 ) ); // Earth radius at location
311 if ( cos_p1 >= +COS_67P5 )
312 height = W / +cos_p1 - rn;
313 else if ( cos_p1 <= -COS_67P5 )
314 height = W / -cos_p1 - rn;
315 else
316 height = z / sin_p1 + rn * ( e2 - 1.0 );
317 if ( hasHeight ) {
318 dstPts[dstOff++] = height;
319 }
320 }
321 srcOff += step;
322 dstOff += step;
323 }
324 }
325
326 /**
327 * Converts geocentric coordinates (x, y, z) to geodetic coordinates
328 * (longitude, latitude, height), according to the current ellipsoid
329 * parameters. The method used here is derived from "An Improved
330 * Algorithm for Geocentric to Geodetic Coordinate Conversion", by
331 * Ralph Toms, Feb 1996.
332 */
333 protected final void inverseTransform( final float[] srcPts, int srcOff, final float[] dstPts,
334 int dstOff, int numPts ) {
335 int step = 0;
336 final int dimSource = getDimSource();
337 final boolean hasHeight = ( dimSource >= 3 );
338 boolean computeHeight = hasHeight;
339 if ( srcPts == dstPts && srcOff < dstOff && srcOff + numPts * dimSource > dstOff ) {
340 step = -dimSource;
341 srcOff -= ( numPts - 1 ) * step;
342 dstOff -= ( numPts - 1 ) * step;
343 }
344 while ( --numPts >= 0 ) {
345 final double x = srcPts[srcOff++];
346 final double y = srcPts[srcOff++];
347 final double z = srcPts[srcOff++];
348
349 // Note: The Java version of 'atan2' work correctly for x==0.
350 // No need for special handling like in the C version.
351 // No special handling neither for latitude. Formulas
352 // below are generic enough, considering that 'atan'
353 // work correctly with infinities (1/0).
354
355 // Note: Variable names follow the notation used in Toms, Feb 1996
356 final double W2 = x * x + y * y; // square of distance from Z axis
357 final double W = Math.sqrt( W2 ); // distance from Z axis
358 final double T0 = z * AD_C; // initial estimate of vertical component
359 final double S0 = Math.sqrt( T0 * T0 + W2 ); // initial estimate of horizontal component
360 final double sin_B0 = T0 / S0; // sin(B0), B0 is estimate of Bowring aux variable
361 final double cos_B0 = W / S0; // cos(B0)
362 final double sin3_B0 = sin_B0 * sin_B0 * sin_B0; // cube of sin(B0)
363 final double T1 = z + b * ep2 * sin3_B0; // corrected estimate of vertical component
364 final double sum = W - a * e2 * ( cos_B0 * cos_B0 * cos_B0 ); // numerator of cos(phi1)
365 final double S1 = Math.sqrt( T1 * T1 + sum * sum ); // corrected estimate of horizontal component
366 final double sin_p1 = T1 / S1; // sin(phi1), phi1 is estimated latitude
367 final double cos_p1 = sum / S1; // cos(phi1)
368
369 final double longitude = Math.toDegrees( Math.atan2( y, x ) );
370 final double latitude = Math.toDegrees( Math.atan( sin_p1 / cos_p1 ) );
371 final double height;
372
373 dstPts[dstOff++] = (float) longitude;
374 dstPts[dstOff++] = (float) latitude;
375 if ( computeHeight ) {
376 final double rn = a / Math.sqrt( 1 - e2 * ( sin_p1 * sin_p1 ) ); // Earth radius at location
377 if ( cos_p1 >= +COS_67P5 )
378 height = W / +cos_p1 - rn;
379 else if ( cos_p1 <= -COS_67P5 )
380 height = W / -cos_p1 - rn;
381 else
382 height = z / sin_p1 + rn * ( e2 - 1.0 );
383 if ( hasHeight ) {
384 dstPts[dstOff++] = (float) height;
385 }
386 }
387 srcOff += step;
388 dstOff += step;
389 }
390 }
391
392 /**
393 * Gets the dimension of input points, which is 2 or 3.
394 */
395 public int getDimSource() {
396 return hasHeight ? 3 : 2;
397 }
398
399 /**
400 * Gets the dimension of output points, which is 3.
401 */
402 public final int getDimTarget() {
403 return 3;
404 }
405
406 /**
407 * Tests whether this transform does not move any points.
408 * This method returns always <code>false</code>.
409 */
410 public final boolean isIdentity() {
411 return false;
412 }
413
414 /**
415 * Returns the inverse of this transform.
416 */
417 public synchronized MathTransform inverse() {
418 if ( inverse == null )
419 inverse = new Inverse();
420 return inverse;
421 }
422
423 /**
424 * Returns a hash value for this transform.
425 */
426 public final int hashCode() {
427 final long code = Double.doubleToLongBits( a )
428 + 37
429 * ( Double.doubleToLongBits( b ) + 37 * ( Double.doubleToLongBits( a2 ) + 37 * ( Double.doubleToLongBits( b2 ) + 37 * ( Double.doubleToLongBits( e2 ) + 37 * ( Double.doubleToLongBits( ep2 ) ) ) ) ) );
430 return (int) code ^ (int) ( code >>> 32 );
431 }
432
433 /**
434 * Compares the specified object with
435 * this math transform for equality.
436 */
437 public final boolean equals( final Object object ) {
438 if ( object == this )
439 return true; // Slight optimization
440 if ( super.equals( object ) ) {
441 final GeocentricTransform that = (GeocentricTransform) object;
442 return Double.doubleToLongBits( this.a ) == Double.doubleToLongBits( that.a )
443 && Double.doubleToLongBits( this.b ) == Double.doubleToLongBits( that.b )
444 && Double.doubleToLongBits( this.a2 ) == Double.doubleToLongBits( that.a2 )
445 && Double.doubleToLongBits( this.b2 ) == Double.doubleToLongBits( that.b2 )
446 && Double.doubleToLongBits( this.e2 ) == Double.doubleToLongBits( that.e2 )
447 && Double.doubleToLongBits( this.ep2 ) == Double.doubleToLongBits( that.ep2 );
448 }
449 return false;
450 }
451
452 /**
453 * Returns the WKT for this math transform.
454 */
455 public final String toString() {
456 return toString( "Ellipsoid_To_Geocentric" );
457 }
458
459 /**
460 * Returns the WKT for this math transform with the
461 * specified classification name. The classification
462 * name should be "Ellipsoid_To_Geocentric" or
463 * "Geocentric_To_Ellipsoid".
464 */
465 final String toString( final String classification ) {
466 final StringBuffer buffer = paramMT( classification );
467 addParameter( buffer, "semi_major", a );
468 addParameter( buffer, "semi_minor", b );
469 buffer.append( ']' );
470 return buffer.toString();
471 }
472
473 /**
474 * Inverse of a geocentric transform.
475 *
476 * @version 1.0
477 * @author Martin Desruisseaux
478 */
479 private final class Inverse extends AbstractMathTransform.Inverse {
480 public Inverse() {
481 GeocentricTransform.this.super();
482 }
483
484 public void transform( final double[] source, final int srcOffset, final double[] dest,
485 final int dstOffset, final int length )
486 throws TransformException {
487 GeocentricTransform.this.inverseTransform( source, srcOffset, dest, dstOffset, length );
488 }
489
490 public void transform( final float[] source, final int srcOffset, final float[] dest,
491 final int dstOffset, final int length )
492 throws TransformException {
493 GeocentricTransform.this.inverseTransform( source, srcOffset, dest, dstOffset, length );
494 }
495
496 public final String toString() {
497 return GeocentricTransform.this.toString( "Geocentric_To_Ellipsoid" );
498 }
499 }
500
501 /**
502 * The provider for {@link GeocentricTransform}.
503 *
504 * @version 1.0
505 * @author Martin Desruisseaux
506 */
507 static final class Provider extends MathTransformProvider {
508 /**
509 * The range of values for the dimension.
510 */
511 private static final Range DIM_RANGE = new Range( Integer.class, new Integer( 2 ),
512 new Integer( 3 ) );
513
514 /**
515 * <code>false</code> for the direct transform,
516 * or <code>true</code> for the inverse transform.
517 */
518 private final boolean inverse;
519
520 /**
521 * Create a provider.
522 *
523 * @param inverse <code>false</code> for the direct transform,
524 * or <code>true</code> for the inverse transform.
525 */
526 public Provider( final boolean inverse ) {
527 super( inverse ? "Geocentric_To_Ellipsoid" : "Ellipsoid_To_Geocentric",
528 0/*ResourceKeys.GEOCENTRIC_TRANSFORM*/, null ); // TODO
529 put( "semi_major", Double.NaN, POSITIVE_RANGE );
530 put( "semi_minor", Double.NaN, POSITIVE_RANGE );
531 putInt( "dim_geoCS", 3, DIM_RANGE ); // Custom parameter: NOT AN OPENGIS SPECIFICATION
532 this.inverse = inverse;
533 }
534
535 /**
536 * Returns a transform for the specified parameters.
537 *
538 * @param parameters The parameter values in standard units.
539 * @return A {@link MathTransform} object of this classification.
540 */
541 public MathTransform create( final ParameterList parameters ) {
542 final double semiMajor = parameters.getDoubleParameter( "semi_major" );
543 final double semiMinor = parameters.getDoubleParameter( "semi_minor" );
544 int dimGeographic = 3;
545 try {
546 dimGeographic = parameters.getIntParameter( "dim_geoCS" );
547 } catch ( IllegalArgumentException exception ) {
548 // the "dim_geoCS" parameter is a custom one required
549 // by our SEAGIS implementation. It is NOT an OpenGIS
550 // one. We can't require clients to know it.
551 }
552 GeocentricTransform transform = new GeocentricTransform( semiMajor, semiMinor,
553 Unit.METRE, dimGeographic != 2 );
554 return ( inverse ) ? transform.inverse() : transform;
555 }
556 }
557 }