Skip to content

Commit

Permalink
Add Z handling to OverlayNG (#645)
Browse files Browse the repository at this point in the history 
* Add Z interpolation and population from ElevationModel
* Disable error when reading Z from CoordinateXY,
to work around mixed-coordinate class issue

Signed-off-by: Martin Davis <[email protected]>
  • Loading branch information
@dr-jts
dr-jts authored Nov 27, 2020
1 parent 5355aef commit 4c88fea
Show file tree
Hide file tree
Showing 8 changed files with 1,053 additions and 54 deletions.
226 changes: 184 additions & 42 deletions modules/core/src/main/java/org/locationtech/jts/algorithm/RobustLineIntersector.java
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@


/*
* Copyright (c) 2016 Vivid Solutions.
* Copyright (c) 2020 Martin Davis.
*
* All rights reserved. This program and the accompanying materials
* are made available under the terms of the Eclipse Public License 2.0
Expand Down Expand Up @@ -73,7 +72,9 @@ protected int computeIntersect(
if ((Qp1>0 && Qp2>0) || (Qp1<0 && Qp2<0)) {
return NO_INTERSECTION;
}

/**
* Intersection is collinear if each endpoint lies on the other line.
*/
boolean collinear = Pq1 == 0
&& Pq2 == 0
&& Qp1 == 0
Expand All @@ -95,6 +96,8 @@ protected int computeIntersect(
* the other line, since at this point we know that the inputLines must
* intersect.
*/
Coordinate p = null;
double z = Double.NaN;
if (Pq1 == 0 || Pq2 == 0 || Qp1 == 0 || Qp2 == 0) {
isProper = false;

Expand All @@ -114,78 +117,111 @@ protected int computeIntersect(
* which used to produce the INCORRECT result: (20.31970698357233, 46.76654261437082, NaN)
*
*/
if (p1.equals2D(q1)
|| p1.equals2D(q2)) {
intPt[0] = p1;
if (p1.equals2D(q1)) {
p = p1;
z = zGet(p1, q1);
}
else if (p2.equals2D(q1)
|| p2.equals2D(q2)) {
intPt[0] = p2;
else if (p1.equals2D(q2)) {
p = p1;
z = zGet(p1, q2);
}
else if (p2.equals2D(q1)) {
p = p2;
z = zGet(p2, q1);
}
else if (p2.equals2D(q2)) {
p = p2;
z = zGet(p2, q2);
}

/**
* Now check to see if any endpoint lies on the interior of the other segment.
*/
else if (Pq1 == 0) {
intPt[0] = new Coordinate(q1);
p = q1;
z = zGetOrInterpolate(q1, p1, p2);
}
else if (Pq2 == 0) {
intPt[0] = new Coordinate(q2);
p = q2;
z = zGetOrInterpolate(q2, p1, p2);
}
else if (Qp1 == 0) {
intPt[0] = new Coordinate(p1);
p = p1;
z = zGetOrInterpolate(p1, q1, q2);
}
else if (Qp2 == 0) {
intPt[0] = new Coordinate(p2);
p = p2;
z = zGetOrInterpolate(p2, q1, q2);
}
}
else {
isProper = true;
intPt[0] = intersection(p1, p2, q1, q2);
p = intersection(p1, p2, q1, q2);
z = zInterpolate(p, p1, p2, q1, q2);
}
intPt[0] = copyWithZ(p, z);
return POINT_INTERSECTION;
}

private int computeCollinearIntersection(Coordinate p1, Coordinate p2,
Coordinate q1, Coordinate q2) {
boolean p1q1p2 = Envelope.intersects(p1, p2, q1);
boolean p1q2p2 = Envelope.intersects(p1, p2, q2);
boolean q1p1q2 = Envelope.intersects(q1, q2, p1);
boolean q1p2q2 = Envelope.intersects(q1, q2, p2);
boolean q1inP = Envelope.intersects(p1, p2, q1);
boolean q2inP = Envelope.intersects(p1, p2, q2);
boolean p1inQ = Envelope.intersects(q1, q2, p1);
boolean p2inQ = Envelope.intersects(q1, q2, p2);

if (p1q1p2 && p1q2p2) {
intPt[0] = q1;
intPt[1] = q2;
if (q1inP && q2inP) {
intPt[0] = copyWithZInterpolate(q1, p1, p2);
intPt[1] = copyWithZInterpolate(q2, p1, p2);
return COLLINEAR_INTERSECTION;
}
if (q1p1q2 && q1p2q2) {
intPt[0] = p1;
intPt[1] = p2;
if (p1inQ && p2inQ) {
intPt[0] = copyWithZInterpolate(p1, q1, q2);
intPt[1] = copyWithZInterpolate(p2, q1, q2);
return COLLINEAR_INTERSECTION;
}
if (p1q1p2 && q1p1q2) {
intPt[0] = q1;
intPt[1] = p1;
return q1.equals(p1) && !p1q2p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
if (q1inP && p1inQ) {
// if pts are equal Z is chosen arbitrarily
intPt[0] = copyWithZInterpolate(q1, p1, p2);
intPt[1] = copyWithZInterpolate(p1, q1, q2);
return q1.equals(p1) && !q2inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
if (p1q1p2 && q1p2q2) {
intPt[0] = q1;
intPt[1] = p2;
return q1.equals(p2) && !p1q2p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
if (q1inP && p2inQ) {
// if pts are equal Z is chosen arbitrarily
intPt[0] = copyWithZInterpolate(q1, p1, p2);
intPt[1] = copyWithZInterpolate(p2, q1, q2);
return q1.equals(p2) && !q2inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
if (p1q2p2 && q1p1q2) {
intPt[0] = q2;
intPt[1] = p1;
return q2.equals(p1) && !p1q1p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
if (q2inP && p1inQ) {
// if pts are equal Z is chosen arbitrarily
intPt[0] = copyWithZInterpolate(q2, p1, p2);
intPt[1] = copyWithZInterpolate(p1, q1, q2);
return q2.equals(p1) && !q1inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
if (p1q2p2 && q1p2q2) {
intPt[0] = q2;
intPt[1] = p2;
return q2.equals(p2) && !p1q1p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
if (q2inP && p2inQ) {
// if pts are equal Z is chosen arbitrarily
intPt[0] = copyWithZInterpolate(q2, p1, p2);
intPt[1] = copyWithZInterpolate(p2, q1, q2);
return q2.equals(p2) && !q1inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
return NO_INTERSECTION;
}

private static Coordinate copyWithZInterpolate(Coordinate p, Coordinate p1, Coordinate p2) {
return copyWithZ(p, zGetOrInterpolate(p, p1, p2));
}

private static Coordinate copyWithZ(Coordinate p, double z) {
Coordinate pCopy = copy(p);
if (! Double.isNaN(z)) {
pCopy.setZ( z );
}
return pCopy;
}

private static Coordinate copy(Coordinate p) {
return new Coordinate(p);
}

/**
* This method computes the actual value of the intersection point.
* To obtain the maximum precision from the intersection calculation,
Expand Down Expand Up @@ -226,7 +262,7 @@ private Coordinate intersection(

// compute a safer result
// copy the coordinate, since it may be rounded later
intPt = new Coordinate(nearestEndpoint(p1, p2, q1, q2));
intPt = copy(nearestEndpoint(p1, p2, q1, q2));
// intPt = CentralEndpointIntersector.getIntersection(p1, p2, q1, q2);

// System.out.println("Segments: " + this);
Expand Down Expand Up @@ -330,5 +366,111 @@ private static Coordinate nearestEndpoint(Coordinate p1, Coordinate p2,
return nearestPt;
}

/**
* Gets the Z value of the first argument if present,
* otherwise the value of the second argument.
*
* @param p a coordinate, possibly with Z
* @param q a coordinate, possibly with Z
* @return the Z value if present
*/
private static double zGet(Coordinate p, Coordinate q) {
double z = p.getZ();
if (Double.isNaN(z)) {
z = q.getZ(); // may be NaN
}
return z;
}

/**
* Gets the Z value of a coordinate if present, or
* interpolates it from the segment it lies on.
* If the segment Z values are not fully populate
* NaN is returned.
*
* @param p a coordinate, possibly with Z
* @param p1 a segment endpoint, possibly with Z
* @param p2 a segment endpoint, possibly with Z
* @return the extracted or interpolated Z value (may be NaN)
*/
private static double zGetOrInterpolate(Coordinate p, Coordinate p1, Coordinate p2) {
double z = p.getZ();
if (! Double.isNaN(z))
return z;
return zInterpolate(p, p1, p2); // may be NaN
}

/**
* Interpolates a Z value for a point along
* a line segment between two points.
* The Z value of the interpolation point (if any) is ignored.
* If either segment point is missing Z,
* returns NaN.
*
* @param p a coordinate
* @param p1 a segment endpoint, possibly with Z
* @param p2 a segment endpoint, possibly with Z
* @return the interpolated Z value (may be NaN)
*/
private static double zInterpolate(Coordinate p, Coordinate p1, Coordinate p2) {
double p1z = p1.getZ();
double p2z = p2.getZ();
if (Double.isNaN(p1z)) {
return p2z; // may be NaN
}
if (Double.isNaN(p2z)) {
return p1z; // may be NaN
}
if (p.equals2D(p1)) {
return p1z; // not NaN
}
if (p.equals2D(p2)) {
return p2z; // not NaN
}
double dz = p2z - p1z;
if (dz == 0.0) {
return p1z;
}
// interpolate Z from distance of p along p1-p2
double dx = (p2.x - p1.x);
double dy = (p2.y - p1.y);
// seg has non-zero length since p1 < p < p2
double seglen = (dx * dx + dy * dy);
double xoff = (p.x - p1.x);
double yoff = (p.y - p1.y);
double plen = (xoff * xoff + yoff * yoff);
double frac = Math.sqrt(plen / seglen);
double zoff = dz * frac;
double zInterpolated = p1z + zoff;
return zInterpolated;
}

/**
* Interpolates a Z value for a point along
* two line segments and computes their average.
* The Z value of the interpolation point (if any) is ignored.
* If one segment point is missing Z that segment is ignored
* if both segments are missing Z, returns NaN.
*
* @param p a coordinate
* @param p1 a segment endpoint, possibly with Z
* @param p2 a segment endpoint, possibly with Z
* @param q1 a segment endpoint, possibly with Z
* @param q2 a segment endpoint, possibly with Z
* @return the averaged interpolated Z value (may be NaN)
*/
private static double zInterpolate(Coordinate p, Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2) {
double zp = zInterpolate(p, p1, p2);
double zq = zInterpolate(p, q1, q2);
if (Double.isNaN(zp)) {
return zq; // may be NaN
}
if (Double.isNaN(zq)) {
return zp; // may be NaN
}
// both Zs have values, so average them
return (zp + zq) / 2.0;
}


}
4 changes: 3 additions & 1 deletion modules/core/src/main/java/org/locationtech/jts/geom/CoordinateXY.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,9 @@ public double getOrdinate(int ordinateIndex) {
case X: return x;
case Y: return y;
}
throw new IllegalArgumentException("Invalid ordinate index: " + ordinateIndex);
return Double.NaN;
// disable for now to avoid regression issues
//throw new IllegalArgumentException("Invalid ordinate index: " + ordinateIndex);
}

@Override
Expand Down
Loading

0 comments on commit 4c88fea

Please sign in to comment.