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