# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1561373665 -7200
# Node ID 879c297c47e950b8ebea6a358bb564d50c447767
# Parent  ec86a71553775c9203eb1dda5a5cae7e066a455d
SR import: Use home-brew concave hull algorithm in multi-beam scan case, too.

diff -r ec86a7155377 -r 879c297c47e9 pkg/imports/sr.go
--- a/pkg/imports/sr.go	Mon Jun 24 11:39:09 2019 +0200
+++ b/pkg/imports/sr.go	Mon Jun 24 12:54:25 2019 +0200
@@ -306,20 +306,22 @@
 	var summary interface{}
 
 	if sr.isSingleBeam() {
-		summary, err = sr.singleBeamScan(
+		summary, err = sr.processScan(
 			ctx,
 			tx,
 			feedback,
+			true,
 			importID,
 			m,
 			xyz,
 			boundary,
 		)
 	} else {
-		summary, err = sr.multiBeamScan(
+		summary, err = sr.processScan(
 			ctx,
 			tx,
 			feedback,
+			false,
 			importID,
 			m,
 			xyz,
@@ -342,24 +344,30 @@
 	return summary, nil
 }
 
-func (sr *SoundingResult) singleBeamScan(
+func (sr *SoundingResult) processScan(
 	ctx context.Context,
 	tx *sql.Tx,
 	feedback Feedback,
+	singleBeam bool,
 	importID int64,
 	m *models.SoundingResultMeta,
 	xyz octree.MultiPointZ,
 	boundary polygonSlice,
 ) (interface{}, error) {
 
-	feedback.Info("Processing as single beam scan.")
+	if singleBeam {
+		feedback.Info("Processing as single beam scan.")
+	} else {
+		feedback.Info("Processing as multi beam scan.")
+	}
+
 	feedback.Info("Reproject XYZ data.")
 
 	start := time.Now()
 
 	xyzWKB := xyz.AsWKB()
 	var reproj []byte
-	var epsg uint
+	var epsg uint32
 
 	if err := tx.QueryRowContext(
 		ctx,
@@ -393,6 +401,7 @@
 		removed                 map[int32]struct{}
 		polygonArea             float64
 		clippingPolygonWKB      []byte
+		tin                     *octree.Tin
 	)
 
 	if boundary == nil {
@@ -451,7 +460,7 @@
 		}
 
 		tin := tri.Tin()
-		tin.EPSG = uint32(epsg)
+		tin.EPSG = epsg
 
 		var str octree.STRTree
 		str.Build(tin)
@@ -459,68 +468,72 @@
 		removed = str.Clip(&clippingPolygon)
 	}
 
-	// Build the first mesh to generate random points on.
+	if singleBeam {
+
+		// Build the first mesh to generate random points on.
+
+		feedback.Info("Build virtual DEM based on original XYZ data.")
+
+		start = time.Now()
 
-	feedback.Info("Build virtual DEM based on original XYZ data.")
+		var tree *octree.Tree
+		{
+			builder := octree.NewBuilder(tri.Tin())
+			builder.Build(removed)
+			tree = builder.Tree()
+		}
+
+		feedback.Info("Building took %v", time.Since(start))
 
-	start = time.Now()
+		feedback.Info("Boundary area: %.2fm²", polygonArea)
+
+		numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
+
+		feedback.Info("Generate %d random points for an average density of ~%d points/m².",
+			numPoints, pointsPerSquareMeter)
+
+		start = time.Now()
+
+		generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0))
 
-	var tree *octree.Tree
-	{
-		builder := octree.NewBuilder(tri.Tin())
-		builder.Build(removed)
-		tree = builder.Tree()
-	}
+		tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
+			generated = append(generated, vertices...)
+		})
+
+		feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
 
-	feedback.Info("Building took %v", time.Since(start))
-
-	feedback.Info("Boundary area: %.2fm²", polygonArea)
+		// Add the boundary to new point cloud.
+		dupes := map[[2]float64]struct{}{}
+		clippingPolygon.Vertices(0, func(x, y float64) {
+			key := [2]float64{x, y}
+			if _, found := dupes[key]; found {
+				return
+			}
+			dupes[key] = struct{}{}
+			if z, ok := tree.Value(x, y); ok {
+				generated = append(generated, octree.Vertex{X: x, Y: y, Z: z})
+			}
+		})
 
-	numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
+		feedback.Info("Triangulate new point cloud.")
+		xyz = octree.MultiPointZ(generated)
+		start = time.Now()
 
-	feedback.Info("Generate %d random points for an average density of ~%d points/m².",
-		numPoints, pointsPerSquareMeter)
+		tri, err = octree.Triangulate(xyz)
+		if err != nil {
+			return nil, err
+		}
+		feedback.Info("Second triangulation took %v.", time.Since(start))
+		feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
+		feedback.Info("Clipping triangles from new mesh.")
+
+	} else { // multi beam
+		// Nothing special
+	}
 
 	start = time.Now()
