0001: /*
0002: * This file is part of the GeOxygene project source files.
0003: *
0004: * GeOxygene aims at providing an open framework which implements OGC/ISO specifications for
0005: * the development and deployment of geographic (GIS) applications. It is a open source
0006: * contribution of the COGIT laboratory at the Institut Géographique National (the French
0007: * National Mapping Agency).
0008: *
0009: * See: http://oxygene-project.sourceforge.net
0010: *
0011: * Copyright (C) 2005 Institut Géographique National
0012: *
0013: * This library is free software; you can redistribute it and/or modify it under the terms
0014: * of the GNU Lesser General Public License as published by the Free Software Foundation;
0015: * either version 2.1 of the License, or any later version.
0016: *
0017: * This library is distributed in the hope that it will be useful, but WITHOUT ANY
0018: * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
0019: * PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
0020: *
0021: * You should have received a copy of the GNU Lesser General Public License along with
0022: * this library (see file LICENSE if present); if not, write to the Free Software
0023: * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
0024: *
0025: */
0026:
0027: package fr.ign.cogit.geoxygene.contrib.geometrie;
0028:
0029: import java.util.ArrayList;
0030: import java.util.Iterator;
0031: import java.util.List;
0032:
0033: import Jama.Matrix;
0034: import fr.ign.cogit.geoxygene.feature.FT_Feature;
0035: import fr.ign.cogit.geoxygene.feature.Population;
0036: import fr.ign.cogit.geoxygene.generalisation.Filtering;
0037: import fr.ign.cogit.geoxygene.spatial.coordgeom.DirectPosition;
0038: import fr.ign.cogit.geoxygene.spatial.coordgeom.DirectPositionList;
0039: import fr.ign.cogit.geoxygene.spatial.coordgeom.GM_LineString;
0040: import fr.ign.cogit.geoxygene.spatial.coordgeom.GM_Polygon;
0041: import fr.ign.cogit.geoxygene.spatial.geomaggr.GM_Aggregate;
0042: import fr.ign.cogit.geoxygene.spatial.geomprim.GM_Point;
0043: import fr.ign.cogit.geoxygene.spatial.geomroot.GM_Object;
0044: import fr.ign.cogit.geoxygene.util.index.Tiling;
0045:
0046: /**
0047: * En vrac un ensemble de méthodes statiques qui manipulent des géométries:
0048: * projections, abscisse curviligne, décalage, orientation...
0049: *
0050: * CONTIENT des méthodes de :
0051: * Projections d'un point
0052: * Manipulation de l'abscisse curviligne d'une ligne
0053: * Mesures sur un polygone
0054: * Offset d'une ligne (décalage)
0055: * Echantillonage d'une ligne
0056: * Regression linéaire
0057: * et beaucoup d'autres choses très diverses
0058: *
0059: * ATTENTION: certaines méthodes n'ont pas été conçues ni testées pour des coordonnées 3D
0060: *
0061: *
0062: * English: Very very diverse set of methods on geometries
0063: *
0064: * @author Mustière / Bonin / Rousseaux / Grosso / Lafragueta
0065: * @version 1.0
0066: */
0067:
0068: public abstract class Operateurs {
0069:
0070: //////////////////////////////////////////////////////////////////////
0071: // Projections d'un point
0072: //////////////////////////////////////////////////////////////////////
0073:
0074: /** Projection de M sur le segment [A,B]
0075: *
0076: * English: Projects M on a [A,B]
0077: * @author Mustière
0078: * */
0079: public static DirectPosition projection(DirectPosition M,
0080: DirectPosition A, DirectPosition B) {
0081: double lambda;
0082: Vecteur uAB, AM;
0083:
0084: if (Distances.distance(A, B) == 0)
0085: return A; // cas ou A et B sont confondus
0086:
0087: uAB = new Vecteur(A, B).vectNorme();
0088: AM = new Vecteur(A, M);
0089: lambda = AM.prodScalaire(uAB);
0090:
0091: if (lambda <= 0)
0092: return A; // Cas ou M se projete en A sur le segment [AB]
0093: if (lambda >= Distances.distance(A, B))
0094: return B; // Cas ou M se projete en B sur le segment [AB]
0095: return translate(A, uAB.multConstante(lambda)); // Cas ou M se projete entre A et B
0096: }
0097:
0098: /** Projection du point sur la polyligne.
0099: * En théorie, il peut y avoir plusieurs points projetés, mais dans ce cas
0100: * cette méthode n'en renvoie qu'un seul (le premier dans le sens de parcours
0101: * de la ligne).
0102: *
0103: * English: Projects M on the lineString
0104: * @author Mustière
0105: */
0106: public static DirectPosition projection(DirectPosition dp,
0107: GM_LineString LS) {
0108: DirectPositionList listePoints = LS.coord();
0109: double d, dmin;
0110: DirectPosition pt, ptmin;
0111:
0112: if (listePoints.size() == 0)
0113: return listePoints.get(0);
0114: ptmin = projection(dp, listePoints.get(0), listePoints.get(1));
0115: dmin = Distances.distance(dp, ptmin);
0116: for (int i = 1; i < listePoints.size() - 1; i++) {
0117: pt = projection(dp, listePoints.get(i), listePoints
0118: .get(i + 1));
0119: d = Distances.distance(dp, pt);
0120: if (d < dmin) {
0121: ptmin = pt;
0122: dmin = d;
0123: }
0124: }
0125: return ptmin;
0126: }
0127:
0128: /** Projection du point sur l'aggregat;
0129: * ATTENTION: ne fonctionne que si l'aggregat ne contient que des GM_Point et GM_LineString.
0130: *
0131: * En théorie, il peut y avoir plusieurs points projetés, mais dans ce cas
0132: * cette méthode n'en renvoie qu'un seul.
0133: *
0134: * English: Projects M on the agregate
0135: * @author Mustière
0136: * */
0137: public static DirectPosition projection(DirectPosition dp,
0138: GM_Aggregate aggr) {
0139: Iterator itComposants = aggr.getList().iterator();
0140: double d = 0, dmin = Double.MAX_VALUE;
0141: DirectPosition pt = null, ptmin = null;
0142: boolean geomOK;
0143:
0144: while (itComposants.hasNext()) {
0145: GM_Object composant = (GM_Object) itComposants.next();
0146: geomOK = false;
0147: if (composant instanceof GM_Point) {
0148: pt = ((GM_Point) composant).getPosition();
0149: d = Distances.distance(pt, dp);
0150: geomOK = true;
0151: }
0152: if (composant instanceof GM_LineString) {
0153: pt = projection(dp, (GM_LineString) composant);
0154: d = Distances.distance(pt, dp);
0155: geomOK = true;
0156: }
0157: if (!geomOK) {
0158: System.out
0159: .println("Projection - Type de géométrie non géré: "
0160: + composant.getClass());
0161: continue;
0162: }
0163: if (d < dmin) {
0164: ptmin = pt;
0165: dmin = d;
0166: }
0167: }
0168: return ptmin;
0169: }
0170:
0171: //////////////////////////////////////////////////////////////////////
0172: // Manipulation de l'abscisse curviligne d'une ligne
0173: //////////////////////////////////////////////////////////////////////
0174:
0175: /** Coordonnées du point situé sur la ligne à l'abscisse curviligne passée en paramètre.
0176: * Renvoie Null si l'abscisse est négative ou plus grande que la longueur de la ligne.
0177: *
0178: * English: Point located at the curvilinear abscisse
0179: * @author Mustière
0180: */
0181: public static DirectPosition pointEnAbscisseCurviligne(
0182: GM_LineString ls, double abscisse) {
0183: int i;
0184: double l = 0;
0185: double d;
0186: DirectPosition pt1, pt2;
0187: Vecteur v1;
0188:
0189: if (abscisse > ls.length() || abscisse < 0)
0190: return null;
0191: pt1 = ls.coord().get(0);
0192: for (i = 1; i < ls.coord().size(); i++) {
0193: pt2 = ls.coord().get(i);
0194: d = Distances.distance(pt1, pt2);
0195: if (d != 0) {
0196: if (l + d > abscisse) {
0197: v1 = new Vecteur(pt1, pt2);
0198: v1 = v1.multConstante((abscisse - l) / d);
0199: return translate(pt1, v1);
0200: }
0201: l = l + d;
0202: pt1 = pt2;
0203: }
0204: }
0205: return ls.coord().get(ls.coord().size() - 1);
0206: }
0207:
0208: /** Abscisse curviligne du ieme point de la ligne ls.
0209: * English: curvilinear abscisse of the ith point
0210: * @author Mustière
0211: */
0212: public static double abscisseCurviligne(GM_LineString ls, int i) {
0213: double abs = 0;
0214: for (int j = 0; j < i; j++) {
0215: abs = abs
0216: + Distances.distance(ls.getControlPoint(j), ls
0217: .getControlPoint(j + 1));
0218: }
0219: return abs;
0220: }
0221:
0222: /** Coordonnées du point situé sur au milieu de la ligne.
0223: * English: Point in the middle of the line
0224: * @author Mustière
0225: */
0226: public static DirectPosition milieu(GM_LineString ls) {
0227: return pointEnAbscisseCurviligne(ls, ls.length() / 2);
0228: }
0229:
0230: /** renvoie le milieu de [A,B].
0231: *
0232: * English: Point in the middle of [A,B]
0233: * @author Mustière
0234: */
0235: public static DirectPosition milieu(DirectPosition A,
0236: DirectPosition B) {
0237: DirectPosition M;
0238: if (!Double.isNaN(A.getZ()) && !Double.isNaN(B.getZ())) {
0239: M = new DirectPosition((A.getX() + B.getX()) / 2,
0240: (A.getY() + B.getY()) / 2,
0241: (A.getZ() + B.getZ()) / 2);
0242: } else {
0243: M = new DirectPosition((A.getX() + B.getX()) / 2,
0244: (A.getY() + B.getY()) / 2, Double.NaN);
0245: }
0246: return M;
0247: }
0248:
0249: /** Premiers points intermédiaires de la ligne ls, situés à moins
0250: * de la longueur curviligne passée en paramètre du point initial.
0251: * Renvoie null si la longueur est négative.
0252: * Renvoie le premier point si et seulement si la longueur est 0.
0253: * Renvoie tous les points si la longueur est supérieure à la longueur de la ligne
0254: * NB: les points sont renvoyés dans l'ordre en partant du premier point.
0255: *
0256: * English: First points of the line.
0257: * @author Mustière
0258: */
0259: public static DirectPositionList premiersPoints(GM_LineString ls,
0260: double longueur) {
0261: int i;
0262: double l = 0;
0263: DirectPositionList listePts = new DirectPositionList();
0264:
0265: if (longueur < 0)
0266: return null;
0267: listePts.add(ls.getControlPoint(0));
0268: for (i = 1; i < ls.coord().size(); i++) {
0269: l = l
0270: + Distances.distance(ls.getControlPoint(i - 1), ls
0271: .getControlPoint(i));
0272: if (l > longueur)
0273: break;
0274: listePts.add(ls.getControlPoint(i));
0275: }
0276: return listePts;
0277: }
0278:
0279: /** Derniers points intermédiaires de la ligne ls, situés à moins
0280: * de la longueur curviligne passée en paramètre du point final.
0281: * Renvoie null si la longueur est négative.
0282: * Renvoie le dernier point seulement si la longueur est 0.
0283: * Renvoie tous les points si la longueur est supérieure à la longueur de la ligne.
0284: * NB: les points sont renvoyés dans l'ordre en partant du dernier point
0285: * (ordre inverse par rapport à la géoémtrie initiale).
0286: *
0287: * English: Last points of the line.
0288: * @author Mustière
0289: */
0290: public static DirectPositionList derniersPoints(GM_LineString ls,
0291: double longueur) {
0292: int i;
0293: double l = 0;
0294: DirectPositionList listePts = new DirectPositionList();
0295: int nbPts = ls.coord().size();
0296:
0297: if (longueur < 0)
0298: return null;
0299: listePts.add(ls.getControlPoint(nbPts - 1));
0300: for (i = nbPts - 2; i >= 0; i--) {
0301: l = l
0302: + Distances.distance(ls.getControlPoint(i), ls
0303: .getControlPoint(i + 1));
0304: if (l > longueur)
0305: break;
0306: listePts.add(ls.getControlPoint(i));
0307: }
0308: return listePts;
0309: }
0310:
0311: //////////////////////////////////////////////////////////////////////
0312: // Mesures sur un polygone
0313: //////////////////////////////////////////////////////////////////////
0314: /** Barycentre 2D (approximatif).
0315: * Il est défini comme le barycentre des points intermédiaires du contour,
0316: * ce qui est très approximatif
0317: *
0318: * English: Center of the points of the polygon.
0319: * @author Mustière
0320: */
0321: public static DirectPosition barycentre2D(GM_Polygon poly) {
0322: DirectPositionList listePoints = poly.coord();
0323: DirectPosition barycentre;
0324: double moyenneX = 0;
0325: double moyenneY = 0;
0326: double sommeX = 0;
0327: double sommeY = 0;
0328:
0329: for (int i = 0; i < listePoints.size() - 1; i++) {
0330: sommeX = sommeX + listePoints.get(i).getX();
0331: sommeY = sommeY + listePoints.get(i).getY();
0332: }
0333: ;
0334: moyenneX = sommeX / (listePoints.size() - 1);
0335: moyenneY = sommeY / (listePoints.size() - 1);
0336:
0337: barycentre = new DirectPosition(moyenneX, moyenneY);
0338: return (barycentre);
0339: }
0340:
0341: //////////////////////////////////////////////////////////////////////
0342: // Offset d'une ligne (décalage)
0343: //////////////////////////////////////////////////////////////////////
0344: /** Calcul d'un offset direct (demi-buffer d'une ligne, ou décalage à gauche).
0345: * Le paramètre offset est la taille du décalage.
0346: *
0347: * English: shift of a line on the left
0348: *
0349: * @author Bonin, Rousseaux.
0350: */
0351: public static GM_LineString directOffset(GM_LineString ls,
0352: double offset) {
0353: DirectPositionList listePoints = ls.coord();
0354: DirectPosition point, pointPrec, pointSuiv, pointRes;
0355: Vecteur u, v, n;
0356: GM_LineString ligneResultat = new GM_LineString();
0357: u = new Vecteur(listePoints.get(0), listePoints.get(1));
0358: u.setZ(0.0);
0359: pointRes = new DirectPosition(listePoints.get(0).getX()
0360: + offset * u.vectNorme().getY(), listePoints.get(0)
0361: .getY()
0362: - offset * u.vectNorme().getX(), listePoints.get(0)
0363: .getZ());
0364: ligneResultat.addControlPoint(pointRes);
0365: for (int j = 1; j < listePoints.size() - 1; j++) {
0366: pointPrec = (DirectPosition) listePoints.get(j - 1);
0367: point = (DirectPosition) listePoints.get(j);
0368: pointSuiv = (DirectPosition) listePoints.get(j + 1);
0369: u = new Vecteur(pointPrec, point);
0370: u.setZ(0);
0371: v = new Vecteur(point, pointSuiv);
0372: v.setZ(0);
0373: n = u.vectNorme().soustrait(v.vectNorme());
0374: if (u.prodVectoriel(v).getZ() > 0) {
0375: pointRes = new DirectPosition(point.getX() + offset
0376: * n.vectNorme().getX(), point.getY() + offset
0377: * n.vectNorme().getY(), point.getZ());
0378: ligneResultat.addControlPoint(pointRes);
0379: } else {
0380: pointRes = new DirectPosition(point.getX() - offset
0381: * n.vectNorme().getX(), point.getY() - offset
0382: * n.vectNorme().getY(), point.getZ());
0383: ligneResultat.addControlPoint(pointRes);
0384: }
0385: }
0386: u = new Vecteur(listePoints.get(listePoints.size() - 2),
0387: listePoints.get(listePoints.size() - 1));
0388: u.setZ(0.0);
0389: pointRes = new DirectPosition(listePoints.get(
0390: listePoints.size() - 1).getX()
0391: + offset * u.vectNorme().getY(), listePoints.get(
0392: listePoints.size() - 1).getY()
0393: - offset * u.vectNorme().getX(), listePoints.get(
0394: listePoints.size() - 1).getZ());
0395: ligneResultat.addControlPoint(pointRes);
0396: return ligneResultat;
0397: }
0398:
0399: /** Calcul d'un offset indirect (demi-buffer d'une ligne, ou décalage à droite).
0400: * Le paramètre offset est la taille du décalage.
0401: *
0402: * English: shift of a line on the right
0403: *
0404: * @author Bonin, Rousseaux.
0405: */
0406: public static GM_LineString indirectOffset(GM_LineString ls,
0407: double offset) {
0408: DirectPositionList listePoints = ls.coord();
0409: DirectPosition point, pointPrec, pointSuiv, pointRes;
0410: Vecteur u, v, n;
0411: GM_LineString ligneResultat = new GM_LineString();
0412: u = new Vecteur(listePoints.get(0), listePoints.get(1));
0413: u.setZ(0);
0414: pointRes = new DirectPosition(listePoints.get(0).getX()
0415: - offset * u.vectNorme().getY(), listePoints.get(0)
0416: .getY()
0417: + offset * u.vectNorme().getX(), listePoints.get(0)
0418: .getZ());
0419: ligneResultat.addControlPoint(pointRes);
0420: for (int j = 1; j < listePoints.size() - 1; j++) {
0421: pointPrec = (DirectPosition) listePoints.get(j - 1);
0422: point = (DirectPosition) listePoints.get(j);
0423: pointSuiv = (DirectPosition) listePoints.get(j + 1);
0424: u = new Vecteur(pointPrec, point);
0425: u.setZ(0);
0426: v = new Vecteur(point, pointSuiv);
0427: v.setZ(0);
0428: n = u.vectNorme().soustrait(v.vectNorme());
0429: if (u.prodVectoriel(v).getZ() < 0) {
0430: pointRes = new DirectPosition(point.getX() + offset
0431: * n.vectNorme().getX(), point.getY() + offset
0432: * n.vectNorme().getY(), point.getZ());
0433: ligneResultat.addControlPoint(pointRes);
0434: } else {
0435: pointRes = new DirectPosition(point.getX() - offset
0436: * n.vectNorme().getX(), point.getY() - offset
0437: * n.vectNorme().getY(), point.getZ());
0438: ligneResultat.addControlPoint(pointRes);
0439: }
0440: }
0441: u = new Vecteur(listePoints.get(listePoints.size() - 2),
0442: listePoints.get(listePoints.size() - 1));
0443: u.setZ(0);
0444: pointRes = new DirectPosition(listePoints.get(
0445: listePoints.size() - 1).getX()
0446: - offset * u.vectNorme().getY(), listePoints.get(
0447: listePoints.size() - 1).getY()
0448: + offset * u.vectNorme().getX(), listePoints.get(
0449: listePoints.size() - 1).getZ());
0450: ligneResultat.addControlPoint(pointRes);
0451: return ligneResultat;
0452: }
0453:
0454: //////////////////////////////////////////////////////////////////////
0455: // Echantillonage
0456: //////////////////////////////////////////////////////////////////////
0457: /** Méthode pour suréchantillonner une GM_LineString.
0458: * Des points intermédiaires écartés du pas sont ajoutés sur chaque segment
0459: * de la ligne ls, à partir du premier point de chaque segment.
0460: * (voir aussi echantillonePasVariable pour une autre méthode )
0461: *
0462: * English: Resampling of a line
0463: *
0464: * @author Bonin, Rousseaux.
0465: */
0466: public static GM_LineString echantillone(GM_LineString ls,
0467: double pas) {
0468: DirectPosition point1, point2, Xins;
0469: Vecteur u1;
0470: int nseg, i;
0471: double longTronc;
0472: Double fseg;
0473: GM_LineString routeSurech = (GM_LineString) ls.clone();
0474: DirectPositionList listePoints = routeSurech.coord();
0475: DirectPositionList listePointsEchant = new DirectPositionList();
0476: for (int j = 1; j < listePoints.size(); j++) {
0477: point1 = (DirectPosition) listePoints.get(j - 1);
0478: listePointsEchant.add(point1);
0479: point2 = (DirectPosition) listePoints.get(j);
0480: longTronc = Distances.distance(point1, point2);
0481: fseg = new Double(longTronc / pas);
0482: nseg = fseg.intValue();
0483: u1 = new Vecteur(point1, point2);
0484: for (i = 0; i < nseg - 1; i++) {
0485: Xins = new DirectPosition(point1.getX() + (i + 1) * pas
0486: * u1.vectNorme().getX(), point1.getY()
0487: + (i + 1) * pas * u1.vectNorme().getY(), point1
0488: .getZ()
0489: + (i + 1) * pas * u1.vectNorme().getZ());
0490: listePointsEchant.add(Xins);
0491: }
0492: ;
0493: }
0494: listePointsEchant.add(listePoints.get(listePoints.size() - 1));
0495: routeSurech = new GM_LineString(listePointsEchant);
0496: return routeSurech;
0497: }
0498:
0499: /** Méthode pour suréchantillonner une GM_LineString.
0500: * A l'inverse de la méthode "echantillone", le pas d'echantillonage
0501: * diffère sur chaque segment de manière à ce que l'on échantillone chaque
0502: * segment en différents mini-segments tous de même longueur.
0503: * Le pas en entrée est le pas maximum autorisé.
0504: *
0505: * English : Resampling of a line
0506: *
0507: * @author Grosso.
0508: */
0509: public static GM_LineString echantillonePasVariable(
0510: GM_LineString ls, double pas) {
0511: DirectPosition point1, point2, Xins;
0512: Vecteur u1;
0513: int nseg, i;
0514: double longTronc;
0515: Double fseg;
0516: GM_LineString routeSurech = (GM_LineString) ls.clone();
0517: DirectPositionList listePoints = routeSurech.coord();
0518: DirectPositionList listePointsEchant = new DirectPositionList();
0519: for (int j = 1; j < listePoints.size(); j++) {
0520: point1 = (DirectPosition) listePoints.get(j - 1);
0521: listePointsEchant.add(point1);
0522: point2 = (DirectPosition) listePoints.get(j);
0523: longTronc = Distances.distance(point1, point2);
0524: fseg = new Double(longTronc / pas);
0525: nseg = fseg.intValue();
0526: double epsilonPas = 0;
0527: if (nseg != 0) {
0528: epsilonPas = longTronc / nseg - pas;
0529: }
0530: double nouveauPas = pas - epsilonPas;
0531: u1 = new Vecteur(point1, point2);
0532: for (i = 0; i < nseg - 1; i++) {
0533: Xins = new DirectPosition(point1.getX() + (i + 1)
0534: * nouveauPas * u1.vectNorme().getX(), point1
0535: .getY()
0536: + (i + 1) * nouveauPas * u1.vectNorme().getY(),
0537: point1.getZ() + (i + 1) * nouveauPas
0538: * u1.vectNorme().getZ());
0539: listePointsEchant.add(Xins);
0540: }
0541: ;
0542: }
0543: listePointsEchant.add(listePoints.get(listePoints.size() - 1));
0544: routeSurech = new GM_LineString(listePointsEchant);
0545: return routeSurech;
0546: }
0547:
0548: /** Renvoie le point translaté de P avec le vecteur V;
0549: * Contrairement au "move" de DirectPosition, on ne deplace pas le point P
0550: *
0551: * English : Shift of a point
0552: */
0553: public static DirectPosition translate(DirectPosition P, Vecteur V) {
0554: if (!Double.isNaN(P.getZ()) && !Double.isNaN(V.getZ())) {
0555: return new DirectPosition(P.getX() + V.getX(), P.getY()
0556: + V.getY(), P.getZ() + V.getZ());
0557: } else {
0558: return new DirectPosition(P.getX() + V.getX(), P.getY()
0559: + V.getY(), Double.NaN);
0560: }
0561: }
0562:
0563: //////////////////////////////////////////////////////////////////////
0564: // (Très) Divers
0565: //////////////////////////////////////////////////////////////////////
0566:
0567: /** Mise bout à bout de plusieurs GM_LineString pour constituer une nouvelle GM_LineString
0568: * La liste en entrée contient des GM_LineString.
0569: * La polyligne créée commence sur l'extrémité libre de la première
0570: * polyligne de la liste.
0571: *
0572: * English: Combination of lines
0573: * @author: Mustière
0574: */
0575: public static GM_LineString compileArcs(List geometries) {
0576: DirectPositionList pointsFinaux = new DirectPositionList();
0577: DirectPosition pointCourant;
0578: GM_LineString LSCourante, LSSuivante, LSCopie;
0579: if (geometries.size() == 0) {
0580: System.out
0581: .println("ATTENTION. Erreur à la compilation de lignes : aucune ligne en entrée");
0582: return null;
0583: }
0584:
0585: LSCourante = (GM_LineString) geometries.get(0);
0586:
0587: if (geometries.size() == 1) {
0588: return LSCourante;
0589: }
0590:
0591: LSSuivante = (GM_LineString) geometries.get(1);
0592:
0593: if (Distances.proche(LSCourante.startPoint(), LSSuivante
0594: .startPoint(), 0)
0595: || Distances.proche(LSCourante.startPoint(), LSSuivante
0596: .endPoint(), 0)) {
0597: //premier point = point finale de la premiere ligne
0598: pointsFinaux.addAll(((GM_LineString) LSCourante.reverse())
0599: .getControlPoint());
0600: pointCourant = LSCourante.startPoint();
0601: } else if (Distances.proche(LSCourante.endPoint(), LSSuivante
0602: .startPoint(), 0)
0603: || Distances.proche(LSCourante.endPoint(), LSSuivante
0604: .endPoint(), 0)) {
0605: //premier point = point initial de la premiere ligne
0606: pointsFinaux.addAll(LSCourante.getControlPoint());
0607: pointCourant = LSCourante.endPoint();
0608: } else {
0609: System.out
0610: .println("ATTENTION. Erreur à la compilation de lignes (Operateurs) : les lignes ne se touchent pas");
0611: return null;
0612: }
0613:
0614: for (int i = 1; i < geometries.size(); i++) {
0615: LSSuivante = (GM_LineString) geometries.get(i);
0616: LSCopie = new GM_LineString(LSSuivante.getControlPoint());
0617: if (Distances.proche(pointCourant, LSSuivante.startPoint(),
0618: 0)) {
0619: //LSSuivante dans le bon sens
0620: LSCopie.removeControlPoint(LSCopie.startPoint());
0621: pointsFinaux.addAll(LSCopie.getControlPoint());
0622:
0623: //quel intéret à cette ligne???
0624: pointCourant = LSCopie.endPoint();
0625: } else if (Distances.proche(pointCourant, LSSuivante
0626: .endPoint(), 0)) {
0627: //LSSuivante dans le bon sens
0628: LSCopie.removeControlPoint(LSCopie.endPoint());
0629: pointsFinaux.addAll(((GM_LineString) LSCopie.reverse())
0630: .getControlPoint());
0631:
0632: //quel intéret à cette ligne???
0633: pointCourant = LSCopie.startPoint();
0634: } else {
0635: System.out
0636: .println("ATTENTION. Erreur à la compilation de lignes (Operateurs) : les lignes ne se touchent pas");
0637: return null;
0638: }
0639: }
0640:
0641: return new GM_LineString(pointsFinaux);
0642: }
0643:
0644: /** Version plus robuste mais aussi potentiellement faussée de l'intersection.
0645: * Si JTS plante au calcul d'intersection, on filtre les surfaces avec Douglas et Peucker,
0646: * progressivement avec 10 seuils entre min et max.
0647: *
0648: * English: Robust intersection of objects (to bypass JTS bugs)
0649: * @author: Mustière
0650: */
0651: public static GM_Object intersectionRobuste(GM_Object A,
0652: GM_Object B, double min, double max) {
0653: GM_Object intersection, Amodif, Bmodif;
0654: double seuilDouglas;
0655: intersection = A.intersection(B);
0656:
0657: if (intersection != null)
0658: return intersection;
0659: for (int i = 0; i < 10; i++) {
0660: seuilDouglas = min + i * (max - min) / 10;
0661: Amodif = Filtering.DouglasPeucker(A, seuilDouglas);
0662: Bmodif = Filtering.DouglasPeucker(B, seuilDouglas);
0663: intersection = Amodif.intersection(Bmodif);
0664: if (intersection != null) {
0665: System.out
0666: .println("Calcul d'intersection fait après filtrage avec Douglas Peucker à "
0667: + seuilDouglas
0668: + "m, pour cause de plantage de JTS");
0669: return intersection;
0670: }
0671: }
0672: System.out
0673: .println("ATTENTION : Plantage du calcul d'intersection, même après nettoyage de la géométrie avec Douglas Peucker");
0674: return null;
0675: }
0676:
0677: /** Version plus robuste mais aussi potentiellement faussée de l'union.
0678: * Si JTS plante au calcul d'union, on filtre les surfaces avec Douglas et Peucker,
0679: * progressivement avec 10 seuils entre min et max.
0680: *
0681: * English: Robust union of objects (to bypass JTS bugs)
0682: * @author: Mustière
0683: */
0684: public static GM_Object unionRobuste(GM_Object A, GM_Object B,
0685: double min, double max) {
0686: GM_Object union, Amodif, Bmodif;
0687: double seuilDouglas;
0688:
0689: union = A.union(B);
0690:
0691: if (union != null)
0692: return union;
0693: for (int i = 0; i < 10; i++) {
0694: seuilDouglas = min + i * (max - min) / 10;
0695: Amodif = Filtering.DouglasPeucker(A, seuilDouglas);
0696: Bmodif = Filtering.DouglasPeucker(B, seuilDouglas);
0697: union = Amodif.union(Bmodif);
0698: if (union != null) {
0699: System.out
0700: .println("Calcul d'union fait après filtrage avec Douglas Peucker à "
0701: + seuilDouglas
0702: + "m, pour cause de plantage de JTS");
0703: return union;
0704: }
0705: }
0706: System.out
0707: .println("ATTENTION : Plantage du calcul d'union, même après nettoyage de la géométrie avec Douglas Peucker");
0708: return null;
0709: }
0710:
0711: // ////////////////////////////////////////////////////////////////////
0712: // Regression linéaire
0713: // ////////////////////////////////////////////////////////////////////
0714: /** Methode qui donne l'angle (radians) par rapport à l'axe des x de la droite passant
0715: * au mieux au milieu d'un nuage de points (regression par moindres carrés).
0716: * Cet angle (défini à pi près) est entre 0 et pi.
0717: *
0718: * English: Linear approximation
0719: * @author: grosso
0720: */
0721: public static Angle directionPrincipale(DirectPositionList listePts) {
0722: Angle ang = new Angle();
0723: double angle;
0724: double x0, y0, x, y, a;
0725: double moyenneX = 0, moyenneY = 0;
0726: Matrix Atrans, A, B, X;
0727: int i;
0728:
0729: // cas où la ligne n'a que 2 pts
0730: if ((listePts.size() == 2)) {
0731: ang = new Angle(listePts.get(0), listePts.get(1));
0732: angle = ang.getAngle();
0733: if (angle >= Math.PI)
0734: return new Angle(angle - Math.PI);
0735: else
0736: return new Angle(angle);
0737: }
0738:
0739: // initialisation des matrices
0740: // On stocke les coordonnées, en se ramenant dans un repère local sur (x0,y0)
0741: A = new Matrix(listePts.size(), 1); // X des points de la ligne
0742: B = new Matrix(listePts.size(), 1); // Y des points de la ligne
0743: x0 = listePts.get(0).getX();
0744: y0 = listePts.get(0).getY();
0745: for (i = 0; i < listePts.size(); i++) {
0746: x = listePts.get(i).getX() - x0;
0747: moyenneX = moyenneX + x;
0748: A.set(i, 0, x);
0749: y = listePts.get(i).getY() - y0;
0750: moyenneY = moyenneY + y;
0751: B.set(i, 0, y);
0752: }
0753: moyenneX = moyenneX / listePts.size();
0754: moyenneY = moyenneY / listePts.size();
0755:
0756: // cas où l'angle est vertical
0757: if (moyenneX == 0)
0758: return new Angle(Math.PI / 2);
0759:
0760: // cas général : on cherche l'angle par régression linéaire
0761: Atrans = A.transpose();
0762: X = (Atrans.times(A)).inverse().times(Atrans.times(B));
0763: a = X.get(0, 0);
0764: angle = Math.atan(a);
0765: // on obtient un angle entre -pi/2 et pi/2 ouvert
0766: // on replace cet angle dans 0 et pi
0767: if (angle < 0)
0768: return new Angle(angle + Math.PI);
0769: else
0770: return new Angle(angle);
0771: }
0772:
0773: /** Methode qui donne l'angle dans [0,2*pi[ par rapport à l'axe des x,
0774: * de la droite orientée passant au mieux au milieu d'un nuage de points ordonnés
0775: * (regression par moindres carrés).
0776: * L'ordre des points en entrée est important, c'est lui qui permet de donner
0777: * l'angle à 2.pi près.
0778: * Exemple: la liste des points peut correspondre à n points d'un arc, l'angle
0779: * représente alors l'orientation générale de ces points, en prenant le premier
0780: * pour point de départ.
0781: *
0782: * English: Linear approximation
0783: * @author: grosso
0784: */
0785: public static Angle directionPrincipaleOrientee(
0786: DirectPositionList listePts) {
0787: double angle;
0788: double x0, y0, x, y, a;
0789: double moyenneX = 0, moyenneY = 0;
0790: Matrix Atrans, A, B, X;
0791: int i;
0792:
0793: // cas où la ligne n'a que 2 pts
0794: if ((listePts.size() == 2)) {
0795: return new Angle(listePts.get(0), listePts.get(1));
0796: }
0797:
0798: // initialisation des matrices
0799: // On stocke les coordonnées, en se ramenant dans un repère local sur (x0,y0)
0800: A = new Matrix(listePts.size(), 1); // X des points de la ligne
0801: B = new Matrix(listePts.size(), 1); // Y des points de la ligne
0802: x0 = listePts.get(0).getX();
0803: y0 = listePts.get(0).getY();
0804: for (i = 0; i < listePts.size(); i++) {
0805: x = listePts.get(i).getX() - x0;
0806: moyenneX = moyenneX + x;
0807: A.set(i, 0, x);
0808: y = listePts.get(i).getY() - y0;
0809: moyenneY = moyenneY + y;
0810: B.set(i, 0, y);
0811: }
0812: moyenneX = moyenneX / listePts.size();
0813: moyenneY = moyenneY / listePts.size();
0814:
0815: // cas où l'angle est vertical
0816: if (moyenneX == 0) {
0817: if (moyenneY < 0)
0818: return new Angle(3 * Math.PI / 2);
0819: if (moyenneY > 0)
0820: return new Angle(Math.PI / 2);
0821: }
0822:
0823: // cas général : on cherche l'angle par régression linéaire
0824: Atrans = A.transpose();
0825: X = (Atrans.times(A)).inverse().times(Atrans.times(B));
0826: a = X.get(0, 0);
0827: // on obtient un angle entre -pi/2 et pi/2 ouvert
0828: angle = Math.atan(a);
0829: // on replace cet angle dans 0 et 2pi
0830: if (moyenneY < 0) {
0831: if (angle >= 0)
0832: return new Angle(angle + Math.PI);
0833: else
0834: return new Angle(angle + 2 * Math.PI);
0835: }
0836: if (angle < 0)
0837: return new Angle(angle + Math.PI);
0838: return new Angle(angle);
0839: }
0840:
0841: //////////////////////////////////////////////////////////////////////
0842: // Divers
0843: //////////////////////////////////////////////////////////////////////
0844: /** Teste si 2 DirectPosition ont les mêmes coordonnées.
0845: *
0846: * English: Tests the equality of geometries
0847: */
0848: public static boolean super poses(DirectPosition pt1,
0849: DirectPosition pt2) {
0850: if (pt1.getX() != pt2.getX())
0851: return false;
0852: if (pt1.getY() != pt2.getY())
0853: return false;
0854: if (pt1.getZ() != Double.NaN) {
0855: if (pt1.getZ() != pt2.getZ())
0856: return false;
0857: }
0858: return true;
0859: }
0860:
0861: /** Teste si 2 GM_Point ont les mêmes coordonnées.
0862: *
0863: * English: Tests the equality of geometries
0864: */
0865: public static boolean super poses(GM_Point pt1, GM_Point pt2) {
0866: return super poses(pt1.getPosition(), pt2.getPosition());
0867: }
0868:
0869: /** Teste la présence d'un DirectPosition (égalité 2D) dans une
0870: * DirectPositionList.
0871: * Renvoie -1 si le directPosition n'est pas dans la liste
0872: *
0873: * English: tests if the line contains the point (in 2D)
0874: */
0875: public static int indice2D(DirectPositionList dpl, DirectPosition dp) {
0876: int i;
0877: DirectPosition dp1;
0878: for (i = 0; i < dpl.size(); i++) {
0879: dp1 = (DirectPosition) dpl.get(i);
0880: if ((dp1.getX() == dp.getX()) && (dp1.getY() == dp.getY()))
0881: return i;
0882: }
0883: return -1;
0884: }
0885:
0886: /** Teste la présence d'un DirectPosition (égalité 3D) dans une
0887: * DirectPositionList.
0888: * Renvoie -1 si le directPosition n'est pas dans la liste
0889: *
0890: * English: tests if the line contains the point (in 3D)
0891: */
0892: public static int indice3D(DirectPositionList dpl, DirectPosition dp) {
0893: int i;
0894: DirectPosition dp1;
0895: for (i = 0; i < dpl.size(); i++) {
0896: dp1 = (DirectPosition) dpl.get(i);
0897: if ((dp1.getX() == dp.getX()) && (dp1.getY() == dp.getY())
0898: && (dp1.getZ() == dp.getZ()))
0899: return i;
0900: }
0901: return -1;
0902: }
0903:
0904: /** Attribue par interpolation un Z aux points d'une ligne en connaissant le Z
0905: * des extrémités.
0906: *
0907: * English: Z interpolation
0908: * @author : Arnaud Lafragueta
0909: */
0910: public static GM_LineString calculeZ(GM_LineString ligne) {
0911:
0912: DirectPosition pointIni = ligne.startPoint();
0913: DirectPosition pointFin = ligne.endPoint();
0914: double z_ini = pointIni.getZ();
0915: double z_fin = pointFin.getZ();
0916: DirectPositionList listePoints = ligne.coord();
0917: double longueur = 0.0;
0918: double zCalc;
0919: DirectPosition pointRoute, pointRoute1;
0920: for (int j = 1; j < listePoints.size() - 1; j++) {
0921: pointRoute = listePoints.get(j);
0922: pointRoute1 = listePoints.get(j - 1);
0923: longueur = longueur
0924: + Distances.distance(pointRoute, pointRoute1);
0925: zCalc = z_ini + (z_fin - z_ini) * longueur / ligne.length();
0926: pointRoute.setZ(zCalc);
0927: ligne.setControlPoint(j, pointRoute);
0928: }
0929:
0930: return ligne;
0931: }
0932:
0933: /** Fusionne les surfaces adjacentes d'une population.
0934: * NB: quand X objets sont fusionnés, un des objets (au hasard) est gardé
0935: * avec ses attributs et sa géoémtrie est remplacée par celle fusionnée.
0936: *
0937: * English: aggregation of surfaces
0938: */
0939: public static void fusionneSurfaces(Population popSurf) {
0940:
0941: Iterator itSurf = popSurf.getElements().iterator();
0942: Iterator itSurfAdjacentes;
0943: List aEnlever = new ArrayList();
0944: GM_Object surfaceAfusionner, surfFusionnee;
0945: FT_Feature objSurf, objAfusionner, objetAEnlever;
0946:
0947: if (!popSurf.hasSpatialIndex())
0948: popSurf.initSpatialIndex(Tiling.class, true);
0949:
0950: while (itSurf.hasNext()) {
0951: objSurf = (FT_Feature) itSurf.next();
0952: if (aEnlever.contains(objSurf))
0953: continue;
0954: List surfAdjacentes = popSurf.select(objSurf.getGeom())
0955: .getElements();
0956: surfAdjacentes.remove(objSurf);
0957: if (surfAdjacentes.size() == 0)
0958: continue;
0959: aEnlever.addAll(surfAdjacentes);
0960: itSurfAdjacentes = surfAdjacentes.iterator();
0961: // ATTENTION: bidouille ci-dessous pour pallier l'absence de "copie" générique de géométrie
0962: surfFusionnee = new GM_Polygon(((GM_Polygon) objSurf
0963: .getGeom()).boundary());
0964: while (itSurfAdjacentes.hasNext()) {
0965: objAfusionner = (FT_Feature) itSurfAdjacentes.next();
0966: surfaceAfusionner = objAfusionner.getGeom();
0967: surfFusionnee = surfFusionnee.union(surfaceAfusionner);
0968: }
0969: objSurf.setGeom(surfFusionnee);
0970: }
0971: Iterator itAEnlever = aEnlever.iterator();
0972: while (itAEnlever.hasNext()) {
0973: objetAEnlever = (FT_Feature) itAEnlever.next();
0974: popSurf.enleveElement(objetAEnlever);
0975: }
0976:
0977: }
0978:
0979: /** Dilate les surfaces de la population.
0980: *
0981: * English: dilates surfaces
0982: */
0983: public static void bufferSurfaces(Population popSurf,
0984: double tailleBuffer) {
0985:
0986: FT_Feature objSurf;
0987: Iterator itSurf = popSurf.getElements().iterator();
0988:
0989: while (itSurf.hasNext()) {
0990: objSurf = (FT_Feature) itSurf.next();
0991: objSurf.setGeom(objSurf.getGeom().buffer(tailleBuffer));
0992: }
0993: }
0994:
0995: /** Surface d'un polygone (trous non gérés).
0996: * Utile pour pallier aux déficiences de JTS qui n'accèpte pas les géométries dégénérées.
0997: *
0998: * Le calcul est effectué dans un repère local centré sur le premier point
0999: * de la surface, ce qui est utile pour minimiser les erreurs de calcul
1000: * si on manipule de grandes coordonnées).
1001: *
1002: * English: surface of a polygon
1003: */
1004: public static double surface(GM_Polygon poly) {
1005: DirectPositionList pts = poly.exteriorCoord();
1006: DirectPosition pt1, pt2;
1007: double ymin;
1008: double surf = 0;
1009:
1010: pt1 = pts.get(0);
1011: ymin = pt1.getY();
1012: for (int i = 1; i < pts.size(); i++) {
1013: pt2 = pts.get(i);
1014: surf = surf
1015: + ((pt2.getX() - pt1.getX()) * (pt2.getY()
1016: + pt1.getY() - 2 * ymin)) / 2;
1017: pt1 = pt2;
1018: }
1019: return Math.abs(surf);
1020: }
1021:
1022: /** Surface d'un polygone (liste de points supposée fermée).
1023: *
1024: * English: surface of a polygon
1025: */
1026: public static double surface(DirectPositionList pts) {
1027: DirectPosition pt1, pt2;
1028: double ymin;
1029: double surf = 0;
1030:
1031: pt1 = pts.get(0);
1032: ymin = pt1.getY();
1033: for (int i = 1; i < pts.size(); i++) {
1034: pt2 = pts.get(i);
1035: surf = surf
1036: + ((pt2.getX() - pt1.getX()) * (pt2.getY()
1037: + pt1.getY() - 2 * ymin)) / 2;
1038: pt1 = pt2;
1039: }
1040: return Math.abs(surf);
1041: }
1042:
1043: /** Détermine si une liste de points tourne dans le sens direct ou non.
1044: * NB : La liste de points est supposée fermée (premier point = dernier point).
1045: * NB : renvoie true pour une surface dégénérée.
1046: *
1047: * English : orientation of a polygon (direct rotation?)
1048: */
1049: public static boolean sensDirect(DirectPositionList pts) {
1050: Iterator itPts = pts.getList().iterator();
1051: DirectPosition pt1, pt2;
1052: double ymin;
1053: double surf = 0;
1054: pt1 = (DirectPosition) itPts.next();
1055: ymin = pt1.getY();
1056: while (itPts.hasNext()) {
1057: pt2 = (DirectPosition) itPts.next();
1058: surf = surf
1059: + ((pt2.getX() - pt1.getX()) * (pt2.getY()
1060: + pt1.getY() - 2 * ymin));
1061: pt1 = pt2;
1062: }
1063: if (surf > 0)
1064: return false;
1065: else
1066: return true;
1067: }
1068: }
|