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