Mercurial > gemma
view pkg/octree/vertex.go @ 4488:bff6c5c1db4f
client: pdf-gen: improve adding bottleneck info to pdf
* Check if the bottleneck is in the current view to add its info to the exported pdf and the pdf filename, this avoid wrong filename and wrong info in pdf in case view has been changed to another location.
* Set the bottleneck to print after moving to it in map.
author | Fadi Abbud <fadi.abbud@intevation.de> |
---|---|
date | Fri, 27 Sep 2019 11:15:02 +0200 |
parents | 4233570de212 |
children | 56f4e8cbfab7 |
line wrap: on
line source
// This is Free Software under GNU Affero General Public License v >= 3.0 // without warranty, see README.md and license for details. // // SPDX-License-Identifier: AGPL-3.0-or-later // License-Filename: LICENSES/AGPL-3.0.txt // // Copyright (C) 2018 by via donau // – Österreichische Wasserstraßen-Gesellschaft mbH // Software engineering by Intevation GmbH // // Author(s): // * Sascha L. Teichmann <sascha.teichmann@intevation.de> package octree import ( "bytes" "encoding/binary" "fmt" "io" "log" "math" "sort" "gemma.intevation.de/gemma/pkg/wkb" ) type ( Point struct { X float64 Y float64 } // Vertex is a 3D vertex. Vertex struct { X float64 Y float64 Z float64 } // Triangle is a triangle consisting of three vertices. Triangle [3]Vertex // Line is a line defined by first vertex on that line // and the second being the direction. Line [2]Vertex // Box is a 3D box. Box [2]Vertex // MultiPointZ is a set of vertices. MultiPointZ []Vertex // LineStringZ is a line string formed of vertices. LineStringZ []Vertex // MultiLineStringZ is a set of line strings. MultiLineStringZ []LineStringZ // Box2D is 2D area from (X1, Y1) to (X2, Y2). Box2D struct { X1 float64 Y1 float64 X2 float64 Y2 float64 } // Plane2D is a 2D plane (a line in 2D space). Plane2D struct { A float64 B float64 C float64 } Plane3D struct { A float64 B float64 C float64 D float64 } ) func (t *Triangle) Plane3D() Plane3D { v0 := t[1].Sub(t[0]) v1 := t[2].Sub(t[0]) n := v0.Cross(v1).Normalize() // x*nx+ y*ny+ z*nz + d = 0 // d = - (x*nx+ y*ny + z*nz) d := -t[0].Dot(n) return Plane3D{ A: n.X, B: n.Y, C: n.Z, D: d, } } func (t *Triangle) BBox() Box2D { minX := math.Min(math.Min(t[0].X, t[1].X), t[2].X) maxX := math.Max(math.Max(t[0].X, t[1].X), t[2].X) minY := math.Min(math.Min(t[0].Y, t[1].Y), t[2].Y) maxY := math.Max(math.Max(t[0].Y, t[1].Y), t[2].Y) return Box2D{ X1: minX, Y1: minY, X2: maxX, Y2: maxY, } } func (p Plane3D) Z(x, y float64) float64 { // p.A*x + p.B*y + p.C*z + p.D = 0 return -(p.A*x + p.B*y + p.D) / p.C } func (p Plane3D) Eval(v Vertex) float64 { return p.A*v.X + p.B*v.Y + p.C*v.Z + p.D } func (v Vertex) Normalize() Vertex { s := 1 / v.Length() return Vertex{ X: s * v.X, Y: s * v.Y, Z: s * v.Z, } } func (v Vertex) Dot(w Vertex) float64 { return v.X*w.X + v.Y*w.Y + v.Z*w.Z } func (v Vertex) Length() float64 { return math.Sqrt(v.Dot(v)) } func area(a, b, c Vertex) float64 { return (b.Y-a.Y)*(c.X-b.X) - (b.X-a.X)*(c.Y-b.Y) } func inCircle(a, b, c, p Vertex) bool { dx := a.X - p.X dy := a.Y - p.Y ex := b.X - p.X ey := b.Y - p.Y fx := c.X - p.X fy := c.Y - p.Y ap := dx*dx + dy*dy bp := ex*ex + ey*ey cp := fx*fx + fy*fy return dx*(ey*cp-bp*fy)-dy*(ex*cp-bp*fx)+ap*(ex*fy-ey*fx) < 0 } func circumradius(a, b, c Vertex) float64 { dx := b.X - a.X dy := b.Y - a.Y ex := c.X - a.X ey := c.Y - a.Y bl := dx*dx + dy*dy cl := ex*ex + ey*ey d := dx*ey - dy*ex x := (ey*bl - dy*cl) * 0.5 / d y := (dx*cl - ex*bl) * 0.5 / d r := x*x + y*y if bl == 0 || cl == 0 || d == 0 || r == 0 { return infinity } return r } func circumcenter(a, b, c Vertex) Vertex { dx := b.X - a.X dy := b.Y - a.Y ex := c.X - a.X ey := c.Y - a.Y bl := dx*dx + dy*dy cl := ex*ex + ey*ey d := dx*ey - dy*ex x := a.X + (ey*bl-dy*cl)*0.5/d y := a.Y + (dx*cl-ex*bl)*0.5/d return Vertex{X: x, Y: y} } func polygonArea(points []Vertex) float64 { var result float64 for i, p := range points { q := points[(i+1)%len(points)] result += (p.X - q.X) * (p.Y + q.Y) } return result / 2 } func polygonPerimeter(points []Vertex) float64 { if len(points) == 0 { return 0 } var result float64 q := points[len(points)-1] for _, p := range points { result += p.Distance2D(q) q = p } return result } func (v Vertex) Distance2D(w Vertex) float64 { return math.Hypot(v.X-w.X, v.Y-w.Y) } func (v Vertex) Distance(w Vertex) float64 { v = v.Sub(w) return math.Sqrt(v.Dot(v)) } // Minimize adjust this vertex v to hold the minimum // values component-wise of v and w. func (v *Vertex) Minimize(w Vertex) { if w.X < v.X { v.X = w.X } if w.Y < v.Y { v.Y = w.Y } if w.Z < v.Z { v.Z = w.Z } } // Maximize adjust this vertex v to hold the maximum // values component-wise of v and w. func (v *Vertex) Maximize(w Vertex) { if w.X > v.X { v.X = w.X } if w.Y > v.Y { v.Y = w.Y } if w.Z > v.Z { v.Z = w.Z } } func (v Vertex) SquaredDistance2D(w Vertex) float64 { dx := v.X - w.X dy := v.Y - w.Y return dx*dx + dy*dy } // Sub2D returns (v - w) component-wise. func (v Vertex) Sub2D(w Vertex) Vertex { return Vertex{ X: v.X - w.X, Y: v.Y - w.Y, } } // Sub returns (v - w) component-wise. func (v Vertex) Sub(w Vertex) Vertex { return Vertex{ v.X - w.X, v.Y - w.Y, v.Z - w.Z, } } // Add returns (v + w) component-wise. func (v Vertex) Add(w Vertex) Vertex { return Vertex{ v.X + w.X, v.Y + w.Y, v.Z + w.Z, } } // Scale returns s*v component-wise. func (v Vertex) Scale(s float64) Vertex { return Vertex{ s * v.X, s * v.Y, s * v.Z, } } // Interpolate returns a function that return s*b[1] + b[0] // component-wise. func (b Box) Interpolate() func(Vertex) Vertex { v1, v2 := b[0], b[1] v2 = v2.Sub(v1) return func(s Vertex) Vertex { return Vertex{ v2.X*s.X + v1.X, v2.Y*s.Y + v1.Y, v2.Z*s.Z + v1.Z, } } } func (b Box) HasX() bool { return math.Abs(b[0].X-b[1].X) > epsPlane } func (b Box) HasY() bool { return math.Abs(b[0].Y-b[1].Y) > epsPlane } func (b Box) HasZ() bool { return math.Abs(b[0].Z-b[1].Z) > epsPlane } // Less returns if one of v component is less than the // corresponing component in w. func (v Vertex) Less(w Vertex) bool { return v.X < w.X || v.Y < w.Y || v.Z < w.Z } // NewLine return a line of point/direction. func NewLine(p1, p2 Vertex) Line { return Line{ p2.Sub(p1), p1, } } // Eval returns the vertex for t*l[0] + l[1]. func (l Line) Eval(t float64) Vertex { return l[0].Scale(t).Add(l[1]) } // IntersectHorizontal returns the intersection point // for a given z value. func (l Line) IntersectHorizontal(h float64) Vertex { t := (h - l[1].Z) / l[0].Z return l.Eval(t) } func side(z, h float64) int { switch { case z < h: return -1 case z > h: return +1 } return 0 } func (v Vertex) Dot2(w Vertex) float64 { return v.X*w.X + v.Y*w.Y } func (t *Triangle) Contains(x, y float64) bool { v0 := t[2].Sub(t[0]) v1 := t[1].Sub(t[0]) v2 := Vertex{X: x, Y: y}.Sub(t[0]) dot00 := v0.Dot2(v0) dot01 := v0.Dot2(v1) dot02 := v0.Dot2(v2) dot11 := v1.Dot2(v1) dot12 := v1.Dot2(v2) // Compute barycentric coordinates invDenom := 1 / (dot00*dot11 - dot01*dot01) u := (dot11*dot02 - dot01*dot12) * invDenom v := (dot00*dot12 - dot01*dot02) * invDenom // Check if point is in triangle return u >= 0 && v >= 0 && u+v < 1 } // IntersectHorizontal calculates the line string that // results when cutting a triangle a a certain height. // Can be empty (on intersection), // one vertex (only touching an vertex) or // two vertices (real intersection). func (t *Triangle) IntersectHorizontal(h float64) LineStringZ { sides := [3]int{ side(t[0].Z, h), side(t[1].Z, h), side(t[2].Z, h), } var points LineStringZ for i := 0; i < 3; i++ { j := (i + 1) % 3 si := sides[i] sj := sides[j] switch { case si == 0 && sj == 0: // both on plane points = append(points, t[i], t[j]) case si == 0 && sj != 0: // first on plane points = append(points, t[i]) case si != 0 && sj == 0: // second on plane points = append(points, t[j]) case si == sj: // both on same side default: // real intersection v := NewLine(t[i], t[j]).IntersectHorizontal(h) points = append(points, v) } } return points } func linearScale(x1, y1, x2, y2 float64) func(Vertex) float64 { dx := x2 - x1 dy := y2 - y1 switch { case dx != 0: return func(v Vertex) float64 { return (v.X - x1) / dx } case dy != 0: return func(v Vertex) float64 { return (v.Y - y1) / dy } default: return func(Vertex) float64 { return 0 } } } func (ls LineStringZ) BBox() Box2D { min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} for _, v := range ls { min.Minimize(v) max.Maximize(v) } return Box2D{ X1: min.X, Y1: min.Y, X2: max.X, Y2: max.Y, } } func (ls LineStringZ) Area() float64 { return polygonArea(ls) } func (ls LineStringZ) Reverse() { for i, j := 0, len(ls)-1; i < j; i, j = i+1, j-1 { ls[i], ls[j] = ls[j], ls[i] } } func (ls LineStringZ) order(position func(Vertex) float64) { type posVertex struct { pos float64 v Vertex } positions := make([]posVertex, len(ls)) for i, v := range ls { positions[i] = posVertex{position(v), v} } sort.Slice(positions, func(i, j int) bool { return positions[i].pos < positions[j].pos }) for i := range positions { ls[i] = positions[i].v } } // EpsEquals returns true if v and w are equal component-wise // with the values within a epsilon range. func (v Vertex) EpsEquals(w Vertex) bool { const eps = 1e-5 return math.Abs(v.X-w.X) < eps && math.Abs(v.Y-w.Y) < eps && math.Abs(v.Z-w.Z) < eps } // JoinOnLine joins the the elements of a given multi line string // under the assumption that the segments are all on the line segment // from (x1, y1) to (x2, y2). func (mls MultiLineStringZ) JoinOnLine(x1, y1, x2, y2 float64) MultiLineStringZ { position := linearScale(x1, y1, x2, y2) type posLineString struct { pos float64 line LineStringZ } positions := make([]posLineString, 0, len(mls)) for _, ls := range mls { if len(ls) == 0 { continue } // order the atoms first ls.order(position) positions = append(positions, posLineString{position(ls[0]), ls}) } sort.Slice(positions, func(i, j int) bool { return positions[i].pos < positions[j].pos }) out := make(MultiLineStringZ, 0, len(positions)) var ignored int for i := range positions { curr := positions[i].line if l := len(out); l > 0 { last := out[l-1] if last[len(last)-1].EpsEquals(curr[0]) { out[l-1] = append(last[:len(last)-1], curr...) continue } if position(last[len(last)-1]) > position(curr[0]) { ignored++ continue } } out = append(out, curr) } log.Printf("info: ignored parts: %d\n", ignored) return out } // Write writes a Vertex as three 64 bit values in little endian order // to the given writer. func (v *Vertex) Write(w io.Writer) error { if err := binary.Write( w, binary.LittleEndian, math.Float64bits(v.X)); err != nil { return err } if err := binary.Write( w, binary.LittleEndian, math.Float64bits(v.Y)); err != nil { return err } return binary.Write( w, binary.LittleEndian, math.Float64bits(v.Z)) } // Read fills this vertex with three 64 bit values stored as // little endian from the given reader. func (v *Vertex) Read(r io.Reader) error { var buf [8]byte b := buf[:] if _, err := io.ReadFull(r, b); err != nil { return nil } v.X = math.Float64frombits(binary.LittleEndian.Uint64(b)) if _, err := io.ReadFull(r, b); err != nil { return nil } v.Y = math.Float64frombits(binary.LittleEndian.Uint64(b)) if _, err := io.ReadFull(r, b); err != nil { return nil } v.Z = math.Float64frombits(binary.LittleEndian.Uint64(b)) return nil } // AsWKB returns the WKB representation of the given multi line string. func (mls MultiLineStringZ) AsWKB() []byte { // pre-calculate size to avoid reallocations. size := 1 + 4 + 4 for _, ml := range mls { size += 1 + 4 + 4 + len(ml)*3*8 } buf := bytes.NewBuffer(make([]byte, 0, size)) binary.Write(buf, binary.LittleEndian, wkb.NDR) binary.Write(buf, binary.LittleEndian, wkb.MultiLineStringZ) binary.Write(buf, binary.LittleEndian, uint32(len(mls))) for _, ml := range mls { binary.Write(buf, binary.LittleEndian, wkb.NDR) binary.Write(buf, binary.LittleEndian, wkb.LineStringZ) binary.Write(buf, binary.LittleEndian, uint32(len(ml))) for _, p := range ml { binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Z)) } } return buf.Bytes() } // AsWKB2D returns the WKB representation of the given multi line string // leaving the z component out. func (mls MultiLineStringZ) AsWKB2D() []byte { // pre-calculate size to avoid reallocations. size := 1 + 4 + 4 for _, ml := range mls { size += 1 + 4 + 4 + len(ml)*2*8 } buf := bytes.NewBuffer(make([]byte, 0, size)) binary.Write(buf, binary.LittleEndian, wkb.NDR) binary.Write(buf, binary.LittleEndian, wkb.MultiLineString) binary.Write(buf, binary.LittleEndian, uint32(len(mls))) for _, ml := range mls { binary.Write(buf, binary.LittleEndian, wkb.NDR) binary.Write(buf, binary.LittleEndian, wkb.LineString) binary.Write(buf, binary.LittleEndian, uint32(len(ml))) for _, p := range ml { binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) } } return buf.Bytes() } // Join joins two lines leaving the first of the second out. func (ls LineStringZ) Join(other LineStringZ) LineStringZ { nline := make(LineStringZ, len(ls)+len(other)-1) copy(nline, ls) copy(nline[len(ls):], other[1:]) return nline } // Merge merges line segments of a given multi line string // by finding common start and end vertices. func (mls MultiLineStringZ) Merge() MultiLineStringZ { var out MultiLineStringZ min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} for _, line := range mls { for _, v := range line { min.Minimize(v) } } type point struct { x int64 y int64 } const precision = 1e7 quant := func(v Vertex) point { return point{ x: int64(math.Round((v.X - min.X) * precision)), y: int64(math.Round((v.Y - min.Y) * precision)), } } heads := make(map[point]*[]LineStringZ) for _, line := range mls { if len(line) < 2 { out = append(out, line) continue } head := quant(line[0]) tail := quant(line[len(line)-1]) if head == tail { // its already a ring out = append(out, line) continue } if hs := heads[tail]; hs != nil { l := len(*hs) last := (*hs)[l-1] if l == 1 { delete(heads, tail) } else { (*hs)[l-1] = nil *hs = (*hs)[:l-1] } line = line.Join(last) if head == quant(line[len(line)-1]) { // its a ring out = append(out, line) continue } } if hs := heads[head]; hs != nil { *hs = append(*hs, line) } else { heads[head] = &[]LineStringZ{line} } } again: for head, lines := range heads { for i, line := range *lines { tail := quant(line[len(line)-1]) for hs := heads[tail]; hs != nil && len(*hs) > 0; hs = heads[tail] { l := len(*hs) last := (*hs)[l-1] (*hs)[l-1] = nil *hs = (*hs)[:l-1] line = line.Join(last) if tail = quant(line[len(line)-1]); head == tail { // its a ring out = append(out, line) // remove from current lines copy((*lines)[i:], (*lines)[i+1:]) (*lines)[len(*lines)-1] = nil *lines = (*lines)[:len(*lines)-1] goto again } // overwrite in current lines (*lines)[i] = line } } } // rings := len(out) // The rest are open line strings. for _, lines := range heads { for _, line := range *lines { out = append(out, line) } } // log.Printf("segments before/after merge: %d/%d (%d rings)\n", // len(mls), len(out), rings) return out } func (a Box2D) Rect(interface{}) ([]float64, []float64) { return []float64{a.X1, a.Y1}, []float64{a.X2, a.Y2} } // Intersects checks if two Box2Ds intersect. func (a Box2D) Intersects(b Box2D) bool { return !(a.X2 < a.X1 || a.X2 < b.X1 || a.Y2 < a.Y1 || a.Y2 < b.Y1) } func (a Box2D) Contains(x, y float64) bool { return a.X1 <= x && x <= a.X2 && a.Y1 <= y && y <= a.Y2 } // Xi returns the i-th x component. func (a Box2D) Xi(i int) float64 { if i == 0 { return a.X1 } return a.X2 } // Yi returns the i-th y component. func (a Box2D) Yi(i int) float64 { if i == 0 { return a.Y1 } return a.Y2 } func (a Box2D) Union(b Box2D) Box2D { return Box2D{ X1: math.Min(a.X1, b.X1), Y1: math.Min(a.Y1, b.Y1), X2: math.Max(a.X2, b.X2), Y2: math.Max(a.Y2, b.Y2), } } func (a Box2D) Area() float64 { return (a.X2 - a.X1) * (a.Y2 - a.Y1) } // NewPlane2D creates a new Plane2D from two given points. func NewPlane2D(x1, y1, x2, y2 float64) Plane2D { b := x2 - x1 a := -(y2 - y1) l := math.Sqrt(a*a + b*b) a /= l b /= l // a*x1 + b*y1 + c = 0 c := -(a*x1 + b*y1) return Plane2D{a, b, c} } // Eval determines the distance of a given point // from the plane. The sign of the result indicates // the sideness. func (p Plane2D) Eval(x, y float64) float64 { return p.A*x + p.B*y + p.C } const epsPlane = 1e-5 func sides(s int, x float64) int { if math.Signbit(x) { return s | 2 } return s | 1 } // IntersectsPlane checks if a Box2D intersects with // a given Plane2D. func (a Box2D) IntersectsPlane(p Plane2D) bool { var s int for i := 0; i < 2; i++ { x := a.Xi(i) for j := 0; j < 2; j++ { y := a.Yi(j) v := p.Eval(x, y) if math.Abs(v) < epsPlane { //log.Printf("on line") return true } if s = sides(s, v); s == 3 { //log.Printf("... on both sides (djt)") return true } } } //log.Printf("side: %d\n", s) return false } // Cross calculates the cross product of two vertices. func (v Vertex) Cross(w Vertex) Vertex { return Vertex{ v.Y*w.Z - v.Z*w.Y, v.Z*w.X - v.X*w.Z, v.X*w.Y - v.Y*w.X, } } // Intersection calcultes the 2D intersection point of // two Plane2Ds. If they do not intersect the returned // bool flags is set to false. func (p Plane2D) Intersection(o Plane2D) (float64, float64, bool) { u1 := Vertex{p.A, p.B, p.C} u2 := Vertex{o.A, o.B, o.C} plane := u1.Cross(u2) if plane.Z == 0 { return 0, 0, false } return plane.X / plane.Z, plane.Y / plane.Z, true } // VerticalLine is a 2D line segment. type VerticalLine struct { x1 float64 y1 float64 x2 float64 y2 float64 line Plane2D } // NewVerticalLine creates a new 2D line segment. 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 relative(v1, v2 Vertex) func(x, y float64) float64 { ls := linearScale(v1.X, v1.Y, v2.X, v2.Y) return func(x, y float64) float64 { return ls(Vertex{x, y, 0}) } } 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{ (v2.X-v1.X)*t + v1.X, (v2.Y-v1.Y)*t + v1.Y, (v2.Z-v1.Z)*t + v1.Z, } } } func inRange(a float64) bool { return 0 <= a && a <= 1 } // Intersection intersects a line segment with a triangle. func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ { var out LineStringZ //defer func() { log.Printf("length out: %d\n", len(out)) }() 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) // log.Printf("s1, s2: %t %t (%f %f)\n", o1, o2, s1, s2) switch { case o1 && o2: pos := relative(t[i], t[j]) t1 := pos(vl.x1, vl.y1) t2 := pos(vl.x2, vl.y2) r1 := inRange(t1) r2 := inRange(t2) switch { case r1 && r2: lin := linear(t[i], t[j]) out = append(out, lin(t1), lin(t2)) return out case !r1 && !r2: // the hole edge out = append(out, t[i], t[j]) return out case !r1: if t1 < 0 { // below first -> clip by first lin := linear(t[i], t[j]) out = append(out, t[i], lin(t2)) } else { // above second -> clip by second lin := linear(t[i], t[j]) out = append(out, lin(t2), t[j]) } return out case !r2: if t2 < 0 { // below first -> clip by first lin := linear(t[i], t[j]) out = append(out, t[i], lin(t1)) } else { // above second -> clip by second lin := linear(t[i], t[j]) out = append(out, lin(t1), t[j]) } return out } case o1: 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 := relative(t[i], t[j])(vl.x2, vl.y2) if !inRange(t2) { continue edges } out = append(out, linear(t[i], t[j])(t2)) default: x, y, intersects := vl.line.Intersection(edge) if !intersects { continue edges } // log.Println("Intersection -----------------------------") t1 := relative(t[i], t[j])(x, y) // log.Printf("relative pos: %f\n", t1) if !inRange(t1) { continue edges } // log.Println("valid ***************") z := interpolate(t[j].Z, t[i].Z)(t1) n := Vertex{x, y, z} if math.Signbit(s1) != math.Signbit(s2) { // log.Println("\ton different sides") // different sides -> vertex on edge. out = append(out, n) } 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(other) if !intersects { continue } t2 := relative(t[k], t[l])(xo, yo) if !inRange(t2) { continue } zo := interpolate(t[k].Z, t[l].Z)(t2) 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 break edges } pos := relative(n, m) tzn := pos(xn, yn) tzf := pos(xf, yf) if !inRange(tzn) { // if near is out of range far is, too. return out } lin := interpolate(n.Z, m.Z) if inRange(tzf) { m.X = xf m.Y = yf m.Z = lin(tzf) } // else m is clipping n.X = xn n.Y = yn n.Z = lin(tzn) out = append(out, n, m) return out } } } } // supress single point touches. if len(out) == 1 { out = out[:0] } return out } // AsWKB returns a WKB representation of the given point cloud. func (mpz MultiPointZ) AsWKB() []byte { size := 1 + 4 + 4 + len(mpz)*(1+4+3*8) buf := bytes.NewBuffer(make([]byte, 0, size)) binary.Write(buf, binary.LittleEndian, wkb.NDR) binary.Write(buf, binary.LittleEndian, wkb.MultiPointZ) binary.Write(buf, binary.LittleEndian, uint32(len(mpz))) perPoint := bytes.NewBuffer(make([]byte, 0, 1+4)) binary.Write(perPoint, binary.LittleEndian, wkb.NDR) binary.Write(perPoint, binary.LittleEndian, wkb.PointZ) hdr := perPoint.Bytes() for _, p := range mpz { buf.Write(hdr) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Z)) } return buf.Bytes() } func (mpz *MultiPointZ) FromWKB(data []byte) error { r := bytes.NewReader(data) endian, err := r.ReadByte() var order binary.ByteOrder switch { case err != nil: return err case endian == wkb.NDR: order = binary.LittleEndian case endian == wkb.XDR: order = binary.BigEndian default: return fmt.Errorf("unknown byte order %x", endian) } var geomType uint32 err = binary.Read(r, order, &geomType) switch { case err != nil: return err case geomType != wkb.MultiPointZ: return fmt.Errorf("unknown geometry type %x", geomType) } var numPoints uint32 if err := binary.Read(r, order, &numPoints); err != nil { return err } points := make(MultiPointZ, numPoints) for i := range points { endian, err = r.ReadByte() switch { case err != nil: return err case endian == wkb.NDR: order = binary.LittleEndian case endian == wkb.XDR: order = binary.BigEndian default: return fmt.Errorf("unknown byte order %x", endian) } err = binary.Read(r, order, &geomType) switch { case err != nil: return err case geomType != wkb.PointZ: return fmt.Errorf("unknown geometry type %x", geomType) } var x, y, z uint64 if err = binary.Read(r, order, &x); err != nil { return err } if err = binary.Read(r, order, &y); err != nil { return err } if err = binary.Read(r, order, &z); err != nil { return err } points[i] = Vertex{ X: math.Float64frombits(x), Y: math.Float64frombits(y), Z: math.Float64frombits(z), } } *mpz = points return nil }