    001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.4_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    }