Bezier Splines

1.3. Defining Edge Splines

Various methods exist to determine the medial axis of an object, such as Voronoi skeletons, distance transform methods and thinning algorithms such as wave propagation via grassfire algorithms. In this study, a Bezier curve approximation technique was used where univariate spline curves (Dierckx 1981, 1993) were fitted to the North and South Canyon edges of VMs via the UnivariateSpline function from Python’s SciPy library (10).

A total of 93 coordinate points were identified on both North and South Canyon edges, with all points evenly spaced and with longitudinal distances preserved within the Cartesian kilometer-based CRS (ESRI:103885).

The edge splines of the main canyon were constructed with parametric cubic spline continuous functions, ensuring first and second derivative continuity and avoiding sharp changes in the curve’s slope or curvature, which are common in lower-degree splines.

Denoting the splines for the North and South edges as N(t) and S(t), where t is the parameter along the spline or arc length, the spline functions for each Canyon edge in a 2D plane can be written as:

(1)

(2)

where:

• B_i,₃ (t) and B_j,₃ (t) are the cubic B-spline basis functions and are defined over a knot vector that is locally supported.
• n and m are the numbers of control points for the North and South edge splines, respectively, determined through the spline fitting process.
• The interval [a, b] represents the domain of the parameter t.
• P_1i and P_2j are the control points for the North and South edge splines and define the initial shape and direction of the spline, which is then adjusted via the optimization process (see below). While control points do not generally lie on the curve, they act as weights in the combination of basis functions that define the curve.

1.4. Spline Optimization

Univariate spline interpolation involves a piecewise-defined cubic polynomial function that fits a set of data points in a single variable, with knots segmenting the dataset and individual polynomial functions fitted within each segment. The placement of knots and the overall smoothness of the curve in the UnivariateSpline algorithm are regulated by a key parameter: the smoothing factor (s).

When s is unspecified, UnivariateSpline defaults to place knots at every data point, effectively interpolating the data. However, specifying a value for s balances the spline’s smoothness with its adherence to the data by a least-squares optimization process, where the spline coefficients are calculated to minimize the sum of squared residuals between the fitted spline and the data points, subject to the smoothness constraint imposed by s.

Given a set of data points (x_i, y_i), where i = 1,2,…,N, and a spline function S(x), the goal is to find the spline that minimizes the following objective function:

(3)

where:

• The sum of the squared residuals, representing the fit of the spline S(x) to the data points (x_i, y_i), (lower values indicate a better fit), is added to the integral of the square of the second derivative of the spline, which acts as a measure of the spline’s smoothness. A smoother curve will have a lower value of this integral.
• λ represents the smoothing factor s in UnivariateSpline and is a parameter that balances the two objectives: smoothness of the spline and proximity to the data points.

This optimization problem is solved in UnivariateSpline via numerical methods, as curve fitting problems typically lack closed-form solutions. With the user-defined s value constraint, the algorithm iteratively adjusts the knots’ positions and the coefficients of the spline to minimize J, capturing the essential trends of the coordinate data while mitigating overfitting.

1.6. Medial Axis Extraction

The medial axis spline is a path equidistant from two edge splines at any point along its length. If the medial axis spline is denoted by M(t), a point M(t_i)=(x_m(t), y_m(t)) on this spline should satisfy the following condition where the distance from any point on the medial axis spline to the nearest point on each edge spline is equal:

‖ M (t) ― N ( t ) ‖ = ‖ M (t) ― S ( t ) ‖

(4)

Finding midpoints between two splines by simply averaging the distance on a vertical bisector so that Mid_y =(North_y + South_y)/2, in the context of VM, is not an optimal method for finding the medial points. The VM system trends ≈20° from northwest to southeast, curving in multiple ways; therefore, line segments joining the Canyon edge splines requires orthogonality to fully capture the curvature and direction of the splines.