-
-	generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0))
-
-	tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
-		generated = append(generated, vertices...)
-	})
-
-	feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
-
-	// Add the boundary to new point cloud.
-	dupes := map[[2]float64]struct{}{}
-	clippingPolygon.Vertices(0, func(x, y float64) {
-		key := [2]float64{x, y}
-		if _, found := dupes[key]; found {
-			return
-		}
-		dupes[key] = struct{}{}
-		if z, ok := tree.Value(x, y); ok {
-			generated = append(generated, octree.Vertex{X: x, Y: y, Z: z})
-		}
-	})
-
-	feedback.Info("Triangulate new point cloud.")
-
-	start = time.Now()
-
-	xyz = octree.MultiPointZ(generated)
-
-	tri, err = octree.Triangulate(xyz)
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Second triangulation took %v.", time.Since(start))
-	feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
-	feedback.Info("Clipping triangles from new mesh.")
-
-	start = time.Now()
-	tin := tri.Tin()
-	tin.EPSG = uint32(epsg)
+	tin = tri.Tin()
+	tin.EPSG = epsg
 
 	var str octree.STRTree
 	str.Build(tin)
@@ -533,8 +546,6 @@
 	feedback.Info("Clipping STR tree took %v.", time.Since(start))
 	feedback.Info("Number of triangles to clip %d.", len(removed))
 
-	start = time.Now()
-
 	feedback.Info("Build final octree index")
 
 	builder := octree.NewBuilder(tin)
@@ -619,160 +630,6 @@
 	return &summary, nil
 }
 
-func (sr *SoundingResult) multiBeamScan(
-	ctx context.Context,
-	tx *sql.Tx,
-	feedback Feedback,
-	importID int64,
-	m *models.SoundingResultMeta,
-	xyz octree.MultiPointZ,
-	boundary polygonSlice,
-) (interface{}, error) {
-
-	feedback.Info("Processing as multi beam scan.")
-	var (
-		id       int64
-		epsg     uint32
-		lat, lon float64
-	)
-	start := time.Now()
-
-	var hull []byte
-
-	xyzWKB := xyz.AsWKB()
-
-	err := tx.QueryRowContext(
-		ctx,
-		insertHullSQL,
-		m.Bottleneck,
-		m.Date.Time,
-		m.DepthReference,
-		xyzWKB,
-		boundary.asWKB(),
-		m.EPSG,
-	).Scan(
-		&id,
-		&lat,
-		&lon,
-		&epsg,
-		&hull,
-	)
-	xyz, boundary = nil, nil // not need from now on.
-	feedback.Info("Calculating hull took %s.", time.Since(start))
-	switch {
-	case err == sql.ErrNoRows:
-		return nil, fmt.Errorf(
-			"No matching bottleneck of given name or time available: %v", err)
-	case err != nil:
-		return nil, err
-	}
-	feedback.Info("Best suited UTM EPSG: %d", epsg)
-
-	start = time.Now()
-
-	var clippingPolygon octree.Polygon
-
-	if err := clippingPolygon.FromWKB(hull); err != nil {
-		return nil, err
-	}
-	clippingPolygon.Indexify()
-	feedback.Info("Building clipping polygon took %v.", time.Since(start))
-
-	start = time.Now()
-
-	var reproj []byte
-
-	if err = tx.QueryRowContext(
-		ctx,
-		reprojectPointsSQL,
-		xyzWKB,
-		m.EPSG,
-		epsg,
-	).Scan(&reproj); err != nil {
-		return nil, err
-	}
-
-	if err := xyz.FromWKB(reproj); err != nil {
-		return nil, err
-	}
-
-	feedback.Info("Reprojecting points took %v.", time.Since(start))
-	feedback.Info("Number of reprojected points: %d", len(xyz))
-
-	start = time.Now()
-
-	tri, err := octree.Triangulate(xyz)
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Triangulation took %v.", time.Since(start))
-
-	start = time.Now()
-
-	tin := tri.Tin()
-
-	var str octree.STRTree
-	str.Build(tin)
-	feedback.Info("Building STR tree took %v", time.Since(start))
-
-	start = time.Now()
-
-	removed := str.Clip(&clippingPolygon)
-	feedback.Info("Clipping STR tree took %v.", time.Since(start))
-	feedback.Info("Number of triangles to clip %d.", len(removed))
-
-	start = time.Now()
-
-	tin.EPSG = epsg
-
-	builder := octree.NewBuilder(tin)
-	builder.Build(removed)
-	octreeIndex, err := builder.Bytes()
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Building octree took %v.", time.Since(start))
-
-	start = time.Now()
-	h := sha1.New()
-	h.Write(octreeIndex)
-	checksum := hex.EncodeToString(h.Sum(nil))
-	_, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Storing octree index took %s.", time.Since(start))
-
-	tree := builder.Tree()
-
-	start = time.Now()
-	err = generateContours(ctx, tx, tree, id)
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Generating and storing contour lines took %s.",
-		time.Since(start))
-
-	// Store for potential later removal.
-	if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
-		return nil, err
-	}
-
-	summary := struct {
-		Bottleneck string      `json:"bottleneck"`
-		Date       models.Date `json:"date"`
-		Lat        float64     `json:"lat"`
-		Lon        float64     `json:"lon"`
-	}{
-		Bottleneck: m.Bottleneck,
-		Date:       m.Date,
-		Lat:        lat,
-		Lon:        lon,
-	}
-
-	return &summary, nil
-}
-
 // CleanUp removes the folder containing the ZIP file with the
 // the sounding result import.
 func (sr *SoundingResult) CleanUp() error {