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.appariement.reseaux;
028:
029: import java.util.ArrayList;
030: import java.util.Iterator;
031: import java.util.List;
032:
033: import fr.ign.cogit.geoxygene.contrib.cartetopo.Arc;
034: import fr.ign.cogit.geoxygene.contrib.cartetopo.CarteTopo;
035: import fr.ign.cogit.geoxygene.contrib.geometrie.Distances;
036: import fr.ign.cogit.geoxygene.contrib.geometrie.Operateurs;
037: import fr.ign.cogit.geoxygene.contrib.geometrie.Vecteur;
038: import fr.ign.cogit.geoxygene.feature.FT_Feature;
039: import fr.ign.cogit.geoxygene.feature.FT_FeatureCollection;
040: import fr.ign.cogit.geoxygene.spatial.coordgeom.DirectPosition;
041: import fr.ign.cogit.geoxygene.spatial.coordgeom.DirectPositionList;
042: import fr.ign.cogit.geoxygene.spatial.coordgeom.GM_LineString;
043: import fr.ign.cogit.geoxygene.util.index.Tiling;
044:
045: /**
046: * Classe supportant les méthodes de comparaison globale de réseaux.
047: *
048: * @author Mustiere - IGN / Laboratoire COGIT
049: * @version 1.0
050: *
051: */
052:
053: public class Comparaison {
054:
055: /** Approximation de l'écart moyen entre deux réseaux.
056: * Cette approximation grossière évalue pour tout point de réseau1 son écart à réseau2
057: * (point le plus proche quelconque), sans réaliser d'appariement de réseau.
058: *
059: * @param reseau1
060: * Un réseau, typiquement un réseau peu précis dont on veut estimer la qualité par rapport à reseau2.
061: *
062: * @param reseau2
063: * Un autre réseau, typiquement un réseau de bonne qualité qui sert de référence.
064: *
065: * @param distanceMax
066: * Sert à éliminer les aberrations dans le calcul de la moyenne:
067: * les arcs de reseau1 situés au moins en un point à plus de distanceMax de reseau2
068: * ne sont pas pris en compte dans le calcul.
069: *
070: * @return
071: * L'écart moyen. Il est calculé comme la moyenne des distances entre les points
072: * des arcs de reseau1 et reseau2. Cette moyenne est pondérée par la longueur
073: * des segments entourant les points en question.
074: */
075: public static double approximationEcartPlaniMoyen(
076: CarteTopo reseau1, CarteTopo reseau2, double distanceMax) {
077:
078: double dtot = 0, ltot = 0, dmax = 0;
079: double ltotArc, dtotArc, dmaxArc, l1, l2, l;
080: int n = 0;
081: Iterator itArcsRef, itArcsComp;
082: FT_FeatureCollection arcsCompProches;
083: FT_FeatureCollection arcsRef = reseau1.getPopArcs();
084: FT_FeatureCollection arcsComp = reseau2.getPopArcs();
085:
086: itArcsRef = arcsRef.getElements().iterator();
087: while (itArcsRef.hasNext()) { //pour chaque arc de this
088: FT_Feature objetRef = (FT_Feature) itArcsRef.next();
089: GM_LineString geomRef = (GM_LineString) objetRef.getGeom();
090: DirectPositionList dp = geomRef.coord();
091: arcsCompProches = arcsComp.select(geomRef, distanceMax);
092:
093: dtotArc = 0;
094: ltotArc = 0;
095: dmaxArc = 0;
096: for (int i = 0; i < dp.size(); i++) { //pour chaque point d'un arc de this
097: DirectPosition ptRef = dp.get(i);
098: if (i == 0)
099: l1 = 0;
100: else
101: l1 = Distances.distance(dp.get(i), dp.get(i - 1));
102: if (i == dp.size() - 1)
103: l2 = 0;
104: else
105: l2 = Distances.distance(dp.get(i), dp.get(i + 1));
106: l = l1 + l2;
107: double dmin = Double.MAX_VALUE;
108: itArcsComp = arcsCompProches.getElements().iterator();
109: while (itArcsComp.hasNext()) {
110: FT_Feature objetComp = (FT_Feature) itArcsComp
111: .next();
112: GM_LineString geomComp = (GM_LineString) objetComp
113: .getGeom();
114: DirectPosition ptProjete = Operateurs.projection(
115: ptRef, geomComp);
116: double d = Distances.distance(ptRef, ptProjete);
117: if (d < dmin)
118: dmin = d;
119: }
120: if (dmin > dmaxArc)
121: dmaxArc = dmin;
122: dtotArc = dtotArc + dmin * l;
123: ltotArc = ltotArc + l;
124: }
125: if (dmaxArc > distanceMax)
126: continue;
127: if (dmaxArc > dmax)
128: dmax = dmaxArc;
129: dtot = dtot + dtotArc;
130: ltot = ltot + ltotArc;
131: n++;
132: }
133:
134: System.out.println("distance moyenne entre les arcs" + dtot
135: / ltot);
136: System.out.println("distance max entre les arcs " + dmax);
137: System.out
138: .println("Nb d'arcs de this pris en compte dans le calcul "
139: + n);
140: System.out
141: .println("Nb d'arcs de this non pris en compte dans le calcul "
142: + (arcsRef.size() - n));
143: return dtot / ltot;
144: }
145:
146: /** Ensemble d'indicateurs évaluant globalement l'écart de position entre le réseau à étudier réseau11
147: * et un réseau de référence réseau2.
148: * Cette évaluation se fait en s'appuyant sur le calcul des écarts entre chaque point
149: * de reseau1 et sa projection au plus près sur le réseau reseau2 (et non sur un réel appariement d'arcs).
150: * Les moyennes sont pondérées par la longueur des segments entourant les points,
151: * pour gommer les effets dus aux pas de découpage variables des lignes.
152: *
153: * @param reseau1
154: * Réseau étudié.
155: *
156: * @param reseau2
157: * Réseau servant de référence.
158: *
159: * @param affichage
160: * Si TRUE alors les résultats sont affichés.
161: *
162: * @param distanceMax
163: * Sert à éliminer les aberrations dans les calculs.
164: * - Les arcs de reseau1 situés en au moins un point à plus de distanceMax
165: * de reseau2 ne sont pas pris en compte dans le calcul des indicateurs sur les arcs.
166: * - Les noeuds de reseau1 situés à plus de distanceMax d'un noeud
167: * de reseau2 ne sont pas pris en compte dans le calcul des indicateurs sur les noeuds
168: *
169: * @return
170: * Liste (de 'double') contenant un ensemble d'indicateurs sur l'écart entre les réseaux :
171: *
172: * ESTIMATEURS SUR LES ARCS
173: * liste(0): longueur des arcs du réseau "this" total
174: * liste(1): longueur des arcs du réseau "this" pris en compte dans les calculs d'évaluation de l'écart
175: * liste(2): longueur des arcs du réseau "reseau"
176: * liste(3): nombre d'arcs du réseau "this" total
177: * liste(4): nombre d'arcs du réseau "this" pris en compte dans les calculs d'évaluation de l'écart
178: * liste(5): nombre d'arcs du réseau "reseau"
179: * liste(6): estimation du biais systématique en X sur les arcs
180: * (valeur en X de la moyenne des vecteurs d'écart entre un point de "this" et son projeté sur "reseau")
181: * liste(7): estimation du biais systématique en Y sur les arcs
182: * (valeur en Y de la moyenne des vecteurs d'écart entre un point de "this" et son projeté sur "reseau")
183: * liste(8): estimation de l'écart moyen sur les arcs
184: * (moyenne des longueurs des vecteurs d'écart entre un point de "this" et son projeté sur "reseau")
185: * liste(9): estimation de l'écart moyen quadratique sur les arcs
186: * (moyenne quadratique des longueurs des vecteurs d'écart entre un point de "this" et son projeté sur "reseau")
187: * liste(10): estimation de l'écart type sur les arcs, i.e. précision une fois le biais corrigé
188: * ( racine(ecart moyen quadratique^2 - biais^2)
189: * liste(11): histogramme de répartition des écarts sur tous les points (en nb de points intermédiaires sur les arcs)
190: *
191: * ESTIMATEURS SUR LES NOEUDS (si ils existent)
192: * liste(12): nombre de noeuds du réseau "this" total
193: * liste(13): nombre de noeuds du réseau "this" pris en compte dans les calculs d'évaluation de l'écart
194: * liste(14): nombre de noeuds du réseau "reseau"
195: * liste(15): estimation du biais systématique en X sur les noeuds
196: * (valeur en X de la moyenne des vecteurs d'écart entre un noeud de "this" et le noeud le plus proche de "reseau")
197: * liste(16): estimation du biais systématique en Y sur les noeuds
198: * (valeur en Y de la moyenne des vecteurs d'écart entre un noeud de "this" et le noeud le plus proche de "reseau")
199: * liste(17): estimation de l'écart moyen sur les noeuds
200: * (moyenne des longueurs des vecteurs d'écart entre un noeud de "this" et le noeud le plus proche de "reseau")
201: * liste(18): estimation de l'écart moyen quadratique sur les arcs
202: * (moyenne quadratique des longueurs des vecteurs d'écart entre un noeud de "this" et le noeud le plus proche de "reseau")
203: * liste(19): estimation de l'écart type sur les noeuds, i.e. précision une fois le biais corrigé
204: * ( racine(ecart moyen quadratique^2 - biais^2)
205: * liste(20): histogramme de répartition des écarts sur tous les noeuds (en nb de noeuds)
206: *
207: */
208: public static List evaluationEcartPosition(CarteTopo reseau1,
209: CarteTopo reseau2, double distanceMax, boolean affichage) {
210:
211: List resultats = new ArrayList();
212: FT_FeatureCollection arcs1 = reseau1.getPopArcs();
213: FT_FeatureCollection arcs2 = reseau2.getPopArcs();
214: Iterator itArcs1, itArcs2;
215: Arc arc1, arc2;
216: GM_LineString geom1, geom2;
217: Vecteur v12, vmin, vPourUnArc1, vTotal;
218: DirectPosition pt1, projete;
219: FT_FeatureCollection arcs2proches;
220: double longArc1, poids, d12, ecartQuadratiquePourUnArc1, l1, l2, poidsPourUnArc1, ecartPourUnArc1, ecartMaxArc1, poidsTotal;
221:
222: // indicateurs finaux
223: double longTotal1, longPrisEnCompte1, longTotal2;
224: double ecartTotal, ecartQuadratiqueTotal, ecartTypeArcs;
225: int nbArcsTotal1, nbArcsPrisEnCompte1, nbArcsTotal2;
226:
227: ///////////////// EVALUATION SUR LES ARCS ////////////////////
228: // indexation des arcs du réseau 2
229: if (!reseau2.getPopArcs().hasSpatialIndex()) {
230: reseau2.getPopArcs().initSpatialIndex(Tiling.class, false,
231: 20);
232: }
233: // parcours du réseau 2 juste pour calculer sa longueur
234: itArcs2 = arcs2.getElements().iterator();
235: longTotal2 = 0;
236: while (itArcs2.hasNext()) {
237: arc2 = (Arc) itArcs2.next();
238: longTotal2 = longTotal2
239: + ((GM_LineString) arc2.getGeom()).length();
240: }
241: nbArcsTotal2 = arcs2.getElements().size();
242: // parcours du réseau 1 pour évaluer les écarts sur les arcs
243: nbArcsTotal1 = arcs1.getElements().size();
244: nbArcsPrisEnCompte1 = 0;
245: poidsTotal = 0;
246: ecartTotal = 0;
247: ecartQuadratiqueTotal = 0;
248: vTotal = new Vecteur(new DirectPosition(0, 0));
249: longTotal1 = 0;
250: longPrisEnCompte1 = 0;
251: itArcs1 = arcs1.getElements().iterator();
252: while (itArcs1.hasNext()) { //pour chaque arc du réseau 1
253: arc1 = (Arc) itArcs1.next();
254: geom1 = (GM_LineString) arc1.getGeom();
255: longArc1 = geom1.length();
256: longTotal1 = longTotal1 + geom1.length();
257:
258: arcs2proches = arcs2.select(geom1, distanceMax);
259: if (arcs2proches.size() == 0)
260: continue;
261: poidsPourUnArc1 = 0;
262: ecartPourUnArc1 = 0;
263: ecartQuadratiquePourUnArc1 = 0;
264: vmin = new Vecteur();
265: vPourUnArc1 = new Vecteur(new DirectPosition(0, 0));
266: ecartMaxArc1 = 0;
267: DirectPositionList dp = geom1.coord();
268: for (int i = 0; i < dp.size(); i++) { //pour chaque point de arc1
269: pt1 = dp.get(i);
270: // calcul du poids de ce point
271: if (i == 0)
272: l1 = 0;
273: else
274: l1 = Distances.distance(dp.get(i), dp.get(i - 1));
275: if (i == dp.size() - 1)
276: l2 = 0;
277: else
278: l2 = Distances.distance(dp.get(i), dp.get(i + 1));
279: poids = l1 + l2;
280: // projection du point sur le réseau2
281: double dmin = Double.MAX_VALUE;
282: itArcs2 = arcs2proches.getElements().iterator();
283: while (itArcs2.hasNext()) { // pour chaque arc du réseau 2
284: arc2 = (Arc) itArcs2.next();
285: geom2 = (GM_LineString) arc2.getGeom();
286: projete = Operateurs.projection(pt1, geom2);
287: v12 = new Vecteur(pt1, projete);
288: d12 = Distances.distance(pt1, projete);
289: if (d12 < dmin) {
290: dmin = d12;
291: vmin = v12;
292: }
293: }
294: // calcul des indicateurs
295: if (dmin > ecartMaxArc1)
296: ecartMaxArc1 = dmin;
297: ecartPourUnArc1 = ecartPourUnArc1 + dmin * poids;
298: ecartQuadratiquePourUnArc1 = ecartPourUnArc1 + poids
299: * Math.pow(dmin, 2);
300: vPourUnArc1 = vPourUnArc1.ajoute(vmin
301: .multConstante(poids));
302: poidsPourUnArc1 = poidsPourUnArc1 + poids;
303: }
304: if (ecartMaxArc1 > distanceMax)
305: continue; // on ne prend pas l'arc en compte
306: longPrisEnCompte1 = longPrisEnCompte1 + longArc1;
307: nbArcsPrisEnCompte1++;
308: poidsTotal = poidsTotal + poidsPourUnArc1;
309: ecartTotal = ecartTotal + ecartPourUnArc1;
310: ecartQuadratiqueTotal = ecartQuadratiqueTotal
311: + ecartQuadratiquePourUnArc1;
312: vTotal = vTotal.ajoute(vPourUnArc1);
313: }
314: vTotal = vTotal.multConstante(1 / poidsTotal);
315: ecartTotal = ecartTotal / poidsTotal;
316: ecartTypeArcs = Math.sqrt((ecartQuadratiqueTotal / poidsTotal)
317: - (Math.pow(vTotal.getX(), 2) + Math.pow(vTotal.getY(),
318: 2)));
319: ecartQuadratiqueTotal = Math.sqrt(ecartQuadratiqueTotal
320: / poidsTotal);
321: // resultats sur les arcs
322: if (affichage)
323: System.out
324: .println("******************* BILAN SUR LES ARCS *******************");
325: if (affichage)
326: System.out
327: .println("** Comparaison globale des réseaux, en nombre et longueur d'arcs :");
328: resultats.add(new Double(longTotal1));
329: if (affichage)
330: System.out
331: .println("** Longueur total des arcs du réseau 1 (km) : "
332: + Math.round(longTotal1 / 1000));
333: resultats.add(new Double(longPrisEnCompte1));
334: if (affichage)
335: System.out
336: .println("** Longueur total des arcs du réseau 1 pris en compte pour le calcul (km) : "
337: + Math.round(longPrisEnCompte1 / 1000));
338: resultats.add(new Double(longTotal2));
339: if (affichage)
340: System.out
341: .println("** Longueur total des arcs du réseau 2(km) : "
342: + Math.round(longTotal2 / 1000));
343: resultats.add(new Double(nbArcsTotal1));
344: if (affichage)
345: System.out.println("** Nombre d'arcs du réseau 1 : "
346: + nbArcsTotal1);
347: resultats.add(new Double(nbArcsPrisEnCompte1));
348: if (affichage)
349: System.out
350: .println("** Nombre d'arcs du réseau 1 pris en compte pour le calcul : "
351: + nbArcsPrisEnCompte1);
352: resultats.add(new Double(nbArcsTotal2));
353: if (affichage)
354: System.out.println("** Nombre d'arcs du réseau 2 : "
355: + nbArcsTotal2);
356:
357: if (affichage)
358: System.out.println("");
359: if (affichage)
360: System.out.println("** Estimateurs d'écart ");
361: resultats.add(new Double(vTotal.getX()));
362: if (affichage)
363: System.out.println("** Biais systématique en X : "
364: + vTotal.getX());
365: resultats.add(new Double(vTotal.getY()));
366: if (affichage)
367: System.out.println("** Biais systématique en Y : "
368: + vTotal.getY());
369: resultats.add(new Double(ecartTotal));
370: if (affichage)
371: System.out.println("** Ecart moyen : " + ecartTotal);
372: resultats.add(new Double(ecartQuadratiqueTotal));
373: if (affichage)
374: System.out.println("** Ecart moyen quadratique : "
375: + ecartQuadratiqueTotal);
376: resultats.add(new Double(ecartTypeArcs));
377: if (affichage)
378: System.out.println("** Ecart type : " + ecartTypeArcs);
379: if (affichage)
380: System.out
381: .println("*********************************************************");
382:
383: return resultats;
384: }
385:
386: }
|