1 files changed, 192 insertions, 0 deletions
diff --git a/climb/profile.go b/climb/profile.go
new file mode 100644
index 0000000..e88b2fa
--- /dev/null
+++ b/climb/profile.go
@@ -0,0 +1,192 @@
+package climb
+
+import "math"
+
+// ProfileSection is one step of the displayed elevation profile: the straight stretch between two
+// bends of the trace.
+type ProfileSection struct {
+	StartDistanceM  float64
+	EndDistanceM    float64
+	StartElevationM float64
+	EndElevationM   float64
+	SlopePercent    float64
+}
+
+// ProfilePolicy tunes how a trace is cut into display sections.
+type ProfilePolicy struct {
+	// MaxSections bounds the drawn steps: the widest bends are kept, the rest is averaged away.
+	MaxSections int
+	// MinBendM is the smallest elevation deviation worth its own step.
+	MinBendM float64
+	// BendPerBaseStep raises MinBendM with the sampling step, so a long climb tolerates
+	// proportionally more smoothing than a short one.
+	BendPerBaseStep float64
+}
+
+// DefaultProfilePolicy is the chart policy the web pages draw with.
+func DefaultProfilePolicy() ProfilePolicy {
+	return ProfilePolicy{
+		MaxSections:     24,
+		MinBendM:        1,
+		BendPerBaseStep: 0.01,
+	}
+}
+
+// profileBaseStepsM is the sampling ladder: the finest step that keeps a trace under
+// profileMaxBaseCells samples, so GPS noise averages out and a ramp still localizes.
+var profileBaseStepsM = []float64{25, 50, 100, 200, 500, 1000}
+
+// profileMaxBaseCells bounds the sampling grid, whatever the trace length.
+const profileMaxBaseCells = 120
+
+type profileNode struct {
+	distanceM  float64
+	elevationM float64
+}
+
+type profileSplit struct {
+	firstIndex int
+	lastIndex  int
+	splitIndex int
+	deviationM float64
+}
+
+// ProfileSections cuts the trace into sections whose width follows its bends: a constant stretch
+// collapses into one long section, a short steep ramp keeps its own, instead of both being averaged
+// into a fixed distance bucket.
+func (route Route) ProfileSections(policy ProfilePolicy) []ProfileSection {
+	if len(route.Distances) != len(route.Elevations) || len(route.Distances) < 2 {
+		return nil
+	}
+	lengthM := route.Distances[len(route.Distances)-1] - route.Distances[0]
+	if lengthM <= 0 {
+		return nil
+	}
+	baseStepM := profileBaseStepM(lengthM)
+	nodes := profileNodes(route, baseStepM)
+	bendM := math.Max(policy.MinBendM, policy.BendPerBaseStep*baseStepM)
+	return sectionsFromNodes(simplifyNodes(nodes, bendM, policy.MaxSections))
+}
+
+// profileBaseStepM picks the finest sampling step that keeps the grid under profileMaxBaseCells.
+func profileBaseStepM(lengthM float64) float64 {
+	for _, stepM := range profileBaseStepsM {
+		if math.Ceil(lengthM/stepM) <= profileMaxBaseCells {
+			return stepM
+		}
+	}
+	return math.Ceil(lengthM / profileMaxBaseCells)
+}
+
+// profileNodes resamples the trace on a uniform grid, interpolating between recorded points.
+func profileNodes(route Route, stepM float64) []profileNode {
+	startM := route.Distances[0]
+	endM := route.Distances[len(route.Distances)-1]
+	elevationAt := elevationInterpolator(route)
+	nodes := make([]profileNode, 0, int((endM-startM)/stepM)+2)
+	for distanceM := startM; distanceM < endM; distanceM += stepM {
+		nodes = append(nodes, profileNode{distanceM: distanceM, elevationM: elevationAt(distanceM)})
+	}
+	return append(nodes, profileNode{distanceM: endM, elevationM: elevationAt(endM)})
+}
+
+// elevationInterpolator returns the trace elevation at a distance. Callers walk it forward, so the
+// point cursor only moves on.
+func elevationInterpolator(route Route) func(distanceM float64) float64 {
+	index := 0
+	lastIndex := len(route.Distances) - 1
+	return func(distanceM float64) float64 {
+		for index < lastIndex-1 && route.Distances[index+1] < distanceM {
+			index++
+		}
+		nextIndex := min(index+1, lastIndex)
+		spanM := route.Distances[nextIndex] - route.Distances[index]
+		if spanM <= 0 {
+			return route.Elevations[index]
+		}
+		fraction := (distanceM - route.Distances[index]) / spanM
+		return route.Elevations[index] + min(max(fraction, 0), 1)*(route.Elevations[nextIndex]-route.Elevations[index])
+	}
+}
+
+// chordDeviationM is how far the node sits above the chord joining the two ends of its segment.
+func chordDeviationM(first, last, node profileNode) float64 {
+	spanM := last.distanceM - first.distanceM
+	if spanM <= 0 {
+		return node.elevationM - first.elevationM
+	}
+	chordElevationM := first.elevationM + ((node.distanceM-first.distanceM)/spanM)*(last.elevationM-first.elevationM)
+	return node.elevationM - chordElevationM
+}
+
+func bestSplit(nodes []profileNode, firstIndex, lastIndex int) profileSplit {
+	split := profileSplit{firstIndex: firstIndex, lastIndex: lastIndex, splitIndex: -1}
+	for index := firstIndex + 1; index < lastIndex; index++ {
+		deviationM := math.Abs(chordDeviationM(nodes[firstIndex], nodes[lastIndex], nodes[index]))
+		if deviationM > split.deviationM {
+			split.deviationM = deviationM
+			split.splitIndex = index
+		}
+	}
+	return split
+}
+
+// simplifyNodes keeps splitting the segment holding the widest bend until the trace fits the section
+// budget or has no bend left above bendM, so the cuts land on the most prominent bends first.
+func simplifyNodes(nodes []profileNode, bendM float64, maxSections int) []profileNode {
+	if len(nodes) <= 2 {
+		return nodes
+	}
+	lastIndex := len(nodes) - 1
+	kept := make([]bool, len(nodes))
+	kept[0] = true
+	kept[lastIndex] = true
+	pending := []profileSplit{bestSplit(nodes, 0, lastIndex)}
+	sectionCount := 1
+	for sectionCount < maxSections {
+		widest := -1
+		for index, split := range pending {
+			if split.splitIndex == -1 {
+				continue
+			}
+			if widest == -1 || split.deviationM > pending[widest].deviationM {
+				widest = index
+			}
+		}
+		if widest == -1 || pending[widest].deviationM <= bendM {
+			break
+		}
+		split := pending[widest]
+		pending = append(pending[:widest], pending[widest+1:]...)
+		kept[split.splitIndex] = true
+		sectionCount++
+		pending = append(pending, bestSplit(nodes, split.firstIndex, split.splitIndex), bestSplit(nodes, split.splitIndex, split.lastIndex))
+	}
+	simplified := make([]profileNode, 0, sectionCount+1)
+	for index, node := range nodes {
+		if kept[index] {
+			simplified = append(simplified, node)
+		}
+	}
+	return simplified
+}
+
+func sectionsFromNodes(nodes []profileNode) []ProfileSection {
+	sections := make([]ProfileSection, 0, len(nodes)-1)
+	for index := 1; index < len(nodes); index++ {
+		start := nodes[index-1]
+		end := nodes[index]
+		spanM := end.distanceM - start.distanceM
+		if spanM <= 0 {
+			continue
+		}
+		sections = append(sections, ProfileSection{
+			StartDistanceM:  start.distanceM,
+			EndDistanceM:    end.distanceM,
+			StartElevationM: start.elevationM,
+			EndElevationM:   end.elevationM,
+			SlopePercent:    (end.elevationM - start.elevationM) / spanM * 100,
+		})
+	}
+	return sections
+}