# HG changeset patch # User Sascha L. Teichmann # Date 1537516246 -7200 # Node ID fe0889460e97bcb8d522bbe1eb00d893b700c9c2 # Parent 18f6f2d89da1f97e1cc786841a04fdc90655c7d7# Parent e9c28c42c9272e2e6708251cb6b76d9d5544b341 Merged default into octree branch. diff -r e9c28c42c927 -r fe0889460e97 3rdpartylibs.sh --- a/3rdpartylibs.sh Fri Sep 21 00:24:41 2018 +0200 +++ b/3rdpartylibs.sh Fri Sep 21 09:50:46 2018 +0200 @@ -8,3 +8,4 @@ go get -u -v gopkg.in/gomail.v2 go get -u -v github.com/rs/cors go get -u -v golang.org/x/net/html/charset +go get -u -v github.com/golang/snappy diff -r e9c28c42c927 -r fe0889460e97 cmd/octree2contour/db.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octree2contour/db.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,39 @@ +package main + +import ( + "database/sql" + "flag" + + "github.com/jackc/pgx" + "github.com/jackc/pgx/stdlib" +) + +var ( + dbhost = flag.String("dbhost", "localhost", "database host") + dbport = flag.Uint("dbport", 5432, "database port") + dbname = flag.String("dbname", "gemma", "database user") + dbuser = flag.String("dbuser", "scott", "database user") + dbpassword = flag.String("dbpw", "tiger", "database password") + dbssl = flag.String("dbssl", "prefer", "database SSL mode") +) + +func run(fn func(*sql.DB) error) error { + + // To ease SSL config ride a bit on parsing. + cc, err := pgx.ParseConnectionString("sslmode=" + *dbssl) + if err != nil { + return err + } + + // Do the rest manually to allow whitespace in user/password. + cc.Host = *dbhost + cc.Port = uint16(*dbport) + cc.User = *dbuser + cc.Password = *dbpassword + cc.Database = *dbname + + db := stdlib.OpenDB(cc) + defer db.Close() + + return fn(db) +} diff -r e9c28c42c927 -r fe0889460e97 cmd/octree2contour/loader.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octree2contour/loader.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,205 @@ +package main + +import ( + "bufio" + "encoding/binary" + "io" + "log" + "math" + "os" + + "github.com/golang/snappy" +) + +type octree struct { + epsg uint32 + + vertices []vertex + triangles [][]int32 + index []int32 + + min vertex + max vertex +} + +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 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 (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) + } + } + } + } +} + +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 e9c28c42c927 -r fe0889460e97 cmd/octree2contour/main.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octree2contour/main.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,106 @@ +package main + +import ( + "flag" + "log" + "runtime" + "sort" + "sync" + "time" +) + +var ( + one = flag.Bool("o", false, "create only a single contour") + step = flag.Float64("s", 0.5, "step with") + max = flag.Float64("m", 10, "max height from lowest point") +) + +func processLevels( + tree *octree, + 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) { + line := t.intersectH(h) + if len(line) > 1 { + lines = append(lines, line) + } + }) + results <- result{h, lines} + } +} + +func process(tree *octree) []result { + + if *one { + var lines multiLineStringZ + tree.horizontal(*step, func(t *triangle) { + line := t.intersectH(*step) + if len(line) > 0 { + lines = append(lines, line) + } + }) + return []result{{*step, lines}} + } + + results := make(chan result) + done := make(chan struct{}) + + var all []result + go func() { + for { + select { + case r := <-results: + all = append(all, r) + case <-done: + return + } + } + }() + + var wg sync.WaitGroup + jobs := make(chan float64) + for i, n := 0, runtime.NumCPU(); i < n; i++ { + wg.Add(1) + go processLevels(tree, jobs, results, &wg) + } + for h := tree.min.z; h <= tree.max.z; h += *step { + jobs <- h + } + close(jobs) + wg.Wait() + done <- struct{}{} + + sort.Slice(all, func(i, j int) bool { + return all[i].h < all[j].h + }) + + return all +} + +func main() { + flag.Parse() + + for _, fname := range flag.Args() { + log.Printf("processing %s\n", fname) + start := time.Now() + tree, err := loadOctree(fname) + if err != nil { + log.Printf("error: %v\n", err) + continue + } + log.Printf("loading took: %v\n", time.Since(start)) + start = time.Now() + all := process(tree) + log.Printf("processing took: %v\n", time.Since(start)) + start = time.Now() + if err = store(all, tree.epsg); err != nil { + log.Printf("error: %v\n", err) + } + log.Printf("storing took: %v\n", time.Since(start)) + } +} diff -r e9c28c42c927 -r fe0889460e97 cmd/octree2contour/store.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octree2contour/store.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,75 @@ +package main + +import ( + "bytes" + "database/sql" + "encoding/binary" + "math" +) + +type multiLineStringZ [][]vertex + +const ( + wkbNDR byte = 1 + wkbPointZ uint32 = 1000 + 1 + wkbLineStringZ uint32 = 1000 + 2 + wkbMultiLineStringZ uint32 = 1000 + 5 +) + +type result struct { + h float64 + lines multiLineStringZ +} + +const insertSQL = ` +INSERT INTO waterway.contour_lines (height, geom) +VALUES ($1, ST_Transform( + ST_SetSRID(ST_GeomFromWKB($2), $3), 4326)::geography) +` + +func store(all []result, epsg uint32) error { + + return run(func(db *sql.DB) error { + + tx, err := db.Begin() + if err != nil { + return err + } + defer tx.Rollback() + + stmt, err := tx.Prepare(insertSQL) + if err != nil { + return err + } + + for _, r := range all { + if _, err := stmt.Exec(r.h, r.lines.asWKB(), epsg); err != nil { + return err + } + } + + 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 e9c28c42c927 -r fe0889460e97 cmd/octree2contour/vertex.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octree2contour/vertex.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,136 @@ +package main + +type vertex struct { + x float64 + y float64 + z float64 +} + +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 +} + +type line [2]vertex + +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) intersectH(h float64) vertex { + t := (h - l[1].z) / l[0].z + return l.eval(t) +} + +type triangle [3]vertex + +func side(z, h float64) int { + switch { + case z < h: + return -1 + case z > h: + return +1 + } + return 0 +} + +func (t *triangle) intersectH(h float64) []vertex { + sides := [3]int{ + side(t[0].z, h), + side(t[1].z, h), + side(t[2].z, h), + } + + var points []vertex + + 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]).intersectH(h) + points = append(points, v) + } + } + + return points +} diff -r e9c28c42c927 -r fe0889460e97 cmd/tin2octree/builder.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/tin2octree/builder.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,154 @@ +package main + +import ( + "encoding/binary" + "io" + "log" +) + +type treeBuilder 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 (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 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 *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 +} diff -r e9c28c42c927 -r fe0889460e97 cmd/tin2octree/main.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/tin2octree/main.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,163 @@ +package main + +import ( + "database/sql" + "flag" + "log" + "math" + "os" + "time" + + "github.com/golang/snappy" + "github.com/jackc/pgx" + "github.com/jackc/pgx/stdlib" +) + +var ( + bottleneck = flag.String("bottleneck", "", "bottleneck id") + date = flag.String("date", "", "date info") + file = flag.String("file", "", "save to file") + + utm = flag.Bool("utm", false, "fetch in matchin UTM zone") + + dbhost = flag.String("dbhost", "localhost", "database host") + dbport = flag.Uint("dbport", 5432, "database port") + dbname = flag.String("dbname", "gemma", "database user") + dbuser = flag.String("dbuser", "scott", "database user") + dbpassword = flag.String("dbpw", "tiger", "database password") + dbssl = flag.String("dbssl", "prefer", "database SSL mode") +) + +func run(fn func(*sql.DB) error) error { + + // To ease SSL config ride a bit on parsing. + cc, err := pgx.ParseConnectionString("sslmode=" + *dbssl) + if err != nil { + return err + } + + // Do the rest manually to allow whitespace in user/password. + cc.Host = *dbhost + cc.Port = uint16(*dbport) + cc.User = *dbuser + cc.Password = *dbpassword + cc.Database = *dbname + + db := stdlib.OpenDB(cc) + defer db.Close() + + return fn(db) +} + +const ( + tinSQL = ` +SELECT ST_AsBinary(ST_DelaunayTriangles(point_cloud::geometry, 0, 2)) +FROM waterway.sounding_results +WHERE bottleneck_id = $1 AND date_info = $2 +` + tinUTMSQL = ` +SELECT ST_AsBinary( + ST_DelaunayTriangles( + ST_Transform(point_cloud::geometry, $3::int), 0, 2)) +FROM waterway.sounding_results +WHERE bottleneck_id = $1 AND date_info = $2 +` + centroidSQL = ` +SELECT ST_X(ST_Centroid(point_cloud::geometry)), ST_Y(ST_Centroid(point_cloud::geometry)) +FROM waterway.sounding_results +WHERE bottleneck_id = $1 AND date_info = $2 +` +) + +func utmZone(x, y float64) int { + var pref int + if y > 0 { + pref = 32600 + } else { + pref = 32700 + } + zone := int(math.Floor((x+180)/6)) + 1 + return zone + pref +} + +func main() { + flag.Parse() + + if *bottleneck == "" || *date == "" { + log.Fatalln("missing bottleneck or date option.") + } + + dateInfo, err := time.Parse("2006-01-02", *date) + if err != nil { + log.Fatalf("error: %v\n", err) + } + + var t tin + + if err := run(func(db *sql.DB) error { + var utmZ int + + if *utm { + var cx, cy float64 + err := db.QueryRow(centroidSQL, *bottleneck, dateInfo).Scan(&cx, &cy) + switch { + case err == sql.ErrNoRows: + return nil + case err != nil: + return err + } + log.Printf("lat/lon: [%f, %f]\n", cx, cy) + utmZ = utmZone(cx, cy) + log.Printf("UTM zone: %d\n", utmZ) + } + + start := time.Now() + var err error + if *utm { + err = db.QueryRow(tinUTMSQL, *bottleneck, dateInfo, utmZ).Scan(&t) + } else { + err = db.QueryRow(tinSQL, *bottleneck, dateInfo).Scan(&t) + } + switch { + case err == sql.ErrNoRows: + return nil + case err != nil: + return err + } + log.Printf("query took: %s\n", time.Since(start)) + + if *utm { + t.epsg = uint32(utmZ) + } + + return nil + }); err != nil { + log.Fatalf("error: %v\n", err) + } + + tb := &treeBuilder{t: &t} + tb.build() + + if *file != "" { + f, err := os.Create(*file) + if err != nil { + log.Printf("error: %v\n", err) + } + out := snappy.NewBufferedWriter(f) + if err := t.Serialize(out); err != nil { + f.Close() + log.Fatalf("error: %v\n", err) + } + if err := tb.Serialize(out); err != nil { + f.Close() + log.Fatalf("error: %v\n", err) + } + if err := out.Flush(); err != nil { + f.Close() + log.Fatalf("error: %v\n", err) + } + if err := f.Close(); err != nil { + log.Fatalf("error: %v\n", err) + } + } +} diff -r e9c28c42c927 -r fe0889460e97 cmd/tin2octree/tin.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/tin2octree/tin.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,245 @@ +package main + +import ( + "bytes" + "encoding/binary" + "errors" + "fmt" + "io" + "log" + "math" +) + +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 []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 (v *vertex) write(w io.Writer) error { + for _, addr := range []*float64{&v.x, &v.y, &v.z} { + if err := binary.Write( + w, binary.LittleEndian, math.Float64bits(*addr)); err != nil { + return err + } + } + return nil +} + +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 +} diff -r e9c28c42c927 -r fe0889460e97 cmd/tin2octree/vertex.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/tin2octree/vertex.go Fri Sep 21 09:50:46 2018 +0200 @@ -0,0 +1,62 @@ +package main + +type vertex struct { + x float64 + y float64 + z float64 +} + +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 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 +}