Mercurial > gemma

changeset 727:41c8dc61f38f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Moved octree loading stuff to octree package.
author Sascha L. Teichmann <teichmann@intevation.de>
date Sat, 22 Sep 2018 21:57:30 +0200
parents 5af9ab39e715
children e9d37300ce99
files cmd/octree2contour/loader.go cmd/octree2contour/main.go cmd/octree2contour/octree.go cmd/octree2contour/store.go cmd/octree2contour/vertex.go pkg/octree/cache.go pkg/octree/loader.go pkg/octree/octree.go pkg/octree/tree.go pkg/octree/vertex.go
diffstat 10 files changed, 410 insertions(+), 409 deletions(-) [+]
line wrap: on
line diff
--- a/cmd/octree2contour/loader.go	Sat Sep 22 21:34:12 2018 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,109 +0,0 @@
-package main
-
-import (
-	"bufio"
-	"encoding/binary"
-	"log"
-	"os"
-
-	"github.com/golang/snappy"
-)
-
-func loadOctreeReader(r *bufio.Reader) (*Octree, error) {
-	tree := new(Octree)
-
-	if err := binary.Read(r, binary.LittleEndian, &tree.EPSG); err != nil {
-		return nil, err
-	}
-
-	log.Printf("EPSG: %d\n", tree.EPSG)
-
-	if err := tree.Min.read(r); err != nil {
-		return nil, err
-	}
-
-	if err := tree.Max.read(r); err != nil {
-		return nil, err
-	}
-
-	log.Printf("BBOX: [[%f, %f, %f], [%f, %f, %f]]\n",
-		tree.Min.x, tree.Min.y, tree.Min.z,
-		tree.Max.x, tree.Max.y, tree.Max.z)
-
-	var numVertices uint32
-	if err := binary.Read(r, binary.LittleEndian, &numVertices); err != nil {
-		return nil, err
-	}
-
-	log.Printf("vertices: %d\n", numVertices)
-
-	vertices := make([]Vertex, numVertices)
-	tree.vertices = vertices
-
-	for i := range vertices {
-		if err := vertices[i].read(r); err != nil {
-			return nil, err
-		}
-	}
-
-	var numTriangles uint32
-	if err := binary.Read(r, binary.LittleEndian, &numTriangles); err != nil {
-		return nil, err
-	}
-
-	log.Printf("triangles: %d\n", numTriangles)
-
-	indices := make([]int32, 3*numTriangles)
-	triangles := make([][]int32, numTriangles)
-	tree.triangles = triangles
-
-	var last int32
-
-	for i := range triangles {
-		tri := indices[:3]
-		indices = indices[3:]
-		triangles[i] = tri
-		for j := range tri {
-			v, err := binary.ReadVarint(r)
-			if err != nil {
-				return nil, err
-			}
-			value := int32(v) + last
-			tri[j] = value
-			last = value
-		}
-	}
-
-	var numNodes uint32
-	if err := binary.Read(r, binary.LittleEndian, &numNodes); err != nil {
-		return nil, err
-	}
-
-	log.Printf("num nodes: %d\n", numNodes)
-
-	tree.index = make([]int32, numNodes)
-	entries := tree.index[1:]
-
-	last = 0
-	for i := range entries {
-		v, err := binary.ReadVarint(r)
-		if err != nil {
-			return nil, err
-		}
-		value := int32(v) + last
-		entries[i] = value
-		last = value
-	}
-
-	return tree, nil
-}
-
-func LoadOctree(fname string) (*Octree, error) {
-
-	f, err := os.Open(fname)
-	if err != nil {
-		return nil, err
-	}
-	defer f.Close()
-	return loadOctreeReader(bufio.NewReader(snappy.NewReader(f)))
-}
--- a/cmd/octree2contour/main.go	Sat Sep 22 21:34:12 2018 +0200
+++ b/cmd/octree2contour/main.go	Sat Sep 22 21:57:30 2018 +0200
@@ -7,6 +7,8 @@
 	"sort"
 	"sync"
 	"time"
