--- a/pkg/imports/isr.go	Sat Jul 10 01:09:57 2021 +0200
+++ b/pkg/imports/isr.go	Sun Jul 11 13:08:49 2021 +0200
@@ -148,6 +148,7 @@
 	if err != nil {
 		return nil, err
 	}
+	heights = heights.Dedup()
 
 	bns, err := fetchBottleneckResults(ctx, conn)
 	if err != nil {
@@ -246,15 +247,7 @@
 			}
 
 			// Re-classify points.
-			classes := heights.Classify(points.All())
-
-			// Should not happen ... Z values over the top.
-			if n := len(classes) - 1; n > 1 && len(classes[n]) > 0 {
-				// Place the over the top values to the class below.
-				classes[n-1] = append(classes[n-1], classes[n]...)
-				classes[n] = nil
-				classes = classes[:n]
-			}
+			classes := heights.Classify(points)
 
 			// Re-insert points
 			for i, class := range classes {

--- a/pkg/imports/sr.go	Sat Jul 10 01:09:57 2021 +0200
+++ b/pkg/imports/sr.go	Sun Jul 11 13:08:49 2021 +0200
@@ -24,7 +24,6 @@
 	"errors"
 	"fmt"
 	"io"
-	"log"
 	"math"
 	"os"
 	"path"
@@ -645,12 +644,10 @@
 	} else {
 		start = time.Now()
 		clippingPolygonBuffered.Indexify()
-		removed = make(map[int32]struct{})
-		for i, v := range xyz {
-			if clippingPolygonBuffered.IntersectionBox2D(v.Box2D()) == mesh.IntersectionOutSide {
-				removed[int32(i)] = struct{}{}
-			}
-		}
+		xyz = xyz.Filter(func(v mesh.Vertex) bool {
+			return clippingPolygonBuffered.IntersectionBox2D(v.Box2D()) ==
+				mesh.IntersectionOutSide
+		})
 		feedback.Info("Clipping took %v.", time.Since(start))
 		feedback.Info("Number of points to clip %d.", len(removed))
 	}
@@ -707,17 +704,19 @@
 			return nil, err
 		}
 		feedback.Info("Storing mesh index took %s.", time.Since(start))
-		if err := generateIsos(ctx, tx, feedback, &final, id); err != nil {
-			return nil, err
-		}
-	} else { // SurveyTypeMarking
-		if err := generateMarkingPoints(
+		err = generateIsoAreas(
 			ctx, tx, feedback,
-			xyz, removed, epsg,
-			id,
-		); err != nil {
-			return nil, err
-		}
+			&final,
+			loadClassBreaks(ctx, tx, feedback, final.Min().Z, final.Max().Z),
+			id)
+	} else { // SurveyTypeMarking
+		err = generateMarkingPoints(
+			ctx, tx, feedback,
+			xyz, epsg,
+			id)
+	}
+	if err != nil {
+		return nil, err
 	}
 
 	// Store for potential later removal.
@@ -937,6 +936,50 @@
 	return heights
 }
 
+func generateMarkingPoints(
+	ctx context.Context,
+	tx *sql.Tx,
+	feedback Feedback,
+	xyz mesh.MultiPointZ,
+	epsg uint32,
+	id int64,
+) error {
+	heights, err := mesh.LoadClassBreaks(
+		ctx, tx,
+		"morphology_classbreaks")
+
+	if err != nil {
+		feedback.Warn("Loading class breaks failed: %v", err)
+		feedback.Info("Using default class breaks")
+		min, max := xyz.MinMax()
+		heights = defaultClassBreaks(min.Z, max.Z)
+	}
+
+	heights = heights.Dedup()
+
+	classes := heights.Classify(xyz)
+
+	stmt, err := tx.PrepareContext(ctx, insertMarkingPointsSQL)
+	if err != nil {
+		return err
+	}
+	defer stmt.Close()
+
+	for i, class := range classes {
+		// Ignore empty classes
+		if len(class) == 0 {
+			continue
+		}
+		if _, err := stmt.ExecContext(
+			ctx, id, heights[i], epsg, class.AsWKB(),
+		); err != nil {
+			return err
+		}
+	}
+
+	return nil
+}
+
 func loadClassBreaks(
 	ctx context.Context,
 	tx *sql.Tx,
@@ -959,68 +1002,6 @@
 	return heights.Dedup()
 }
 
