view cmd/isoareas/main.go @ 4532:769f723c2581 iso-areas

Cut triangles against class breaks.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 24 Sep 2019 18:03:43 +0200
parents c9b6be8d81c8
children 3998a9ab69c6
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) 2019 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 (
	"bufio"
	"io"
	"log"
	"math"
	"os"
	"strconv"
	"strings"
	"time"

	"gemma.intevation.de/gemma/pkg/octree"
)

const classBreaks = `1:#ff00dd,1.5,1.7,1.9,2.1,2.3,2.5:#f25f20,2.7,2.9,3.1:#f7e40e,3.3,3.5,4:#8ad51a,4.5,5,5.5,6,6.5,7:#1414ff`

func loadXYZ(r io.Reader) (octree.MultiPointZ, error) {

	scanner := bufio.NewScanner(r)

	points := make(octree.MultiPointZ, 0, 2000000)

	var x, y, z float64
	var err error

	for scanner.Scan() {
		line := strings.TrimSpace(scanner.Text())
		if len(line) == 0 || strings.HasPrefix(line, "#") {
			continue
		}
		parts := strings.SplitN(line, " ", 3)
		if len(parts) != 3 {
			continue
		}

		if x, err = strconv.ParseFloat(parts[0], 64); err != nil {
			return nil, err
		}
		if y, err = strconv.ParseFloat(parts[1], 64); err != nil {
			return nil, err
		}
		if z, err = strconv.ParseFloat(parts[2], 64); err != nil {
			return nil, err
		}
		points = append(points, octree.Vertex{X: x, Y: y, Z: z})
	}

	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)
	}
}

func intersectTriangles(tri *octree.Triangulation, heights []float64) {

	type indexedCut struct {
		cut   octree.LineStringZ
		index int
	}

	cuts := make([][]indexedCut, len(heights))

	for i := 0; i < len(tri.Triangles); i += 3 {
		ti := tri.Triangles[i : i+3]
		t := octree.Triangle{
			tri.Points[ti[0]],
			tri.Points[ti[1]],
			tri.Points[ti[2]],
		}
		min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z))
		max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z))

		for j, h := range heights {
			if h < min {
				continue
			}
			if h > max {
				break
			}
			cut := t.IntersectHorizontal(h)
			if len(cut) < 2 {
				continue
			}
			cuts[j] = append(cuts[j], indexedCut{cut, i})
		}
	}

	log.Println("cuts")
	for i := range cuts {
		log.Printf("%.3f: %d\n", heights[i], len(cuts[i]))
	}

	// TODO: sew segments together.

}

func main() {

	heights, err := octree.ParseClassBreaks(classBreaks)
	check(err)

	xyz, err := loadXYZ(os.Stdin)
	check(err)

	log.Printf("num vertices: %d\n", len(xyz))

	min, max := minMax(xyz)

	log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
		min.X, min.Y, min.Z,
		max.X, max.Y, max.Z)

	heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z)

	log.Printf("classes: %d\n", len(heights))

	tri, err := octree.Triangulate(xyz)
	check(err)

	start := time.Now()
	intersectTriangles(tri, heights)
	log.Printf("intersecting triangles took %v.\n", time.Since(start))
}