+
+	"gemma.intevation.de/gemma/pkg/octree"
 )
 
 var (
@@ -16,15 +18,15 @@
 )
 
 func processLevels(
-	tree *Octree,
+	tree *octree.Tree,
 	jobs <-chan float64,
 	results chan<- result,
 	wg *sync.WaitGroup,
 ) {
 	defer wg.Done()
 	for h := range jobs {
-		var lines MultiLineStringZ
-		tree.Horizontal(h, func(t *Triangle) {
+		var lines octree.MultiLineStringZ
+		tree.Horizontal(h, func(t *octree.Triangle) {
 			line := t.IntersectHorizontal(h)
 			if len(line) > 1 {
 				lines = append(lines, line)
@@ -34,11 +36,11 @@
 	}
 }
 
-func process(tree *Octree) []result {
+func process(tree *octree.Tree) []result {
 
 	if *one {
-		var lines MultiLineStringZ
-		tree.Horizontal(*step, func(t *Triangle) {
+		var lines octree.MultiLineStringZ
+		tree.Horizontal(*step, func(t *octree.Triangle) {
 			line := t.IntersectHorizontal(*step)
 			if len(line) > 0 {
 				lines = append(lines, line)
@@ -68,7 +70,7 @@
 		wg.Add(1)
 		go processLevels(tree, jobs, results, &wg)
 	}
-	for h := tree.Min.z; h <= tree.Max.z; h += *step {
+	for h := tree.Min.Z; h <= tree.Max.Z; h += *step {
 		jobs <- h
 	}
 	close(jobs)
@@ -88,7 +90,7 @@
 	for _, fname := range flag.Args() {
 		log.Printf("processing %s\n", fname)
 		start := time.Now()
-		tree, err := LoadOctree(fname)
+		tree, err := octree.LoadTree(fname)
 		if err != nil {
 			log.Printf("error: %v\n", err)
 			continue
--- a/cmd/octree2contour/octree.go	Sat Sep 22 21:34:12 2018 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-package main
-
-import (
-	"math"
-)
-
-type Octree struct {
-	EPSG uint32
-
-	vertices  []Vertex
-	triangles [][]int32
-	index     []int32
-
-	Min Vertex
-	Max Vertex
-}
-
-func (ot *Octree) Horizontal(h float64, fn func(*Triangle)) {
-
-	type frame struct {
-		pos int32
-		min float64
-		max float64
-	}
-
-	if h < ot.Min.z || ot.Max.z < h {
-		return
-	}
-
-	stack := []frame{{1, ot.Min.z, ot.Max.z}}
-
-	dupes := map[int32]struct{}{}
-
-	for len(stack) > 0 {
-		top := stack[len(stack)-1]
-		stack = stack[:len(stack)-1]
-
-		pos := top.pos
-		if pos == 0 {
-			continue
-		}
-		min, max := top.min, top.max
-
-		if pos > 0 { // node
-			if mid := (max-min)*0.5 + min; h >= mid {
-				pos += 4 // nodes with z-bit set
-				min = mid
-			} else {
-				max = mid
-			}
-			stack = append(stack,
-				frame{ot.index[pos+0], min, max},
-				frame{ot.index[pos+1], min, max},
-				frame{ot.index[pos+2], min, max},
-				frame{ot.index[pos+3], min, max})
-		} else { // leaf
-			pos = -pos - 1
-			n := ot.index[pos]
-			//log.Printf("%d %d %d\n", pos, n, len(ot.index))
-			indices := ot.index[pos+1 : pos+1+n]
-
-			for _, idx := range indices {
-				if _, found := dupes[idx]; found {
-					continue
-				}
-				tri := ot.triangles[idx]
-				t := Triangle{
-					ot.vertices[tri[0]],
-					ot.vertices[tri[1]],
-					ot.vertices[tri[2]],
-				}
-
-				if !(math.Min(t[0].z, math.Min(t[1].z, t[2].z)) > h ||
-					math.Max(t[0].z, math.Max(t[1].z, t[2].z)) < h) {
-					dupes[idx] = struct{}{}
-					fn(&t)
-				}
-			}
-		}
-	}
-}
--- a/cmd/octree2contour/store.go	Sat Sep 22 21:34:12 2018 +0200
+++ b/cmd/octree2contour/store.go	Sat Sep 22 21:57:30 2018 +0200
@@ -1,22 +1,14 @@
 package main
 
 import (
-	"bytes"
 	"database/sql"
-	"encoding/binary"
-	"math"
-)
 
-const (
-	wkbNDR              byte   = 1
-	wkbPointZ           uint32 = 1000 + 1
-	wkbLineStringZ      uint32 = 1000 + 2
-	wkbMultiLineStringZ uint32 = 1000 + 5
+	"gemma.intevation.de/gemma/pkg/octree"
 )
 
 type result struct {
 	h     float64
-	lines MultiLineStringZ
+	lines octree.MultiLineStringZ
 }
 
 const insertSQL = `
@@ -41,7 +33,7 @@
 		}
 
 		for _, r := range all {
-			if _, err := stmt.Exec(r.h, r.lines.asWKB(), epsg); err != nil {
+			if _, err := stmt.Exec(r.h, r.lines.AsWKB(), epsg); err != nil {
 				return err
 			}
 		}
@@ -49,25 +41,3 @@
 		return tx.Commit()
 	})
 }
-
-func (mls MultiLineStringZ) asWKB() []byte {
-
-	var buf bytes.Buffer
-
-	binary.Write(&buf, binary.LittleEndian, wkbNDR)
-	binary.Write(&buf, binary.LittleEndian, wkbMultiLineStringZ)
-	binary.Write(&buf, binary.LittleEndian, uint32(len(mls)))
-
-	for _, ml := range mls {
-		binary.Write(&buf, binary.LittleEndian, wkbNDR)
-		binary.Write(&buf, binary.LittleEndian, wkbLineStringZ)
-		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()
-}
--- a/cmd/octree2contour/vertex.go	Sat Sep 22 21:34:12 2018 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,165 +0,0 @@
-package main
-
-import (
-	"encoding/binary"
-	"io"
-	"math"
-)
-
-type (
-	Vertex struct {
-		x float64
-		y float64
-		z float64
-	}
-
-	Triangle [3]Vertex
-
-	Line [2]Vertex
-
-	LineStringZ      []Vertex
-	MultiLineStringZ []LineStringZ
-)
-
-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
-	}
-}
-
-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) Sub(w Vertex) Vertex {
-	return Vertex{
-		v.x - w.x,
-		v.y - w.y,
-		v.z - w.z,
-	}
-}
-
-func (v Vertex) Add(w Vertex) Vertex {
-	return Vertex{
-		v.x + w.x,
-		v.y + w.y,
-		v.z + w.z,
-	}
-}
-
-func (v Vertex) scale(s float64) Vertex {
-	return Vertex{
-		s * v.x,
-		s * v.y,
-		s * v.z,
-	}
-}
-
-func Interpolate(v1, v2 Vertex) func(Vertex) Vertex {
-	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 (a Vertex) Less(b Vertex) bool {
-	return a.x < b.x || a.y < b.y || a.z < b.z
-}
-
-func NewLine(p1, p2 Vertex) Line {
-	return Line{
-		p2.Sub(p1),
-		p1,
-	}
-}
-
-func (l Line) Eval(t float64) Vertex {
-	return l[0].scale(t).Add(l[1])
-}
-
-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 (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 (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
-}
--- a/pkg/octree/cache.go	Sat Sep 22 21:34:12 2018 +0200
+++ b/pkg/octree/cache.go	Sat Sep 22 21:57:30 2018 +0200
@@ -15,7 +15,7 @@
 
 	cacheEntry struct {
 		checksum string
-		tree     *Octree
+		tree     *Tree
 		access   time.Time
 	}
 	OctreeCache struct {
@@ -77,7 +77,7 @@
 func (oc *OctreeCache) Get(
 	bottleneck string, date time.Time,
 	conn *sql.Conn, ctx context.Context,
-) (*Octree, error) {
+) (*Tree, error) {
 	oc.Lock()
 	defer oc.Unlock()
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/loader.go	Sat Sep 22 21:57:30 2018 +0200
@@ -0,0 +1,119 @@
+package octree
+
+import (
+	"bufio"
+	"bytes"
+	"encoding/binary"
+	"log"
+	"os"
+
+	"github.com/golang/snappy"
+)
+
+func loadReader(r *bufio.Reader) (*Tree, error) {
+	tree := new(Tree)
+
+	if err := binary.Read(r, binary.LittleEndian, &tree.EPSG); err != nil {
+		return nil, err
+	}
+
+	log.Printf("EPSG: %d\n", tree.EPSG)
+
+	if err := tree.Min.read(r); err != nil {
+		return nil, err
+	}
+
+	if err := tree.Max.read(r); err != nil {
+		return nil, err
+	}
+
+	log.Printf("BBOX: [[%f, %f, %f], [%f, %f, %f]]\n",
+		tree.Min.X, tree.Min.Y, tree.Min.Z,
+		tree.Max.X, tree.Max.Y, tree.Max.Z)
+
+	var numVertices uint32
+	if err := binary.Read(r, binary.LittleEndian, &numVertices); err != nil {
+		return nil, err
+	}
+
+	log.Printf("vertices: %d\n", numVertices)
+
+	vertices := make([]Vertex, numVertices)
+	tree.vertices = vertices
+
+	for i := range vertices {
+		if err := vertices[i].read(r); err != nil {
+			return nil, err
+		}
+	}
+
+	var numTriangles uint32
+	if err := binary.Read(r, binary.LittleEndian, &numTriangles); err != nil {
+		return nil, err
+	}
+
+	log.Printf("triangles: %d\n", numTriangles)
+
+	indices := make([]int32, 3*numTriangles)
+	triangles := make([][]int32, numTriangles)
+	tree.triangles = triangles
+
+	var last int32
+
+	for i := range triangles {
+		tri := indices[:3]
+		indices = indices[3:]
+		triangles[i] = tri
+		for j := range tri {
+			v, err := binary.ReadVarint(r)
+			if err != nil {
+				return nil, err
+			}
+			value := int32(v) + last
+			tri[j] = value
+			last = value
+		}
+	}
+
+	var numNodes uint32
+	if err := binary.Read(r, binary.LittleEndian, &numNodes); err != nil {
+		return nil, err
+	}
+
+	log.Printf("num nodes: %d\n", numNodes)
+
+	tree.index = make([]int32, numNodes)
+	entries := tree.index[1:]
+
+	last = 0
+	for i := range entries {
+		v, err := binary.ReadVarint(r)
+		if err != nil {
+			return nil, err
+		}
+		value := int32(v) + last
+		entries[i] = value
+		last = value
+	}
+
+	return tree, nil
+}
+
+func LoadTree(fname string) (*Tree, error) {
+
+	f, err := os.Open(fname)
+	if err != nil {
+		return nil, err
+	}
+	defer f.Close()
+	return loadReader(
+		bufio.NewReader(
+			snappy.NewReader(f)))
+}
+
+func Deserialize(data []byte) (*Tree, error) {
+	return loadReader(
+		bufio.NewReader(
+			snappy.NewReader(
+				bytes.NewReader(data))))
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/octree.go	Sat Sep 22 21:57:30 2018 +0200
@@ -0,0 +1,81 @@
+package octree
+
+import (
+	"math"
+)
+
+type Tree struct {
+	EPSG uint32
+
+	vertices  []Vertex
+	triangles [][]int32
+	index     []int32
+
+	Min Vertex
+	Max Vertex
+}
+
+func (ot *Tree) Horizontal(h float64, fn func(*Triangle)) {
+
+	type frame struct {
+		pos int32
+		min float64
+		max float64
+	}
+
+	if h < ot.Min.Z || ot.Max.Z < h {
+		return
+	}
+
+	stack := []frame{{1, ot.Min.Z, ot.Max.Z}}
+
+	dupes := map[int32]struct{}{}
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		pos := top.pos
+		if pos == 0 {
+			continue
+		}
+		min, max := top.min, top.max
+
+		if pos > 0 { // node
+			if mid := (max-min)*0.5 + min; h >= mid {
+				pos += 4 // nodes with z-bit set
+				min = mid
+			} else {
+				max = mid
+			}
+			stack = append(stack,
+				frame{ot.index[pos+0], min, max},
+				frame{ot.index[pos+1], min, max},
+				frame{ot.index[pos+2], min, max},
+				frame{ot.index[pos+3], min, max})
+		} else { // leaf
+			pos = -pos - 1
+			n := ot.index[pos]
+			//log.Printf("%d %d %d\n", pos, n, len(ot.index))
+			indices := ot.index[pos+1 : pos+1+n]
+
+			for _, idx := range indices {
+				if _, found := dupes[idx]; found {
+					continue
+				}
+				tri := ot.triangles[idx]
+				t := Triangle{
+					ot.vertices[tri[0]],
+					ot.vertices[tri[1]],
+					ot.vertices[tri[2]],
+				}
+
+				if !(math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z)) > h ||
+					math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) < h) {
+					dupes[idx] = struct{}{}
+					fn(&t)
+				}
+			}
+		}
+	}
+}
--- a/pkg/octree/tree.go	Sat Sep 22 21:34:12 2018 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,11 +0,0 @@
-package octree
-
-import "fmt"
-
-type Octree struct {
-	// TODO: Implement me!
-}
-
-func Deserialize(data []byte) (*Octree, error) {
-	return nil, fmt.Errorf("Not implemented, yet.")
-}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/vertex.go	Sat Sep 22 21:57:30 2018 +0200
@@ -0,0 +1,195 @@
+package octree
+
+import (
+	"bytes"
+	"encoding/binary"
+	"io"
+	"math"
+)
+
+type (
+	Vertex struct {
+		X float64
+		Y float64
+		Z float64
+	}
+
+	Triangle [3]Vertex
+
+	Line [2]Vertex
+
+	LineStringZ      []Vertex
+	MultiLineStringZ []LineStringZ
+)
+
+const (
+	wkbNDR              byte   = 1
+	wkbPointZ           uint32 = 1000 + 1
+	wkbLineStringZ      uint32 = 1000 + 2
+	wkbMultiLineStringZ uint32 = 1000 + 5
+)
+
+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
+	}
+}
+
+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) Sub(w Vertex) Vertex {
+	return Vertex{
+		v.X - w.X,
+		v.Y - w.Y,
+		v.Z - w.Z,
+	}
+}
+
+func (v Vertex) Add(w Vertex) Vertex {
+	return Vertex{
+		v.X + w.X,
+		v.Y + w.Y,
+		v.Z + w.Z,
+	}
+}
+
+func (v Vertex) scale(s float64) Vertex {
+	return Vertex{
+		s * v.X,
+		s * v.Y,
+		s * v.Z,
+	}
+}
+
+func Interpolate(v1, v2 Vertex) func(Vertex) Vertex {
+	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 (a Vertex) Less(b Vertex) bool {
+	return a.X < b.X || a.Y < b.Y || a.Z < b.Z
+}
+
+func NewLine(p1, p2 Vertex) Line {
+	return Line{
+		p2.Sub(p1),
+		p1,
+	}
+}
+
+func (l Line) Eval(t float64) Vertex {
+	return l[0].scale(t).Add(l[1])
+}
+
+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 (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 (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
+}
+
+func (mls MultiLineStringZ) AsWKB() []byte {
+
+	var buf bytes.Buffer
+
+	binary.Write(&buf, binary.LittleEndian, wkbNDR)
+	binary.Write(&buf, binary.LittleEndian, wkbMultiLineStringZ)
+	binary.Write(&buf, binary.LittleEndian, uint32(len(mls)))
+
+	for _, ml := range mls {
+		binary.Write(&buf, binary.LittleEndian, wkbNDR)
+		binary.Write(&buf, binary.LittleEndian, wkbLineStringZ)
+		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()
+}