-func generateMarkingPoints(
-	ctx context.Context,
-	tx *sql.Tx,
-	feedback Feedback,
-	xyz mesh.MultiPointZ,
-	removed map[int32]struct{},
-	epsg uint32,
-	id int64,
-) error {
-	log.Printf("debug: generateMarkingPoints")
-
-	min, max := mesh.MinMaxVertex(xyz.FilterRemoved(removed))
-
-	log.Printf("debug: min/max %.2f/%.2f\n", min.Z, max.Z)
-
-	heights := loadClassBreaks(ctx, tx, feedback, min.Z, max.Z)
-
-	classes := heights.Classify(xyz.FilterRemoved(removed))
-
-	// Should not happen ... Z values over the top.
-	if n := len(classes) - 1; n > 1 && len(classes[n]) > 0 {
-		// Place the over the top values to the class below.
-		classes[n-1] = append(classes[n-1], classes[n]...)
-		classes[n] = nil
-		classes = classes[:n]
-	}
-
-	stmt, err := tx.PrepareContext(ctx, insertMarkingPointsSQL)
-	if err != nil {
-		return err
-	}
-	defer stmt.Close()
-
-	for i, class := range classes {
-		// Ignore empty classes
-		if len(class) == 0 {
-			continue
-		}
-		log.Printf("debug: class %d: %d\n", i, len(class))
-		if _, err := stmt.ExecContext(
-			ctx, id, heights[i], epsg, class.AsWKB(),
-		); err != nil {
-			return err
-		}
-	}
-
-	return nil
-}
-
-func generateIsos(
-	ctx context.Context,
-	tx *sql.Tx,
-	feedback Feedback,
-	tree *mesh.STRTree,
-	id int64,
-) error {
-
-	heights := loadClassBreaks(ctx, tx, feedback, tree.Min().Z, tree.Max().Z)
-
-	return generateIsoAreas(ctx, tx, feedback, tree, heights, id)
-}
-
 func generateIsoAreas(
 	ctx context.Context,
 	tx *sql.Tx,

--- a/pkg/mesh/classbreaks.go	Sat Jul 10 01:09:57 2021 +0200
+++ b/pkg/mesh/classbreaks.go	Sun Jul 11 13:08:49 2021 +0200
@@ -157,10 +157,14 @@
 	return ClassBreaks(common.DedupFloat64s(cbs))
 }
 
-func (cbs ClassBreaks) Classify(points func() (Vertex, bool)) []MultiPointZ {
-	classes := make([]MultiPointZ, len(cbs)+1)
-	for v, ok := points(); ok; v, ok = points() {
-		idx := len(cbs)
+func (cbs ClassBreaks) Classify(points MultiPointZ) []MultiPointZ {
+	if len(cbs) == 0 {
+		return nil
+	}
+	classes := make([]MultiPointZ, len(cbs))
+	for _, v := range points {
+		// Place in last class if greater than all.
+		idx := len(cbs) - 1
 		for i, cb := range cbs {
 			if v.Z <= cb {
 				idx = i

--- a/pkg/mesh/vertex.go	Sat Jul 10 01:09:57 2021 +0200
+++ b/pkg/mesh/vertex.go	Sun Jul 11 13:08:49 2021 +0200
@@ -1144,47 +1144,23 @@
 	return out
 }
 
-// All returns all points as an iterator.
-func (mpz MultiPointZ) All() func() (Vertex, bool) {
-	var idx int
-	return func() (v Vertex, ok bool) {
-		if idx >= len(mpz) {
-			ok = false
-			return
+// Filter returns a copy removed the vertices which
+// don't pass the filter test.
+func (mpz MultiPointZ) Filter(filter func(Vertex) bool) MultiPointZ {
+	n := make(MultiPointZ, 0, len(mpz))
+	for _, v := range mpz {
+		if filter(v) {
+			n = append(n, v)
 		}
-		v, ok = mpz[idx], true
-		idx++
-		return
 	}
+	return n
 }
 
-// FilterRemoved returns an iterator that only delivers the vertices
-// which indices are not marked as removed.
-func (mpz MultiPointZ) FilterRemoved(removed map[int32]struct{}) func() (Vertex, bool) {
-	var idx int32
-	return func() (v Vertex, ok bool) {
-		for {
-			if idx >= int32(len(mpz)) {
-				ok = false
-				return
-			}
-			if _, rm := removed[idx]; rm {
-				idx++
-				continue
-			}
-			break
-		}
-		v, ok = mpz[idx], true
-		idx++
-		return
-	}
-}
-
-// MinMaxVertex runs over a point iterator and figures out its boundary.
-func MinMaxVertex(points func() (Vertex, bool)) (Vertex, Vertex) {
+// MinMaxVertex returns the extend of the point set.
+func (mpz MultiPointZ) MinMax() (Vertex, Vertex) {
 	min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64}
 	max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64}
-	for v, ok := points(); ok; v, ok = points() {
+	for _, v := range mpz {
 		min.Minimize(v)
 		max.Maximize(v)
 	}