001 //$HeadURL: svn+ssh://rbezema@svn.wald.intevation.org/deegree/base/tags/2.1/src/org/deegree/model/csct/ct/AbridgedMolodenskiTransform.java $
002 /*---------------- FILE HEADER ------------------------------------------
003
004 This file is part of deegree.
005 Copyright (C) 2001 by:
006 EXSE, Department of Geography, University of Bonn
007 http://www.giub.uni-bonn.de/exse/
008 lat/lon GmbH
009 http://www.lat-lon.de
010
011 It has been implemented within SEAGIS - An OpenSource implementation of OpenGIS specification
012 (C) 2001, Institut de Recherche pour le D�veloppement (http://sourceforge.net/projects/seagis/)
013 SEAGIS Contacts: Surveillance de l'Environnement Assist�e par Satellite
014 Institut de Recherche pour le D�veloppement / US-Espace
015 mailto:seasnet@teledetection.fr
016
017
018 This library is free software; you can redistribute it and/or
019 modify it under the terms of the GNU Lesser General Public
020 License as published by the Free Software Foundation; either
021 version 2.1 of the License, or (at your option) any later version.
022
023 This library is distributed in the hope that it will be useful,
024 but WITHOUT ANY WARRANTY; without even the implied warranty of
025 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
026 Lesser General Public License for more details.
027
028 You should have received a copy of the GNU Lesser General Public
029 License along with this library; if not, write to the Free Software
030 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
031
032 Contact:
033
034 Andreas Poth
035 lat/lon GmbH
036 Aennchenstr. 19
037 53115 Bonn
038 Germany
039 E-Mail: poth@lat-lon.de
040
041 Klaus Greve
042 Department of Geography
043 University of Bonn
044 Meckenheimer Allee 166
045 53115 Bonn
046 Germany
047 E-Mail: klaus.greve@uni-bonn.de
048
049
050 ---------------------------------------------------------------------------*/
051 package org.deegree.model.csct.ct;
052
053 // OpenGIS dependencies (SEAGIS)
054 import java.io.Serializable;
055
056 import javax.media.jai.ParameterList;
057
058 import org.deegree.model.csct.cs.Ellipsoid;
059 import org.deegree.model.csct.cs.HorizontalDatum;
060 import org.deegree.model.csct.cs.WGS84ConversionInfo;
061 import org.deegree.model.csct.resources.css.ResourceKeys;
062
063
064 /**
065 * Transforms a three dimensional geographic points using
066 * abridged versions of formulas derived by Molodenski.
067 *
068 * @version 1.00
069 * @author OpenGIS (www.opengis.org)
070 * @author Martin Desruisseaux
071 */
072 class AbridgedMolodenskiTransform extends AbstractMathTransform implements Serializable
073 {
074 /**
075 * Serial number for interoperability with different versions.
076 */
077 //private static final long serialVersionUID = ?;
078
079 /**
080 * X,Y,Z shift in meters
081 */
082 private final double dx, dy, dz;
083
084 /**
085 * Source equatorial radius in meters.
086 */
087 private final double a;
088
089 /**
090 * Source polar radius in meters
091 */
092 private final double b;
093
094 /**
095 * Source flattening factor.
096 */
097 private final double f;
098
099 /**
100 * Difference in the semi-major axes (a1 - a2)
101 * of the first and second ellipsoids.
102 */
103 private final double da;
104
105 /**
106 * Difference in the flattening of the two ellipsoids.
107 */
108 private final double df;
109
110 /**
111 * Square of the eccentricity of the ellipsoid.
112 */
113 private final double e2;
114
115 /**
116 * Defined as <code>(a*df) + (f*da)</code>.
117 */
118 private final double adf;
119
120 /**
121 * Construct a transform.
122 */
123 protected AbridgedMolodenskiTransform(final HorizontalDatum source, final HorizontalDatum target)
124 {
125 final WGS84ConversionInfo srcInfo = source.getWGS84Parameters();
126 final WGS84ConversionInfo tgtInfo = source.getWGS84Parameters();
127 final Ellipsoid srcEllipsoid = source.getEllipsoid();
128 final Ellipsoid tgtEllipsoid = target.getEllipsoid();
129 dx = srcInfo.dx - tgtInfo.dx;
130 dy = srcInfo.dy - tgtInfo.dy;
131 dz = srcInfo.dz - tgtInfo.dz;
132 a = srcEllipsoid.getSemiMajorAxis();
133 b = srcEllipsoid.getSemiMinorAxis();
134 f = 1 / srcEllipsoid.getInverseFlattening();
135 da = a - tgtEllipsoid.getSemiMajorAxis();
136 df = f - 1/tgtEllipsoid.getInverseFlattening();
137 e2 = 1 - (b*b)/(a*a);
138 adf = (a*df) + (f*da);
139 }
140
141 /**
142 * Transforms a list of coordinate point ordinal values.
143 */
144 public void transform(final double[] srcPts, int srcOff, final double[] dstPts, int dstOff, int numPts)
145 {
146 int step = 0;
147 if (srcPts==dstPts && srcOff<dstOff && srcOff+numPts*getDimSource()>dstOff)
148 {
149 step = -getDimSource();
150 srcOff -= (numPts-1)*step;
151 dstOff -= (numPts-1)*step;
152 }
153 final double rho = Double.NaN; // TODO: Definition???
154 while (--numPts >= 0)
155 {
156 double x = Math.toRadians(srcPts[srcOff++]);
157 double y = Math.toRadians(srcPts[srcOff++]);
158 double z = srcPts[srcOff++];
159 final double sinX = Math.sin(x);
160 final double cosX = Math.cos(x);
161 final double sinY = Math.sin(y);
162 final double cosY = Math.cos(y);
163 final double sin2Y = sinY*sinY;
164 final double nu = a / Math.sqrt(1 - e2*sin2Y);
165
166 // Note: Computation of 'x' and 'y' ommit the division by sin(1"), because
167 // 1/sin(1") / (60*60*180/PI) = 1.0000000000039174050898603898692...
168 // (60*60 is for converting the final result from seconds to degrees,
169 // and 180/PI is for converting degrees to radians). This is an error
170 // of about 8E-7 arc seconds, probably close to rounding errors anyway.
171 y += (dz*cosY - sinY*(dy*sinX + dx*cosX) + adf*Math.sin(2*y)) / rho;
172 x += (dy*cosX - dx*sinX) / (nu*cosY);
173 z += dx*cosY*cosX + dy*cosY*sinX + dz*sinY + adf*sin2Y - da;
174
175 dstPts[dstOff++] = Math.toDegrees(x);
176 dstPts[dstOff++] = Math.toDegrees(y);
177 dstPts[dstOff++] = z;
178 srcOff += step;
179 dstOff += step;
180 }
181 }
182
183 /**
184 * Transforms a list of coordinate point ordinal values.
185 */
186 public void transform(final float[] srcPts, final int srcOff, final float[] dstPts, final int dstOff, int numPts)
187 {
188 // TODO: Copy the implementation from 'transform(double[]...)'.
189 try
190 {
191 super.transform(srcPts, srcOff, dstPts, dstOff, numPts);
192 }
193 catch (TransformException exception)
194 {
195 // Should not happen.
196 }
197 }
198
199 /**
200 * Gets the dimension of input points, which is 3.
201 */
202 public int getDimSource()
203 {return 3;}
204
205 /**
206 * Gets the dimension of output points, which
207 * is the same than {@link #getDimSource()}.
208 */
209 public final int getDimTarget()
210 {return getDimSource();}
211
212 /**
213 * Tests whether this transform does not move any points.
214 * This method returns always <code>false</code>.
215 */
216 public final boolean isIdentity()
217 {return false;}
218
219 /**
220 * Returns a hash value for this transform.
221 */
222 public final int hashCode()
223 {
224 final long code = Double.doubleToLongBits(dx) +
225 37*(Double.doubleToLongBits(dy) +
226 37*(Double.doubleToLongBits(dz) +
227 37*(Double.doubleToLongBits(a ) +
228 37*(Double.doubleToLongBits(b ) +
229 37*(Double.doubleToLongBits(da) +
230 37*(Double.doubleToLongBits(df)))))));
231 return (int) code ^ (int) (code >>> 32);
232 }
233
234 /**
235 * Compares the specified object with
236 * this math transform for equality.
237 */
238 public final boolean equals(final Object object)
239 {
240 if (object==this) return true; // Slight optimization
241 if (super.equals(object))
242 {
243 final AbridgedMolodenskiTransform that = (AbridgedMolodenskiTransform) object;
244 return Double.doubleToLongBits(this.dx) == Double.doubleToLongBits(that.dx) &&
245 Double.doubleToLongBits(this.dy) == Double.doubleToLongBits(that.dy) &&
246 Double.doubleToLongBits(this.dz) == Double.doubleToLongBits(that.dz) &&
247 Double.doubleToLongBits(this.a ) == Double.doubleToLongBits(that.a ) &&
248 Double.doubleToLongBits(this.b ) == Double.doubleToLongBits(that.b ) &&
249 Double.doubleToLongBits(this.da) == Double.doubleToLongBits(that.da) &&
250 Double.doubleToLongBits(this.df) == Double.doubleToLongBits(that.df);
251 }
252 return false;
253 }
254
255 /**
256 * Returns the WKT for this math transform.
257 */
258 public final String toString()
259 {
260 final StringBuffer buffer = paramMT("Abridged_Molodenski");
261 addParameter(buffer, "dim", getDimSource());
262 addParameter(buffer, "dx", dx);
263 addParameter(buffer, "dy", dy);
264 addParameter(buffer, "src_semi_major", a);
265 addParameter(buffer, "src_semi_minor", b);
266 addParameter(buffer, "tgt_semi_major", a-da);
267 // addParameter(buffer, "tgt_semi_minor", b); // TODO
268 buffer.append(']');
269 return buffer.toString();
270 }
271
272 /**
273 * The provider for {@link AbridgedMolodenskiTransform}.
274 *
275 * @version 1.0
276 * @author Martin Desruisseaux
277 */
278 static final class Provider extends MathTransformProvider
279 {
280 /**
281 * Create a provider.
282 */
283 public Provider()
284 {
285 super("Abridged_Molodenski", ResourceKeys.ABRIDGED_MOLODENSKI_TRANSFORM, null);
286 put("dim", Double.NaN, POSITIVE_RANGE);
287 put("dx", Double.NaN, null);
288 put("dy", Double.NaN, null);
289 put("src_semi_major", Double.NaN, POSITIVE_RANGE);
290 put("src_semi_minor", Double.NaN, POSITIVE_RANGE);
291 put("tgt_semi_major", Double.NaN, POSITIVE_RANGE);
292 put("tgt_semi_minor", Double.NaN, POSITIVE_RANGE);
293 }
294
295 /**
296 * Returns a transform for the specified parameters.
297 *
298 * @param parameters The parameter values in standard units.
299 * @return A {@link MathTransform} object of this classification.
300 */
301 public MathTransform create(final ParameterList parameters)
302 {return null;} // TODO
303 }
304 }