001 //$HeadURL$
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.model.spatialschema;
038
039 import javax.vecmath.Point3d;
040 import javax.vecmath.Vector2d;
041 import javax.vecmath.Vector3d;
042
043 import org.deegree.framework.log.ILogger;
044 import org.deegree.framework.log.LoggerFactory;
045
046 /**
047 * Utility class for the linearization of arcs and circles.
048 *
049 * @author <a href="mailto:elmasry@lat-lon.de">Moataz Elmasry</a>
050 * @author last edited by: $Author: elmasri$
051 *
052 * @version $Revision: $, $Date: 9 May 2008 13:09:29$
053 */
054 public class LinearizationUtil {
055
056 private static ILogger LOG = LoggerFactory.getLogger( LinearizationUtil.class );
057
058 private static double PI2 = Math.PI * 2;
059
060 private static double EPSILON = 1E-11;
061
062 /**
063 * Takes in a three Position of an arc and computes a linearized Arc (linestring) using 150 positions/vertices
064 *
065 * @param p0
066 * @param p1
067 * @param p2
068 * @return An array of linearized Positions
069 */
070 public static Position[] linearizeArc( Position p0, Position p1, Position p2 ) {
071 return linearizeArc( p0, p1, p2, 150 );
072 }
073
074 /**
075 * Takes in a three Positions of an Arc and computes a linearized Arc (linestring).
076 *
077 * @param p0
078 * @param p1
079 * @param p2
080 *
081 * @param numberOfOpPositions
082 * the number of linearized Positions to produce
083 * @return An array of linearized Positions
084 */
085 public static Position[] linearizeArc( Position p0, Position p1, Position p2, int numberOfOpPositions ) {
086
087 // shift the points down (to reduce the occurrence of floating point errors), independently on the x and y axes
088 double shiftx = findShiftX( p0, p1, p2 );
089 double shifty = findShiftY( p0, p1, p2 );
090 // if the points are already shifted, this does no harm!
091 Position p0Shifted = new PositionImpl( p0.getX() - shiftx, p0.getY() - shifty );
092 Position p1Shifted = new PositionImpl( p1.getX() - shiftx, p1.getY() - shifty );
093 Position p2Shifted = new PositionImpl( p2.getX() - shiftx, p2.getY() - shifty );
094
095 if ( areCollinear( p0Shifted, p1Shifted, p2Shifted ) ) {
096 Position[] Positions = new Position[3];
097 Positions[0] = p0;
098 Positions[1] = p1;
099 Positions[2] = p2;
100 return Positions;
101 }
102
103 int j = numberOfOpPositions - 1;
104 if ( j < 1 ) {
105 return null;
106 }
107
108 double[] positionsDouble = computeArc( p0Shifted.getX(), p0Shifted.getY(), p1Shifted.getX(), p1Shifted.getY(),
109 p2Shifted.getX(), p2Shifted.getY() );
110 double d11 = computeArcOrientation( p0Shifted.getX(), p0Shifted.getY(), p1Shifted.getX(), p1Shifted.getY(),
111 p2Shifted.getX(), p2Shifted.getY() );
112 double d12 = positionsDouble[3];
113 double d13 = positionsDouble[5];
114 if ( d11 > 0 && d13 < d12 ) {
115 d13 += PI2;
116 } else if ( d11 < 0 && d13 > d12 ) {
117 d13 -= PI2;
118 }
119 Position outPositions[] = new Position[numberOfOpPositions];
120 for ( int k = 0; k <= j; k++ ) {
121 double d15 = d12 + ( ( d13 - d12 ) * k ) / j;
122 double xCoord = positionsDouble[0] + positionsDouble[2] * Math.cos( d15 );
123 double yCoord = positionsDouble[1] + positionsDouble[2] * Math.sin( d15 );
124 outPositions[k] = GeometryFactory.createPosition( xCoord, yCoord );
125 }
126
127 // ensure numerical stability for start and end points
128 outPositions[0] = GeometryFactory.createPosition( p0Shifted.getX(), p0Shifted.getY() );
129 outPositions[outPositions.length - 1] = GeometryFactory.createPosition( p2Shifted.getX(), p2Shifted.getY() );
130
131 // shift the points back up
132 for ( int i = 0; i < outPositions.length; i++ ) {
133 outPositions[i] = new PositionImpl( outPositions[i].getX() + shiftx, outPositions[i].getY() + shifty );
134 }
135
136 return outPositions;
137 }
138
139 private static double findShiftY( Position p0, Position p1, Position p2 ) {
140 double miny = p0.getY();
141 if ( p1.getY() < miny ) {
142 miny = p1.getY();
143 }
144 if ( p2.getY() < miny ) {
145 miny = p2.getY();
146 }
147
148 double maxy = p0.getY();
149 if ( p1.getX() > maxy ) {
150 maxy = p1.getY();
151 }
152 if ( p2.getY() > maxy ) {
153 maxy = p2.getY();
154 }
155 return ( miny + maxy ) / 2;
156 }
157
158 private static double findShiftX( Position p0, Position p1, Position p2 ) {
159 double minx = p0.getX();
160 if ( p1.getX() < minx ) {
161 minx = p1.getX();
162 }
163 if ( p2.getX() < minx ) {
164 minx = p2.getX();
165 }
166
167 double maxx = p0.getX();
168 if ( p1.getX() > maxx ) {
169 maxx = p1.getX();
170 }
171 if ( p2.getX() > maxx ) {
172 maxx = p2.getX();
173 }
174 return ( minx + maxx ) / 2;
175 }
176
177 /**
178 * Calculates a sequence of {@link Position}s on the arc of a circle.
179 * <p>
180 * The circle is constructed from the input points: All three points belong to the arc of the circle. They must be
181 * distinct and non-colinear. To form a complete circle, the arc is extended past the third control point until the
182 * first control point is encountered.
183 *
184 * @param p0
185 * start point on the arc of the circle
186 * @param p1
187 * second point on the arc of the circle
188 * @param p2
189 * third point on the arc of the circle
190 * @param numberOfPositions
191 * number of interpolation points, the returned array contains numberOfPositions+1 entries, to ensure
192 * that the first point is identical to the last
193 * @return Position array with numberOfPositions+1 entries, first entry and last entry are identical to p0
194 * @throws IllegalArgumentException
195 * if no order can be determined, because the points are identical or co-linear
196 */
197 public static Position[] linearizeCircle( Position p0, Position p1, Position p2, int numberOfPositions ) {
198
199 // shift the points down (to reduce the occurrence of floating point errors), independently on the x and y axes
200 double shiftx = findShiftX( p0, p1, p2 );
201 double shifty = findShiftY( p0, p1, p2 );
202 // if the points are already shifted, this does no harm!
203 Position p0Shifted = new PositionImpl( p0.getX() - shiftx, p0.getY() - shifty );
204 Position p1Shifted = new PositionImpl( p1.getX() - shiftx, p1.getY() - shifty );
205 Position p2Shifted = new PositionImpl( p2.getX() - shiftx, p2.getY() - shifty );
206
207 Position[] positions = new Position[numberOfPositions + 1];
208
209 Position center = findCircleCenter( p0Shifted, p1Shifted, p2Shifted );
210
211 double centerX = center.getX();
212 double centerY = center.getY();
213 double dx = p0Shifted.getX() - centerX;
214 double dy = p0Shifted.getY() - centerY;
215 double radius = Math.sqrt( dx * dx + dy * dy );
216
217 double angleStep = 0.0;
218 if ( isClockwise( p0Shifted, p1Shifted, p2Shifted ) ) {
219 angleStep = -Math.PI * 2 / numberOfPositions;
220 } else {
221 angleStep = Math.PI * 2 / numberOfPositions;
222 }
223 double angle = Math.atan2( dy, dx );
224
225 for ( int i = 0; i < numberOfPositions; i++ ) {
226 double x = centerX + Math.cos( angle ) * radius;
227 double y = centerY + Math.sin( angle ) * radius;
228 positions[i] = GeometryFactory.createPosition( x, y );
229 angle += angleStep;
230 }
231 positions[numberOfPositions] = positions[0];
232
233 // shift the points back up
234 for ( int i = 0; i <= numberOfPositions; i++ ) {
235 positions[i] = new PositionImpl( positions[i].getX() + shiftx, positions[i].getY() + shifty );
236 }
237
238 return positions;
239 }
240
241 /**
242 * Returns whether the order of the given three points is clockwise or counterclockwise. Uses the (signed) area of a
243 * planar triangle to get to know about the order of the points.
244 *
245 * @param p0
246 * first point
247 * @param p1
248 * second point
249 * @param p2
250 * third point
251 * @return true, if order is clockwise, otherwise false
252 * @throws IllegalArgumentException
253 * if no order can be determined, because the points are identical or co-linear
254 */
255 static final boolean isClockwise( Position p0, Position p1, Position p2 )
256 throws IllegalArgumentException {
257
258 if ( areCollinear( p0, p1, p2 ) ) {
259 throw new IllegalArgumentException( "Cannot evaluate isClockwise(). The three points are colinear." );
260 }
261 double res = ( p2.getX() - p0.getX() ) * ( ( p2.getY() + p0.getY() ) / 2 ) + ( p1.getX() - p2.getX() )
262 * ( ( p1.getY() + p2.getY() ) / 2 ) + ( p0.getX() - p1.getX() ) * ( ( p0.getY() + p1.getY() ) / 2 );
263 return res < 0.0 ? true : false;
264 }
265
266 /**
267 * Tests if the given three points are collinear.
268 * <p>
269 * NOTE: Only this method should be used throughout the whole linearization process for testing collinearity to
270 * avoid inconsistent results (the necessary EPSILON would differ).
271 *
272 * @param p0
273 * @param p1
274 * @param p2
275 * @return true if the points are collinear, false otherwise
276 */
277 public static boolean areCollinear( Position p0, Position p1, Position p2 ) {
278 double res = ( p2.getX() - p0.getX() ) * ( ( p2.getY() + p0.getY() ) / 2 ) + ( p1.getX() - p2.getX() )
279 * ( ( p1.getY() + p2.getY() ) / 2 ) + ( p0.getX() - p1.getX() ) * ( ( p0.getY() + p1.getY() ) / 2 );
280 return Math.abs( res ) < EPSILON;
281 }
282
283 /**
284 * Computes an Arc from of three points (wich must not be collinear).
285 *
286 * @return array of points representing the computed arc
287 */
288 private static final double[] computeArc( double d1, double d2, double d3, double d4, double d5, double d6 ) {
289
290 double d13 = d1 - d3;
291 double d14 = d3 - d5;
292 double d15 = d4 - d6;
293 double d16 = d2 - d4;
294 double d17 = d1 + d3;
295 double d18 = d3 + d5;
296 double d19 = d2 + d4;
297 double d20 = d4 + d6;
298 double d21 = d13 * d15 - d14 * d16;
299
300 double d7 = ( ( d15 * d17 * d13 - d16 * d18 * d14 ) + d16 * d15 * ( d16 + d15 ) ) / d21;
301 double d8;
302 if ( Math.abs( d4 - d6 ) < 5E-008D ) {
303 d8 = ( d19 * d16 + ( d17 - d7 ) * d13 ) / d16;
304 } else {
305 d8 = ( d20 * d15 + ( d18 - d7 ) * d14 ) / d15;
306 }
307 d7 *= 0.5;
308 d8 *= 0.5;
309 double d9 = Math.sqrt( ( d7 - d1 ) * ( d7 - d1 ) + ( d8 - d2 ) * ( d8 - d2 ) );
310 double d10 = Math.atan2( d2 - d8, d1 - d7 );
311 if ( d10 < 0 ) {
312 d10 += PI2;
313 }
314 double d11 = Math.atan2( d4 - d8, d3 - d7 );
315 if ( d11 < 0 ) {
316 d11 += PI2;
317 }
318 double d12 = Math.atan2( d6 - d8, d5 - d7 );
319 if ( d12 < 0 ) {
320 d12 += PI2;
321 }
322
323 double ad[] = new double[6];
324 ad[0] = d7;
325 ad[1] = d8;
326 ad[2] = d9;
327 ad[3] = d10;
328 ad[4] = d11;
329 ad[5] = d12;
330 return ad;
331
332 }
333
334 private static final double computeArcOrientation( double d, double d1, double d2, double d3, double d4, double d5 ) {
335 return ( d * d3 + d2 * d5 + d4 * d1 ) - ( d4 * d3 + d2 * d1 + d * d5 );
336 }
337
338 /**
339 * Transforms an array of Position to an array of double
340 *
341 * @param positions
342 * @return array of double
343 */
344 protected static double[] positionToDouble( Position[] positions ) {
345 double[] ad = new double[positions.length * 2];
346 for ( int i = 0; i < positions.length; i++ ) {
347 ad[i] = positions[i].getX();
348 ad[i + 1] = positions[i].getY();
349 }
350 return ad;
351 }
352
353 /**
354 * Transforms an array of double to an array of Position
355 *
356 * @param doubles
357 * @return generated array of Position or null if doubles%2 != 0
358 */
359 protected static Position[] doubleToPosition( double[] doubles ) {
360 // The number of doubles has to be even
361 if ( doubles.length % 2 != 0 ) {
362 return null;
363 }
364 Position[] positions = new Position[doubles.length / 2];
365 for ( int i = 0; i < positions.length; i++ ) {
366 positions[i] = GeometryFactory.createPosition( doubles[i], doubles[i + 1] );
367 }
368 return positions;
369 }
370
371 /**
372 * Finds the center of a cirlce/arc using three points that lie on the circle. Thanks to wikipedia:
373 * http://en.wikipedia.org/wiki/Circumradius#Coordinates_of_circumcenter.
374 *
375 * @param p0
376 * @param p1
377 * @param p2
378 * @return center of the circle
379 * @throws IllegalArgumentException
380 * if the points are co linear, e.g. on a line.
381 */
382 static Position findCircleCenter( Position p0, Position p1, Position p2 )
383 throws IllegalArgumentException {
384
385 // shift the points down (to reduce the occurrence of floating point errors), independently on the x and y axes
386 double shiftx = findShiftX( p0, p1, p2 );
387 double shifty = findShiftY( p0, p1, p2 );
388 // if the points are already shifted, this does no harm!
389 Position p0Shifted = new PositionImpl( p0.getX() - shiftx, p0.getY() - shifty );
390 Position p1Shifted = new PositionImpl( p1.getX() - shiftx, p1.getY() - shifty );
391 Position p2Shifted = new PositionImpl( p2.getX() - shiftx, p2.getY() - shifty );
392
393 if ( areCollinear( p0Shifted, p1Shifted, p2Shifted ) ) {
394 throw new IllegalArgumentException( "The given points are co linear, no circum center can be calculated." );
395 }
396
397 Vector3d a = new Vector3d( p0Shifted.getAsPoint3d() );
398 Vector3d b = new Vector3d( p1Shifted.getAsPoint3d() );
399 Vector3d c = new Vector3d( p2Shifted.getAsPoint3d() );
400
401 if ( Double.isNaN( a.z ) ) {
402 a.z = 0;
403 }
404 if ( Double.isNaN( b.z ) ) {
405 b.z = 0;
406 }
407 if ( Double.isNaN( c.z ) ) {
408 c.z = 0;
409 }
410
411 Vector3d ab = new Vector3d( a );
412 Vector3d ac = new Vector3d( a );
413 Vector3d bc = new Vector3d( b );
414 Vector3d ba = new Vector3d( b );
415 Vector3d ca = new Vector3d( c );
416 Vector3d cb = new Vector3d( c );
417
418 ab.sub( b );
419 ac.sub( c );
420 bc.sub( c );
421 ba.sub( a );
422 ca.sub( a );
423 cb.sub( b );
424
425 Vector3d cros = new Vector3d();
426 cros.cross( ab, bc );
427 double crosSquare = 2 * cros.length() * cros.length();
428
429 if ( LOG.isDebug() ) {
430 double radius = ( ab.length() * bc.length() * ca.length() ) / ( 2 * cros.length() );
431 LOG.logDebug( "Radius: " + radius );
432 }
433
434 a.scale( ( ( bc.length() * bc.length() ) * ab.dot( ac ) ) / crosSquare );
435 b.scale( ( ( ac.length() * ac.length() ) * ba.dot( bc ) ) / crosSquare );
436 c.scale( ( ( ab.length() * ab.length() ) * ca.dot( cb ) ) / crosSquare );
437
438 Point3d circle = new Point3d( a );
439 circle.add( b );
440 circle.add( c );
441
442 // shift the center circle back up
443 circle.x += shiftx;
444 circle.y += shifty;
445
446 return GeometryFactory.createPosition( circle );
447 }
448
449 /**
450 * Finds the angle between three points, p0,p1 and a center, where the angle to be measured. In other words, its the
451 * angle between the segments p0-center and p1-center.
452 *
453 * @param p0
454 * @param p1
455 * @param centre
456 * @return calculated angle in degree
457 */
458 private static double findAngle( Position p0, Position p1, Position centre ) {
459 Vector2d vec1 = calcVector( centre, p0 );
460 Vector2d vec2 = calcVector( centre, p1 );
461 double cosTheta = ( vec1.dot( vec2 ) ) / ( vec1.length() * vec2.length() );
462 int round = 1000;
463 // Using the 1000th decimal to round the cosTheta.
464 cosTheta = Math.ceil( cosTheta * round ) / round;
465 // convert from radian to degree
466 return Math.acos( cosTheta ) / Math.PI * 180;
467
468 }
469
470 /**
471 * Finds the vector between two points. which is the difference between p1 and p0 and the direction is from p0 to
472 * p1.
473 *
474 * @param p0
475 * @param p1
476 * @return calculated vector
477 */
478 protected static Vector2d calcVector( Position p0, Position p1 ) {
479 return new Vector2d( p1.getX() - p0.getX(), p1.getY() - p0.getY() );
480 }
481 }