# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1538028928 -7200
# Node ID 9f3a4a60dc04640f63ad310fa5971621c28ac153
# Parent  3d927e06b92cf3e2cf12b1ae3b6d9bc2f0628197
Cleaned up interpolation mess in vertical triangle intersection.

diff -r 3d927e06b92c -r 9f3a4a60dc04 pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Wed Sep 26 18:36:12 2018 +0200
+++ b/pkg/octree/vertex.go	Thu Sep 27 08:15:28 2018 +0200
@@ -428,7 +428,7 @@
 	return math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))
 }
 
-func interpolate(v1, v2 Vertex) func(x, y float64) float64 {
+func relative(v1, v2 Vertex) func(x, y float64) float64 {
 
 	length := dist(v1.X, v1.Y, v2.X, v2.Y)
 
@@ -446,6 +446,12 @@
 	}
 }
 
+func interpolate(a, b float64) func(float64) float64 {
+	return func(x float64) float64 {
+		return (b-a)*x + a
+	}
+}
+
 func linear(v1, v2 Vertex) func(t float64) Vertex {
 	return func(t float64) Vertex {
 		return Vertex{
@@ -456,19 +462,12 @@
 	}
 }
 
-func inRange(a, b, c float64) bool {
-	if a > b {
-		a, b = b, a
-	}
-	return c >= a && c <= b
-}
+func inRange(a float64) bool { return 0 <= a && a <= 1 }
 
 func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ {
 
 	var out LineStringZ
 
-	length := dist(vl.x1, vl.y1, vl.x2, vl.y2)
-
 edges:
 	for i := 0; i < 3 && len(out) < 2; i++ {
 		j := (i + 1) % 3
@@ -482,19 +481,17 @@
 
 		switch {
 		case o1 && o2:
-			// XXX: Broken!
-			pos := interpolate(t[i], t[j])
+			pos := relative(t[i], t[j])
 			t1 := pos(vl.x1, vl.y1)
 			t2 := pos(vl.x2, vl.y2)
 
-			r1 := inRange(0, length, t1)
-			r2 := inRange(0, length, t2)
+			r1 := inRange(t1)
+			r2 := inRange(t2)
 
 			switch {
 			case r1 && r2:
-				l := linear(t[i], t[j])
-				// XXX: Broken!
-				out = append(out, l(t1/length), l(t2/length))
+				lin := linear(t[i], t[j])
+				out = append(out, lin(t1), lin(t2))
 				return out
 
 			case !r1 && !r2: // the hole edge
@@ -507,15 +504,15 @@
 			}
 
 		case o1:
-			t1 := interpolate(t[i], t[j])(vl.x1, vl.y1)
-			if !inRange(0, length, t1) {
+			t1 := relative(t[i], t[j])(vl.x1, vl.y1)
+			if !inRange(t1) {
 				continue edges
 			}
 			out = append(out, linear(t[i], t[j])(t1))
 
 		case o2:
-			t2 := interpolate(t[i], t[j])(vl.x2, vl.y2)
-			if !inRange(0, length, t2) {
+			t2 := relative(t[i], t[j])(vl.x2, vl.y2)
+			if !inRange(t2) {
 				continue edges
 			}
 			out = append(out, linear(t[i], t[j])(t2))
@@ -525,18 +522,17 @@
 			if !intersects {
 				continue edges
 			}
-			t1 := interpolate(t[i], t[j])(x, y)
-			if !inRange(0, length, t1) {
+			t1 := relative(t[i], t[j])(x, y)
+			if !inRange(t1) {
 				continue edges
 			}
 
-			// XXX: Ugly!
-			z := (t[j].Z-t[i].Z)*(t1/length) + t[1].Z
+			z := interpolate(t[j].Z, t[i].Z)(t1)
 			n := Vertex{x, y, z}
 
 			if math.Signbit(s1) != math.Signbit(s2) {
 				// different sides -> vertex on edge.
-				out = append(out)
+				out = append(out, n)
 			} else { // same side -> inside. Find peer
 				if len(out) > 0 { // we already have our peer
 					out = append(out, n)
@@ -553,11 +549,12 @@
 					if !intersects {
 						continue
 					}
-					// XXX: Broken!!!
-					zo := interpolate(t[k], t[j])(xo, yo)
-					if !inRange(t[k].Z, t[k].Z, zo) {
+					t2 := relative(t[k], t[j])(xo, yo)
+					if !inRange(t2) {
 						continue
 					}
+					zo := interpolate(t[k].Z, t[j].Z)(t2)
+
 					m := Vertex{xo, yo, zo}
 
 					var xn, yn, xf, yf float64
@@ -574,29 +571,30 @@
 					if onPlane(other.Eval(xn, yn)) {
 						// triangle intersect in only point
 						// treat as no intersection
-						return out
+						break edges
 					}
 
-					// XXX: Broken!!!
-					zint := interpolate(n, m)
+					pos := relative(n, m)
 
-					zn := zint(xn, yn)
-					zf := zint(xf, yf)
+					tzn := pos(xn, yn)
+					tzf := pos(xf, yf)
 
-					if !inRange(n.Z, m.Z, zn) {
+					if !inRange(tzn) {
 						// if near is out of range far is, too.
 						return out
 					}
 
-					if inRange(n.Z, m.Z, zf) {
+					lin := interpolate(n.Z, m.Z)
+
+					if inRange(tzf) {
 						m.X = xf
 						m.Y = yf
-						m.Z = zf
+						m.Z = lin(tzf)
 					} // else m is clipping
 
 					n.X = xn
 					n.Y = yn
-					n.Z = zn
+					n.Z = lin(tzn)
 
 					out = append(out, n, m)
 					return out