The interpolation overshoots the next point and reverses direction. ParametricPlot[{line[[1]][tt], line[[2]][tt]}, {tt, 0.2, 0.3}, Epilog -> {Point[points[[All, 1 ;; 2]]]}] ![Mathematica graphics](https://i.sstatic.net/rYjZZ.png) You can reduce the interpolation order to `1` or use a centripetal parametrization `parametrizeCurve` from [J.M.'s answer](http://mathematica.stackexchange.com/questions/10273/higher-order-periodic-interpolation-curve-fitting/10277#10277). parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] := FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; MatrixQ[pts, NumericQ] tvals = parametrizeCurve[points]; line = Interpolation[Transpose[{tvals, #}]] & /@ Transpose@points; s[t_?NumericQ] := NIntegrate[Norm[{line[[1]][tt], line[[2]][tt]}], Evaluate @ DeleteDuplicates @ Flatten[{tt, 0, Select[Flatten[p["Grid"]], 0 < # < t &], t}]]; ParametricPlot[Evaluate@{s[ss], line[[3]][ss]}, {ss, 0, 1}, PlotRange -> {{0, All}, {0, All}}, AspectRatio -> 1/2, Epilog -> {Point@Table[{s[ss], line[[3]][ss]}, {ss, 0, 1, 0.1}]}] ![Mathematica graphics](https://i.sstatic.net/zxKlv.png) Including the grid points in `NIntegrate`, Evaluate@DeleteDuplicates@Flatten[{tt, 0, Select[Flatten[p["Grid"]], 0 < # < t &], t}] speeds up the integration. For a really fast implementation use s = NDSolveValue[{ss'[t] == Norm[{line[[1]][t], line[[2]][t]}], ss[0] == 0}, ss, {t, 0, 1}] which constructs an `InterpolatingFunction` for the arc length.