Menu
▾
▴
[Fudaa-svn-commit] SF.net SVN: fudaa:[4468] branches/FudaaModeleur_TC1Bis/fudaa_devel
|
From: <emm...@us...> - 2009-02-13 15:33:09
|
Revision: 4468
http://fudaa.svn.sourceforge.net/fudaa/?rev=4468&view=rev
Author: emmanuel_martin
Date: 2009-02-13 15:33:04 +0000 (Fri, 13 Feb 2009)
Log Message:
-----------
Modeleur 1d : la trace de profil lors de l'exportation somporte jusqu'?\195?\160 4 points.
Modified Paths:
--------------
branches/FudaaModeleur_TC1Bis/fudaa_devel/dodico/src/org/fudaa/dodico/mascaret/io/MascaretWriter.java
branches/FudaaModeleur_TC1Bis/fudaa_devel/fudaa/src/org/fudaa/fudaa/modeleur/modeleur1d/model/UtilsProfil1d.java
Modified: branches/FudaaModeleur_TC1Bis/fudaa_devel/dodico/src/org/fudaa/dodico/mascaret/io/MascaretWriter.java
===================================================================
--- branches/FudaaModeleur_TC1Bis/fudaa_devel/dodico/src/org/fudaa/dodico/mascaret/io/MascaretWriter.java 2009-02-13 15:28:57 UTC (rev 4467)
+++ branches/FudaaModeleur_TC1Bis/fudaa_devel/dodico/src/org/fudaa/dodico/mascaret/io/MascaretWriter.java 2009-02-13 15:33:04 UTC (rev 4468)
@@ -20,8 +20,7 @@
import org.fudaa.ctulu.fileformat.FortranInterface;
import org.fudaa.ctulu.gis.GISAttributeConstants;
import org.fudaa.ctulu.gis.GISAttributeModelDoubleInterface;
-import org.fudaa.ctulu.gis.GISGeometryFactory;
-import org.fudaa.ctulu.gis.GISPoint;
+import org.fudaa.ctulu.gis.GISLib;
import org.fudaa.ctulu.gis.GISPolyligne;
import org.fudaa.ctulu.gis.GISZoneCollection;
import org.fudaa.dodico.commun.DodicoLib;
@@ -218,11 +217,20 @@
throw new MascaretDataError("Un profil \xE0 plusieurs ou aucun point(s) d'intersection(s) avec l'axe hydraulique.");
profilAbs.coordIntersectionAxeHydraulique=new double[]{intersection.getCoordinate().x, intersection.getCoordinate().y};
// Abscisse curviligne axe hydraulique
+ int idxAttrCommHydraulique=zoneProfil.getIndiceOf(GISAttributeConstants.COMMENTAIRE_HYDRO);
+ if(idxAttrCommHydraulique==-1)
+ throw new MascaretDataError("L'attribut COMMENTAIRE_HYDRO doit \xEAtre pr\xE9sent.");
+ String comm=(String) zoneProfil.getValue(idxAttrCommHydraulique, idxProfil);
+ if(!GISLib.isHydroCommentValued(comm, "PK"))
+ throw new MascaretDataError("La propri\xE9t\xE9 PK de COMMENTAIRE_HYDRO doit \xEAtre valu\xE9e.");
+ profilAbs.abscisseCurvilingeAxeHydraulique=GISLib.getHydroCommentDouble(comm, "PK");
CoordinateSequence coordSeq=axeHydraulique.getCoordinateSequence();
- profilAbs.abscisseCurvilingeAxeHydraulique=absCurv(coordSeq, intersection.getCoordinate());
// Trace profil
coordSeq=profil.getCoordinateSequence();
+ Coordinate coord[]=getPointsRuptures(coordSeq);
profilAbs.coordXYTraceProfil.add(new double[]{coordSeq.getX(0), coordSeq.getY(0)});
+ for(int j=0;j<coord.length;j++)
+ profilAbs.coordXYTraceProfil.add(new double[]{coord[j].x, coord[j].y});
profilAbs.coordXYTraceProfil.add(new double[]{coordSeq.getX(coordSeq.size()-1), coordSeq.getY(coordSeq.size()-1)});
// Information atomiques sur le profil \\
@@ -325,20 +333,124 @@
}
/**
- * Calcul l'abscisse curviligne du point indiqu\xE9 en param\xE8tre
+ * Retourne la liste des points de ruptures.
*/
- private double absCurv(CoordinateSequence _coordSeq, Coordinate _point){
- int k=idxIntersection(_coordSeq, _point);
- // Calcule de l'abscisse curviligne
- double valueCurviligne=0;
- for (int j=1; j<=k; j++)
- valueCurviligne+=_coordSeq.getCoordinate(j).distance(_coordSeq.getCoordinate(j-1));
- // Mise \xE0 jour de l'abscisse curviligne
- valueCurviligne+=_coordSeq.getCoordinate(k).distance(_point);
- return valueCurviligne;
+ private Coordinate[] getPointsRuptures(CoordinateSequence _seq) {
+ List<Coordinate> coord=new ArrayList<Coordinate>();
+ int idx1=0;
+ int idx2=1;
+ int idx3=2;
+ int idxPointRupturePrecedent=-1;
+ while(idx3<_seq.size()) {
+ if(idx1==idx2)
+ idx2++;
+ else if(idx2==idx3)
+ idx3++;
+ else if(_seq.getCoordinate(idx1).equals(_seq.getCoordinate(idx2)))
+ idx2++;
+ else if(_seq.getCoordinate(idx2).equals(_seq.getCoordinate(idx3)))
+ idx3++;
+ else {
+ if(getAngle(_seq.getCoordinate(idx1), _seq.getCoordinate(idx2), _seq.getCoordinate(idx3))<(Math.PI-0.0001))
+ if(idxPointRupturePrecedent<0||!_seq.getCoordinate(idxPointRupturePrecedent).equals(_seq.getCoordinate(idx2)))
+ coord.add(_seq.getCoordinate(idx2));
+ if(idx2+1<idx3)
+ idx2++;
+ else
+ idx1=idx2;
+ }
+ }
+ return coord.toArray(new Coordinate[0]);
}
/**
+ * Retourne l'angle form\xE9 par les segments _a_b et _b_c. _a, _b et _c doivent
+ * \xEAtre non confondus sous peine d'avoir un {@link IllegalArgumentException}.
+ * L'angle est toujours positif.
+ * L'ordonn\xE9 z n'est pas prise en compte.
+ * L'angle retourn\xE9 est en radian.
+ */
+ private double getAngle(Coordinate _a, Coordinate _b, Coordinate _c) {
+ Coordinate a=new Coordinate(_a.x, _a.y, 0);
+ Coordinate b=new Coordinate(_b.x, _b.y, 0);
+ Coordinate c=new Coordinate(_c.x, _c.y, 0);
+ if(egal(a, b)||egal(b, c)||egal(a, c))
+ throw new IllegalArgumentException("Les trois points doivent \xEAtre non confonus.");
+ // Projection de b sur ac
+ Coordinate b2=proj(b, a, c);
+ // Calcul des angles interm\xE9diaires
+ double angleCBB2=Math.asin(b2.distance(c)/c.distance(b));
+ double angleB2BA=Math.asin(b2.distance(a)/a.distance(b));
+ if(angleCBB2==Double.NaN||angleB2BA==Double.NaN)
+ return Math.PI;
+ // Selon que B2 appartient ou non \xE0 ac, la d\xE9termination de l'angle final diff\xE8re.
+ double angle;
+ int position=getPositionXY(b2, a, c);
+ if(position==-1)
+ angle=angleCBB2-angleB2BA;
+ else if(position==1)
+ angle=angleB2BA-angleCBB2;
+ else
+ angle=angleCBB2+angleB2BA;
+ return angle;
+ }
+
+ /**
+ * Soit _a, _b, _c align\xE9s et non confondus. Retourne -1 si _a avant _b_c, 0
+ * si _a est entre _b_c (inclu), +1 sinon. Si _a, _b et _c ne sont pas
+ * align\xE9s, {@link IllegalArgumentException} est lev\xE9.
+ */
+ private int getPositionXY(Coordinate _a, Coordinate _b, Coordinate _c) {
+ if(!egal(proj(_a, _b, _c), _a)||egal(_a, _b)||egal(_b, _c)||egal(_a, _c))
+ throw new IllegalArgumentException("Les trois points doivent \xEAtre align\xE9s et non confonus.");
+ Coordinate vecBC=vec(_b, _c);
+ double k;
+ if(vecBC.x!=0)
+ k=(_a.x-_b.x)/vecBC.x;
+ else
+ k=(_a.y-_b.y)/vecBC.y;
+ if(k<0)
+ return -1;
+ else if(k>1)
+ return 1;
+ else
+ return 0;
+ }
+
+ /**
+ * Retourne la projection de _a sur la droite _b_c.
+ * Attention : Il s'agit bien d'une projection sur une
+ * droite et non un segment de droite : le r\xE9sultat
+ * peut \xEAtre en dehors de [_b, _c]
+ */
+ private Coordinate proj(Coordinate _a, Coordinate _b, Coordinate _c){
+ Coordinate vBC=vec(_b, _c);
+ Coordinate vBA=vec(_b, _a);
+ double normeBC=Math.sqrt(vBC.x*vBC.x+vBC.y*vBC.y+vBC.z*vBC.z);
+ double produitScalaireBCBA=vBC.x*vBA.x+vBC.y*vBA.y+vBC.z*vBA.z;
+ double rapportSurBC=produitScalaireBCBA/(normeBC*normeBC);
+ return new Coordinate(vBC.x*rapportSurBC+_b.x, vBC.y*rapportSurBC+_b.y, vBC.z*rapportSurBC+_b.z);
+ }
+
+ /**
+ * Compars les deux vecteurs avec la tol\xE9rance donn\xE9e. Le dernier param\xE8tre
+ * indique si un normalisation des deux vecteurs doit \xEAtre faite avec la
+ * comparaison. Les vecteurs ne sont pas modif\xE9s dans l'op\xE9ration. _testZ doit
+ * \xEAtre mit \xE0 vrai pour que la comparaison tienne compte du z.
+ */
+ private boolean egal(Coordinate _a, Coordinate _b) {
+ Coordinate diff=vec(_a, _b);
+ return Math.abs(diff.x)<0.0001&&Math.abs(diff.y)<0.0001;
+ }
+
+ /**
+ * retourne le vecteur _a, _b sous forme de Coordinate.
+ */
+ private Coordinate vec(Coordinate _a, Coordinate _b){
+ return new Coordinate(_b.x-_a.x, _b.y-_a.y, _b.z-_a.z);
+ }
+
+ /**
* Calcul l'abscisse curviligne du point indiqu\xE9 en param\xE8tre
*/
private double absCurv(CoordinateSequence _coordSeq, int _idxPoint){
@@ -349,19 +461,4 @@
valueCurviligne+=_coordSeq.getCoordinate(j).distance(_coordSeq.getCoordinate(j-1));
return valueCurviligne;
}
-
- /**
- * Retourne l'indice pr\xE9c\xE9dent d'un point d'intersectin avec une geometrie.
- */
- private int idxIntersection(CoordinateSequence _coordSeq, Coordinate _point){
- /* Recherche de l'index du dernier point de la polyligne \xE0 prendre en
- * compte. */
- boolean fini=false;
- int k=-1;
- GISPoint point= new GISPoint(_point);
- while (!fini&&++k<_coordSeq.size()-1)
- fini=(GISGeometryFactory.INSTANCE
- .createLineString(new Coordinate[]{_coordSeq.getCoordinate(k), _coordSeq.getCoordinate(k+1)})).intersects(point);
- return k;
- }
}
Modified: branches/FudaaModeleur_TC1Bis/fudaa_devel/fudaa/src/org/fudaa/fudaa/modeleur/modeleur1d/model/UtilsProfil1d.java
===================================================================
--- branches/FudaaModeleur_TC1Bis/fudaa_devel/fudaa/src/org/fudaa/fudaa/modeleur/modeleur1d/model/UtilsProfil1d.java 2009-02-13 15:28:57 UTC (rev 4467)
+++ branches/FudaaModeleur_TC1Bis/fudaa_devel/fudaa/src/org/fudaa/fudaa/modeleur/modeleur1d/model/UtilsProfil1d.java 2009-02-13 15:33:04 UTC (rev 4468)
@@ -52,7 +52,7 @@
GISGeometryFactory.INSTANCE.createLineString(new Coordinate[]{seq.getCoordinate(i), seq.getCoordinate(i+1)}))
.getNumPoints()>1)
ok=false;
- // Aucun point d'intersectin avec les axes qui suivent
+ // Aucun point d'intersection avec les axes qui suivent
int j=i;
while (ok&&++j+1<seq.size())
ok=!mainAxe.intersects(GISGeometryFactory.INSTANCE.createLineString(new Coordinate[]{seq.getCoordinate(j),
@@ -65,49 +65,30 @@
// / D\xE9termination des points de ruptures et verification qu'il y en a
// au maximun deux. \\\
int nbPointsRupture=0;
- double upsilon=0.0001; // Marge d'erreur des calculs sur double
- double coefDirecteur=0;
- // coefDirX : Indique si le coefficient directeur est calcul\xE9 sur les x
- // ou sur les y (au d\xE9nominateur).
- boolean coefDirX;
- // Initialisation : calcul du premier cofficient directeur \\
- int idx=1;
- // Calcul du coefficient directeur
- if (seq.getX(idx)!=seq.getX(idx-1)) {
- // Calcul sur les x
- coefDirX=true;
- coefDirecteur=(seq.getY(idx)-seq.getY(idx-1))/(seq.getX(idx)-seq.getX(idx-1));
- }
- else {
- // Calcul sur les y
- coefDirX=false;
- coefDirecteur=0;
- }
- // Coeur : Calcul de tous les coefficients directeur et m\xE9morisation des
- // ruptures \\
- for (idx=2; idx<seq.size(); idx++) {
- double newCoefDirecteur;
- boolean newCoefDirX;
- // Calcul du coefficient directeur
- if (seq.getX(idx)!=seq.getX(idx-1)) {
- // Calcul sur les x
- newCoefDirX=true;
- newCoefDirecteur=(seq.getY(idx)-seq.getY(idx-1))/(seq.getX(idx)-seq.getX(idx-1));
- }
+ int idx1=0;
+ int idx2=1;
+ int idx3=2;
+ int idxPointRupturePrecedent=-1;
+ while(idx3<seq.size()) {
+ if(idx1==idx2)
+ idx2++;
+ else if(idx2==idx3)
+ idx3++;
+ else if(seq.getCoordinate(idx1).equals(seq.getCoordinate(idx2)))
+ idx2++;
+ else if(seq.getCoordinate(idx2).equals(seq.getCoordinate(idx3)))
+ idx3++;
else {
- // Calcul sur les y
- newCoefDirX=false;
- newCoefDirecteur=0;
- }
- // D\xE9tection des ruptures
- if (coefDirX!=newCoefDirX||Math.abs(coefDirecteur-newCoefDirecteur)>upsilon) {
- nbPointsRupture++;
+ if(UtilsProfil1d.getAngle(seq.getCoordinate(idx1), seq.getCoordinate(idx2), seq.getCoordinate(idx3))<(Math.PI-0.0001))
+ if(idxPointRupturePrecedent<0||!seq.getCoordinate(idxPointRupturePrecedent).equals(seq.getCoordinate(idx2)))
+ nbPointsRupture++;
if (nbPointsRupture>2)
return false;
+ if(idx2+1<idx3)
+ idx2++;
+ else
+ idx1=idx2;
}
- // Pr\xE9paration de l'it\xE9ration suivante
- coefDirecteur=newCoefDirecteur;
- coefDirX=newCoefDirX;
}
return true;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
