Mercurial > gemma
changeset 968:a4fe07a21ba7
Moved octree builder into octree package to be reusable by the sounding result import job.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 17 Oct 2018 16:00:49 +0200 |
parents | 2025074ad835 |
children | a4e003ba0074 |
files | cmd/tin2octree/builder.go cmd/tin2octree/main.go cmd/tin2octree/tin.go pkg/octree/builder.go pkg/octree/tin.go pkg/octree/vertex.go pkg/octree/wkb.go |
diffstat | 7 files changed, 421 insertions(+), 422 deletions(-) [+] |
line wrap: on
line diff
--- a/cmd/tin2octree/builder.go Wed Oct 17 15:51:48 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,156 +0,0 @@ -package main - -import ( - "encoding/binary" - "io" - "log" - - "gemma.intevation.de/gemma/pkg/octree" -) - -type treeBuilder struct { - t *tin - nodes int - leaves int - index []int32 -} - -var cubes = [8][2]octree.Vertex{ - makeCube(0), - makeCube(1), - makeCube(2), - makeCube(3), - makeCube(4), - makeCube(5), - makeCube(6), - makeCube(7), -} - -func makeCube(i int) [2]octree.Vertex { - var d octree.Vertex - if i&1 == 1 { - d.X = 0.5 - } - if i&2 == 2 { - d.Y = 0.5 - } - if i&4 == 4 { - d.Z = 0.5 - } - return [2]octree.Vertex{ - octree.Vertex{0.0, 0.0, 0.0}.Add(d), - octree.Vertex{0.5, 0.5, 0.5}.Add(d), - } -} - -func (tb *treeBuilder) build() { - - triangles := make([]int32, len(tb.t.triangles)) - for i := range triangles { - triangles[i] = int32(i) - } - - tb.index = append(tb.index, 0) - - tb.buildRecursive(triangles, tb.t.min, tb.t.max, 0) - tb.index[0] = int32(len(tb.index)) - log.Printf("num nodes: %d\n", tb.index[0]) - - log.Printf("nodes: %d leaves: %d index %d\n", - tb.nodes, tb.leaves, tb.index[0]) -} - -func (tb *treeBuilder) buildRecursive( - triangles []int32, - min, max octree.Vertex, - depth int, -) int32 { - if len(triangles) <= 16 || depth > 8 { - pos := len(tb.index) - tb.index = append(tb.index, int32(len(triangles))) - tb.index = append(tb.index, triangles...) - //log.Printf("leaf entries: %d (%d)\n", len(triangles), depth) - tb.leaves++ - return int32(-(pos + 1)) - } - - pos := len(tb.index) - tb.index = append(tb.index, - 0, 0, 0, 0, - 0, 0, 0, 0) - - bbox := octree.Interpolate(min, max) - - bboxes := make([][2]octree.Vertex, len(cubes)) - - for i := range cubes { - bboxes[i] = [2]octree.Vertex{ - bbox(cubes[i][0]), - bbox(cubes[i][1]), - } - } - - var quandrants [8][]int32 - - for _, tri := range triangles { - triangle := tb.t.triangles[tri] - v0 := tb.t.vertices[triangle[0]] - v1 := tb.t.vertices[triangle[1]] - v2 := tb.t.vertices[triangle[2]] - - l := v0 - l.Minimize(v1) - l.Minimize(v2) - - h := v0 - h.Maximize(v1) - h.Maximize(v2) - - for i := range bboxes { - if !(h.Less(bboxes[i][0]) || bboxes[i][1].Less(l)) { - quandrants[i] = append(quandrants[i], tri) - } - } - } - - for i := range quandrants { - if len(quandrants[i]) > 0 { - child := tb.buildRecursive( - quandrants[i], - bboxes[i][0], bboxes[i][1], - depth+1) - tb.index[pos+i] = child - } - } - tb.nodes++ - return int32(pos) -} - -func (tb *treeBuilder) Serialize(w io.Writer) error { - var buf [binary.MaxVarintLen32]byte - - if err := binary.Write(w, binary.LittleEndian, tb.index[0]); err != nil { - return err - } - - var last int32 - var written int - - for _, x := range tb.index[1:] { - delta := x - last - n := binary.PutVarint(buf[:], int64(delta)) - for p := buf[:n]; len(p) > 0; p = p[n:] { - var err error - if n, err = w.Write(p); err != nil { - return err - } - written += n - } - - last = x - } - log.Printf("compressed octree index in bytes: %d (%d)\n", - written, 4*len(tb.index)) - - return nil -}
--- a/cmd/tin2octree/main.go Wed Oct 17 15:51:48 2018 +0200 +++ b/cmd/tin2octree/main.go Wed Oct 17 16:00:49 2018 +0200 @@ -13,7 +13,7 @@ "os" "time" - "github.com/golang/snappy" + "gemma.intevation.de/gemma/pkg/octree" "github.com/jackc/pgx" "github.com/jackc/pgx/stdlib" ) @@ -97,7 +97,7 @@ log.Fatalf("error: %v\n", err) } - var t tin + var t octree.Tin if err := run(func(db *sql.DB) error { var utmZ int @@ -132,7 +132,7 @@ log.Printf("query took: %s\n", time.Since(start)) if *utm { - t.epsg = uint32(utmZ) + t.EPSG = uint32(utmZ) } return nil @@ -140,8 +140,8 @@ log.Fatalf("error: %v\n", err) } - tb := &treeBuilder{t: &t} - tb.build() + tb := octree.NewBuilder(&t) + tb.Build() if *insert { var w io.Writer @@ -157,7 +157,7 @@ } var buf bytes.Buffer - if err := write(&buf, tb); err != nil { + if err := tb.WriteTo(&buf); err != nil { log.Fatalf("error: %v\n", err) } data := buf.String() @@ -188,7 +188,7 @@ if err != nil { log.Printf("error: %v\n", err) } - err = write(f, tb) + err = tb.WriteTo(f) if err2 := f.Close(); err == nil { if err != nil { err = err2 @@ -199,14 +199,3 @@ } } } - -func write(w io.Writer, tb *treeBuilder) error { - out := snappy.NewBufferedWriter(w) - if err := tb.t.Serialize(out); err != nil { - return err - } - if err := tb.Serialize(out); err != nil { - return err - } - return out.Flush() -}
--- a/cmd/tin2octree/tin.go Wed Oct 17 15:51:48 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,237 +0,0 @@ -package main - -import ( - "bytes" - "encoding/binary" - "errors" - "fmt" - "io" - "log" - "math" - - "gemma.intevation.de/gemma/pkg/octree" -) - -var ( - errNoByteSlice = errors.New("Not a byte slice") - errTooLessPoints = errors.New("Too less points") -) - -const ( - wkbXDR byte = 0 - wkbNDR byte = 1 - wkbTinZ uint32 = 1000 + 16 - wkbTriangleZ uint32 = 1000 + 17 -) - -const wgs84 = 4326 - -type tin struct { - epsg uint32 - vertices []octree.Vertex - triangles [][]int32 - - min octree.Vertex - max octree.Vertex -} - -func (t *tin) FromWKB(data []byte) error { - log.Printf("data length %d\n", len(data)) - - r := bytes.NewReader(data) - - endian, err := r.ReadByte() - - var order binary.ByteOrder - - switch { - case err != nil: - return err - case endian == wkbNDR: - order = binary.LittleEndian - case endian == wkbXDR: - 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 != wkbTinZ: - return fmt.Errorf("unknown geometry type %x", geomType) - } - - var num uint32 - if err = binary.Read(r, order, &num); err != nil { - return err - } - - vertices := make([]octree.Vertex, 0, 100000) - - var v octree.Vertex - - v2i := make(map[octree.Vertex]int32, 100000) - - var indexPool []int32 - - allocIndices := func() []int32 { - if len(indexPool) == 0 { - indexPool = make([]int32, 3*8*1024) - } - ids := indexPool[:3] - indexPool = indexPool[3:] - return ids - } - - var triangles [][]int32 - - min := octree.Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} - max := octree.Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} - - for i := uint32(0); i < num; i++ { - - endian, err = r.ReadByte() - switch { - case err != nil: - return err - case endian == wkbNDR: - order = binary.LittleEndian - case endian == wkbXDR: - 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 != wkbTriangleZ: - return fmt.Errorf("unknown geometry type %d", geomType) - } - - var rings uint32 - if err = binary.Read(r, order, &rings); err != nil { - return err - } - triangle := allocIndices() - - for ring := uint32(0); ring < rings; ring++ { - var npoints uint32 - if err = binary.Read(r, order, &npoints); err != nil { - return err - } - - if npoints < 3 { - return errTooLessPoints - } - - for p := uint32(0); p < npoints; p++ { - var x, y, z uint64 - for _, addr := range []*uint64{&x, &y, &z} { - if err = binary.Read(r, order, addr); err != nil { - return err - } - } - if p >= 3 || ring >= 1 { - // Don't store the forth point. - continue - } - // Do this conversion later to spare reflect calls - // and allocs in binary.Read. - v.X = math.Float64frombits(x) - v.Y = math.Float64frombits(y) - v.Z = math.Float64frombits(z) - idx, found := v2i[v] - if !found { - idx = int32(len(vertices)) - v2i[v] = idx - vertices = append(vertices, v) - min.Minimize(v) - max.Maximize(v) - } - triangle[p] = idx - } - } - triangles = append(triangles, triangle) - } - - log.Printf("bbox: [[%f, %f], [%f, %f]]\n", - min.X, min.Y, max.X, max.Y) - - *t = tin{ - epsg: wgs84, - vertices: vertices, - triangles: triangles, - min: min, - max: max, - } - - return nil -} - -func (t *tin) Scan(raw interface{}) error { - - data, ok := raw.([]byte) - if !ok { - return errNoByteSlice - } - return t.FromWKB(data) -} - -func (t *tin) Serialize(w io.Writer) error { - - if err := binary.Write(w, binary.LittleEndian, t.epsg); err != nil { - return err - } - - if err := t.min.Write(w); err != nil { - return err - } - if err := t.max.Write(w); err != nil { - return err - } - - if err := binary.Write( - w, binary.LittleEndian, uint32(len(t.vertices))); err != nil { - return err - } - - for _, v := range t.vertices { - if err := v.Write(w); err != nil { - return err - } - } - log.Printf("vertices %d (%d)\n", len(t.vertices), len(t.vertices)*3*8) - - if err := binary.Write( - w, binary.LittleEndian, uint32(len(t.triangles))); err != nil { - return err - } - - var buf [binary.MaxVarintLen32]byte - var written int - var last int32 - for _, triangle := range t.triangles { - for _, idx := range triangle { - value := idx - last - n := binary.PutVarint(buf[:], int64(value)) - for p := buf[:n]; len(p) > 0; p = p[n:] { - var err error - if n, err = w.Write(p); err != nil { - return err - } - written += n - } - last = idx - } - } - log.Printf("compressed tin indices in bytes: %d (%d)\n", - written, 3*4*len(t.triangles)) - - return nil -}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/builder.go Wed Oct 17 16:00:49 2018 +0200 @@ -0,0 +1,171 @@ +package octree + +import ( + "encoding/binary" + "io" + "log" + + "github.com/golang/snappy" +) + +type Builder struct { + t *Tin + nodes int + leaves int + index []int32 +} + +var cubes = [8][2]Vertex{ + makeCube(0), + makeCube(1), + makeCube(2), + makeCube(3), + makeCube(4), + makeCube(5), + makeCube(6), + makeCube(7), +} + +func makeCube(i int) [2]Vertex { + var d Vertex + if i&1 == 1 { + d.X = 0.5 + } + if i&2 == 2 { + d.Y = 0.5 + } + if i&4 == 4 { + d.Z = 0.5 + } + return [2]Vertex{ + Vertex{0.0, 0.0, 0.0}.Add(d), + Vertex{0.5, 0.5, 0.5}.Add(d), + } +} + +func NewBuilder(t *Tin) *Builder { + return &Builder{t: t} +} + +func (tb *Builder) Build() { + + triangles := make([]int32, len(tb.t.Triangles)) + for i := range triangles { + triangles[i] = int32(i) + } + + tb.index = append(tb.index, 0) + + tb.buildRecursive(triangles, tb.t.Min, tb.t.Max, 0) + tb.index[0] = int32(len(tb.index)) + log.Printf("num nodes: %d\n", tb.index[0]) + + log.Printf("nodes: %d leaves: %d index %d\n", + tb.nodes, tb.leaves, tb.index[0]) +} + +func (tb *Builder) buildRecursive( + triangles []int32, + min, max Vertex, + depth int, +) int32 { + if len(triangles) <= 16 || depth > 8 { + pos := len(tb.index) + tb.index = append(tb.index, int32(len(triangles))) + tb.index = append(tb.index, triangles...) + //log.Printf("leaf entries: %d (%d)\n", len(triangles), depth) + tb.leaves++ + return int32(-(pos + 1)) + } + + pos := len(tb.index) + tb.index = append(tb.index, + 0, 0, 0, 0, + 0, 0, 0, 0) + + bbox := Interpolate(min, max) + + bboxes := make([][2]Vertex, len(cubes)) + + for i := range cubes { + bboxes[i] = [2]Vertex{ + bbox(cubes[i][0]), + bbox(cubes[i][1]), + } + } + + var quandrants [8][]int32 + + for _, tri := range triangles { + triangle := tb.t.Triangles[tri] + v0 := tb.t.Vertices[triangle[0]] + v1 := tb.t.Vertices[triangle[1]] + v2 := tb.t.Vertices[triangle[2]] + + l := v0 + l.Minimize(v1) + l.Minimize(v2) + + h := v0 + h.Maximize(v1) + h.Maximize(v2) + + for i := range bboxes { + if !(h.Less(bboxes[i][0]) || bboxes[i][1].Less(l)) { + quandrants[i] = append(quandrants[i], tri) + } + } + } + + for i := range quandrants { + if len(quandrants[i]) > 0 { + child := tb.buildRecursive( + quandrants[i], + bboxes[i][0], bboxes[i][1], + depth+1) + tb.index[pos+i] = child + } + } + tb.nodes++ + return int32(pos) +} + +func (tb *Builder) Serialize(w io.Writer) error { + var buf [binary.MaxVarintLen32]byte + + if err := binary.Write(w, binary.LittleEndian, tb.index[0]); err != nil { + return err + } + + var last int32 + var written int + + for _, x := range tb.index[1:] { + delta := x - last + n := binary.PutVarint(buf[:], int64(delta)) + for p := buf[:n]; len(p) > 0; p = p[n:] { + var err error + if n, err = w.Write(p); err != nil { + return err + } + written += n + } + + last = x + } + log.Printf("compressed octree index in bytes: %d (%d)\n", + written, 4*len(tb.index)) + + return nil +} + +func (tb *Builder) WriteTo(w io.Writer) error { + out := snappy.NewBufferedWriter(w) + if err := tb.t.Serialize(out); err != nil { + return err + } + if err := tb.Serialize(out); err != nil { + return err + } + return out.Flush() +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/tin.go Wed Oct 17 16:00:49 2018 +0200 @@ -0,0 +1,228 @@ +package octree + +import ( + "bytes" + "encoding/binary" + "errors" + "fmt" + "io" + "log" + "math" +) + +var ( + errNoByteSlice = errors.New("Not a byte slice") + errTooLessPoints = errors.New("Too less points") +) + +const wgs84 = 4326 + +type Tin struct { + EPSG uint32 + Vertices []Vertex + Triangles [][]int32 + + Min Vertex + Max Vertex +} + +func (t *Tin) FromWKB(data []byte) error { + log.Printf("data length %d\n", len(data)) + + r := bytes.NewReader(data) + + endian, err := r.ReadByte() + + var order binary.ByteOrder + + switch { + case err != nil: + return err + case endian == wkbNDR: + order = binary.LittleEndian + case endian == wkbXDR: + 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 != wkbTinZ: + return fmt.Errorf("unknown geometry type %x", geomType) + } + + var num uint32 + if err = binary.Read(r, order, &num); err != nil { + return err + } + + vertices := make([]Vertex, 0, 100000) + + var v Vertex + + v2i := make(map[Vertex]int32, 100000) + + var indexPool []int32 + + allocIndices := func() []int32 { + if len(indexPool) == 0 { + indexPool = make([]int32, 3*8*1024) + } + ids := indexPool[:3] + indexPool = indexPool[3:] + return ids + } + + var triangles [][]int32 + + min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} + max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} + + for i := uint32(0); i < num; i++ { + + endian, err = r.ReadByte() + switch { + case err != nil: + return err + case endian == wkbNDR: + order = binary.LittleEndian + case endian == wkbXDR: + 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 != wkbTriangleZ: + return fmt.Errorf("unknown geometry type %d", geomType) + } + + var rings uint32 + if err = binary.Read(r, order, &rings); err != nil { + return err + } + triangle := allocIndices() + + for ring := uint32(0); ring < rings; ring++ { + var npoints uint32 + if err = binary.Read(r, order, &npoints); err != nil { + return err + } + + if npoints < 3 { + return errTooLessPoints + } + + for p := uint32(0); p < npoints; p++ { + var x, y, z uint64 + for _, addr := range []*uint64{&x, &y, &z} { + if err = binary.Read(r, order, addr); err != nil { + return err + } + } + if p >= 3 || ring >= 1 { + // Don't store the forth point. + continue + } + // Do this conversion later to spare reflect calls + // and allocs in binary.Read. + v.X = math.Float64frombits(x) + v.Y = math.Float64frombits(y) + v.Z = math.Float64frombits(z) + idx, found := v2i[v] + if !found { + idx = int32(len(vertices)) + v2i[v] = idx + vertices = append(vertices, v) + min.Minimize(v) + max.Maximize(v) + } + triangle[p] = idx + } + } + triangles = append(triangles, triangle) + } + + log.Printf("bbox: [[%f, %f], [%f, %f]]\n", + min.X, min.Y, max.X, max.Y) + + *t = Tin{ + EPSG: wgs84, + Vertices: vertices, + Triangles: triangles, + Min: min, + Max: max, + } + + return nil +} + +func (t *Tin) Scan(raw interface{}) error { + + data, ok := raw.([]byte) + if !ok { + return errNoByteSlice + } + return t.FromWKB(data) +} + +func (t *Tin) Serialize(w io.Writer) error { + + if err := binary.Write(w, binary.LittleEndian, t.EPSG); err != nil { + return err + } + + if err := t.Min.Write(w); err != nil { + return err + } + if err := t.Max.Write(w); err != nil { + return err + } + + if err := binary.Write( + w, binary.LittleEndian, uint32(len(t.Vertices))); err != nil { + return err + } + + for _, v := range t.Vertices { + if err := v.Write(w); err != nil { + return err + } + } + log.Printf("vertices %d (%d)\n", len(t.Vertices), len(t.Vertices)*3*8) + + if err := binary.Write( + w, binary.LittleEndian, uint32(len(t.Triangles))); err != nil { + return err + } + + var buf [binary.MaxVarintLen32]byte + var written int + var last int32 + for _, triangle := range t.Triangles { + for _, idx := range triangle { + value := idx - last + n := binary.PutVarint(buf[:], int64(value)) + for p := buf[:n]; len(p) > 0; p = p[n:] { + var err error + if n, err = w.Write(p); err != nil { + return err + } + written += n + } + last = idx + } + } + log.Printf("compressed tin indices in bytes: %d (%d)\n", + written, 3*4*len(t.Triangles)) + + return nil +}
--- a/pkg/octree/vertex.go Wed Oct 17 15:51:48 2018 +0200 +++ b/pkg/octree/vertex.go Wed Oct 17 16:00:49 2018 +0200 @@ -38,17 +38,6 @@ } ) -const ( - wkbNDR byte = 1 - - wkbLineString uint32 = 2 - wkbMultiLineString uint32 = 5 - wkbPointZ uint32 = 1000 + 1 - wkbLineStringZ uint32 = 1000 + 2 - wkbMultiPointZ uint32 = 1000 + 4 - wkbMultiLineStringZ uint32 = 1000 + 5 -) - func (v *Vertex) Minimize(w Vertex) { if w.X < v.X { v.X = w.X
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/wkb.go Wed Oct 17 16:00:49 2018 +0200 @@ -0,0 +1,15 @@ +package octree + +const ( + wkbXDR byte = 0 + wkbNDR byte = 1 + + wkbLineString uint32 = 2 + wkbMultiLineString uint32 = 5 + wkbPointZ uint32 = 1000 + 1 + wkbLineStringZ uint32 = 1000 + 2 + wkbMultiPointZ uint32 = 1000 + 4 + wkbMultiLineStringZ uint32 = 1000 + 5 + wkbTinZ uint32 = 1000 + 16 + wkbTriangleZ uint32 = 1000 + 17 +)