--- a/pkg/controllers/diff.go	Sat Oct 19 23:07:59 2019 +0200
+++ b/pkg/controllers/diff.go	Mon Oct 21 02:01:56 2019 +0200
@@ -21,7 +21,6 @@
 	"log"
 	"net/http"
 	"runtime"
-	"sync"
 	"time"
 
 	"gemma.intevation.de/gemma/pkg/auth"
@@ -92,6 +91,22 @@
     4326
   )
 `
+	commonDiffBBoxSQL = `
+WITH joined AS (
+  SELECT
+    sr.area      AS area,
+    sr.date_info AS date_info
+  FROM waterway.sounding_results sr
+  WHERE sr.bottleneck_id = $1
+),
+bbox AS (
+  SELECT ST_Extent(ST_intersection(
+    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
+    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
+  )) AS area
+)
+SELECT ST_XMin(area), ST_YMin(area), ST_XMax(area), ST_YMax(area) FROM bbox
+`
 )
 
 type (
@@ -162,6 +177,38 @@
 			dci.Minuend.Format(common.DateFormat))
 	}
 
+	epsg := minuendTree.EPSG()
+
+	var box octree.Box2D
+
+	switch err := conn.QueryRowContext(
+		ctx,
+		commonDiffBBoxSQL,
+		dci.Bottleneck,
+		epsg,
+		dci.Minuend.Time,
+		dci.Subtrahend.Time,
+	).Scan(&box.X1, &box.Y1, &box.X2, &box.Y2); {
+	case err == sql.ErrNoRows:
+		return 0, errors.New("No such intersection")
+	case err != nil:
+		return 0, err
+	}
+
+	if box.Empty() {
+		return 0, errors.New("Intersection is empty")
+	}
+
+	log.Printf("info: bbox of intersection: (%.2f, %.2f) - (%.2f, %.2f)\n",
+		box.X1, box.Y1, box.X2, box.Y2)
+
+	start = time.Now()
+	raster := octree.NewRaster(box, isoCellSize)
+	raster.Rasterize(minuendTree.Value)
+	log.Printf("info: rasterizing minuend took %v\n", time.Since(start))
+
+	minuendTree = nil
+
 	start = time.Now()
 
 	subtrahendTree, err := octree.FromCache(
@@ -180,85 +227,15 @@
 	}
 
 	// We need a slow path implementation for this.
-	epsg := minuendTree.EPSG()
 	if epsg != subtrahendTree.EPSG() {
 		return 0, errors.New("Calculating differences between two different " +
 			"EPSG code meshes are not supported, yet.")
 	}
 
 	start = time.Now()
-	points := minuendTree.Diff(subtrahendTree)
+	raster.Diff(subtrahendTree.Value)
 	log.Printf("info: A - B took %v\n", time.Since(start))
-
-	minuendTree, subtrahendTree = nil, nil
-
-	// The Triangulation and the loading of the clipping
-	// polygon can be done concurrently.
-
-	jobs := make(chan func())
-
-	wg := new(sync.WaitGroup)
-	for i := 0; i < 2; i++ {
-		wg.Add(1)
-		go func() {
-			defer wg.Done()
-			for job := range jobs {
-				job()
-			}
-		}()
-	}
-
-	var (
-		tri     *octree.Triangulation
-		triErr  error
-		clip    *octree.Polygon
-		clipErr error
-	)
-
-	jobs <- func() {
-		start := time.Now()
-		tri, triErr = points.Triangulate()
-		log.Printf("info: triangulation took %v\n", time.Since(start))
-	}
-
-	jobs <- func() {
-		start := time.Now()
-		clip, clipErr = octree.LoadClippingPolygon(
-			ctx, conn,
-			epsg,
-			dci.Bottleneck,
-			dci.Minuend.Time,
-			dci.Subtrahend.Time)
-		log.Printf("info: loading clipping polygon took %v\n", time.Since(start))
-	}
-	close(jobs)
-	wg.Wait()
-
-	switch {
-	case triErr != nil && clipErr != nil:
-		return 0, fmt.Errorf("%v %v", triErr, clipErr)
-	case triErr != nil:
-		return 0, triErr
-	case clipErr != nil:
-		return 0, clipErr
-	}
-
-	start = time.Now()
-	tin := tri.Tin()
-	removed := tin.Clip(clip)
-	clip = nil
-	log.Printf("info: clipping TIN took %v\n", time.Since(start))
-
-	log.Printf("info: Number of triangles to clip: %d\n", len(removed))
-
-	start = time.Now()
-	var tree octree.STRTree
-
-	tree.BuildWithout(tin, removed)
-
-	log.Printf("info: Building final mesh took: %v\n", time.Since(start))
-
-	start = time.Now()
+	subtrahendTree = nil
 
 	// XXX: Maybe we should start this transaction earlier!?
 	var tx *sql.Tx
@@ -267,6 +244,13 @@
 	}
 	defer tx.Rollback()
 
+	zMin, zMax, ok := raster.ZExtent()
+	if !ok {
+		return 0, errors.New("Scans do not have common points")
+	}
+
+	log.Printf("info: z range: %.3f - %.3f\n", zMin, zMax)
+
 	var heights []float64
 
 	heights, err = octree.LoadClassBreaks(
@@ -275,12 +259,12 @@
 	if err != nil {
 		log.Printf("warn: Loading class breaks failed: %v\n", err)
 		err = nil
-		heights = octree.SampleDiffHeights(tin.Min.Z, tin.Max.Z, contourStep)
+		heights = octree.SampleDiffHeights(zMin, zMax, contourStep)
 	} else {
-		heights = octree.ExtrapolateClassBreaks(heights, tin.Min.Z, tin.Max.Z)
+		heights = octree.ExtrapolateClassBreaks(heights, zMin, zMax)
 	}
 
-	log.Printf("info: z range: %.3f - %.3f\n", tin.Min.Z, tin.Max.Z)
+	heights = common.DedupFloat64s(heights)
 
 	log.Printf("info: num heights: %d\n", len(heights))
 
@@ -302,9 +286,9 @@
 		return 0, err
 	}
 
-	heights = common.DedupFloat64s(heights)
+	areas := raster.Trace(heights)
 
-	areas := octree.TraceAreas(heights, isoCellSize, tin.Min, tin.Max, tree.Value)
+	raster = nil
 
 	var size int
 
@@ -327,7 +311,7 @@
 	log.Printf("info: Transferred WKB size: %.2fMB.\n",
 		float64(size)/(1024*1024))
 
-	log.Printf("info: calculating and storing iso lines took %v\n",
+	log.Printf("info: calculating and storing iso areas took %v\n",
 		time.Since(start))
 
 	if err = tx.Commit(); err != nil {