Mercurial > gemma
changeset 787:3d927e06b92c
Triangle intersection. WIP. Currently interpolation is messed up.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 26 Sep 2018 18:36:12 +0200 |
parents | 1bee00039973 |
children | 9f3a4a60dc04 |
files | pkg/controllers/octreecross.go pkg/octree/tree.go pkg/octree/vertex.go |
diffstat | 3 files changed, 214 insertions(+), 5 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/controllers/octreecross.go Wed Sep 26 17:04:47 2018 +0200 +++ b/pkg/controllers/octreecross.go Wed Sep 26 18:36:12 2018 +0200 @@ -103,10 +103,11 @@ c1 := &coords[i] c2 := &coords[i+1] + verticalLine := octree.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon) + var line octree.MultiLineStringZ tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *octree.Triangle) { - - if ls := t.IntersectVertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon); len(ls) > 0 { + if ls := verticalLine.Intersection(t); len(ls) > 0 { line = append(line) } })
--- a/pkg/octree/tree.go Wed Sep 26 17:04:47 2018 +0200 +++ b/pkg/octree/tree.go Wed Sep 26 18:36:12 2018 +0200 @@ -100,9 +100,7 @@ v1 := line.Eval(t[1].X, t[1].Y) v2 := line.Eval(t[2].X, t[2].Y) - if math.Abs(v0) < epsPlane || - math.Abs(v1) < epsPlane || - math.Abs(v2) < epsPlane || + if onPlane(v0) || onPlane(v1) || onPlane(v2) || sides(sides(sides(0, v0), v1), v2) == 3 { fn(&t) }
--- a/pkg/octree/vertex.go Wed Sep 26 17:04:47 2018 +0200 +++ b/pkg/octree/vertex.go Wed Sep 26 18:36:12 2018 +0200 @@ -402,3 +402,213 @@ return p.X / p.Z, p.Y / p.Z, true } + +type VerticalLine struct { + x1 float64 + y1 float64 + x2 float64 + y2 float64 + + line Plane2D +} + +func NewVerticalLine(x1, y1, x2, y2 float64) *VerticalLine { + return &VerticalLine{ + x1: x1, + y1: y1, + x2: x2, + y2: y2, + line: NewPlane2D(x1, y1, x2, y2), + } +} + +func onPlane(x float64) bool { return math.Abs(x) < epsPlane } + +func dist(x1, y1, x2, y2 float64) float64 { + return math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)) +} + +func interpolate(v1, v2 Vertex) func(x, y float64) float64 { + + length := dist(v1.X, v1.Y, v2.X, v2.Y) + + // f(0) = v1.Z + // f(length) = v2.Z + // v1.Z = 0*m + b + // v2.Z = length*m + b + + b := v1.Z + m := (v2.Z - v1.Z) / length + + return func(x, y float64) float64 { + l := dist(v1.X, v1.Y, x, y) + return m*l + b + } +} + +func linear(v1, v2 Vertex) func(t float64) Vertex { + return func(t float64) Vertex { + return Vertex{ + (v2.X-v1.X)*t + v1.X, + (v2.Y-v1.Y)*t + v1.Y, + (v2.Z-v1.Z)*t + v1.Z, + } + } +} + +func inRange(a, b, c float64) bool { + if a > b { + a, b = b, a + } + return c >= a && c <= b +} + +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 + edge := NewPlane2D(t[i].X, t[i].Y, t[j].X, t[j].Y) + + s1 := edge.Eval(vl.x1, vl.y1) + s2 := edge.Eval(vl.x2, vl.y2) + + o1 := onPlane(s1) + o2 := onPlane(s2) + + switch { + case o1 && o2: + // XXX: Broken! + pos := interpolate(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) + + switch { + case r1 && r2: + l := linear(t[i], t[j]) + // XXX: Broken! + out = append(out, l(t1/length), l(t2/length)) + return out + + case !r1 && !r2: // the hole edge + out = append(out, t[i], t[j]) + return out + case !r1: + // TODO + case !r2: + // TODO + } + + case o1: + t1 := interpolate(t[i], t[j])(vl.x1, vl.y1) + if !inRange(0, length, 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) { + continue edges + } + out = append(out, linear(t[i], t[j])(t2)) + + default: + x, y, intersects := vl.line.Intersection(edge) + if !intersects { + continue edges + } + t1 := interpolate(t[i], t[j])(x, y) + if !inRange(0, length, t1) { + continue edges + } + + // XXX: Ugly! + z := (t[j].Z-t[i].Z)*(t1/length) + t[1].Z + n := Vertex{x, y, z} + + if math.Signbit(s1) != math.Signbit(s2) { + // different sides -> vertex on edge. + out = append(out) + } else { // same side -> inside. Find peer + if len(out) > 0 { // we already have our peer + out = append(out, n) + return out + } + + for k := 0; k < 3; k++ { + if k == i { + continue + } + l := (k + 1) % 3 + other := NewPlane2D(t[k].X, t[k].Y, t[l].X, t[l].Y) + xo, yo, intersects := vl.line.Intersection(edge) + if !intersects { + continue + } + // XXX: Broken!!! + zo := interpolate(t[k], t[j])(xo, yo) + if !inRange(t[k].Z, t[k].Z, zo) { + continue + } + m := Vertex{xo, yo, zo} + + var xn, yn, xf, yf float64 + + // Which is nearer to current edge? + if math.Abs(s1) < math.Abs(s2) { + xn, yn = vl.x1, vl.y1 + xf, yf = vl.x2, vl.y2 + } else { + xn, yn = vl.x2, vl.y2 + xf, yf = vl.x1, vl.y1 + } + + if onPlane(other.Eval(xn, yn)) { + // triangle intersect in only point + // treat as no intersection + return out + } + + // XXX: Broken!!! + zint := interpolate(n, m) + + zn := zint(xn, yn) + zf := zint(xf, yf) + + if !inRange(n.Z, m.Z, zn) { + // if near is out of range far is, too. + return out + } + + if inRange(n.Z, m.Z, zf) { + m.X = xf + m.Y = yf + m.Z = zf + } // else m is clipping + + n.X = xn + n.Y = yn + n.Z = zn + + out = append(out, n, m) + return out + } + } + } + } + + // supress single point touches. + if len(out) == 1 { + out = out[:0] + } + + return out +}