Mercurial > gemma
view cmd/octreediff/main.go @ 2476:efe332e985b9 octree-diff
Calculate the iso lines of the difference.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Thu, 28 Feb 2019 11:51:31 +0100 |
parents | 19beb7d17337 |
children | c85b16db8a02 |
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 main import ( "bytes" "context" "database/sql" "encoding/binary" "errors" "flag" "fmt" "log" "math" "runtime" "sync" "time" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/octree" "gemma.intevation.de/gemma/pkg/wkb" ) var ( bottleneck = flag.String("bottleneck", "", "name of the bottleneck") first = flag.String("first", "", "date of the first sounding result") second = flag.String("second", "", "date of the second sounding result") ) const ( loadOctreeSQL = ` SELECT sr.octree_index FROM waterway.sounding_results sr JOIN waterway.bottlenecks bn ON sr.bottleneck_id = bn.id WHERE bn.bottleneck_id = $1 AND sr.date_info = $2::date AND sr.octree_index IS NOT NULL` triangulateSQL = ` SELECT ST_AsBinary(ST_DelaunayTriangles(ST_GeomFromWKB($1, $2::int), 0, 2)) ` ) func check(err error) { if err != nil { log.Fatalf("error: %v\n", err) } } type point struct { x float64 y float64 } type pointMap map[point]float64 func (pm pointMap) asWKB() []byte { size := 1 + 4 + 4 + len(pm)*(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(pm))) 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, z := range pm { 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(z)) } return buf.Bytes() } func sliceWork( vs []octree.Vertex, dst pointMap, fn func([]octree.Vertex, func([]octree.Vertex) []octree.Vertex), ) { n := runtime.NumCPU() wg := new(sync.WaitGroup) slices := make(chan []octree.Vertex) out := make(chan []octree.Vertex) pool := make(chan []octree.Vertex, n) const pageSize = 2048 turn := func(p []octree.Vertex) []octree.Vertex { if p != nil { out <- p } select { case p = <-pool: default: p = make([]octree.Vertex, 0, pageSize) } return p } for i := 0; i < n; i++ { wg.Add(1) go func() { defer wg.Done() for slice := range slices { fn(slice, turn) } }() } done := make(chan struct{}) go func() { defer close(done) for s := range out { for i := range s { v := &s[i] key := point{x: v.X, y: v.Y} if z, found := dst[key]; found { dst[key] = (z + v.Z) * 0.5 } else { dst[key] = v.Z } } select { case pool <- s[:0:pageSize]: default: } } }() size := len(vs)/n + 1 for len(vs) > 0 { var l int if len(vs) < size { l = len(vs) } else { l = size } slices <- vs[:l] vs = vs[l:] } close(slices) wg.Wait() close(out) <-done } func process(bottleneck string, firstDate, secondDate time.Time) error { start := time.Now() defer func() { log.Printf("processing took %v\n", time.Since(start)) }() ctx := context.Background() return run(func(db *sql.DB) error { conn, err := db.Conn(ctx) if err != nil { return err } defer conn.Close() type load struct { date time.Time data []byte err *error dst **octree.Tree } out := make(chan *load) wg := new(sync.WaitGroup) n := runtime.NumCPU() if n > 2 { n = 2 } for i := 0; i < n; i++ { wg.Add(1) go func() { defer wg.Done() for l := range out { if *l.err == nil { *l.dst, *l.err = octree.Deserialize(l.data) } } }() } var firstErr, secondErr error var first, second *octree.Tree for _, l := range []*load{ {date: firstDate, dst: &first, err: &firstErr}, {date: secondDate, dst: &second, err: &secondErr}, } { var data []byte if err := conn.QueryRowContext( ctx, loadOctreeSQL, bottleneck, l.date, ).Scan(&data); err != nil { *l.err = err } else { l.data = data } out <- l } close(out) wg.Wait() if firstErr != nil || secondErr != nil { if firstErr != nil && secondErr != nil { return fmt.Errorf("%v, %v", firstErr, secondErr) } if firstErr != nil { return firstErr } return secondErr } if first.EPSG != second.EPSG { return errors.New("EPSG codes mismatch. Needs transformation slow pass.") } now := time.Now() log.Printf("loading took %v\n", now.Sub(start)) last := now firstVs, secondVs := first.Vertices(), second.Vertices() result := make(pointMap, len(firstVs)+len(secondVs)) sliceWork( firstVs, result, func( slice []octree.Vertex, turn func([]octree.Vertex) []octree.Vertex, ) { p := turn(nil) for i := range slice { v := &slice[i] if z, found := second.Value(v.X, v.Y); found { p = append(p, octree.Vertex{v.X, v.Y, v.Z - z}) if len(p) == cap(p) { p = turn(p) } } } if len(p) > 0 { turn(p) } }) sliceWork( secondVs, result, func( slice []octree.Vertex, turn func([]octree.Vertex) []octree.Vertex, ) { p := turn(nil) for i := range slice { v := &slice[i] if z, found := first.Value(v.X, v.Y); found { p = append(p, octree.Vertex{v.X, v.Y, z - v.Z}) if len(p) == cap(p) { p = turn(p) } } } if len(p) > 0 { turn(p) } }) now = time.Now() log.Printf("setting in took %v\n", now.Sub(last)) last = now log.Printf("num points: %d\n", len(result)) data := result.asWKB() now = time.Now() log.Printf("turing into WKB took %v\n", now.Sub(last)) last = now log.Printf("WKB size %.3fMB\n", float64(len(data))/(1024*1024)) var tin octree.Tin if err := conn.QueryRowContext( ctx, triangulateSQL, data, first.EPSG, ).Scan(&tin); err != nil { return err } now = time.Now() log.Printf("triangulation took %v\n", now.Sub(last)) last = now builder := octree.NewBuilder(&tin) builder.Build() now = time.Now() log.Printf("building octree took %v\n", now.Sub(last)) last = now tree := builder.Tree() log.Printf("min/max: %f %f\n", tree.Min.Z, tree.Max.Z) var heights []float64 switch { case tree.Min.Z >= 0: // All values positive. for v := 0.0; v <= tree.Max.Z; v += 0.1 { if v >= tree.Min.Z { heights = append(heights, v) } } case tree.Max.Z <= 0: // All values negative. for v := 0.0; v >= tree.Min.Z; v -= 0.1 { if v <= tree.Max.Z { heights = append(heights, v) } } default: // Positive and negative. for v := 0.1; v <= tree.Max.Z; v += 0.1 { heights = append(heights, v) } for i, j := 0, len(heights)-1; i < j; i, j = i+1, j-1 { heights[i], heights[j] = heights[j], heights[i] } for v := 0.0; v >= tree.Min.Z; v -= 0.1 { heights = append(heights, v) } } var dataSize int octree.DoContours(tree, heights, func(res *octree.ContourResult) { // TODO: Store them. log.Printf("%f: segments: %d\n", res.Height, len(res.Lines)) dataSize += len(res.Lines.AsWKB2D()) }) now = time.Now() log.Printf("Number of iso lines: %d\n", len(heights)) log.Printf("Total WKB size: %.2fMB\n", float64(dataSize)/(1024*1024)) log.Printf("generating iso lines took %v\n", now.Sub(last)) last = now return nil }) } func main() { flag.Parse() firstDate, err := time.Parse(common.DateFormat, *first) check(err) secondDate, err := time.Parse(common.DateFormat, *second) check(err) if *bottleneck == "" { log.Fatalln("Missing bottleneck name") } check(process(*bottleneck, firstDate, secondDate)) }