The goal is to find an equal number of points along each of the Canyon edge splines where lines joining pairs of points are perpendicular (normal) to both splines at their respective points of tangency and then find the midpoints M(t_i) of these lines. This multidimensional optimization problem has a 2D parameter space that is required to optimize two variables (tnorth) and (tsouth), which are parameters along North and South splines at intersection points with a line segment between both splines (11). The optimization sequence can be written as:

• Smoothing: Apply s = 5000 to both North and South splines.
• Point Pair Identification: Generate a sequence of point pairs {N(t_i, S(t_i)} for i=1,2,…,n along the entire length of the North and South splines.
• Slope threshold: Calculate the slope at each point on the southern spline as the first derivative of the spline at each point, denoted as Slope_S(t_i). Define a threshold value for the slope as slope_threshold
• Conditional processing: For each point pair {N(t_i, S(t_i)}, check the slope at S(t_i):
  If: Slope_S(t_i)< slope_threshold, proceed with subsequent steps (5--10).
  Else: Skip to the next point pair {N(t_i, S(t_i)} and proceed with subsequent steps (5 to 10)
• Vector: Define L as the vector connecting a point on the North spline at (tnorth) to a corresponding point on the South spline at (tsouth) expressed as L= S(tsouth)−N(tnorth).
• Tangent Vectors: Calculate the tangent vectors at tnorth and tsouth as TN (tnorth) and TS (tsouth), respectively.
• The perpendicularity function f (tnorth, tsouth) is:

  f(tnorth, tsouth) = [L ⋅ TN(tnorth)]² + [ L ⋅ TS(tsouth)]²
• This function quantifies the squared sum of the dot products of the vector L with the vectors at TN (tnorth) and TS (tsouth).
• Multidimensional Optimization: adjust tnorth and tsouth to minimize f (tnorth, tsouth) to zero and therefore the perpendicularity of L with both TN (tnorth) and TS (tsouth).
• Midpoint Calculation: Use the perpendicular bisector to calculate the midpoint M(t_i) of L, which represents the optimized medial axis point for each t_i
• Iteration: Repeat the process for 1,000 points to compile the set {M(ti)}.

(This sequence forms the basis of the logic to produce Python code, which, in this study, was written to determine the Medial Axis with the SciPy-optimize function from Python’s SciPy library (11)).

Normal line equivalence at points (tnorth)_and(tsouth) requires the tangents of the curve at each point to be parallel. For the original dataset of the Canyon edge coordinate points, this condition is potentially difficult to meet, as there may be significant direction changes in the Canyon edges.

However, the smoothing factor used in the first step of the above procedure helps to regularize the data by reducing noise and minor fluctuations and facilitates the multidimensional optimization process by making the behavior of the spline (especially its tangents) more consistent and less prone to abrupt changes. Too much variability in the coordinate data might otherwise mislead an optimization algorithm to either prevent convergence or alternatively produce a local solution that is not globally optimal. For consistency in this study, a smoothing factor of 5,000 was used for both the North and South splines.

A specific conditionality function is also integrated into the above sequence (steps 3 & 4) to address the anomaly at the midpoint of the South Spline, where the region of Melas Chasma imposes pronounced southward semicircular extrusions along the canyon edge and spline representation (see Figure 3, above). The steeper tangent slopes in this region potentially lead to complexities in the algorithm’s ability to resolve accurately. To avoid this, the algorithm employs a threshold-based decision process where if the calculated slope at any given point exceeds a predetermined angle – set for this study at ±25 degrees – the algorithm omits processing for that particular point pair, advancing instead to the subsequent pair. Given the localized nature of the Melas Chasma steep tangent slopes, this adjustment impacts only a negligible portion of the overall spline representation.

The determination of each coordinate point M(t_i) in the saved array {M(ti)} can be described more formally by optimizing the objective function J₁, which enforces equidistance and orthogonality constraints for points on the medial axis spline:

(6)

where:

• d₁(t_i) and d₂(t_i) are the distances from M(t_i) to N(t_i) and S(t_i), respectively.
• O₁(t_i) and O₂(t_i) are measures of the deviation from orthogonality at M(t_i) to N(t_i) and S(t_i), respectively. This is computed as the sum of the squares of the dot product of the vector L with the tangent vectors T_N and T_S.
• λ₁ and λ₂ are weighting parameters for the importance of equidistance and orthogonality in the objective function.
• t_i represents the parameter values at which the spline and constraints are evaluated.
• N is the number of points along the spline where calculations are made.

The array of medial axis coordinate points {M(t_i)} was used to construct a medial axis spline in a 2D plane via the same method as the Canyon edge splines were created via the UnivariateSpline algorithm, formally expressed as:

(7)

where:

• M(t) represents the medial axis spline.
• B_k,₃(t) are the cubic B-spline basis functions for the medial axis spline.
• P_mk are the control points of the medial axis spline.
• p is the number of control points for the medial axis spline, which is determined through the spline fitting process.
• The interval [c, d] represents the domain of the parameter t for the medial axis spline.

For the medial axis spline optimization the objective function J_M balances the dual goals of fitting the medial axis spline closely to the determined medial axis points while also maintaining a smooth curve controlled by the choice of the smoothing parameter:

(8)

where:

• (x_i, y_i) are the coordinate points that the medial axis spline is intended to fit.
• q is the number of data points.
• The sum of the squared residuals, representing the fit of the medial axis spline to the data points (x_i, y_i), is added to ∫ [M’′(x)]²dx, the integral of the square of the second derivative of the spline, and acts as a measure of the spline’s smoothness.
• λ is the smoothing parameter that balances the trade-off between the fit to the data points and the smoothness of the spline (represented by the smoothing factor s in UnivariateSpline).

1.7. Medial Axis Spline Charts

Figure 4. Medial axis spline fitted to 1,000 median axis points. The smoothing factor of 5,000, resulted in 12 knots and an SSR of 5,000.

Figure 4. Medial axis spline fitted to 1,000 median axis points. The smoothing factor of 5,000, resulted in 12 knots and an SSR of 5,000.

Figure 5. Positions relative to a THEMIS daytime IR mosaic composite image of Valles Marineris of the fitted splines (North, South & Medial). The shaded area represents the extent of the VM Main Canyon.

Figure 5. Positions relative to a THEMIS daytime IR mosaic composite image of Valles Marineris of the fitted splines (North, South & Medial). The shaded area represents the extent of the VM Main Canyon.

Regression Analysis

2.2. Model: Cubic Polynomial

Cubic polynomial models are widely used in situations where a balance is needed between model complexity and the ability to accurately represent the underlying trends of the data without the risk of overfitting associated with higher degree models. The loss function for n data points can be defined as:

(10)

This loss function was minimized by a least squares algorithm implemented in Python by scipy.optimize.curve, and the fitted cubic polynomial was determined for 1000 coordinate points as:

y = 0.0000x³ + -0.0035x² + 59.0218x + -332751.1942

(11)

2.5. Model Comparison - 100 Coordinate Points

Figure 7. Chart showing 4 Models fitted to medial axis coordinates (100).

Figure 7. Chart showing 4 Models fitted to medial axis coordinates (100).

Table 4. Comparison of regression metrics.

Table 4. Comparison of regression metrics.

	Log-Likelihood	AIC	BIC	SD of Residuals	R²	Adjusted R²	RMSE
Model
Cubic Polynomial Model	-425.196	858.392	868.812	17.082	0.990	0.990	16.996
Sinusoidal Model	-429.416	866.832	877.253	17.818	0.989	0.989	17.729
Linear Regression Model	-519.145	1042.290	1047.500	43.707	0.935	0.935	43.488
Geodesic Path	-552.370	1108.740	1113.950	35.240	0.876	0.873	60.319

The plot shown in Figure 7 above, is not at a sufficient scale to distinguish separate sinusoidal and cubic polynomial fitted curves - there is considerable overlap. However, a clearer way of discriminating between each model can be seen in Table 4 above, where a dataset of 100 medial axis points was examined according to a variety of standard metrics. The cubic polynomial model’s ability to adapt to more complex patterns, although numerically showing the best performance, does not provide a significant advantage, suggesting that the complexity it adds may not be necessary to capture the primary dynamics of the data, thereby tending to validate the simpler sinusoidal model.

Cubic polynomial models have more coefficients for each degree of the polynomial, allowing a fit that captures more fluctuations in the data than a sinusoidal model, which is characterized by a single frequency, amplitude, and phase. The Akaike information criterion (AIC) and the Bayesian information criterion (BIC) are standard techniques for comparing models with different complexities, with the AIC penalizing complexity by adding twice the number of parameters (2k) and the BIC further multiplying the number of parameters by the logarithm of the data points (log(n) *k). However, in both cases, while complexity is penalized, a good fit is always rewarded, and if the polynomial model’s complexity leads to a better fit, this can suggest improved performance despite a greater number of parameters.

4. Discussion

5. Limitations

6. Conclusion

References

1. Klokočník J, Kletetschka G, Kostelecký J, Bezděk A, 2023. Gravity aspects for Mars, Icarus, Volume 406, 2023, 115729, ISSN 0019-1035. [CrossRef]
2. Hynek, B. M., Beach, M., Hoke M. R. 2010. Updated global map of Martian valley networks and implications for climate and hydrologic processes. Journal of Geophysical Research: Planets, 115(E9). [CrossRef]
3. Hansjoerg J. Seybold et al. 2018. Branching geometry of valley networks on Mars and Earth and its implications for early Martian climate. Sci. Adv.4, eaar6692. [CrossRef]
4. Brustel, C., Flahaut, J., Hauber, E., Fueten, F., Quantin, C., Stesky, R. Davies, G. R., 2017. Valles Marineris tectonic and volcanic history inferred from dikes in eastern Coprates Chasma, Journal of Geophysical Research: Planets, 10.1002/2016 JE005231, 122, 6, (1353-1371). [CrossRef]
5. Fueten, F., Novakovic, N., Stesky, R., Flahaut, J., Hauber, E., 2017. The Evolution of Juventae Chasma, Valles Marineris, Mars: Progressive Collapse and Sedimentation. Journal of Geophysical Research. Planets, 122 (11), pp.2223-2249. ⟨10.1002/2017JE005334⟩. ⟨hal-02130274⟩. [CrossRef]
6. Andrews-Hanna, J. C. 2012. The formation of Valles Marineris: 1. Tectonic architecture and the relative roles of extension and subsidence. Journal of Geophysical Research: Planets, 117(E3). [CrossRef]
7. Gurgurewicz, J., Mège, D., Schmidt, F. et al. Megashears and hydrothermalism at the Martian crustal dichotomy in Valles Marineris. Commun Earth Environ 3, 282 (2022). [CrossRef]
8. https://doi.org/10.1038/s43247-022-00612. [CrossRef]
9. Wang, S., Liu, C., Panozzo, D., Zorin, D., & Jacobson, A. (2023). Bézier Spline Simplification Using Locally Integrated Error Metrics. SIGGRAPH Asia.
10. Dierckx, P.,1981. An improved algorithm for curve fitting with spline functions, report tw54, Dept. Computer Science, K.U. Leuven,.
11. Dierckx, P., 1993. Curve and surface fitting with splines, Monographs on Numerical Analysis, Oxford University Press, 1993.
12. scipy.interpolate.UnivariateSpline. SciPy v1.11.4 Manual https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.
13. Optimization and root finding. SciPy v1.11.4 Manual https://docs.scipy.org/doc/scipy/reference/optimize.html.
14. Ghosh, A., Holt, W. E., & Bahadori, A. (2019). Role of Large-Scale Tectonic Forces in Intraplate Earthquakes of Central and Eastern North America. Geochemistry, Geophysics, Geosystems, 20(4), 2134-2156. [CrossRef]