view cmd/isoareas/main.go @ 4533:3998a9ab69c6 iso-areas

Find out which triangles are fully inside classes.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 25 Sep 2019 13:04:08 +0200
parents 769f723c2581
children b7d31a449dd2
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 int32
	}

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

	classes := make([][]int32, len(heights)+1)

all:
	for i := int32(0); i < int32(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))

		if min > heights[len(heights)-1] {
			classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
			continue
		}
		max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z))

		for j, h := range heights {

			var l float64
			if j > 0 {
				l = heights[j-1]
			} else {
				l = -math.MaxFloat64
			}

			if l < min && max < h {
				classes[j] = append(classes[j], i/3)
				continue all
			}
			if min > h || max < h {
				continue
			}
			cut := t.IntersectHorizontal(h)
			if len(cut) >= 2 {
				cuts[j] = append(cuts[j], indexedCut{cut, i / 3})
			}
		}
	}

	log.Println("inside classes:")
	for _, c := range classes {
		log.Printf("\t%d\n", len(c))
	}

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

	// TODO: sew segments together.

}

func main() {

	heights, err := octree.ParseClassBreaks(classBreaks)
	check(err)
	log.Printf("num classes: %d\n", len(heights))

	start := time.Now()
	xyz, err := loadXYZ(os.Stdin)
	check(err)
	log.Printf("loading took %v\n", time.Since(start))

	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)

	start = time.Now()
	tri, err := octree.Triangulate(xyz)
	check(err)
	log.Printf("triangulation took %v\n", time.Since(time.Now()))

	log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)

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