# HG changeset patch # User Sascha L. Teichmann # Date 1537646250 -7200 # Node ID 41c8dc61f38f4f46bed96c8dd2e232ab224391a1 # Parent 5af9ab39e71575fd65218b28ced5ff7a3bc5618c Moved octree loading stuff to octree package. diff -r 5af9ab39e715 -r 41c8dc61f38f cmd/octree2contour/loader.go --- 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))) -} diff -r 5af9ab39e715 -r 41c8dc61f38f cmd/octree2contour/main.go --- 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 diff -r 5af9ab39e715 -r 41c8dc61f38f cmd/octree2contour/octree.go --- 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) - } - } - } - } -} diff -r 5af9ab39e715 -r 41c8dc61f38f cmd/octree2contour/store.go --- 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() -} diff -r 5af9ab39e715 -r 41c8dc61f38f cmd/octree2contour/vertex.go --- 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 -} diff -r 5af9ab39e715 -r 41c8dc61f38f pkg/octree/cache.go --- 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() diff -r 5af9ab39e715 -r 41c8dc61f38f pkg/octree/loader.go --- /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)))) +} diff -r 5af9ab39e715 -r 41c8dc61f38f pkg/octree/octree.go --- /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) + } + } + } + } +} diff -r 5af9ab39e715 -r 41c8dc61f38f pkg/octree/tree.go --- 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.") -} diff -r 5af9ab39e715 -r 41c8dc61f38f pkg/octree/vertex.go --- /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() +}