001    //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_testing/src/org/deegree/crs/transformations/helmert/Helmert.java $
002    /*----------------------------------------------------------------------------
003     This file is part of deegree, http://deegree.org/
004     Copyright (C) 2001-2009 by:
005       Department of Geography, University of Bonn
006     and
007       lat/lon GmbH
008    
009     This library is free software; you can redistribute it and/or modify it under
010     the terms of the GNU Lesser General Public License as published by the Free
011     Software Foundation; either version 2.1 of the License, or (at your option)
012     any later version.
013     This library is distributed in the hope that it will be useful, but WITHOUT
014     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
015     FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
016     details.
017     You should have received a copy of the GNU Lesser General Public License
018     along with this library; if not, write to the Free Software Foundation, Inc.,
019     59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
020    
021     Contact information:
022    
023     lat/lon GmbH
024     Aennchenstr. 19, 53177 Bonn
025     Germany
026     http://lat-lon.de/
027    
028     Department of Geography, University of Bonn
029     Prof. Dr. Klaus Greve
030     Postfach 1147, 53001 Bonn
031     Germany
032     http://www.geographie.uni-bonn.de/deegree/
033    
034     e-mail: info@deegree.org
035    ----------------------------------------------------------------------------*/
036    package org.deegree.crs.transformations.helmert;
037    
038    import static org.deegree.crs.projections.ProjectionUtils.EPS11;
039    
040    import java.util.List;
041    
042    import javax.vecmath.Matrix4d;
043    import javax.vecmath.Point3d;
044    
045    import org.deegree.crs.Identifiable;
046    import org.deegree.crs.coordinatesystems.CoordinateSystem;
047    import org.deegree.crs.exceptions.TransformationException;
048    import org.deegree.crs.transformations.Transformation;
049    
050    /**
051     * Parameters for a geographic transformation into another datum. The Bursa Wolf parameters should be applied to
052     * geocentric coordinates, where the X axis points towards the Greenwich Prime Meridian, the Y axis points East, and the
053     * Z axis points North.
054     *
055     *
056     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
057     *
058     * @author last edited by: $Author: mschneider $
059     *
060     * @version $Revision: 18195 $, $Date: 2009-06-18 17:55:39 +0200 (Do, 18 Jun 2009) $
061     *
062     */
063    public class Helmert extends Transformation {
064    
065        private static final long serialVersionUID = -1490341500541671907L;
066    
067        /** Bursa Wolf shift in meters. */
068        public double dx;
069    
070        /** Bursa Wolf shift in meters. */
071        public double dy;
072    
073        /** Bursa Wolf shift in meters. */
074        public double dz;
075    
076        /** Bursa Wolf rotation in arc seconds, which is 1/3600 of a degree. */
077        public double ex;
078    
079        /** Bursa Wolf rotation in arc seconds. */
080        public double ey;
081    
082        /** Bursa Wolf rotation in arc seconds. */
083        public double ez;
084    
085        /** Bursa Wolf scaling in parts per million. */
086        public double ppm;
087    
088        private Matrix4d transformMatrix = null;
089    
090        private Matrix4d inverseMatrix = null;
091    
092        private boolean rotationInRadians = false;
093    
094        /**
095         * @param dx
096         *            Bursa Wolf shift in meters.
097         * @param dy
098         *            Bursa Wolf shift in meters.
099         * @param dz
100         *            Bursa Wolf shift in meters.
101         * @param ex
102         *            Bursa Wolf rotation in arc seconds or in radians (by setting the flag).
103         * @param ey
104         *            Bursa Wolf rotation in arc seconds or in radians (by setting the flag).
105         * @param ez
106         *            Bursa Wolf rotation in arc seconds or in radians (by setting the flag).
107         * @param ppm
108         *            Bursa Wolf scaling in parts per million.
109         * @param sourceCRS
110         *            of this helmert transformation
111         * @param targetCRS
112         *            of this helmert transformation
113         * @param identifiable
114         *            object containing all relevant id.
115         * @param inRadians
116         *            true if the rotation parameters are in radians
117         */
118        public Helmert( double dx, double dy, double dz, double ex, double ey, double ez, double ppm,
119                        CoordinateSystem sourceCRS, CoordinateSystem targetCRS, Identifiable identifiable, boolean inRadians ) {
120            super( sourceCRS, targetCRS, identifiable );
121            this.dx = dx;
122            this.dy = dy;
123            this.dz = dz;
124            this.ex = ex;
125            this.ey = ey;
126            this.ez = ez;
127            this.ppm = ppm;
128            rotationInRadians = inRadians;
129        }
130    
131        /**
132         * @param dx
133         *            Bursa Wolf shift in meters.
134         * @param dy
135         *            Bursa Wolf shift in meters.
136         * @param dz
137         *            Bursa Wolf shift in meters.
138         * @param ex
139         *            Bursa Wolf rotation in arc seconds.
140         * @param ey
141         *            Bursa Wolf rotation in arc seconds.
142         * @param ez
143         *            Bursa Wolf rotation in arc seconds.
144         * @param ppm
145         *            Bursa Wolf scaling in parts per million.
146         * @param sourceCRS
147         *            of this helmert transformation
148         * @param targetCRS
149         *            of this helmert transformation
150         * @param identifiable
151         *            object containing all relevant id.
152         */
153        public Helmert( double dx, double dy, double dz, double ex, double ey, double ez, double ppm,
154                        CoordinateSystem sourceCRS, CoordinateSystem targetCRS, Identifiable identifiable ) {
155            this( dx, dy, dz, ex, ey, ez, ppm, sourceCRS, targetCRS, identifiable, false );
156        }
157    
158        /**
159         * Construct a conversion info with all parameters set to 0;
160         *
161         * @param sourceCRS
162         *            of this helmert transformation
163         * @param targetCRS
164         *            of this helmert transformation
165         *
166         * @param identifiers
167         * @param names
168         * @param versions
169         * @param descriptions
170         * @param areasOfUse
171         */
172        public Helmert( CoordinateSystem sourceCRS, CoordinateSystem targetCRS, String[] identifiers, String[] names,
173                        String[] versions, String[] descriptions, String[] areasOfUse ) {
174            this( 0, 0, 0, 0, 0, 0, 0, sourceCRS, targetCRS, new Identifiable( identifiers, names, versions, descriptions,
175                                                                               areasOfUse ) );
176        }
177    
178        /**
179         * Construct a conversion info with all parameters set to 0;
180         *
181         * @param sourceCRS
182         *            of this helmert transformation
183         * @param targetCRS
184         *            of this helmert transformation
185         * @param identifier
186         */
187        public Helmert( CoordinateSystem sourceCRS, CoordinateSystem targetCRS, String identifier ) {
188            this( sourceCRS, targetCRS, new String[] { identifier } );
189        }
190    
191        /**
192         * Construct a conversion info with all parameters set to 0;
193         *
194         * @param sourceCRS
195         *            of this helmert transformation
196         * @param targetCRS
197         *            of this helmert transformation
198         * @param identifiers
199         */
200        public Helmert( CoordinateSystem sourceCRS, CoordinateSystem targetCRS, String[] identifiers ) {
201            this( sourceCRS, targetCRS, identifiers, null, null, null, null );
202        }
203    
204        /**
205         * @param dx
206         *            Bursa Wolf shift in meters.
207         * @param dy
208         *            Bursa Wolf shift in meters.
209         * @param dz
210         *            Bursa Wolf shift in meters.
211         * @param ex
212         *            Bursa Wolf rotation in arc seconds.
213         * @param ey
214         *            Bursa Wolf rotation in arc seconds.
215         * @param ez
216         *            Bursa Wolf rotation in arc seconds.
217         * @param ppm
218         *            Bursa Wolf scaling in parts per million.
219         * @param sourceCRS
220         *            of this helmert transformation
221         * @param targetCRS
222         *            of this helmert transformation
223         * @param identifiers
224         * @param names
225         * @param versions
226         * @param descriptions
227         * @param areaOfUses
228         */
229        public Helmert( double dx, double dy, double dz, double ex, double ey, double ez, double ppm,
230                        CoordinateSystem sourceCRS, CoordinateSystem targetCRS, String[] identifiers, String[] names,
231                        String[] versions, String[] descriptions, String[] areaOfUses ) {
232            this( dx, dy, dz, ex, ey, ez, ppm, sourceCRS, targetCRS, new Identifiable( identifiers, names, versions,
233                                                                                       descriptions, areaOfUses ) );
234        }
235    
236        /**
237         * @param dx
238         *            Bursa Wolf shift in meters.
239         * @param dy
240         *            Bursa Wolf shift in meters.
241         * @param dz
242         *            Bursa Wolf shift in meters.
243         * @param ex
244         *            Bursa Wolf rotation in arc seconds.
245         * @param ey
246         *            Bursa Wolf rotation in arc seconds.
247         * @param ez
248         *            Bursa Wolf rotation in arc seconds.
249         * @param ppm
250         *            Bursa Wolf scaling in parts per million.
251         * @param sourceCRS
252         *            of this helmert transformation
253         * @param targetCRS
254         *            of this helmert transformation
255         * @param identifier
256         * @param name
257         * @param version
258         * @param description
259         * @param areaOfUse
260         */
261        public Helmert( double dx, double dy, double dz, double ex, double ey, double ez, double ppm,
262                        CoordinateSystem sourceCRS, CoordinateSystem targetCRS, String identifier, String name,
263                        String version, String description, String areaOfUse ) {
264            this( dx, dy, dz, ex, ey, ez, ppm, sourceCRS, targetCRS, new String[] { identifier }, new String[] { name },
265                  new String[] { version }, new String[] { description }, new String[] { areaOfUse } );
266        }
267    
268        /**
269         * @param dx
270         *            Bursa Wolf shift in meters.
271         * @param dy
272         *            Bursa Wolf shift in meters.
273         * @param dz
274         *            Bursa Wolf shift in meters.
275         * @param ex
276         *            Bursa Wolf rotation in arc seconds.
277         * @param ey
278         *            Bursa Wolf rotation in arc seconds.
279         * @param ez
280         *            Bursa Wolf rotation in arc seconds.
281         * @param ppm
282         *            Bursa Wolf scaling in parts per million.
283         * @param sourceCRS
284         *            of this helmert transformation
285         * @param targetCRS
286         *            of this helmert transformation
287         * @param identifiers
288         */
289        public Helmert( double dx, double dy, double dz, double ex, double ey, double ez, double ppm,
290                        CoordinateSystem sourceCRS, CoordinateSystem targetCRS, String[] identifiers ) {
291            this( dx, dy, dz, ex, ey, ez, ppm, sourceCRS, targetCRS, identifiers, null, null, null, null );
292        }
293    
294        /**
295         * @param dx
296         *            Bursa Wolf shift in meters.
297         * @param dy
298         *            Bursa Wolf shift in meters.
299         * @param dz
300         *            Bursa Wolf shift in meters.
301         * @param ex
302         *            Bursa Wolf rotation in arc seconds.
303         * @param ey
304         *            Bursa Wolf rotation in arc seconds.
305         * @param ez
306         *            Bursa Wolf rotation in arc seconds.
307         * @param ppm
308         *            Bursa Wolf scaling in parts per million.
309         * @param sourceCRS
310         *            of this helmert transformation
311         * @param targetCRS
312         *            of this helmert transformation
313         * @param identifier
314         */
315        public Helmert( double dx, double dy, double dz, double ex, double ey, double ez, double ppm,
316                        CoordinateSystem sourceCRS, CoordinateSystem targetCRS, String identifier ) {
317            this( dx, dy, dz, ex, ey, ez, ppm, sourceCRS, targetCRS, new String[] { identifier } );
318        }
319    
320        /**
321         * Returns an affine transformation also known as the "Helmert" transformation. The matrix representation of this
322         * transformation (also known as "Bursa Wolf" formula) is as follows:
323         *
324         * <blockquote>
325         *
326         * <pre>
327         *       S = 1 + {@link #ppm}*1E-6
328         *
329         *       [ X ]     [ S          -{@link #ez}*S  +{@link #ey}*S   {@link #dx} ]  [ X ]
330         *       [ Y ]  = [ +{@link #ez}*S  S          -{@link #ex}*S   {@link #dy} ]  [ Y ]
331         *       [ Z ]     [ -{@link #ey}*S   +{@link #ex}*S   S         {@link #dz} ]  [ Z ]
332         *       [ 1 ]     [ 0           0           0           1 ]  [ 1 ]
333         * </pre>
334         *
335         * </blockquote>
336         *
337         * This affine transform can be applied to transform <code>geocentric</code> coordinates from one datum into
338         * <code>geocentric</code> coordinates of an other datum. see <a
339         * href="http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs35.html#CS3523_helmert">
340         * http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs35.html</a> for more information.
341         *
342         * @return the affine "Helmert" transformation as a Matrix4d.
343         */
344        public Matrix4d getAsAffineTransform() {
345            if ( transformMatrix != null ) {
346                return new Matrix4d( transformMatrix );
347            }
348            // Note: (ex, ey, ez) is a rotation in arc seconds. We need to convert it into radians (the
349            // R factor in RS).
350            final double S = 1 + ( ppm * 1E-6 );
351            double arcToRad = .00000484813681109535; // Math.PI / ( 180. * 3600. )
352            if ( rotationInRadians ) {
353                arcToRad = 1;
354            }
355            final double RS = arcToRad * S;
356    
357            transformMatrix = new Matrix4d( S, -ez * RS, +ey * RS, dx, +ez * RS, S, -ex * RS, dy, -ey * RS, +ex * RS, S,
358                                            dz, 0, 0, 0, 1. );
359            return new Matrix4d( transformMatrix );
360        }
361    
362        /**
363         * @return true if any of the helmert parameters were set.
364         */
365        public boolean hasValues() {
366            return !( ex == 0 && ey == 0 && ez == 0 && dx == 0 && dy == 0 && dz == 0 && ppm == 0 );
367        }
368    
369        @Override
370        public boolean equals( final Object other ) {
371            if ( other != null && other instanceof Helmert ) {
372                final Helmert that = (Helmert) other;
373                return ( Math.abs( this.dx - that.dx ) < EPS11 ) && ( Math.abs( this.dy - that.dy ) < EPS11 )
374                       && ( Math.abs( this.dz - that.dz ) < EPS11 ) && ( Math.abs( this.ex - that.ex ) < EPS11 )
375                       && ( Math.abs( this.ey - that.ey ) < EPS11 ) && ( Math.abs( this.ez - that.ez ) < EPS11 )
376                       && ( Math.abs( this.ppm - that.ppm ) < EPS11 ) && super.equals( that );
377    
378            }
379            return false;
380        }
381    
382        /**
383         * Returns the Well Know Text (WKT) for this object. The WKT is part of OpenGIS's specification and looks like
384         * <code>TOWGS84[dx, dy, dz, ex, ey, ez, ppm]</code>.
385         *
386         * @return the Well Know Text (WKT) for this object.
387         */
388        @Override
389        public String toString() {
390            final StringBuffer buffer = new StringBuffer( super.getIdAndName() );
391            buffer.append( "\n[\"" );
392            buffer.append( dx );
393            buffer.append( ", " );
394            buffer.append( dy );
395            buffer.append( ", " );
396            buffer.append( dz );
397            buffer.append( ", " );
398            buffer.append( ex );
399            buffer.append( ", " );
400            buffer.append( ey );
401            buffer.append( ", " );
402            buffer.append( ez );
403            buffer.append( ", " );
404            buffer.append( ppm );
405            buffer.append( ']' );
406            return buffer.toString();
407        }
408    
409        /**
410         * Implementation as proposed by Joshua Block in Effective Java (Addison-Wesley 2001), which supplies an even
411         * distribution and is relatively fast. It is created from field <b>f</b> as follows:
412         * <ul>
413         * <li>boolean -- code = (f ? 0 : 1)</li>
414         * <li>byte, char, short, int -- code = (int)f</li>
415         * <li>long -- code = (int)(f ^ (f &gt;&gt;&gt;32))</li>
416         * <li>float -- code = Float.floatToIntBits(f);</li>
417         * <li>double -- long l = Double.doubleToLongBits(f); code = (int)(l ^ (l &gt;&gt;&gt; 32))</li>
418         * <li>all Objects, (where equals(&nbsp;) calls equals(&nbsp;) for this field) -- code = f.hashCode(&nbsp;)</li>
419         * <li>Array -- Apply above rules to each element</li>
420         * </ul>
421         * <p>
422         * Combining the hash code(s) computed above: result = 37 * result + code;
423         * </p>
424         *
425         * @return (int) ( result >>> 32 ) ^ (int) result;
426         *
427         * @see java.lang.Object#hashCode()
428         */
429        @Override
430        public int hashCode() {
431            // the 2nd millionth prime, :-)
432            long code = 32452843;
433            long tmp = Double.doubleToLongBits( dx );
434            code = code * 37 + (int) ( tmp ^ ( tmp >>> 32 ) );
435    
436            tmp = Double.doubleToLongBits( dy );
437            code = code * 37 + (int) ( tmp ^ ( tmp >>> 32 ) );
438    
439            tmp = Double.doubleToLongBits( dz );
440            code = code * 37 + (int) ( tmp ^ ( tmp >>> 32 ) );
441    
442            tmp = Double.doubleToLongBits( ex );
443            code = code * 37 + (int) ( tmp ^ ( tmp >>> 32 ) );
444    
445            tmp = Double.doubleToLongBits( ey );
446            code = code * 37 + (int) ( tmp ^ ( tmp >>> 32 ) );
447    
448            tmp = Double.doubleToLongBits( ez );
449            code = code * 37 + (int) ( tmp ^ ( tmp >>> 32 ) );
450    
451            tmp = Double.doubleToLongBits( ppm );
452            code = code * 37 + (int) ( tmp ^ ( tmp >>> 32 ) );
453            return (int) ( code >>> 32 ) ^ (int) code;
454        }
455    
456        @Override
457        public synchronized List<Point3d> doTransform( List<Point3d> srcPts )
458                                throws TransformationException {
459    
460            if ( srcPts == null || srcPts.size() == 0 ) {
461                return srcPts;
462            }
463    
464            // lazy instantiation
465            if ( transformMatrix == null ) {
466                transformMatrix = getAsAffineTransform();
467            }
468            Matrix4d matrix = transformMatrix;
469    
470            // create the inverse matrix
471            if ( isInverseTransform() ) {
472                if ( inverseMatrix == null ) {
473                    transformMatrix.invert( inverseMatrix );
474                }
475                matrix = inverseMatrix;
476            }
477            for ( Point3d p : srcPts ) {
478                matrix.transform( p );
479            }
480    
481            return srcPts;
482        }
483    
484        @Override
485        public String getImplementationName() {
486            return "Helmert";
487        }
488    
489        @Override
490        public boolean isIdentity() {
491            return hasValues();
492        }
493    }