Mercurial > gemma
changeset 4552:7ed5a4a94efc iso-areas
Used fogleman's contourmap as a tracing alternative.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 01 Oct 2019 13:18:16 +0200 |
parents | b5b23b6d76f1 |
children | 4c476d65d1bb |
files | cmd/isoareas/main.go go.mod go.sum |
diffstat | 3 files changed, 209 insertions(+), 23 deletions(-) [+] |
line wrap: on
line diff
--- a/cmd/isoareas/main.go Tue Oct 01 11:07:33 2019 +0200 +++ b/cmd/isoareas/main.go Tue Oct 01 13:18:16 2019 +0200 @@ -15,17 +15,21 @@ import ( "bufio" + "flag" "fmt" "io" "log" "math" "math/rand" "os" + "runtime" "strconv" "strings" + "sync" "time" svg "github.com/ajstarks/svgo" + "github.com/fogleman/contourmap" "gemma.intevation.de/gemma/pkg/octree" ) @@ -66,21 +70,6 @@ return points, nil } -func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) { - if len(points) == 0 { - return octree.Vertex{}, octree.Vertex{} - } - - min, max := points[0], points[0] - - for _, v := range points { - min.Minimize(v) - max.Maximize(v) - } - - return min, max -} - func check(err error) { if err != nil { log.Fatalf("error: %v\n", err) @@ -108,6 +97,77 @@ } } +func dumpContoursSVG( + w io.Writer, + contours [][]contourmap.Contour, +) (err error) { + + var ( + minX = math.MaxFloat64 + minY = math.MaxFloat64 + maxX = -math.MaxFloat64 + maxY = -math.MaxFloat64 + ) + + for _, c := range contours { + for _, r := range c { + for _, p := range r { + if p.X < minX { + minX = p.X + } + if p.Y < minY { + minY = p.Y + } + if p.X > maxX { + maxX = p.X + } + if p.Y > maxY { + maxY = p.Y + } + } + } + } + ratio := (maxX - minX) / (maxY - minY) + + log.Printf("ratio: %.2f\n", ratio) + + const width = 50000 + height := int(math.Ceil(width * ratio)) + + px := linear(minX, 0, maxX, width) + py := linear(minY, float64(height), maxY, 0) + + out := bufio.NewWriter(w) + defer func() { err = out.Flush() }() + + canvas := svg.New(out) + canvas.Start(width, height) + + rnd := rand.New(rand.NewSource(42)) + + for _, c := range contours { + + r := byte(rnd.Intn(256)) + g := byte(rnd.Intn(256)) + b := byte(rnd.Intn(256)) + style := fmt.Sprintf("fill:#%02x%02x%02x", r, g, b) + + for _, r := range c { + x := make([]int, len(r)) + y := make([]int, len(r)) + for i, p := range r { + x[i] = int(math.Round(px(p.X))) + y[i] = int(math.Round(py(p.Y))) + } + + canvas.Polygon(x, y, style) + } + } + + canvas.End() + return +} + func dumpPolygonsSVG( w io.Writer, min, max octree.Vertex, @@ -218,6 +278,12 @@ func main() { + cellSize := float64(0.5) + + flag.Float64Var(&cellSize, "cellsize", 0.5, "width/height of raster cell [m]") + + flag.Parse() + heights, err := octree.ParseClassBreaks(classBreaks) check(err) log.Printf("num classes: %d\n", len(heights)) @@ -229,22 +295,139 @@ log.Printf("num vertices: %d\n", len(xyz)) - min, max := minMax(xyz) + start = time.Now() + tri, err := octree.Triangulate(xyz) + check(err) + log.Printf("triangulation took %v\n", time.Since(start)) + tooLongEdge := tri.EstimateTooLong() + log.Printf("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) + + start = time.Now() + polygon, removed := tri.ConcaveHull(tooLongEdge) + + polygonArea := polygon.Area() + if polygonArea < 0.0 { // counter clockwise + polygonArea = -polygonArea + polygon.Reverse() + } + + var clippingPolygon octree.Polygon + clippingPolygon.FromLineStringZ(polygon) + clippingPolygon.Indexify() + + tin := tri.Tin() + // tin.EPSG = epsg + + var str octree.STRTree + str.Build(tin) + + removed = str.Clip(&clippingPolygon) + + builder := octree.NewBuilder(tin) + builder.Build(removed) + + tree := builder.Tree() + + min, max := tin.Min, tin.Max log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n", min.X, min.Y, min.Z, max.X, max.Y, max.Z) + width := max.X - min.X + height := max.Y - min.Y + + log.Printf("width/height: %.2f / %.2f\n", width, height) + + xcells := int(math.Ceil(width / cellSize)) + ycells := int(math.Ceil(height / cellSize)) + + log.Printf("raster size: (%d, %d)\n", xcells, ycells) + + raster := make([]float64, xcells*ycells) + + const closed = -math.MaxFloat64 + for i := range raster { + raster[i] = closed + } + start = time.Now() - tri, err := octree.Triangulate(xyz) - check(err) - log.Printf("triangulation took %v\n", time.Since(start)) + + var wg sync.WaitGroup + + rows := make(chan int) - log.Printf("number of triangles: %d\n", len(tri.Triangles)/3) + rasterRow := func() { + defer wg.Done() + for i := range rows { + pos := i * xcells + row := raster[pos : pos+xcells] + //log.Printf("len: %d\n", len(row)) + py := min.Y + float64(i)*cellSize + cellSize/2 + px := min.X + cellSize/2 + for j := range row { + v, ok := tree.Value(px, py) + if ok { + row[j] = v + } + px += cellSize + } + } + } start = time.Now() - polygons := intersectTriangles(tri, heights, min, max) - log.Printf("intersecting triangles took %v.\n", time.Since(start)) + + for n := runtime.NumCPU(); n >= 1; n-- { + wg.Add(1) + go rasterRow() + } + + for i := 0; i < ycells; i++ { + rows <- i + } + close(rows) + + wg.Wait() + + log.Printf("Rasterizing took %v.\n", time.Since(start)) + + start = time.Now() + cm := contourmap.FromFloat64s(xcells, ycells, raster).Closed() + + start = time.Now() + + contours := make([][]contourmap.Contour, len(heights)) + + cnts := make(chan int) - check(dumpPolygonsSVG(os.Stdout, min, max, polygons)) + doContours := func() { + defer wg.Done() + for i := range cnts { + contours[i] = cm.Contours(heights[i]) + } + } + + for n := runtime.NumCPU(); n >= 1; n-- { + wg.Add(1) + go doContours() + } + + for i := range heights { + cnts <- i + } + close(cnts) + + wg.Wait() + + log.Printf("Calculating contours took %v.\n", time.Since(start)) + check(dumpContoursSVG(os.Stdout, contours)) + + /* + + start = time.Now() + polygons := intersectTriangles(tri, heights, min, max) + log.Printf("intersecting triangles took %v.\n", time.Since(start)) + + check(dumpPolygonsSVG(os.Stdout, min, max, polygons)) + */ }
--- a/go.mod Tue Oct 01 11:07:33 2019 +0200 +++ b/go.mod Tue Oct 01 13:18:16 2019 +0200 @@ -6,6 +6,7 @@ github.com/ajstarks/svgo v0.0.0-20180226025133-644b8db467af github.com/cockroachdb/apd v1.1.0 // indirect github.com/etcd-io/bbolt v1.3.3 + github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199 github.com/gofrs/uuid v3.2.0+incompatible // indirect github.com/golang/snappy v0.0.1 github.com/gorilla/mux v1.7.3
--- a/go.sum Tue Oct 01 11:07:33 2019 +0200 +++ b/go.sum Tue Oct 01 13:18:16 2019 +0200 @@ -26,6 +26,8 @@ github.com/dgryski/go-sip13 v0.0.0-20181026042036-e10d5fee7954/go.mod h1:vAd38F8PWV+bWy6jNmig1y/TA+kYO4g3RSRF0IAv0no= github.com/etcd-io/bbolt v1.3.3 h1:gSJmxrs37LgTqR/oyJBWok6k6SvXEUerFTbltIhXkBM= github.com/etcd-io/bbolt v1.3.3/go.mod h1:ZF2nL25h33cCyBtcyWeZ2/I3HQOfTP+0PIEvHjkjCrw= +github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199 h1:kufr0u0RIG5ACpjFsPRbbuHa0FhMWsS3tnSFZ2hf07s= +github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199/go.mod h1:mqaaaP4j7nTF8T/hx5OCljA7BYWHmrH2uh+Q023OchE= github.com/fogleman/gg v1.2.1-0.20190220221249-0403632d5b90/go.mod h1:R/bRT+9gY/C5z7JzPU0zXsXHKM4/ayA+zqcVNZzPa1k= github.com/fsnotify/fsnotify v1.4.7 h1:IXs+QLmnXW2CcXuY+8Mzv/fWEsPGWxqefPtCP5CnV9I= github.com/fsnotify/fsnotify v1.4.7/go.mod h1:jwhsz4b93w/PPRr/qN1Yymfu8t87LnFCMoQvtojpjFo=