001    //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_testing/src/org/deegree/crs/transformations/polynomial/LeastSquareApproximation.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    
037    package org.deegree.crs.transformations.polynomial;
038    
039    import static org.deegree.crs.projections.ProjectionUtils.EPS11;
040    
041    import java.awt.geom.Point2D;
042    import java.util.ArrayList;
043    import java.util.List;
044    
045    import javax.media.jai.WarpCubic;
046    import javax.media.jai.WarpGeneralPolynomial;
047    import javax.media.jai.WarpPolynomial;
048    import javax.media.jai.WarpQuadratic;
049    import javax.vecmath.Point3d;
050    
051    import org.deegree.crs.Identifiable;
052    import org.deegree.crs.coordinatesystems.CoordinateSystem;
053    import org.deegree.crs.exceptions.TransformationException;
054    import org.deegree.framework.log.ILogger;
055    import org.deegree.framework.log.LoggerFactory;
056    import org.deegree.i18n.Messages;
057    
058    /**
059     * <code>LeastSquareApproximation</code> is a polynomial transformation which uses the least square method to
060     * approximate a function given by some measured values.
061     *
062     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
063     *
064     * @author last edited by: $Author: mschneider $
065     *
066     * @version $Revision: 18195 $, $Date: 2009-06-18 17:55:39 +0200 (Do, 18 Jun 2009) $
067     *
068     */
069    public class LeastSquareApproximation extends PolynomialTransformation {
070    
071        private static final long serialVersionUID = -8228847133699923749L;
072    
073        private static ILogger LOG = LoggerFactory.getLogger( LeastSquareApproximation.class );
074    
075        private WarpPolynomial leastSquarePolynomial;
076    
077        private final int order;
078    
079        private final float scaleX;
080    
081        private final float scaleY;
082    
083        /**
084         * @param firstParameters
085         *            of the polynomial
086         * @param secondParameters
087         *            of the polynomial
088         * @param sourceCRS
089         *            of this transformation
090         * @param targetCRS
091         *            of this transformation
092         * @param scaleX
093         *            to apply to incoming data's x value, if 1 (or 0) no scale will be applied.
094         * @param scaleY
095         *            to apply to incoming data's y value, if 1 (or 0) no scale will be applied.
096         * @param id
097         *            an identifiable instance containing information about this transformation
098         */
099        public LeastSquareApproximation( List<Double> firstParameters, List<Double> secondParameters,
100                                         CoordinateSystem sourceCRS, CoordinateSystem targetCRS, float scaleX,
101                                         float scaleY, Identifiable id ) {
102            super( firstParameters, secondParameters, sourceCRS, targetCRS, id );
103            if ( getSecondParams().size() != getFirstParams().size() ) {
104                throw new IllegalArgumentException( "The given parameter lists do not have equal length" );
105            }
106            // calc from (n+1)*(n+2) = 2*size;
107            order = (int) Math.floor( ( -3 + Math.sqrt( 9 + ( 4 * ( getFirstParams().size() * 2 ) - 2 ) ) ) * 0.5 );
108            float[] aParams = new float[getFirstParams().size()];
109            for ( int i = 0; i < firstParameters.size(); ++i ) {
110                aParams[i] = firstParameters.get( i ).floatValue();
111            }
112            float[] bParams = new float[getSecondParams().size()];
113            for ( int i = 0; i < secondParameters.size(); ++i ) {
114                bParams[i] = secondParameters.get( i ).floatValue();
115            }
116            if ( Float.isNaN( scaleX ) || Math.abs( scaleX ) < EPS11 ) {
117                scaleX = 1;
118            }
119            this.scaleX = scaleX;
120    
121            if ( Float.isNaN( scaleY ) || Math.abs( scaleY ) < EPS11 ) {
122                scaleY = 1;
123            }
124            this.scaleY = scaleY;
125            switch ( order ) {
126            case 2:
127                leastSquarePolynomial = new WarpQuadratic( aParams, bParams, this.scaleX, this.scaleY, 1f / this.scaleX,
128                                                           1f / this.scaleY );
129                break;
130            case 3:
131    
132                leastSquarePolynomial = new WarpCubic( aParams, bParams, this.scaleX, this.scaleY, 1f / this.scaleX,
133                                                       1f / this.scaleY );
134                break;
135            default:
136                leastSquarePolynomial = new WarpGeneralPolynomial( aParams, bParams, this.scaleX, this.scaleY,
137                                                                   1f / this.scaleX, 1f / this.scaleY );
138                break;
139            }
140        }
141    
142        /**
143         * Sets the id to EPSG::9645 ( General polynomial of degree 2 ).
144         *
145         * @param firstParameters
146         *            of the polynomial
147         * @param secondParameters
148         *            of the polynomial
149         * @param sourceCRS
150         *            of this transformation
151         * @param targetCRS
152         *            of this transformation
153         * @param scaleX
154         *            to apply to incoming data's x value, if 1 (or 0) no scale will be applied.
155         * @param scaleY
156         *            to apply to incoming data's y value, if 1 (or 0) no scale will be applied.
157         */
158        public LeastSquareApproximation( List<Double> firstParameters, List<Double> secondParameters,
159                                         CoordinateSystem sourceCRS, CoordinateSystem targetCRS, float scaleX, float scaleY ) {
160            this( firstParameters, secondParameters, sourceCRS, targetCRS, scaleX, scaleY, new Identifiable( "EPSG::9645" ) );
161        }
162    
163        @Override
164        public List<Point3d> applyPolynomial( List<Point3d> srcPts )
165                                throws TransformationException {
166            if ( srcPts == null || srcPts.size() == 0 ) {
167                return srcPts;
168            }
169            List<Point3d> result = new ArrayList<Point3d>( srcPts.size() );
170            for ( Point3d p : srcPts ) {
171                Point2D r = leastSquarePolynomial.mapDestPoint( new Point2D.Double( p.x, p.y ) );
172                if ( r != null ) {
173                    result.add( new Point3d( r.getX(), r.getY(), p.z ) );
174                } else {
175                    throw new TransformationException( Messages.getMessage( "CRS_POLYNOMIAL_TRANSFORM_ERROR", p.toString() ) );
176                }
177            }
178            return result;
179        }
180    
181        @Override
182        public String getImplementationName() {
183            return "leastsquare";
184        }
185    
186        @Override
187        public float[][] createVariables( List<Point3d> originalPoints, List<Point3d> projectedPoints, int polynomalOrder ) {
188            float[] sourceCoords = new float[originalPoints.size() * 2];
189            int count = 0;
190            double minX = Double.MAX_VALUE;
191            double minY = Double.MAX_VALUE;
192    
193            double maxX = Double.MIN_VALUE;
194            double maxY = Double.MIN_VALUE;
195    
196            for ( Point3d coord : originalPoints ) {
197                if ( minX > coord.x ) {
198                    minX = coord.x;
199                }
200                if ( maxX < coord.x ) {
201                    maxX = coord.x;
202                }
203                if ( minY > coord.y ) {
204                    minY = coord.y;
205                }
206                if ( maxY < coord.y ) {
207                    maxY = coord.y;
208                }
209                sourceCoords[count++] = (float) coord.x;
210                sourceCoords[count++] = (float) coord.y;
211            }
212    
213            float sX = Math.abs( ( maxX - minX ) ) > EPS11 ? (float) ( 1. / ( maxX - minX ) ) : 1;
214            float sY = Math.abs( ( maxY - minY ) ) > EPS11 ? (float) ( 1. / ( maxY - minY ) ) : 1;
215            float[] targetCoords = new float[projectedPoints.size() * 2];
216    
217            count = 0;
218            for ( Point3d coord : projectedPoints ) {
219                targetCoords[count++] = (float) coord.x;
220                targetCoords[count++] = (float) coord.y;
221            }
222    
223            StringBuilder sb = new StringBuilder( "\nCalculated scales are:\n" );
224            sb.append( "<crs:scaleX>" ).append( sX ).append( "</crs:scaleX>\n" );
225            sb.append( "<crs:scaleY>" ).append( sY ).append( "</crs:scaleY>\n" );
226            LOG.logInfo( sb.toString() );
227    
228            /**
229             * create warp object from reference points and desired interpolation, Because jai only implements the
230             * mapDestPoint function, we must calculate the polynomial variables by using the target coordinates as the
231             * source.
232             */
233            WarpPolynomial warp = WarpPolynomial.createWarp( targetCoords, 0, sourceCoords, 0, sourceCoords.length, sX, sY,
234                                                             1f / sX, 1f / sY, polynomalOrder );
235            return warp.getCoeffs();
236        }
237    
238        @Override
239        public int getOrder() {
240            return order;
241        }
242    
243        /**
244         * @return the scale which will be applied to the x value of all incoming data
245         */
246        public final float getScaleX() {
247            return scaleX;
248        }
249    
250        /**
251         * @return the scale which will be applied to the y value of all incoming data
252         */
253        public final float getScaleY() {
254            return scaleY;
255        }
256    
257    }