001: /*
002: * This file is part of the GeOxygene project source files.
003: *
004: * GeOxygene aims at providing an open framework which implements OGC/ISO specifications for
005: * the development and deployment of geographic (GIS) applications. It is a open source
006: * contribution of the COGIT laboratory at the Institut Géographique National (the French
007: * National Mapping Agency).
008: *
009: * See: http://oxygene-project.sourceforge.net
010: *
011: * Copyright (C) 2005 Institut Géographique National
012: *
013: * This library is free software; you can redistribute it and/or modify it under the terms
014: * of the GNU Lesser General Public License as published by the Free Software Foundation;
015: * either version 2.1 of the License, or any later version.
016: *
017: * This library is distributed in the hope that it will be useful, but WITHOUT ANY
018: * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
019: * PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
020: *
021: * You should have received a copy of the GNU Lesser General Public License along with
022: * this library (see file LICENSE if present); if not, write to the Free Software
023: * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
024: *
025: */
026:
027: package fr.ign.cogit.geoxygene.contrib.geometrie;
028:
029: import java.util.Iterator;
030:
031: import fr.ign.cogit.geoxygene.spatial.coordgeom.DirectPosition;
032: import fr.ign.cogit.geoxygene.spatial.coordgeom.DirectPositionList;
033: import fr.ign.cogit.geoxygene.spatial.coordgeom.GM_LineString;
034: import fr.ign.cogit.geoxygene.spatial.coordgeom.GM_Polygon;
035: import fr.ign.cogit.geoxygene.spatial.geomaggr.GM_MultiSurface;
036: import fr.ign.cogit.geoxygene.spatial.geomprim.GM_Point;
037: import fr.ign.cogit.geoxygene.spatial.geomprim.GM_Surface;
038: import fr.ign.cogit.geoxygene.spatial.geomroot.GM_Object;
039:
040: /** Méthodes statiques de calcul de distance.
041: *
042: * English: Computation of distances (static methods)
043: *
044: * @author Mustière/Bonin
045: * @version 1.0
046: */
047:
048: public abstract class Distances {
049:
050: // Organisation du code:
051: // - Distances entre points
052: // - Distances entre un point et un autre type de géoémtrie
053: // - Distances entre lignes
054: // - Distances entre surfaces
055:
056: //////////////////////////////////////////////////////////////
057: // //
058: // Distances entre points //
059: // //
060: //////////////////////////////////////////////////////////////
061:
062: /** Distance euclidienne entre 2 points (en 2D ou 3D si les points ont un Z). */
063: public static double distance(DirectPosition dp1, DirectPosition dp2) {
064: if (!Double.isNaN(dp1.getZ()) && !Double.isNaN(dp2.getZ())) {
065: return Math.sqrt(Math.pow(dp1.getX() - dp2.getX(), 2)
066: + Math.pow(dp1.getY() - dp2.getY(), 2)
067: + Math.pow(dp1.getZ() - dp2.getZ(), 2));
068: } else {
069: return Math.sqrt(Math.pow(dp1.getX() - dp2.getX(), 2)
070: + Math.pow(dp1.getY() - dp2.getY(), 2));
071: }
072: }
073:
074: /** Distance euclidienne calculée en 2 dimensions XY, même sur des objets 3D. */
075: public static double distance2D(DirectPosition dp1,
076: DirectPosition dp2) {
077: return Math.sqrt(Math.pow(dp1.getX() - dp2.getX(), 2)
078: + Math.pow(dp1.getY() - dp2.getY(), 2));
079: }
080:
081: /** Est-ce que les deux points sont distants de moins du seuil passé en paramètre ?
082: * Méthode optimisée pour accélérer les requêtes spatiales. */
083: public static boolean proche(DirectPosition dp1,
084: DirectPosition dp2, double distance) {
085: if (Math.abs(dp1.getX() - dp2.getX()) > distance)
086: return false;
087: if (Math.abs(dp1.getY() - dp2.getY()) > distance)
088: return false;
089: if (distance(dp1, dp2) > distance)
090: return false;
091: return true;
092: }
093:
094: //////////////////////////////////////////////////////////////
095: // //
096: // Distances entre un point et un autre type de géométrie //
097: // //
098: //////////////////////////////////////////////////////////////
099:
100: /** Distance euclidienne du point M au segment [A,B] */
101: public static double distancePointSegment(DirectPosition M,
102: DirectPosition A, DirectPosition B) {
103: return distance(M, Operateurs.projection(M, A, B));
104: }
105:
106: /** Distance euclidienne d'un point P à une ligne. */
107: public static double distance(DirectPosition M, GM_LineString L) {
108: DirectPositionList listePoints = L.coord();
109: double distmin, dist;
110:
111: distmin = distance(listePoints.get(0), M);
112: for (int i = 0; i < listePoints.size() - 1; i++) {
113: dist = distancePointSegment(M, listePoints.get(i),
114: listePoints.get(i + 1));
115: if (dist < distmin)
116: distmin = dist;
117: }
118: ;
119: return (distmin);
120: }
121:
122: ///////////////////////////////////////////////////////////
123: // //
124: // Distances entre lignes //
125: // //
126: ///////////////////////////////////////////////////////////
127:
128: /** Approximation de la première composante de Hausdorff d'une ligne vers une autre.
129: * Elle est calculee comme le maximum des distances des points intermédiaires
130: * de la première ligne L1 à l'autre ligne L2.
131: */
132: public static double premiereComposanteHausdorff(GM_LineString L1,
133: GM_LineString L2) {
134: DirectPositionList listePoints = L1.coord();
135: double dist, distmax = 0;
136:
137: for (int i = 0; i < listePoints.size(); i++) {
138: dist = distance(listePoints.get(i), L2);
139: if (dist > distmax)
140: distmax = dist;
141: }
142: ;
143: return distmax;
144: }
145:
146: /** Approximation (très proche) de la distance de Hausdorff entre deux lignes.
147: * Elle est calculee comme le maximum des distances d'un point intermediaire
148: * d'une des lignes a l'autre ligne. Dans certains cas cette definition
149: * diffère de la définition theorique pure car la distance de Hausdorff ne se
150: * realise pas necessairement sur un point intermediaire. Mais cela est rare
151: * sur des données réelles. Cette implementation est un bon compromis entre
152: * simplicité et précision.
153: */
154: public static double hausdorff(GM_LineString L1, GM_LineString L2) {
155: return Math.max(premiereComposanteHausdorff(L1, L2),
156: premiereComposanteHausdorff(L2, L1));
157: }
158:
159: /** Distance de Hausdorff entre un point P et une ligne L.
160: * C'est-à-dire distance au point P du point intermédiaire de
161: * la ligne L le plus éloigné du point P.
162: */
163: public static double hausdorff(GM_LineString L, GM_Point P) {
164: Iterator itPts = L.coord().getList().iterator();
165: DirectPosition point;
166: double distmax = 0, dist;
167:
168: while (itPts.hasNext()) {
169: point = (DirectPosition) itPts.next();
170: dist = distance(point, P.getPosition());
171: if (dist > distmax)
172: distmax = dist;
173: }
174: ;
175: return distmax;
176: }
177:
178: /** Distance moyenne entre deux polylignes,
179: * définie comme le rapport de l'aire séparant deux polylignes
180: * sur la moyenne de leurs longueurs.
181: *
182: * IMPORTANT: la méthode suppose que les lignes sont orientées globalement
183: * dans le même sens.
184: */
185: public static double distanceMoyenne(GM_LineString L1,
186: GM_LineString L2) {
187: GM_Polygon poly;
188: GM_LineString perimetre;
189: Iterator itPts;
190:
191: //fabrication de la surface delimitée par les lignes
192: perimetre = new GM_LineString();
193: itPts = L1.coord().getList().iterator();
194: while (itPts.hasNext()) {
195: DirectPosition pt = (DirectPosition) itPts.next();
196: perimetre.addControlPoint(0, pt);
197: }
198: itPts = L2.coord().getList().iterator();
199: while (itPts.hasNext()) {
200: DirectPosition pt = (DirectPosition) itPts.next();
201: perimetre.addControlPoint(0, pt);
202: }
203: perimetre.addControlPoint(L1.endPoint());
204: poly = new GM_Polygon(perimetre);
205:
206: return 2 * poly.area() / (L1.length() + L2.length());
207: }
208:
209: /** Mesure d'écart entre deux polylignes, défini comme une approximation de la
210: * surface séparant les polylignes.
211: * Plus précisément, cet écart est égal à la somme, pour chaque point P de L1,
212: * de (distance de P à L2) * (moyenne des longueurs des segments autour de P)
213: *
214: * NB: Ce n'est pas une distance au sens mathématique du terme,
215: * et en particulier cet écart n'est pas symétrique: ecart(L1,L2) != ecart(L2,L1)
216: */
217: public static double ecartSurface(GM_LineString L1, GM_LineString L2) {
218: double ecartTotal = 0, distPt, long1, long2;
219: DirectPositionList pts = L1.coord();
220: for (int i = 0; i < pts.size(); i++) {
221: distPt = distance(pts.get(i), L2);
222: if (i == 0)
223: long1 = 0;
224: else
225: long1 = distance(pts.get(i), pts.get(i - 1));
226: if (i == pts.size() - 1)
227: long2 = 0;
228: else
229: long2 = distance(pts.get(i), pts.get(i + 1));
230: ecartTotal = ecartTotal + distPt * (long1 + long2) / 2;
231: }
232: return ecartTotal;
233:
234: }
235:
236: ////////////////////////////////////////////////////////////
237: // //
238: // Distances entre surfaces //
239: // //
240: ////////////////////////////////////////////////////////////
241:
242: /** Distance surfacique entre deux GM_Polygon.
243: *
244: * Définition : 1 - surface(intersection)/surface(union)
245: * Ref [Vauglin 97]
246: *
247: * NB: renvoie 2 en cas de problème lors du calcul d'intersection avec JTS
248: * (bug en particulier si les surfaces sont dégénérées ou trop complexes).
249: */
250: public static double distanceSurfacique(GM_Polygon A, GM_Polygon B) {
251: GM_Object inter = A.intersection(B);
252: if (inter == null)
253: return 2;
254: GM_Object union = A.union(B);
255: if (union == null)
256: return 1;
257: return 1 - inter.area() / union.area();
258: }
259:
260: /** Distance surfacique entre deux GM_MultiSurface.
261: *
262: * Définition : 1 - surface(intersection)/surface(union)
263: * Ref [Vauglin 97]
264: *
265: * NB: renvoie 2 en cas de problème lors du calcul d'intersection avec JTS
266: * (bug en particulier si les surfaces sont dégénérées ou trop complexes).
267: */
268: public static double distanceSurfacique(GM_MultiSurface A,
269: GM_MultiSurface B) {
270: GM_Object inter = A.intersection(B);
271: //en cas de problème d'intersection avec JTS, la méthode retourne 2
272: if (inter == null)
273: return 2;
274: GM_Object union = A.union(B);
275: if (union == null)
276: return 1;
277: return 1 - inter.area() / union.area();
278: }
279:
280: /** Distance surfacique "robuste" entre deux polygones.
281: *
282: * Il s'agit ici d'une pure bidouille pour contourner certains bugs de JTS:
283: * Si JTS plante au calcul d'intersection, on filtre les surfaces avec Douglas et Peucker,
284: * progressivement avec 10 seuils entre min et max. Min et Max doivent être fixer donc de
285: * l'ordre de grandeur de la précision des données sinon le calcul risque d'être trop faussé.
286: *
287: * Définition : 1 - surface(intersection)/surface(union)
288: * Ref [Vauglin 97]
289: *
290: * NB: renvoie 2 en cas de problème lors du calcul d'intersection avec JTS
291: * (bug en particulier si les surfaces sont dégénérées ou trop complexes).
292: * */
293: public static double distanceSurfaciqueRobuste(GM_Polygon A,
294: GM_Polygon B, double min, double max) {
295: GM_Object inter = Operateurs
296: .intersectionRobuste(A, B, min, max);
297: //en cas de problème d'intersection avec JTS, la méthode retourne 2
298: if (inter == null)
299: return 2;
300: GM_Object union = A.union(B);
301: if (union == null)
302: return 1;
303: return 1 - inter.area() / union.area();
304: }
305:
306: /** Distance surfacique entre deux GM_MultiSurface.
307: *
308: * Cette méthode contourne des bugs de JTS, qui sont trop nombreux sur les agrégats.
309: * En contrepartie, cette méthode n'est valable que si les GM_Polygon composant A [resp. B]
310: * ne s'intersectent pas entre elles.
311: *
312: * Définition : 1 - surface(intersection)/surface(union)
313: * Ref [Vauglin 97]
314: *
315: * NB: renvoie 2 en cas de problème résiduer lors du calcul d'intersection avec JTS
316: * (bug en particulier si les surfaces sont dégénérées ou trop complexes).
317: */
318: public static double distanceSurfaciqueRobuste(GM_MultiSurface A,
319: GM_MultiSurface B) {
320: double inter = surfaceIntersection(A, B);
321: if (inter == -1) {
322: System.out
323: .println("Plantage JTS, renvoi 2 à la distance surfacique de deux multi_surfaces");
324: return 2;
325: }
326: return 1 - inter / (A.area() + B.area() - inter);
327: }
328:
329: /** Surface de l'intersection.
330: *
331: * Cette méthode contourne des bugs de JTS, qui sont trop nombreux sur les agrégats.
332: * En contrepartie, cette méthode n'est valable que si les GM_Polygon composant A [resp. B]
333: * ne s'intersectent pas entre elles.
334: *
335: * NB: renvoie -1 en cas de problème résiduer lors du calcul d'intersection avec JTS
336: * (bug en particulier si les surfaces sont dégénérées ou trop complexes).
337: */
338: public static double surfaceIntersection(GM_MultiSurface A,
339: GM_MultiSurface B) {
340: Iterator itA = A.getList().iterator();
341: Iterator itB;
342: double inter = 0;
343:
344: while (itA.hasNext()) {
345: GM_Surface surfA = (GM_Surface) itA.next();
346: itB = B.getList().iterator();
347: while (itB.hasNext()) {
348: GM_Surface surfB = (GM_Surface) itB.next();
349: if (surfB.intersection(surfA) == null) {
350: System.out
351: .println("Plantage JTS, renvoi -1 à l'intersection de deux multi_surfaces");
352: return -1;
353: }
354: inter = inter + surfB.intersection(surfA).area();
355: }
356: }
357: return inter;
358: }
359:
360: /** Surface de l'union.
361: *
362: * Cette méthode contourne des bugs de JTS, qui sont trop nombreux sur les agrégats.
363: * En contrepartie, cette méthode n'est valable que si les GM_Polygon composant A [resp. B]
364: * ne s'intersectent pas entre elles.
365: *
366: * NB: renvoie -1 en cas de problème résiduer lors du calcul d'intersection avec JTS
367: * (bug en particulier si les surfaces sont dégénérées ou trop complexes).
368: */
369: public static double surfaceUnion(GM_MultiSurface A,
370: GM_MultiSurface B) {
371: double inter = surfaceIntersection(A, B);
372: if (inter == -1) {
373: System.out
374: .println("Plantage JTS, renvoi -1 à l'union de deux 2 multi_surfaces");
375: return -1;
376: }
377: return A.area() + B.area() - inter;
378: }
379:
380: /** Mesure dite "Exactitude" entre 2 surfaces.
381: * Ref : [Bel Hadj Ali 2001]
382: *
383: * Définition : Surface(A inter B) / Surface(A)
384: */
385: public static double exactitude(GM_Polygon A, GM_Polygon B) {
386: GM_Object inter = A.intersection(B);
387: if (inter == null)
388: return 0;
389: return inter.area() / A.area();
390: }
391:
392: /** Mesure dite "Complétude" entre 2 surfaces.
393: * Ref : [Bel Hadj Ali 2001]
394: *
395: * Définition : Surface(A inter B) / Surface(B)
396: */
397: public static double completude(GM_Polygon A, GM_Polygon B) {
398: return exactitude(B, A);
399: }
400:
401: /** Mesure dite "Exactitude" entre 2 GM_MultiSurface.
402: * Ref : [Bel Hadj Ali 2001]
403: * Définition : Surface(A inter B) / Surface(A)
404: */
405: public static double exactitude(GM_MultiSurface A, GM_MultiSurface B) {
406: GM_Object inter = A.intersection(B);
407: if (inter == null)
408: return 0;
409: return inter.area() / A.area();
410: }
411:
412: /** Mesure dite "Complétude" entre 2 GM_MultiSurface.
413: *
414: * Ref : [Bel Hadj Ali 2001]
415: * Définition : Surface(A inter B) / Surface(B)
416: */
417: public static double completude(GM_MultiSurface A, GM_MultiSurface B) {
418: return exactitude(B, A);
419: }
420:
421: /** Mesure d'association entre deux surfaces (cf. [Bel Hadj Ali 2001]).
422: * <BR> <STRONG> Definition : </STRONG> associationSurfaces(A,B) = vrai si
423: * <UL>
424: * <LI> Surface(intersection) > min (min etant la resolution minimum des deux bases) </LI>
425: * <LI> ET (Surface(intersection) > surface(A) * coeff </LI>
426: * <LI> OU Surface(intersection) > surface(B) * coeff ) </LI>
427: * </UL>
428: * <BR> associationSurfaces(A,B) = faux sinon.
429: *
430: */
431: public static boolean associationSurfaces(GM_Object A, GM_Object B,
432: double min, double coeff) {
433: GM_Object inter = A.intersection(B);
434: if (inter == null)
435: return false;
436: double interArea = inter.area();
437: if (interArea < min)
438: return false;
439: if (interArea > A.area() * coeff)
440: return true;
441: if (interArea > B.area() * coeff)
442: return true;
443: return false;
444: }
445:
446: /** Test d'association "robuste" entre deux surfaces (cf. [Bel Hadj Ali 2001]).
447: *
448: * Il s'agit ici d'une pure bidouille pour contourner certains bugs de JTS:
449: * Si JTS plante au calcul , on filtre les surfaces avec Douglas et Peucker,
450: * progressivement avec 10 seuils entre min et max. Min et Max doivent être fixer donc de
451: * l'ordre de grandeur de la précision des données sinon le calcul risque d'être trop faussé.
452: *
453: * <BR> <STRONG> Definition : </STRONG> associationSurfaces(A,B) = vrai si
454: * <UL>
455: * <LI> Surface(intersection) > min (min etant la resolution minimum des deux bases) </LI>
456: * <LI> ET (Surface(intersection) > surface(A) * coeff </LI>
457: * <LI> OU Surface(intersection) > surface(B) * coeff ) </LI>
458: * </UL>
459: * <BR> associationSurfaces(A,B) = faux sinon.
460: *
461: */
462: public static boolean associationSurfacesRobuste(GM_Object A,
463: GM_Object B, double min, double coeff, double minDouglas,
464: double maxDouglas) {
465: GM_Object inter = Operateurs.intersectionRobuste(A, B,
466: minDouglas, maxDouglas);
467: if (inter == null)
468: return false;
469: double interArea = inter.area();
470: if (interArea < min)
471: return false;
472: if (interArea > A.area() * coeff)
473: return true;
474: if (interArea > B.area() * coeff)
475: return true;
476: return false;
477: }
478:
479: }
|