--- a/pkg/controllers/cross.go	Thu Sep 27 12:49:45 2018 +0200
+++ b/pkg/controllers/cross.go	Thu Sep 27 13:18:50 2018 +0200
@@ -1,86 +1,141 @@
 package controllers
 
 import (
+	"context"
 	"database/sql"
 	"log"
 	"net/http"
 	"time"
 
 	"gemma.intevation.de/gemma/pkg/models"
+	"gemma.intevation.de/gemma/pkg/octree"
 )
 
-const crossSQL = `
-WITH line AS (
-SELECT ST_LineMerge(ST_Union(ST_3DIntersection(
-  ST_Translate(
-    ST_Extrude(
-      ST_GeomFromWKB($1, 4326),
-    0, 0, 1000),
-   0, 0, -500),
-   geom))) AS geom
-FROM waterway.meshes m JOIN waterway.sounding_results sr ON m.sounding_result_id = sr.id
-WHERE ST_Intersects(geom, ST_GeomFromWKB($1, 4326)) AND
-  sr.bottleneck_id = $2 AND sr.date_info = $3
-)
-SELECT ST_AsBinary((ST_Dump(ST_Intersection(line.geom, sr.area::geometry))).geom)
-  FROM line, waterway.sounding_results sr
-  WHERE sr.bottleneck_id = $2 AND sr.date_info = $3
-`
+const WGS84 = 4326
+
+func reproject(
+	rp *models.Reprojector,
+	src models.GeoJSONLineCoordinates,
+	ctx context.Context,
+) (models.GeoJSONLineCoordinates, error) {
+
+	dst := make(models.GeoJSONLineCoordinates, len(src))
+	for i, s := range src {
+		var err error
+		if dst[i].Lat, dst[i].Lon, err = rp.Reproject(
+			s.Lat, s.Lon, ctx,
+		); err != nil {
+			return nil, err
+		}
+	}
+	return dst, nil
+}
+
+func projectBack(
+	rp *models.Reprojector,
+	lines octree.MultiLineStringZ,
+	ctx context.Context,
+) (models.GeoJSONMultiLineCoordinatesZ, error) {
+
+	out := make(models.GeoJSONMultiLineCoordinatesZ, len(lines))
+
+	for i, segment := range lines {
+
+		coords := make(models.GeoJSONLineCoordinatesZ, len(segment))
+		out[i] = coords
+		for j, v := range segment {
+			lat, lon, err := rp.Reproject(v.X, v.Y, ctx)
+			if err != nil {
+				return nil, err
+			}
+			coords[j] = models.GeoJSONCoordinateZ{Lat: lat, Lon: lon, Z: v.Z}
+		}
+	}
+	return out, nil
+}
 
 func crossSection(
 	input interface{},
 	req *http.Request,
-	db *sql.Conn,
+	conn *sql.Conn,
 ) (jr JSONResult, err error) {
 
 	csi := input.(*models.CrossSectionInput)
 
 	start := time.Now()
 
-	var rows *sql.Rows
-	if rows, err = db.QueryContext(
-		req.Context(),
-		crossSQL,
-		csi.Geometry.Coordinates.AsWKB(),
-		csi.Properties.Bottleneck,
-		csi.Properties.Date.Time,
+	tree, err := octree.Cache.Get(
+		csi.Properties.Bottleneck, csi.Properties.Date.Time,
+		conn, req.Context())
+
+	log.Printf("loading octree took: %s\n", time.Since(start))
+	if err != nil {
+		return
+	}
+
+	// The coordinate system of the octree is an UTM projection.
+	// The input coordinates are in WGS84.
+	// So we need to reproject them.
+
+	start = time.Now()
+
+	var rp *models.Reprojector
+	if rp, err = models.NewReprojector(
+		conn, req.Context(),
+		WGS84, tree.EPSG,
 	); err != nil {
 		return
 	}
-	defer rows.Close()
-
-	var segments models.GeoJSONMultiLineCoordinatesZ
+	defer rp.Close()
 
-	for rows.Next() {
-		var segment models.GeoJSONLineCoordinatesZ
-		if err = rows.Scan(&segment); err != nil {
-			return
-		}
-		segments = append(segments, segment)
-	}
+	coords, err := reproject(rp, csi.Geometry.Coordinates, req.Context())
 
-	if err = rows.Err(); err != nil {
+	log.Printf("transforming input coords took: %s\n", time.Since(start))
+	if err != nil {
 		return
 	}
 
-	log.Printf("query took: %v\n", time.Since(start))
-
 	start = time.Now()
 
-	joined := segments.Join()
+	var segments octree.MultiLineStringZ
+
+	for i := 0; i < len(coords)-1; i++ {
+		c1 := &coords[i]
+		c2 := &coords[i+1]
 
-	log.Printf("joining took: %v\n", time.Since(start))
-	log.Printf("segments before/after: %d %d\n", len(segments), len(joined))
+		verticalLine := octree.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon)
+
+		var line octree.MultiLineStringZ
+		tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *octree.Triangle) {
+			if ls := verticalLine.Intersection(t); len(ls) > 0 {
+				line = append(line, ls)
+			}
+		})
 
-	/*
-		if err2 := dumpProfile(
-			csi.Geometry.Coordinates[0],
-			csi.Geometry.Coordinates[len(csi.Geometry.Coordinates)-1],
-			joined,
-		); err2 != nil {
-			log.Printf("error: %v\n", err2)
+		if len(line) > 0 {
+			log.Printf("line length: %d\n", len(line))
+			// They are all on the segment (c1.Lat, c1.Lon) - (c2.Lat, c2.Lon).
+			// Sort them by project them on this line.
+			joined := line.JoinOnLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon)
+			log.Printf("joined length: %d\n", len(joined))
+			segments = append(segments, joined...)
 		}
-	*/
+
+	}
+	log.Printf("octree traversal took: %s\n", time.Since(start))
+
+	// The result have to be WGS84. So project the result back.
+	start = time.Now()
+
+	rp.FromEPSG, rp.ToEPSG = tree.EPSG, WGS84
+
+	var joined models.GeoJSONMultiLineCoordinatesZ
+	joined, err = projectBack(rp, segments, req.Context())
+
+	log.Printf("projecting back took: %s\n", time.Since(start))
+	if err != nil {
+		return
+	}
 
 	jr = JSONResult{
 		Result: &models.CrossSectionOutput{
@@ -95,75 +150,3 @@
 
 	return
 }
-
-/*
-func dumpProfile(
-	start models.GeoJSONCoordinate,
-	stop models.GeoJSONCoordinate,
-	segments models.GeoJSONMultiLineCoordinatesZ,
-) error {
-
-	w, err := os.Create("/tmp/cross.txt")
-	if err != nil {
-		return err
-	}
-	defer w.Close()
-
-	begin := models.GeoJSONCoordinateZ{
-		Lat: start.Lat,
-		Lon: start.Lon,
-	}
-	end := models.GeoJSONCoordinateZ{
-		Lat: stop.Lat,
-		Lon: stop.Lon,
-	}
-
-	last := begin
-
-	var total float64
-
-	minz := 10000.0
-
-	for _, line := range segments {
-		for _, coord := range line {
-			if coord.Z < minz {
-				minz = coord.Z
-			}
-			total += last.Distance(coord)
-			last = coord
-		}
-	}
-
-	log.Printf("length: %.3f minz:  %f\n", total, minz)
-
-	var pos float64
-
-	for i, line := range segments {
-		for j, coord := range line {
-			if i == 0 && j == 0 {
-				if !begin.Equals(coord) {
-					pos = begin.Distance(coord)
-					fmt.Fprintf(w, "%.3f %f\n", 0.0, 200.0)
-					fmt.Fprintf(w, "%.3f %f\n", pos-0.0001, 200.0)
-					fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz)
-					continue
-				}
-			} else if j == 0 {
-				fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0)
-				fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(coord)-0.0001, 200.0)
-				continue
-			}
-			pos += last.Distance(coord)
-			fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz)
-			last = coord
-		}
-	}
-
-	if !last.Equals(end) {
-		fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0)
-		fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(end), 200.0)
-	}
-
-	return nil
-}
-*/