001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.3_testing/src/org/deegree/crs/transformations/coordinate/MatrixTransform.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.coordinate;
037    
038    import static org.deegree.crs.projections.ProjectionUtils.EPS11;
039    
040    import java.util.List;
041    
042    import javax.media.jai.PerspectiveTransform;
043    import javax.vecmath.GMatrix;
044    import javax.vecmath.Matrix3d;
045    import javax.vecmath.Matrix4d;
046    import javax.vecmath.Point3d;
047    
048    import org.deegree.crs.Identifiable;
049    import org.deegree.crs.coordinatesystems.CoordinateSystem;
050    
051    /**
052     *
053     * The <code>MatrixTransform</code> class allows transformations using matrices. Although technically n &times; m
054     * matrices are possible, at the moment only 2 &times; 2, 3 &times; 3 and 4 &times; 4 matrices are supported.
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 MatrixTransform extends CRSTransformation {
064    
065        /**
066         * Serial number for interoperability with different versions.
067         */
068        private static final long serialVersionUID = -2104496465933824935L;
069    
070        /**
071         * the number of rows.
072         */
073        private final int numRow;
074    
075        /**
076         * the number of columns.
077         */
078        private final int numCol;
079    
080        private GMatrix matrix = null;
081    
082        private GMatrix invertMatrix = null;
083    
084        private Matrix3d matrix3D = null;
085    
086        private Matrix3d invertMatrix3D = null;
087    
088        private Matrix4d matrix4D = null;
089    
090        private Matrix4d invertMatrix4D = null;
091    
092        private String transformationName = "Matrix-Transform";
093    
094        /**
095         * Set super values and numRow,numCol.
096         *
097         * @param source
098         * @param target
099         * @param numRow
100         * @param numCol
101         */
102        private MatrixTransform( CoordinateSystem source, CoordinateSystem target, int numRow, int numCol, Identifiable id ) {
103            super( source, target, id );
104            this.numCol = numCol;
105            this.numRow = numRow;
106        }
107    
108        /**
109         * Construct a transform.
110         *
111         * @param source
112         *            the source coordinate system
113         * @param target
114         *            the target coordinate system.
115         *
116         * @param matrix
117         * @param id
118         *            an identifiable instance containing information about this transformation
119         */
120        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, final GMatrix matrix, Identifiable id ) {
121            this( source, target, matrix.getNumRow(), matrix.getNumCol(), id );
122            if ( numCol == numRow ) {
123                if ( numCol == 3 ) {
124                    this.matrix3D = new Matrix3d();
125                    matrix.get( this.matrix3D );
126                    invertMatrix3D = new Matrix3d();
127                    invertMatrix3D.invert( matrix3D );
128                }
129                if ( numCol == 4 ) {
130                    this.matrix4D = new Matrix4d();
131                    matrix.get( this.matrix4D );
132                    invertMatrix4D = new Matrix4d();
133                    invertMatrix4D.invert( matrix4D );
134                }
135            } else {
136                this.matrix = new GMatrix( matrix );
137                invertMatrix = new GMatrix( matrix );
138                invertMatrix.invert();
139            }
140    
141        }
142    
143        /**
144         * Construct a transform.
145         *
146         * @param source
147         *            the source coordinate system
148         * @param target
149         *            the target coordinate system.
150         *
151         * @param matrix
152         */
153        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, final GMatrix matrix ) {
154            this( source, target, matrix, new Identifiable( createFromTo( source.getIdentifier(), target.getIdentifier() ) ) );
155        }
156    
157        /**
158         * Construct a 3d transform.
159         *
160         * @param source
161         *            the source coordinate system
162         * @param target
163         *            the target coordinate system.
164         *
165         * @param matrix
166         * @param id
167         *            an identifiable instance containing information about this transformation
168         */
169        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, final Matrix3d matrix, Identifiable id ) {
170            this( source, target, 3, 3, id );
171            this.matrix3D = new Matrix3d( matrix );
172            invertMatrix3D = new Matrix3d();
173            invertMatrix3D.invert( matrix3D );
174        }
175    
176        /**
177         * Construct a 3d transform.
178         *
179         * @param source
180         *            the source coordinate system
181         * @param target
182         *            the target coordinate system.
183         *
184         * @param matrix
185         */
186        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, final Matrix3d matrix ) {
187            this( source, target, matrix, new Identifiable( createFromTo( source.getIdentifier(), target.getIdentifier() ) ) );
188        }
189    
190        /**
191         * Construct a 4d transform.
192         *
193         * @param source
194         *            the source coordinate system
195         * @param target
196         *            the target coordinate system.
197         *
198         * @param matrix
199         * @param id
200         *            an identifiable instance containing information about this transformation
201         */
202        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, Matrix4d matrix, Identifiable id ) {
203            this( source, target, 4, 4, id );
204            matrix4D = new Matrix4d( matrix );
205            invertMatrix4D = new Matrix4d();
206            invertMatrix4D.invert( matrix4D );
207        }
208    
209        /**
210         * Construct a 4d transform.
211         *
212         * @param source
213         *            the source coordinate system
214         * @param target
215         *            the target coordinate system.
216         *
217         * @param matrix
218         */
219        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, Matrix4d matrix ) {
220            this( source, target, matrix, new Identifiable( createFromTo( source.getIdentifier(), target.getIdentifier() ) ) );
221        }
222    
223        /**
224         * Construct a 4d transform.
225         *
226         * @param source
227         *            the source coordinate system
228         * @param target
229         *            the target coordinate system.
230         *
231         * @param matrix
232         * @param transformationName
233         *            the 'optional' name of the transformation, which is useful to specify the 'helmert' transformation.
234         * @param id
235         *            an identifiable instance containing information about this transformation
236         */
237        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, Matrix4d matrix,
238                                String transformationName, Identifiable id ) {
239            this( source, target, matrix, id );
240            if ( transformationName != null ) {
241                this.transformationName = transformationName;
242            }
243        }
244    
245        /**
246         * Construct a 4d transform.
247         *
248         * @param source
249         *            the source coordinate system
250         * @param target
251         *            the target coordinate system.
252         *
253         * @param matrix
254         * @param transformationName
255         *            the 'optional' name of the transformation, which is useful to specify the 'helmert' transformation.
256         */
257        public MatrixTransform( CoordinateSystem source, CoordinateSystem target, Matrix4d matrix, String transformationName ) {
258            this( source, target, matrix, new Identifiable( createFromTo( source.getIdentifier(), target.getIdentifier() ) ) );
259        }
260    
261        @Override
262        public List<Point3d> doTransform( List<Point3d> srcPts ) {
263            if ( isIdentity() ) {
264                return srcPts;
265            }
266            // List<Point3d> results = new ArrayList<Point3d>( srcPts );
267            if ( isInverseTransform() ) {
268                if ( matrix3D != null ) {
269                    transform( invertMatrix3D, srcPts );
270                } else if ( matrix4D != null ) {
271                    transform( invertMatrix4D, srcPts );
272                } else {
273                    transform( invertMatrix, srcPts );
274                }
275            } else {
276                if ( matrix3D != null ) {
277                    transform( matrix3D, srcPts );
278                } else if ( matrix4D != null ) {
279                    transform( matrix4D, srcPts );
280                } else {
281                    transform( matrix, srcPts );
282                }
283            }
284            return srcPts;
285        }
286    
287        /**
288         * @return the dimension of input points.
289         */
290        public int getDimSource() {
291            return numCol - 1;
292        }
293    
294        /**
295         * @return the dimension of output points.
296         */
297        public int getDimTarget() {
298            return numRow - 1;
299        }
300    
301        /**
302         * @return true if this transformation holds an identity matrix (e.g. doesn't transform at all).
303         */
304        @Override
305        public boolean isIdentity() {
306            if ( numRow != numCol ) {
307                return false;
308            }
309            for ( int row = 0; row < numRow; row++ ) {
310                for ( int col = 0; col < numCol; col++ ) {
311                    double value = ( matrix3D != null ) ? matrix3D.getElement( row, col )
312                                                       : ( ( matrix4D != null ) ? matrix4D.getElement( row, col )
313                                                                               : matrix.getElement( row, col ) );
314                    if ( Math.abs( value - ( col == row ? 1 : 0 ) ) > EPS11 ) {
315                        return false;
316                    }
317                }
318            }
319            return true;
320        }
321    
322        @Override
323        public boolean equals( final Object object ) {
324            if ( object == this ) {
325                return true; // Slight optimization
326            }
327            if ( object != null && super.equals( object ) ) {
328                return matrix.equals( ( (MatrixTransform) object ).matrix );
329            }
330            return false;
331        }
332    
333        /**
334         * Transforms an array of floating point coordinates by this matrix. Because of the usage of a point3d some
335         * assumptions are made,
336         * <ol>
337         * <li>The result point will have the dimension of the number of rows-1.</li>
338         * <li>Only the first number of Columns -1 (=input dimension) are used from the incoming point to calculate the
339         * result.</li>
340         * <li>An output dimension > 3 (e.g. the number of rows > 4 ) results in an IllegalArgumentException.</li>
341         * <li>An input dimension > 3 (e.g. the number of columns > 4 ) results in an IllegalArgumentException.</li>
342         * <li>Only inputDimension values of the incoming point will be taken into account.</li>
343         * <li>If the input dimension == 4 the fourth value is assumed to be 1.</li>
344         * </ol>
345         *
346         * <p>
347         * For example, a square matrix of size 4&times;4. Is a three-dimensional transformation for incoming and outgoing
348         * coordinates. The transformed points <code>(x',y',z')</code> are computed as below (note that this computation is similar to
349         * {@link PerspectiveTransform}):
350         *
351         * <blockquote> <code>
352         * <pre>
353         *    [ a ]     [ m&lt;sub&gt;00&lt;/sub&gt;  m&lt;sub&gt;01&lt;/sub&gt;  m&lt;sub&gt;02&lt;/sub&gt;  m&lt;sub&gt;03&lt;/sub&gt; ] [ x ]
354         *    [ b ]  =  [ m&lt;sub&gt;10&lt;/sub&gt;  m&lt;sub&gt;11&lt;/sub&gt;  m&lt;sub&gt;12&lt;/sub&gt;  m&lt;sub&gt;13&lt;/sub&gt; ] [ y ]
355         *    [ c ]     [ m&lt;sub&gt;20&lt;/sub&gt;  m&lt;sub&gt;21&lt;/sub&gt;  m&lt;sub&gt;22&lt;/sub&gt;  m&lt;sub&gt;23&lt;/sub&gt; ] [ z ]
356         *    [ w ]     [ m&lt;sub&gt;30&lt;/sub&gt;  m&lt;sub&gt;31&lt;/sub&gt;  m&lt;sub&gt;32&lt;/sub&gt;  m&lt;sub&gt;33&lt;/sub&gt; ] [ 1 ]
357         *
358         *     x' = a/w
359         *     y' = b/w
360         *     z' = c/w
361         * </pre>
362         * </code> </blockquote>
363         *
364         * @param srcPts
365         *            list containing the source point coordinates.
366         */
367        private void transform( GMatrix gm, List<Point3d> srcPts ) {
368            final int inputDimension = numCol - 1;
369            final int outputDimension = numRow - 1;
370            if ( inputDimension > 3 ) {
371                throw new IllegalArgumentException(
372                                                    "Number of collumns: "
373                                                                            + numCol
374                                                                            + " of the given matrix exceed the maximum dimension (3) supported by this Transformation" );
375            }
376            if ( outputDimension > 3 ) {
377                throw new IllegalArgumentException(
378                                                    "Number of rows: "
379                                                                            + numRow
380                                                                            + " of the given matrix exceed the maximum dimension (3) supported by this Transformation" );
381            }
382    
383            final double[] tmpPoint = new double[numRow];
384            for ( Point3d p : srcPts ) {
385                for ( int row = 0; row < numRow; ++row ) {
386                    tmpPoint[row] = gm.getElement( row, 0 ) * p.x;
387                    if ( numCol >= 2 ) {
388                        tmpPoint[row] += gm.getElement( row, 1 ) * p.y;
389                        if ( numCol >= 3 ) {
390                            tmpPoint[row] += gm.getElement( row, 2 )
391                                             * ( ( !Double.isNaN( p.z ) && !Double.isInfinite( p.z ) ) ? p.z : 1 );
392                            if ( numCol == 4 ) { // assume 1
393                                tmpPoint[row] += gm.getElement( row, 3 );
394                            }
395                        }
396    
397                    }
398                }
399                final double w = tmpPoint[outputDimension];
400                if ( outputDimension >= 1 ) {
401                    p.x = tmpPoint[0] / w;
402                    if ( outputDimension >= 2 ) {
403                        p.y = tmpPoint[1] / w;
404                        if ( outputDimension == 3 ) {
405                            p.z = tmpPoint[2] / w;
406                        }
407                    }
408                }
409            }
410    
411        }
412    
413        /**
414         * Use the given GMatrix to transform the given points inplace.
415         *
416         * @param m4d
417         *            the matrix to use (e.g. the inverse matrix or the forward matrix.
418         * @param srcPts
419         *            The array containing the source point coordinates.
420         */
421        private void transform( Matrix4d m4d, List<Point3d> srcPts ) {
422            for ( Point3d p : srcPts ) {
423                m4d.transform( p );
424            }
425        }
426    
427        /**
428         * Use the given GMatrix to transform the given points in-place.
429         *
430         * @param m3d
431         *            the matrix to use (e.g. the inverse matrix or the forward matrix).
432         * @param srcPts
433         *            The array containing the source point coordinates.
434         */
435        private void transform( Matrix3d m3d, List<Point3d> srcPts ) {
436            for ( Point3d p : srcPts ) {
437    
438                boolean zIsNaN = Double.isNaN( p.z );
439                if ( zIsNaN ) {
440                    p.z = 1;
441                }
442                m3d.transform( p );
443                if ( zIsNaN ) {
444                    p.z = Double.NaN;
445                }
446            }
447        }
448    
449        /**
450         * @return the matrix.
451         */
452        public final GMatrix getMatrix() {
453            if ( matrix != null ) {
454                return isInverseTransform() ? invertMatrix : matrix;
455            }
456            GMatrix result = new GMatrix( numRow, numCol );
457            if ( matrix3D != null ) {
458                if ( isInverseTransform() ) {
459                    result.set( invertMatrix3D );
460                } else {
461                    result.set( matrix3D );
462                }
463            }
464            if ( matrix4D != null ) {
465                if ( isInverseTransform() ) {
466                    result.set( invertMatrix4D );
467                } else {
468                    result.set( matrix4D );
469                }
470            }
471            return result;
472        }
473    
474        @Override
475        public String getImplementationName() {
476            return transformationName;
477        }
478    
479    }