Hybrid surface fitting framework using local parametric optimization and adaptive truncated hierarchical refinement

Abstract

Accurate reconstruction of freeform engineering surfaces from measured point clouds is essential in modern CAD, CAM, CAE and industrial reverse-engineering pipelines. Conventional B-spline fitting often struggles to balance geometric fidelity, parameterization quality and computational efficiency, particularly when input data contain both smooth regions and localized irregularities. Although truncated hierarchical B-splines enable localized refinement, classical hierarchical fitting frequently applies refinement even when errors are predominantly caused by distorted parameterization rather than insufficient geometric resolution. This work presents a hybrid fitting framework that integrates local parametric optimization with adaptive truncated hierarchical refinement to achieve high accuracy while maintaining a compact surface representation. A global least-squares B-spline approximation establishes the baseline geometry, and an error-driven analysis identifies regions requiring improvement. Each region is first corrected through local optimization of parametric values to reduce mapping distortion, while refinement is introduced only when optimization alone provides insufficient error reduction. Truncation confines refinement to the selected region, limiting the recomputation to the associated control points. The framework is evaluated on several engineering geometries and consistently reduces refinement depth and degrees of freedom while achieving accuracy comparable to or better than classical hierarchical fitting. The results indicate that coupling parametric optimization with localized hierarchical refinement provides a computationally efficient and robust strategy for engineering-oriented freeform surface reconstruction.

Keywords

B-Spline surface least squares fitting parametric values optimization THB-Spline refinement hybrid fitting

1 Introduction

1.1 Motivation

Accurate reconstruction of freeform surfaces from measured or simulated data is essential in many engineering tasks, including hull-form modeling, geometric inspection, reverse engineering and the generation of watertight CAD and CAE models. Such reconstructions are also increasingly important in modern digital engineering environments, where accurate geometric representations serve as the foundation for digital twins, simulation-driven design and data-integrated product lifecycle workflows. Modern point clouds commonly include smoothly varying regions together with isolated or densely distributed geometric features. Effective reconstruction requires high geometric fidelity, adaptive refinement capabilities, computational efficiency and a compact surface representation.

The reconstruction process typically begins with the acquisition of geometric data through 3D scanning or simulation. The resulting point cloud represents a set of sampled points describing the object's geometry. Depending on the scanning technology and acquisition resolution, the point cloud may contain dense sampling, noise and other artefacts that require filtering or segmentation.^29,43 Although a point cloud provides geometric coordinates, it does not constitute a usable numerical representation for advanced computational procedures and therefore must be converted into an appropriate parametric model to obtain a functional and accurate freeform surface representation.⁹ The resulting set of shape parameters should be compact yet sufficiently informative to accurately represent the original geometry. Transforming a point cloud into a parametric model typically requires solving complex and often nonlinear numerical problems.

The reconstruction process proceeds with the selection of an appropriate parametric model that can describe the input data and support subsequent shape optimization. Model selection is often guided by practical considerations, since different parametric formulations vary in flexibility and their ability to represent complex geometries. Inadequate model choice may lead to unnecessarily high dimensionality, increased computational cost and reduced numerical stability.⁹

For engineering applications, B-splines and NURBS are among the most commonly used parametric surface models.

For geometries requiring local refinement, adaptive spline models such as T-splines and truncated hierarchical B-splines (THB-Spline) have become common, since they allow modification of selected subregions without affecting the global surface.

Industrial geometries often contain large smooth regions alongside small but critical features such as fillets, bolt seats or curvature transitions. Global spline refinement results in a higher number of degrees of freedom and may negatively affect downstream simulation and design tasks. In many cases, fitting errors arise not from inadequate geometric resolution but from parameterization distortion introduced by the scanning process or by the surface topology. These considerations motivate methods that improve parameterization when possible and introduce geometric refinement only when reparameterization cannot achieve the desired accuracy.

This paper proposes a hybrid fitting method that integrates gradient-based parametric optimization with hierarchical spline refinement. The method relies on the interaction between global parametric optimization and local truncated hierarchical refinement. In the first stage, gradient-based optimization is applied to the initial B-spline surface to improve its global topology and establish an adaptive parametric distribution. This stage not only enhances the quality of the initial fit but also provides valuable information for determining which regions exhibit higher geometric deviation and therefore require further refinement. The second stage represents the core of the proposed method. When a subregion of the THB-Spline surface is selected for refinement, gradient-based parametric optimization is first applied locally within that region. This allows the fitting to adapt more effectively to the geometric features without immediately increasing the number of hierarchical levels. If the local optimization achieves sufficient improvement, additional THB refinement can be avoided, resulting in lower computational cost and fewer shape parameters. Conversely, if no significant improvement is achieved through parametric optimization, a new THB refinement level is introduced to increase local resolution. In this way, the hybrid procedure intelligently balances global and local fitting accuracy, achieving robustness, adaptability and computational efficiency while minimizing unnecessary refinement.

The paper is organized as follows. Section 2 describes theory behind B-Splines and THB-Splines, Section 3 describes procedure behind B-Spline surface construction via least squares approximation and reducing the problem to linear system. Section 4 and 5 narrates developed methodologies upon which proposed hybrid fitting method is developed and proposed fitting method, respectively. In Section 6, several numerical examples are given, with Section 7 concluding the paper.

1.2 Literature review

NURBS and B-splines have become the de facto standard parametric models for surface fitting and CAD representation. This work focuses on fitting methods that incorporate parameterization of data acquired through 3D scanning. The first step toward constructing a CAD model is to parameterize the input point cloud into a two-dimensional parametric domain, typically the unit square in variables u and v. Early studies applied uniform, chord-length and centripetal parameterization to one-dimensional fitting problems.^18,46 For three-dimensional data, more sophisticated parameterization methods are required. These approaches include mapping techniques that compute a parameterization by minimizing deformation energy,⁴¹ as well as other energy-based strategies.^19,20 More recent work includes neural network–based parameterization methods,^27,48 in which models are trained to reproduce established mapping techniques.⁵⁷ provided a comprehensive review of point-cloud parameterization methods. For single-patch NURBS fitting, extensive research has focused on optimizing knot placement, control points, along with the associated weights. The control points are typically computed using a least-squares approximation.

Many B-spline fitting methods focus on optimizing knot values. Early studies employed uniform, chord-length and centripetal strategies for knot placement,^13,14 while subsequent work introduced adaptive methods based on local data deviation,³⁰ dominant-point strategies⁴⁴ and knot insertion and removal algorithms.^31,54 Metaheuristic approaches, including genetic algorithms and particle swarm optimization, have been applied to optimize knot positions, parametric values and control points jointly.^22,24,38 Ćurković et al. proposed dynamic reparameterization procedures combining gradient-based optimization of the parametric grid with auxiliary functions¹⁰ and principal component analysis with Gaussian convolution.¹¹ Neural network models have recently been applied to predict parametric values and knot locations.^49,55

NURBS fitting extends B-spline fitting through the inclusion of weight parameters. Gradient-based methods addressing the global NURBS fitting problem are presented in,^7,37 while metaheuristic strategies have been employed to optimize knots, weights and control points simultaneously.^8,23,53 Neural network architectures incorporating B-spline basis functions have also been proposed for NURBS parameter identification.^16,52

In cases involving more complex topologies, single-patch B-spline and NURBS representations often fail to achieve satisfactory fitting accuracy. To enable more precise and localized control over shape parameters in specific regions, adaptive parametric models with local refinement capabilities have been developed. Prominent examples include T-splines,⁵⁰ Locally Refined (LR) B-splines¹⁵ and hierarchical B-splines.²¹ T-splines are constructed on T-meshes, which facilitate local refinement through the introduction of T-junctions. However, T-spline basis functions are not guaranteed to be linearly independent.⁶ To address this issue,³⁹ proposed a class of T-splines that maintain key mathematical properties such as non-negativity, linear independence and the partition of unity. Similarly, LR B-splines achieve linear independence of basis functions under specific conditions. Building on this,⁴ and² extended the formulation of LR B-splines to ensure partition of unity, local support and linear independence. While both T-splines and LR B-splines can be viewed as B-splines defined over locally refined knot vectors, hierarchical B-splines adopt a fundamentally different approach. They are based on a multilevel construction within the tensor-product space of basis functions, enabling a more robust and computationally efficient form of local refinement.²⁶ Nonetheless, standard hierarchical B-splines lack the partition of unity property, which limits their direct applicability in certain contexts. To overcome this limitation,²⁵ introduced truncated basis functions, giving rise to THB-splines, which preserve the partition of unity property. This development paved the way for the efficient use of hierarchical B-spline spaces in geometric modeling. The first surface fitting approach employing THB-splines for CAD model construction utilized least-squares approximation by²⁵ and.³² Subsequently, quasi-interpolation schemes were introduced to compute spline coefficients.^5,51 More recently,³ proposed a hybrid approach that combines least-squares fitting with quasi-interpolation: in the first stage, a least-squares solution identifies regions requiring refinement, and in the second stage, quasi-interpolation defines the coefficients of the locally refined spline functions.

Recent studies in computer-aided engineering further highlight the growing role of data-driven modeling and advanced computational techniques in engineering applications. For example,³⁶ investigated automated surface defect detection in industrial manufacturing, while¹⁷ introduced a parametric and feature-based CAD dataset supporting human–computer interaction in 3D shape learning. Computational modeling and simulation approaches have also been explored in broader engineering contexts, including regional seismic damage assessment of building curtain walls⁴⁵ and gravity-constrained SLAM methods for improving mapping accuracy in large-scale environments.⁵⁶ In robotics,⁴² examined the role of joint range-of-motion development in improving learning performance. These studies illustrate the increasing integration of data-driven and optimization-based techniques across modern computer-aided engineering workflows.

1.3 Relation to previous work and contribution

Our earlier research primarily investigated global parameter optimization techniques for single-patch B-spline and NURBS surfaces.¹⁰ Proposed nonlinear and evolutionary optimization frameworks aimed at improving the global parameterization of single B-spline surfaces, whereas¹¹ derived analytical knot-sensitivity formulations and compared parameter and knot optimization within the broader context of shape optimization. Further studies from our group explored parameter selection, knot insertion and sensitivity analysis, yet these investigations remained confined to the paradigm of a single surface patch applied to relatively simple or moderately complex geometries.

A central limitation of these earlier methods lies in their underlying assumption that the entire geometry can be adequately represented by a single tensor-product B-spline or NURBS patch. While effective for smooth, compact and moderately sized surfaces, such global representations become increasingly unsuitable as geometric complexity grows. Complex engineering surfaces often exhibit strong curvature variations, disconnected regions of local detail and spatial scales that render global parameterizations unstable or prone to excessive distortion.

The present study addresses this fundamental constraint by extending our previous framework toward locally adaptive, multilevel surface representations. Specifically, we employ THB-Splines to achieve localized refinement while maintaining essential mathematical properties such as partition of unity and basis linear independence. Moreover, local parametric optimization is introduced to enhance parameter quality at the subdomain level, thereby mitigating mapping distortions without increasing the number of global control variables. This integration of adaptive refinement and local optimization establishes a hybrid modeling framework capable of reconstructing large, geometrically intricate surfaces that cannot be represented effectively by a single global B-spline or NURBS patch.

The principal contributions of this work are summarized as follows:

A hybrid surface fitting framework that integrates local parametric optimization with adaptive truncated hierarchical refinement.

A patchwise parameter-correction mechanism that mitigates mapping distortions without introducing additional global degrees of freedom.

A truncation-consistent refinement strategy that ensures strictly localized hierarchical enrichment.

A formal complexity analysis demonstrating that the proposed hybrid scheme requires fewer refinements and fewer unknowns than classical hierarchical fitting.

A topology-sensitive evaluation that quantifies the interplay between parameter correction and refinement intensity.

Applications to three representative engineering geometries that demonstrate the method's robustness and practical relevance.

2 B-Spline and THB-spline parameterization

2.1 B-Spline parameterization

B-Spline surface is defined as a product of bidirectional net of control points $Q_{i_{o}, i}$ and univariate B-Spline basis functions tensor product⁴⁶

\begin{aligned} S (u, v) = \sum_{i_{0} = 0}^{n_{0} - 1} \sum_{i_{1} = 0}^{n_{1} - 1} N_{i_{0}, p_{0}} (u) N_{i_{1}, p_{1}} (v) Q_{i_{0}, i_{1}}, \end{aligned}

(1)

where

n_{0}

and

n_{1}

are number of control points,

N_{i_{0}, p_{0}} (u)

and

N_{i_{1}, p_{1}} (v)

are the B-Spline basis functions of polynomial degree

p_{0}, p_{1}

defined per u and v parametric axis, respectively. Generally speaking,

0

index will define parameters which are corresponding to the u parametric axis, and index

1

to the v parametric axis parameters. B-Spline basis functions are evaluated by Cox-de Boor recursive formula

\begin{aligned} \begin{aligned} N_{i, 0} (u) & = {\begin{array}{ll} 1, & if {\bar{u}}_{i_{0}} \leq u \leq {\bar{u}}_{i_{0} + 1}, \\ 0, & otherwise, \end{array} \\ N_{i_{0}, p} (u) & = \frac{u - {\bar{u}}_{i_{0}}}{{\bar{u}}_{i_{0} + p} - {\bar{u}}_{i_{0}}} N_{i_{0}, p - 1} (u) \\ + \frac{{\bar{u}}_{i_{0} + p + 1} - u}{{\bar{u}}_{i_{0} + p + 1} - {\bar{u}}_{i_{0} + 1}} N_{i_{0} + 1, p - 1} (u) . \end{aligned} \end{aligned}

(2)

Knots $\bar{u}$ define local influence of B-Spline basis functions and are usually set as

\begin{aligned} \bar{u} & = {\underset{p_{0} + 1}{\underset{⏟}{0, \dots, 0}}, \frac{1}{n_{0} - p_{0}}, \frac{2}{n_{0} - p_{0}}, \dots, \\ \frac{n_{0} - p_{0} - 1}{n_{0} - p_{0}}, \underset{p_{0} + 1}{\underset{⏟}{1, \dots, 1}}} . \end{aligned}

(3)

Basis functions $N_{i_{1}, p_{1}} (v)$ and knots $\bar{v}$ per v axis are evaluated analogously. B-Spline basis functions of degree p are piecewise polynomials of degree p. Their smoothness increases with the degree: at a knot of multiplicity m, the functions are $C^{p - m}$ continuous. As the degree increases, smoothness is preserved as a combination of smooth functions. However, the problematic points in terms of smoothness are always where polynomial pieces meet. Due to the appropriate combination, at the knot the function is smoothed so that the smoothness increases by one degree. As the degree increases, the support of each basis function covers a larger number of knot spans so additional knots with multiplicity $p_{0} + 1$ are inserted at the ends of the knot vector. In the interior, we have no knots of higher multiplicity. The only seemingly higher multiplicity occurs at the boundaries. If splines were considered outside the domain, they would indeed have reduced smoothness. However, within the domain itself, smoothness is preserved, and formally this is a matter of left-hand and right-hand derivatives. In the multivariate case, the foregoing equality applies analogously along each parametric axis, with the tensor-product structure of B-Splines ensuring that the extension to higher dimensions is direct. Matrix of control points $Q_{i_{0}, i_{1}}$ is given as

\begin{aligned} Q = [\begin{array}{ccc} Q_{0, 0} & \dots & Q_{0, n_{1} - 1} \\ ⋮ & ⋱ & ⋮ \\ Q_{n_{0} - 1, 0} & \dots & Q_{n_{0} - 1, n_{1} - 1} \end{array}] \in R^{3 (n_{0}) \times (n_{1})} . \end{aligned}

(4)

2.2 THB-Spline parameterization

Further refinement of a B-spline parametric model is achieved through knot refinement, yielding a new refined global knot vector.²¹ To restrict refinement to regions of interest,^34,35 introduced the hierarchical B-spline (HB-Spline) space, defined via tensor product B-spline basis functions with nested spline spaces

\begin{aligned} V_{0} \subset V_{1} \subset V_{2} \dots \end{aligned}

(5)

The HB-spline space supports adaptive local refinement while preserving linear independence, non-negativity and the nested structure of spline spaces; however, it does not satisfy the partition-of-unity property. This limitation was addressed by²⁵ through the introduction of the truncation operator, which removes the contributions of refined basis functions at coarser levels, thereby restoring the partition-of-unity property and yielding the THB-Spline basis. The truncation mechanism exploits the fact that any basis function $β \in B^{l}$ at level l can be expressed as a linear combination of basis functions at the finer level $l + 1$ ²⁵

\begin{aligned} τ = \sum_{i = 0}^{p + 1} c_{i}^{l + 1} β_{i}^{l + 1}, \end{aligned}

(6)

where

τ

is B-Spline basis function denoted in THB refined domain

Ω^{l + 1}

and

c_{i}^{l + 1}

are refinement coefficients.^1,40 Besides preserving the partition-of-unity property, truncation improves stability and restricts the fine-level basis functions to their intended regions, thereby enhancing locality. Under dyadic refinement, which is commonly used, a B-spline basis function is refined into

p + 2

finer basis functions. Truncation of

τ

with respect to

B^{l + 1}

and

Ω^{l + 1}

is defined as

\begin{aligned} t r u n c^{l + 1} τ = \sum_{β \in B^{l + 1}, s u p p β ⊊ Ω^{l + 1}} c_{β}^{l + 1} (τ) β \end{aligned}

(7)

With truncation operator defined THB-Splines are defined as

Definition 1. The THB-Spline basis space $T$ is recursively defined 1.

Initialization:

\begin{aligned} T^{0} = B^{0} \end{aligned}

Recursive construction of $T^{l + 1}$ :

\begin{aligned} T^{l + 1} = T_{A}^{l + 1} \cup T_{B}^{l + 1}, for l = 0, \dots, N - 2, \end{aligned}

where

\begin{aligned} T_{A}^{l + 1} & = {t r u n c^{l + 1} τ : τ \in T^{l} a n d s u p p τ ⊈ Ω^{l + 1}}, \end{aligned}

\begin{aligned} T_{B}^{l + 1} & = {β \in B^{l + 1} : s u p p β \subseteq Ω^{l + 1}} = H_{B}^{l + 1} . \end{aligned}

Result:

\begin{aligned} T = T^{N - 1} \end{aligned}

The THB basis preserves the key properties of the HB basis, including linear independence, non-negativity and the nested structure of spline spaces, while also satisfying the partition-of-unity property. These features make it well suited for fitting procedures and numerical analysis such as isogeometric analysis. A simple example of THB basis construction is presented on Figures 1 and 2, in accordance with Definition 1.

Figure 1.

Hierarchical mesh (a) and nested sequence of domains (b) with hierarchy $Ω^{l} \supseteq Ω^{l + 1} for\; l = 0, \dots, 2$ .

Figure 2.

Starting from initial B-Spline basis $B^{0}$ (a), first level refinement is achieved with truncated basis from previous level $T_{A}^{1}$ (b) and standard basis $T_{B}^{1}$ (c) on refined level $Ω^{l}$ . Following the recursive construction from Definition 1, second level THB basis space is defined with truncated basis from previous levels $T_{A}^{2}$ (c) and with standard basis function $T_{B}^{2}$ (d) supported on $Ω^{2}$ . With two level refinement, final THB basis $T$ is established (f).

Presented in Figure 1, starting from initial grid $Ω^{0} = (\bar{u} \times \bar{v}) = [0, 1] \times [0, 1]$ , parametric domain is refined with two level dyadic refinement $Ω^{1} = [0.4, 1] \times [0.2, 1]$ and $Ω^{2} = [0.6, 1] \times [0.5, 0.8]$ , respectively. Recursive definition of THB basis space for a mesh hierarchy presented in Figure 1 is displayed on Figure 2, with every level THB basis defined. Starting from initial B-Spline basis space $B^{0}$ , every new level refined basis are displayed (Figure 2(b)-2(c)) until construction of final THB basis space (Figure 2(f)).

3 Construction of a B-spline surface using the least squares formulation

This section explains the process used in this paper to obtain an initial B-spline surface from point-cloud geometry. From the input point-cloud geometry, a matrix-based data set is obtained through harmonic mapping and interpolation. The geometry in matrix form is associated with an initial, equidistant grid of parametric values, which is subsequently optimized. The baseline parametric values and their role in B-spline surface construction, which plays a key role in the proposed fitting procedure, are explained in detail in this section. Finally, after defining the matrix data set and the grid of parametric values, the B-spline surface is constructed using a least-squares linear system. Crucially, the initially constructed B-spline surface serves as the only globally computed solution within the proposed procedure.

3.1 Geometry projection into 2D parametric space

The approach for constructing a B-Spline surface, used throughout this paper, requires a structured matrix data set $P$ , consisting of $m_{0} \times m_{1}$ points, derived from the point cloud of the 3D-scanned geometry. The point cloud is projected onto the parametric domain $(u, v) \in [0, 1] \times [0, 1]$ by means of harmonic mapping,⁴⁷ which determines the projection by solving the Laplace equation $Δ f = 0$ subject to Dirichlet boundary conditions, reducing to two decoupled linear systems $Ku = 0$ and $Kv = 0$ , where the boundary values $u_{k}, v_{k}$ are prescribed along the boundary contour of the geometry. The projected parametric coordinates are subsequently used to interpolate the point cloud data over an equidistant $(m_{0} \times m_{1})$ grid, yielding the matrix data set

\begin{aligned} P = [\begin{array}{ccc} P_{0, 0} & \dots & P_{0, m_{1} - 1} \\ ⋮ & ⋱ & ⋮ \\ P_{m_{0} - 1, 0} & \dots & P_{m_{0} - 1, m_{1} - 1} \end{array}] \in R^{3 (m_{0}) \times (m_{1})} . \end{aligned}

(8)

which establishes the discrete correspondence between the physical and parametric domains and serves as the basis for the subsequent fitting procedure. The equidistant grid simultaneously defines the initial distribution of parametric values

U

and

V

. The input point cloud, its projection onto the parametric domain, and the resulting matrix data set

P

are shown in Figure 3.

Figure 3.

Triangulated point cloud data of boat hull (a), its harmonic projection (b) and interpolated matrix data set (c).

Figure 4.

Negative quadratic surface with Gaussian function bump.

Figure 5.

Example of topological parameters optimization on quadratic function with asymmetrical Gaussian bump (Figure 4). B-Spline least squares approximation with equidistant cparametric grid in the first row (a-e) with sum of squares error $S S^{0} = 88.818$ , and with optimal parametric grid in the second row (b-f), with $S S^{o p t} = 2.334$ . Red dots on surface plots (c, d) and on contour plot (e, f) represent control points obtained with least squares linear system from Eq. (11.).

The resulting matrix dataset serves as the reference geometry for subsequent surface fitting and error evaluation procedures. It is important to emphasize that the matrix representation significantly reduces algorithmic complexity and computational cost during B-spline surface construction, as it provides a structured and uniform data layout. However, this representation is not suitable for capturing large or highly intricate point cloud geometries, since the regular matrix discretization may fail to preserve the detailed spatial characteristics of the original surface.

3.2 Least squares formulation for fitting B-Spline surface

When fitting a B-spline surface to a geometry defined by a matrix dataset, each parametric point $(u, v) \in [0, 1] \times [0, 1]$ , associated with an interpolated value in the matrix dataset $P$ , corresponds to a surface value $S (u_{i}, v_{j})$ on the B-Spline surface. The objective of the fitting procedure is to determine the control coefficients such that $S (u_{i}, v_{j})$ approximates the reference geometry $P$ in the least-squares sense. The parametric values are initially defined through two uniform parametric vectors, $u \in [0, 1]^{m_{0}}$ and $v \in [0, 1]^{m_{1}}$ , whose components are distributed equidistantly as ${k / (m_{0, 1} - 1}_{k = 0}^{m_{0, 1} - 1}$ .

Fitting procedure of B-Spline surface relies on minimizing error defined by standard least squares formulation¹⁰

\begin{aligned} E (U, V, Q) = \frac{1}{2} \sum_{i = 0}^{m_{0} - 1} \sum_{j = 0}^{m_{1} - 1} ‖ S (u_{i}, v_{j}) - P_{i, j} ‖^{2} \end{aligned}

(9)

The B-spline surface parameters, namely the polynomial degrees $p_{0}$ and $p_{1}$ , together with the control point grid dimensions $n_{0} \times n_{1}$ , are defined prior to the fitting procedure. Based on these parameters, the corresponding knot vectors are evaluated according to Eq. (3). Once the parametric configuration and knot vectors are established, the least-squares fitting problem formulated in Eq. (9) can be expressed as the minimization of the surface fitting error function. This objective function quantifies the deviation between the interpolated data points in $P$ and the evaluated B-Spline surface $S (u_{i}, v_{j})$ across the defined parametric grid.

\begin{aligned} E (Q) = \frac{1}{2} \sum_{i = 0}^{m_{0} - 1} \sum_{j = 0}^{m_{1} - 1} ‖ S (u_{i}, v_{j}) - P_{i, j} ‖^{2} \end{aligned}

(10)

which yields the linear system

\begin{aligned} \frac{\partial E}{\partial Q} = A^{T} A Q B^{T} B - A^{T} P B = 0 \end{aligned}

(11)

where

A

and

B

are matrices of basis function values. The least-squares formulation yields the control point matrix

Q

, as determined by Eq. (11). Once

Q

is computed, the B-Spline surface is fully defined and can be directly evaluated over the prescribed parametric domain. This final step completes the reconstruction process, establishing a continuous and analytically smooth surface representation that approximates the reference geometry

P

within the least-squares tolerance.

4 Methodologies for defining hybrid fitting method

The proposed hybrid fitting framework integrates two complementary methodological approaches. The first approach builds upon the work of,^10,11 where reparameterization of the parametric values is achieved through a gradient-based optimization procedure. This optimization ensures improved parameter uniformity and reduced mapping distortion across the parametric domain. The second approach employs THB-spline refinement by²⁵ and,^33,32 enabling adaptive local refinement within regions exhibiting elevated geometric approximation error. By combining these two strategies, the hybrid method achieves both enhanced parametric quality and locally controlled surface fidelity.

4.1 Optimization of parametric values

When applying the fitting procedure defined in Eq. (10) each point from the geometry matrix $P$ is associated with a corresponding parametric pair $(u_{i_{0}, j_{0}}, v_{i_{1}, j_{1}})$ , which form parametric grid. As previously discussed, the matrices $U$ and $V$ are initialized with equidistantly distributed parameter values, providing a uniform sampling of the parametric domain. This uniform parametric structure is illustrated in Figure 5(a).

The influence of parametric values optimization can be interpreted through B-Spline surface theory.¹⁰ Each control point $Q_{i_{0}, i_{1}}$ affects a local subdomain of the surface defined by the knot span $[{\bar{u}}_{i_{0}}, {\bar{u}}_{i_{0} + p_{0} + 1}) \times [{\bar{v}}_{i_{1}}, {\bar{v}}_{i_{1} + p_{1} + 1})$ . Consequently, for each geometry point $P_{i, j}$ , there exists a corresponding surface point $S (u_{i}, v_{j}) = P_{i, j}$ . By adjusting the parametric values $u_{i} \in [{\bar{u}}_{i_{0}}, {\bar{u}}_{i_{0} + p_{0} + 1})$ and $v_{j} \in [{\bar{v}}_{i_{1}}, {\bar{v}}_{i_{1} + p_{1} + 1})$ , the influence regions of control points are modified, effectively adding or removing parametric samples from their respective knot spans. This manipulation alters the spatial weighting of basis functions and thus adjusts the geometric contribution of each control point within the parametric domain. The optimization is carried out using a gradient-based procedure, specifically the Levenberg–Marquardt algorithm, where the parametric vectors $u \in [0, 1]^{m_{0}}$ and $v \in [0, 1]^{m_{1}}$ , which define the matrices $U$ and $V$ , serve as the optimization variables in a nonlinear least-squares formulation. During each iteration, updated control points are computed by solving the local least-squares system defined in Eq. (11). The effect of this parameter optimization is demonstrated using a test surface consisting of a negatively curved quadratic base with an asymmetrically positioned anisotropic Gaussian bump

\begin{aligned} f (x, y) = - x^{2} - y^{2} + 2 \exp (- [\frac{{(x - 0.3)}^{2}}{{0.1}^{2}} + \frac{y^{2}}{{0.25}^{2}}]), \end{aligned}

as shown in Figure 4.

The presented test function, defined on a grid of size $m_{0} \times m_{1} = 100 \times 100$ , is approximated using a cubic B-Spline surface with a control point grid of $n_{0} \times n_{1} = 9 \times 9$ . The solution of the least-squares linear system described by Eq. (11), computed using the initial equidistant parameter grid (Figure 5(a)), is shown in Figure 5(c) and 5(e). Following the nonlinear gradient-based optimization of the global parametric vectors, the results in Figure 5(b) illustrate the redistribution of parametric values in the vicinity of the Gaussian bump, the most dominant geometric feature of the surface. This localized dispersion of parametric samples increases the number of active knot intervals in that region, consequently concentrating more control points near the feature and thereby enhancing its local fitting accuracy. Figure 5(f) depicts the relocation of control points from an initially uniform distribution (Figure 5(e) toward a denser configuration around the geometric feature. The optimization of parametric values thus introduces an additional degree of freedom, enabling adaptive relocation of control points in accordance with the underlying surface geometry, ultimately improving the fitting performance. Further methodological details, including a comprehensive analysis of potential challenges and method adaptations, are available in work by.^11,10

4.2 Refinement subdomains selection and THB-Spline parameterization

Following global optimization of the parametric values, the fitting procedure proceeds with THB-Spline refinement applied to subdomains of pronounced geometric error. The fitting error is quantified as the Euclidean distance $D (P, S) = ‖ P - S ‖_{2},$ evaluated at each node of the $(m_{0} \times m_{1})$ parametric grid, where $P$ is the reference matrix data set and $S (u_{i}, v_{j})$ is the B-Spline surface evaluated at the optimized parametric configuration. Grid points satisfying the refinement condition

\begin{aligned} D_{i, j} > C_{e r r} max (D), \end{aligned}

(12)

where

C_{e r r} \in [0, 1)

is a user-defined error threshold, are selected as candidate regions for hierarchical refinement. The influence of

C_{e r r}

on the refinement behavior is summarized in Table 17 in Appendix B. The parametric coordinates of the selected points are extracted into vectors

u^{e r r}

and

v^{e r r}

and grouped into clusters using a parametric distance threshold

Δ u_{r e f}

: two points are assigned to the same cluster if their Euclidean distance in the parametric domain is less than

Δ u_{r e f}

; otherwise, a new cluster is initiated. The parametric distance threshold

Δ u_{r e f}

is a user-defined parameter whose value is determined based on the geometric characteristics of the input surface, and is therefore subject to the user's judgment. For each cluster, the refinement subdomain

Ω_{r, k}

is defined by expanding the cluster center by a knot-span offset

Δ \bar{u}

Δ \bar{v}

sufficient to encompass all cluster members, with subdomain boundaries snapped to the nearest knot values of the current knot vectors to ensure parametric consistency. It should be noted that the clustering procedure is not universally robust and may require manual intervention in certain configurations. An illustrative example of the refinement subdomain for a negative quadratic surface with a Gaussian bump is shown in Figure 6. For each subdomain

Ω_{r, k}

, the corresponding THB-Spline control points at the refined level are computed by solving the local least-squares system defined in Eq. (11).^25,32 The local subsystem is extracted from the globally assembled least-squares formulation and solved independently for each subdomain. Recursive application of THB-Spline refinement is admissible provided the control point grid dimension at the subsequent level satisfies

\begin{aligned} (q_{0} \times q_{1})^{l + 1} \leq m_{0} \times m_{1} . \end{aligned}

(13)

Figure 6.

Example of refinement subdomain definition: extracted error points $u^{e r r}$ , $v^{e r r}$ (blue crosses), cluster centroid (green square), and subdomain boundaries (red dashed line).

If this condition is violated, the least-squares system becomes underdetermined and further refinement is infeasible, indicating that the recursive process is inherently constrained by the resolution of the input data grid.

5 Proposed hybrid fitting procedure

At the beginning of the hybrid fitting procedure, the input geometry is provided in matrix form. This matrix representation inherently contains a predefined parametric grid initialized with equidistant parametric values. The same matrix-based geometry is employed as the reference surface for evaluating fitting errors between the constructed and target geometries. For an arbitrary B-Spline polynomial degree $(p_{0}, p_{1})$ and control point grid size $(n_{0} \times n_{1})$ , the initial B-Spline surface is generated by solving the least-squares problem defined in Eq. (11). The first stage of the proposed hybrid fitting framework involves global optimization of the parametric values using a gradient-based method, as detailed in Section 4.1. The optimization is performed with the Levenberg–Marquardt algorithm implemented in the ALGLIB C++ library. Following the optimization, refinement subdomains are identified according to the procedure described in Section 4.2. Once the subdomains are established, the hybrid fitting framework proceeds with three distinct fitting strategies:

Local parametric values optimization with THB-Spline parameterization.

THB-Spline recursive parameterization with optimal parametric values.

THB-Spline recursive parameterization with initial parametric values.

The flowchart presented in Figure 7 summarizes the proposed hybrid fitting framework. After the refinement regions are defined, the three hybrid fitting strategies are executed independently and their respective results are evaluated to identify the configuration that yields the most accurate surface reconstruction. The first strategy applies gradient-based optimization locally within each refinement subdomain, combined with a single level of THB-Spline parameterization. If no significant improvement in fitting accuracy is achieved or if the optimization process fails to converge, the second strategy is employed. In this approach, recursive THB-Spline parameterization is performed using the geometry associated with the globally optimized parametric values. The third strategy applies THB-Spline parameterization over the previously defined refinement subdomain, but utilizes the initial equidistant parametric configuration instead of the optimized one. Upon completing all three fitting procedures, the results are compared and analyzed by the user to determine the optimal solution. The selection criteria depend on the specific application requirements and may prioritize minimal fitting error, reduced computational time, or a smaller number of active shape parameters.

Figure 7.

Flowchart of proposed hybrid fitting method.

The primary objective of the proposed hybrid fitting framework is to enhance THB-Spline parameterization accuracy through optimization of parametric values. The global optimization stage improves the initial B-Spline surface fit to a degree determined by the geometric complexity of the input geometry, while simultaneously identifying regions that require local refinement. Consequently, the subsequent THB-Spline parameterization achieves superior fitting accuracy with adaptive control over precision, computational cost, and the extent of local enrichment. The workflow and outcomes of each stage are demonstrated on the reference example of a negative quadratic surface with an embedded Gaussian bump. Since the proposed procedure exhibits topology-dependent behavior, the effectiveness of each fitting step varies according to the geometric and topological characteristics of the input surface. A comparative analysis is conducted across all fitting stages, evaluating the reduction in fitting error, computational time, and the number of active shape parameters, thereby identifying the fitting strategy that yields the optimal balance between geometric accuracy, computational efficiency, and model compactness.

5.1 Local parametric values optimization with THB-Spline parameterization

Each refinement subdomain is parameterized using a single-level THB-Spline model and new control points are computed through local least-squares approximation. The local parametric vectors $u_{r, k}$ and $v_{r, k}$ , corresponding to the rows and columns of the refinement subdomain, serve as optimization variables. The optimization procedure follows the same approach described in Section 4.1, utilizing the Levenberg–Marquardt algorithm implemented in the ALGLIB C++ library. During each optimization iteration, the control points of the single-level THB-Spline surface are recalculated using the local least-squares approximation scheme. Figure 8 illustrates the results of the local gradient-based optimization performed on the negative quadratic surface with the Gaussian bump. After applying a single level of THB-Spline parameterization, a local control point grid of size $(q_{0} \times q_{1})^{l = 1} = 5 \times 8$ is introduced, reducing the sum of squared errors to $S S = 1.98$ . Further optimization of the local parametric values yields an additional reduction in the fitting error, achieving a final value of $S S = 1.7335$ . If the resulting reduction in fitting error exceeds a predefined threshold, further refinement is not performed, thereby preserving the current number of control points.

Figure 8.

Results of local parametric values optimization on negative quadratic surface with Gaussian bump. Blue points represent solution of global parametric values optimization and red points are locally optimized parametric points (a), with constructed THB-Spline surface (b) where red dots are control points from initial B-Spline with optimized parametric values and green dots are control points in first level of THB parameterization.

5.2 THB-Spline recursive parameterization with optimal and initial parametric values

If the reparameterization does not yield a sufficient reduction in fitting error, truncated hierarchical refinement is applied within the subdomain $Ω_{r, k}$ in a dyadic manner. At this stage, new basis functions from the subsequent hierarchical level are activated, while the corresponding coarse-level basis functions are truncated to preserve locality and ensure partition of unity. The control points associated with each new THB-Spline level that are active within $Ω_{r, k}$ are recalculated using the local least-squares formulation described in Section 4.2, leaving the remainder of the surface unchanged. The subdomains are recursively refined up to a prescribed maximum number of hierarchical levels, which is limited by the resolution of the matrix data grid. Two refinement strategies are considered: one that employs the parametric values obtained from global gradient-based optimization, and another that utilizes the initial equidistant parametric configuration corresponding to the classical THB-Spline parameterization. Following the same reference example of the negative quadratic surface with the Gaussian bump, both parametric configurations (optimized and initial) are refined up to three THB-Spline levels, resulting in a final local control grid of $(q_{0} \times q_{1})^{l = 3} = 29 \times 32$ . The final sum of squared errors obtained for the initial parametric configuration is $S S = 1.58$ , whereas the optimized parametric configuration achieves a significantly lower error of $S S = 0.68$ . The final reconstructed THB-Spline representations of the negative quadratic surface with the Gaussian bump are illustrated in Figure 9.

Figure 9.

Three level refined THB-Spline example function with optimal parametric values (a) and initial parametric values (b) where green dots represent last refined level control points.

5.3 Engineering application and complexity analysis

The hybrid method is designed for practical engineering workflows where accuracy, compact representation and computational efficiency are essential. Its ability to adaptively balance global and local fitting procedures enables efficient reconstruction of complex geometries while minimizing the number of control points. By integrating parametric optimization with THB-Spline refinement, the method achieves high geometric fidelity with reduced computational cost, making it suitable for CAD model reconstruction, reverse engineering and surface design tasks that demand both precision and scalability. Its ability to:

perform refinement only where necessary,

correct parameterization distortions,

preserve compatibility with CAD and CAE systems,

maintain local support and stable refinement through truncation,

makes the proposed hybrid framework particularly suitable for practical engineering applications such as reverse engineering, quality inspection and simulation preprocessing. These characteristics ensure that the method not only delivers accurate and computationally efficient surface representations but also integrates seamlessly into existing industrial modeling and analysis workflows.

To formalize the selection among the three fitting strategies introduced in Section 5, a composite decision functional can be defined as

\begin{aligned} J = α \frac{R M S E}{R M S E_{0}} + (1 - α) \frac{D o F}{D o F_{max}}, \end{aligned}

(14)

where

α \in [0, 1]

controls the trade-off between approximation accuracy and model complexity,

R M S E_{0}

denotes the error of the initial B-Spline fit, and

D o F_{max}

represents the maximum number of active degrees of freedom admissible within the considered refinement setting. Minimizing J yields an automatic strategy-selection criterion that remains consistent with the user-driven comparison adopted in the present work, while preserving flexibility for application-specific priorities. Larger values of

α

favor accuracy, whereas smaller values promote model compactness.

The algorithm performs two types of operations: (i) local parametric optimization and (ii) local truncated hierarchical refinement. We examine best, average and worst-case behaviors.

Cost of Local Parametric Optimization:

Let $m_{r, k}$ denote the number of points in each subdomain $Ω_{r, k}$ and p the spline degree. Basis evaluation per each subdomain costs:

\begin{aligned} O (m_{r, k} (p + 1)^{2}) . \end{aligned}

The Levenberg-Marquardt update solves a system of size $m_{r, k}$ , giving:

\begin{aligned} O (m_{r, k}^{3}) . \end{aligned}

Since $m_{r, k}$ is relatively small, this step is significantly cheaper than refinement and does not alter the dimension of the spline space. 2.

Cost of Local Refinement:

Let $n_{r, k}$ denote the total number of basis functions overlapping in $Ω_{r, k}$ . The local least-squares update for refinement costs:

\begin{aligned} O (n_{r, k}^{3}) . \end{aligned}

Typically $n_{r, k} > m_{r, k}$ , making refinement considerably more expensive than optimization. 3.

Average-Case Behavior:

In many engineering surfaces, errors arise primarily from parameterization distortion. Thus:

\begin{aligned} cos t_{hybrid} = O (m_{r, k}^{3}), cos t_{standard\; THB} = O (n_{r, k}^{3}) . \end{aligned}

The numerical experiments indicate that parametric optimization can resolve a significant portion of high-error regions before additional refinement becomes necessary for geometries with localized dominant features. However, its effectiveness decreases for surfaces characterized by dense and uniformly distributed geometric variations, where hierarchical refinement remains the primary mechanism for accuracy improvement. This confirms that the contribution of parametric optimization is inherently topology-dependent. 4.

Worst-Case Behavior:

The worst case occurs when parameter optimization is attempted for every patch but refinement is eventually required for all:

\begin{aligned} cos t_{worst} = O (m_{r, k}^{3}) + O (n_{r, k}^{3}) . \end{aligned}

Because $m_{r, k} ≪ n_{r, k}$ , this overhead is negligible. Thus:

\begin{aligned} cos t_{worst} (hybrid) \approx cost (standard\; THB) . \end{aligned}

Best-Case Behavior:

If all high-error patches are caused by parametric distortion, no refinement occurs:

\begin{aligned} cos t_{best} (hybrid) = \sum_{r} O (m_{r, k}^{3}), \end{aligned}

while standard THB always incurs refinement cost

\sum_{r} O (n_{r, k}^{3})

5.4 Implementation details and reproducibility

The proposed framework is implemented using hierarchical B-Spline (THB-Spline) representations with fixed polynomial degree p. All experiments are initialized using an equidistant parameterization of the input data, which serves as the baseline configuration for both the classical THB-Spline method and the proposed hybrid approach.

The refinement process is guided by a user-defined error threshold $C_{e r r}$ , which controls the identification of high-error regions. In the experiments, $C_{e r r}$ is selected within the range $[0.3, 0.7]$ , with $C_{e r r} = 0.5$ used as the default configuration to balance approximation accuracy and model complexity, as further analyzed in Appendix B. Refinement is performed hierarchically using admissible subdivision rules, ensuring nested spline spaces and controlled growth of degrees of freedom.

The iterative fitting procedure is terminated when no significant reduction in the approximation error is observed or when the maximum admissible refinement level $L_{max}$ is reached. Local parametric optimization is applied selectively in regions where parameterization distortion is identified as a dominant source of error. These adjustments are performed through localized perturbations of parametric coordinates within bounded local neighborhoods, without altering the underlying hierarchical structure.

All comparative experiments are conducted using identical spline degrees, refinement criteria, and admissible refinement depths for both the classical and hybrid approaches, ensuring a fair and consistent evaluation.

Additional sensitivity analyses, including the influence of $C_{e r r}$ and refinement depth, are reported in Appendix B. The structured matrix datasets and implementation details used to generate the reported results are available from the authors upon reasonable request.

6 Results and discussion

All numerical experiments employed a data set of size $m_{0} \times m_{1} = 100 \times 100$ and a cubic B-Spline surface, with the control point grid selected individually for each test case based on the geometric characteristics of the surface. Global gradient-based optimization of the parametric values was performed using the ALGLIB C++ library with the following configuration: convergence criterion $ϵ = 10^{- 4}$ , differentiation step $Δ x = 10^{- 5}$ , and maximum step length $s t e p_{m a x} = 0.5$ . Each optimization problem comprised 196 free variables, as the boundary elements of the parametric vectors were fixed to $u_{0} = v_{0} = 0$ and $u_{m_{0} - 1} = v_{m_{1} - 1} = 1$ , with all remaining variables initialized as equidistant parameters and constrained within $[0.000001, 0.999999]$ .

Local parametric optimization within each refinement subdomain $Ω_{r, k} = ([u_{min}, u_{max}] \times [v_{min}, v_{max}])_{r, k}$ followed the same gradient-based procedure with a modified differentiation step $Δ x = 10^{- 3}$ . Initial values were taken from the preceding global optimization solution, with a small perturbation $δ = 10^{- 6}$ added to each variable to avoid premature convergence at previously attained local minima; the influence of this perturbation is examined in Table 18 in Appendix B.

An error threshold of $C_{e r r} = 0.5$ was applied uniformly across all examples, while the parametric distance threshold $Δ u_{r e f}$ was assigned individually for each test case based on the geometric characteristics of the respective surface.

Fitting accuracy was assessed using the root mean squared error

\begin{aligned} R M S E = \sqrt{\frac{1}{m_{0} m_{1}} \sum_{i = 0}^{m_{0} - 1} \sum_{j = 0}^{m_{1} - 1} {[S (u_{i}, v_{j}) - P_{i, j}]}^{2}}, \end{aligned}

(15)

and the improvement relative to the initial B-Spline fit

R M S E_{0}

is quantified by the relative change

\begin{aligned} R C = \frac{R M S E_{0} - R M S E_{i}}{R M S E_{0}} [\%], \end{aligned}

(16)

where higher values of

R C

indicate greater improvement. Spatial error distributions are represented as contour plots of the Euclidean distance field

D (P, S)

. Degrees of freedom (DoF) will be denoted as cumulative sum of all active control point in all active THB levels. Reported computational times are decomposed into two components: (i) construction and solution of the local least-squares system for THB-Spline fitting, and (ii) parametric optimization.

All procedures were implemented in Cython with C/C++ extensions and executed on a Linux-based workstation (Debian 12) equipped with a Quad-Core Intel ^® i7-1165G7 processor at 4.70 GHz and 16 GB of RAM. The numerical examples presented in the following sections were selected for their distinct topological characteristics, enabling a systematic demonstration of the flexibility, robustness, and accuracy of the proposed hybrid fitting framework.

To evaluate the performance of the proposed framework, a comparative analysis with classical THB-Spline fitting is presented in Table 1. The classical baseline corresponds to recursive THB-Spline parameterization with the initial equidistant parametric values, whereas the hybrid result corresponds to the best-performing configuration obtained after introducing parametric optimization. Since both approaches are evaluated on identical matrix datasets, using the same spline degree, refinement subdomains and admissible refinement depth, the comparison isolates the effect of parametric correction on fitting accuracy and model complexity.

Table 1.

Benchmark comparison between classical THB-Spline fitting and the proposed hybrid framework on identical datasets.

	Classical THB-Spline				Proposed hybrid framework
Geometry	RMSE	RC [%]	Level	DoF	RMSE	RC [%]	Level	DoF
Boat hull	2.640	46.172	3	1860	2.184	55.473	3	1860
Clutch housing	0.935	38.187	2	852	0.791	47.710	2	852
Square plate	1.044	51.041	2	1868	1.054	50.604	2	1868
Square plate with elevation	0.916	26.930	2	1740	0.919	26.660	2	1740
Three peaks function	0.0053	84.565	4	2089	0.0055	83.894	4	2089

The results demonstrate that the hybrid framework outperforms the classical THB-Spline baseline on geometries where parameterization distortion is a dominant source of error, such as the boat hull and clutch housing. In contrast, for geometries characterized by dense local variations, both methods converge to comparable solutions once the same refinement depth is reached. This behavior highlights the topology-dependent interaction between parameter correction and hierarchical refinement.

The key solver settings, stopping criteria, perturbation parameters, and refinement thresholds used in the experiments are reported in Section 5.4 and detailed in Appendix B. The structured matrix datasets and implementation details required to reproduce the reported results are available from the authors upon reasonable request.

6.1 Numerical examples

The first numerical example considers the reconstruction of a boat hull geometry, representing a moderately complex surface characterized by a dominant bow feature and several longitudinal ridges converging at the bow. The input point cloud comprises 4111 data points, and the initial control point grid was set to $(q_{0} \times q_{1})^{0} = 15 \times 12$ . The input point cloud, the interpolated matrix data set $P$ , and the initial B-Spline surface are shown in Figure 10. The initial fit yielded $R M S E_{0} = 4.905$ . Global parametric optimization results are reported in Table 2 and illustrated in Figure 11. The optimization produced a noticeable redistribution of parametric values in the vicinity of the bow and along the longitudinal ridges, consistent with the dominant influence these features exert on the overall surface shape. The resulting concentration of control points in these regions reduced the $R M S E$ by $27.024 \%$ .

Figure 10.

Boat hull geometry: input data and initial B-Spline fit.

Figure 11.

Global parametric optimization results for the boat hull geometry.

Table 2.

Global parametric optimization results for the boat hull geometry.

Time [min:sec]	Iterations	RMSE	RC [%]
01:20.5	1969	3.579	27.024

Based on the error distribution and with parametric distance threshold $Δ u_{r e f} = 0.3$ , the refinement subdomain was defined as $Ω_{r, 0} = [0.25, 0.833333] \times [0.555556, 1.0]$ , as depicted in Figure 12. The selected subdomain contained a $4 \times 4$ active control point grid from the initial solution; after applying a single level of THB-Spline parameterization, the local control grid expanded to $(q_{0} \times q_{1})_{r, 0}^{l = 1} = 11 \times 8$ . The results of THB-Spline parameterization and subsequent local parametric optimization are summarized in Table 3 and shown in Figure 13. Single-level THB-Spline parameterization achieved an $R M S E$ reduction of $45.69 \%$ , while the additional local parametric optimization yielded only a marginal further improvement of $0.1 \%$ .

Figure 12.

Definition of the THB-Spline refinement subdomain for the boat hull geometry.

Figure 13.

THB-Spline surface with initial B-Spline control points (red) and first-level refinement control points after local optimization (green).

Table 3.

THB-Spline parameterization and local parametric optimization results for the boat hull geometry.

	Subdomains	$(n_{0} \times n_{1})_{r, k}$	Time [min:sec]	Iterations	RMSE	RC [%]
THB parameterization	All	/	00:00.002	/	2.664	45.69
Local optimization	$Ω_{r, 0}$	$57 \times 45$	01:31.19	209	2.663	45.703

Recursive\enlargethispage{-18pt} THB-Spline parameterization was subsequently performed independently for both the initial and globally optimized parametric configurations. Results up to the maximum admissible refinement level, determined by the data grid resolution and subdomain size, are reported in Table 4 and illustrated in Figures 14 and 15. Three refinement levels were attainable within the given subdomain. The optimized parametric configuration consistently outperformed the initial one at each level, achieving $R C = 55.473 \%$ at level 3 compared to $46.172 \%$ for the initial parametric values. Notably, the combination of global parametric optimization with single-level THB-Spline parameterization yields a comparable $R M S E$ reduction to three-level recursive THB-Spline parameterization with initial parametric values, demonstrating the effectiveness of the proposed hybrid strategy.

Figure 14.

$RC$ metric per refinement level for recursive THB-Spline parameterization with initial $(u \times v)^{0}$ and optimized $(u \times v)^{opt}$ parametric grids.

Figure 15.

Recursive THB-Spline surfaces at maximum refinement level; red dots denote initial B-Spline control points, green dots denote final refinement level control points.

Table 4.

Recursive THB-Spline parameterization results for the boat hull geometry with initial and optimized parametric values.

	Domain	Level	Time [min:s]	RMSE	RC [%]	$(q_{0} \times q_{1})_{r, k}^{l + 1}$	DoF
Initial parametric values	$Ω_{r, 0}$	1	00:00.002	3.330	32.107	$11 \times 8$	252
		2	00:00.005	2.773	43.467	$25 \times 16$	564
		3	00:00.018	2.640	46.172	$53 \times 32$	1860
Optimized parametric values	$Ω_{r, 0}$	1	00:00.002	2.664	45.690	$11 \times 8$	252
		2	00:00.005	2.272	53.668	$25 \times 16$	564
		3	00:00.017	2.184	55.473	$53 \times$	1860

The spatial distribution of the Euclidean distance error at each stage of the fitting procedure is presented in Figure 16, providing a visual assessment of the progressive reduction in geometric deviation achieved by the proposed hybrid framework.

Figure 16.

Euclidean distance error distribution at each stage of the hybrid fitting procedure for the boat hull geometry.

The second numerical example considers the reconstruction of an asymmetric half of a clutch housing geometry. Clutch housings are typically cylindrical components with smooth surface transitions; the half containing three asymmetrically positioned bolt features was selected, as these constitute dominant geometric irregularities inducing a nonuniform surface distribution. The input point cloud comprises 34,249 data points, and the initial control point grid was set to $(q_{0} \times q_{1})^{0} = 16 \times 12$ . The input point cloud, the interpolated matrix data set $P$ , and the initial B-Spline surface are shown in Figure 17. It is noted that the interpolated matrix data set $P$ , shown in Figure 17(b), exhibits a noticeable geometric loss relative to the original point cloud as a result of the projection and interpolation process. Although this loss could be reduced through alternative projection techniques,¹² it is disregarded in the present analysis, which focuses on verifying the performance of the proposed hybrid fitting framework in the context of parametric value optimization and THB-Spline parameterization. The initial fit yielded $R M S E_{0} = 1.5133$ . Global parametric optimization results are reported in Table 5 and illustrated in Figure 18. The optimization produced a redistribution of parametric values along the housing edges, reducing the $R M S E$ by $25.22 \%$ . Despite this improvement, the bolt regions remain insufficiently approximated.

Figure 17.

Clutch housing geometry: input data and initial B-Spline fit.

Figure 18.

Global parametric optimization results for the clutch housing geometry.

Table 5.

Global parametric optimization results for the clutch housing geometry.

Time [min:sec]	Iterations	RMSE	RC [%]
00:24.685	1190	1.132	25.22

Based on the error distribution and $Δ u_{r e f} = 0.3$ , two refinement subdomains were defined as $Ω_{r, 0} = [0.538462, 1.0] \times [0.0, 0.444444]$ , $Ω_{r, 1} = [0.0, 0.384615] \times [0.0, 0.444444]$ , as depicted in Figure 19. The selected subdomains encompass the bisected housing edge and two of the three bolt features; the central bolt was intentionally excluded from refinement and retained as a visual reference for assessing fitting quality in the neighboring regions. The initial control point grids within the refinement subdomains were $6 \times 4$ for $Ω_{r, 0}$ and $5 \times 4$ for $Ω_{r, 1}$ . After applying a single level of THB-Spline parameterization, the local control grids expanded to $(q_{0} \times q_{1})_{r, 0}^{l = 1} = 12 \times 8$ and $(q_{0} \times q_{1})_{r, 1}^{l = 1} = 10 \times 8$ , respectively. The results of THB Spline parameterization and subsequent local parametric optimization are summarized in Table 6 and shown in Figure 20.

Figure 19.

Definition of the THB-Spline refinement subdomain for the clutch housing geometry.

Figure 20.

Local parametric optimization results within the THB-Spline refinement subdomain.

Table 6.

THB-Spline parameterization and local parametric optimization results for the clutch housing geometry.

	Subdomains	$(n_{0} \times n_{1})_{r, k}$	Time [min:sec]	Iterations	RMSE	RC [%]
THB parameterization	All	/	00:00.002	/	0.908	40.01
Local optimization	$Ω_{r, 0}$	$44 \times 46$	03:47.65	367	0.901	40.47
Local optimization	$Ω_{r, 1}$	$37 \times 46$	03:29.95	340	0.890	41.186

Single-level THB-Spline parameterization reduced the RMSE by 40.01%. Local optimization was performed sequentially, first for $Ω_{r, 0}$ and subsequently for $Ω_{r, 1}$ yielding an additional reduction of 1.185% and a total cumulative RC of 41.186% relative to the initial B-Spline surface.

Recursive THB-Spline parameterization was likewise performed sequentially, beginning with $Ω_{r, 0}$ refined to the maximum admissible level, followed by $Ω_{r, 1}$ . Results for both parametric configurations are reported in Table 7 and illustrated in Figures 21 and 22. A maximum of two refinement levels was attainable within each subdomain. The optimized parametric configuration achieved a final RC of 47.710% compared to 38.187% for the initial parametric values, demonstrating a consistent and significant advantage of the hybrid strategy across both subdomains.

Figure 21.

$R C$ metric per refinement level for recursive THB-Spline parameterization with initial $(u \times v)^{0}$ and optimized $(u \times v)^{o p t}$ parametric grids.

Figure 22.

Recursive THB-Spline surfaces at maximum refinement level; red dots denote initial B-Spline control points, green dots denote final refinement level control points.

Table 7.

Recursive THB-Spline parameterization results for the clutch housing geometry with initial and optimized parametric values.

	Domain	Level	Time [min:s]	RMSE	RC [%]	$(q_{0} \times q_{1})_{r, k}^{l + 1}$	DoF
Initial parametric values	$Ω_{r, 0}$	1	00:00.002	1.232	18.607	$12 \times 8$	244
	$Ω_{r, 0}$	2	00:00.005	1.145	24.328	$24 \times 16$	532
	$Ω_{r, 1}$	1	00:00.002	1.026	32.224	$10 \times 8$	612
	$Ω_{r, 1}$	2	00:00.005	0.935	38.187	$20 \times 16$	852
Optimized parametric values	$Ω_{r, 0}$	1	00:00.002	0.987	34.773	$12 \times 8$	244
	$Ω_{r, 0}$	2	00:00.005	0.914	39.595	$24 \times 16$	532
	$Ω_{r, 1}$	1	00:00.002	0.855	43.483	$10 \times 8$	612
	$Ω_{r, 1}$	2	00:00.005	0.791	47.710	$20 \times 16$	852

The spatial distribution of the Euclidean distance error at each stage of the fitting procedure is presented in Figure 23, demonstrating the progressive reduction in geometric deviation across all fitting stages, with particular improvement evident along the housing edges and around the bolt regions.

Figure 23.

Euclidean distance error distribution at each stage of the hybrid fitting procedure for the clutch housing geometry.

6.2 Discussion and future research

In the proposed fitting procedure, first thing that requires more attention is the mapping and interpolation procedure. The implemented mapping and interpolation procedure demonstrably introduces non-negligible information loss, as evidenced by systematic discrepancies observed across multiple test cases. Presented interpolation is done linearly on the uniform grid size, without any focus on the geometric feature, where the influential loss of data occurs. This can be improved by applying different interpolation procedure which should focus on geometric features.¹²

The recursive THB-Spline parameterization using the initial parametric configuration serves as a baseline representing purely hierarchical refinement without parametric optimization. Using the same refinement levels allows a direct comparison with the proposed hybrids trategy and isolates the influence of parametric optimization on the fitting performance. The classical THB-Spline parameterization using the initial parametric configuration serves as a baseline representing purely hierarchical refinement without parameter optimization. This enables a direct comparison and isolates the influence of parametric optimization on fitting performance. For further reference, Appendix A provides a benchmark approximation from standard THB works, reporting computational results in terms of error reduction and degrees of freedom, and distinguishing the advantages and limitations of the proposed hybrid method relative to existing approaches.

The sensitivity of the proposed method to geometric topology is consistently observed across multiple test cases, most notably during the global optimization of parametric values (Appendix A). For geometries exhibiting highly distinct geometric features, the global optimization stage fails to yield a significant reduction in approximation error (RC < 10%), rendering the framework predominantly reliant on classical THB-Spline refinement in such configurations. Two contributing factors to this limitation are identified. First, the Levenberg–Marquardt implementation employed in the proposed framework is sourced from an external library, whereby the gradient of the objective function is evaluated numerically rather than analytically; this contrasts with the original formulation¹² in which analytical derivatives with respect to the knot vector are used, yielding a more accurate and efficient gradient computation. Second, further tuning of the Levenberg–Marquardt hyperparameters, or the adoption of alternative optimization strategies better suited to complex geometric topologies, represents a viable direction for future improvement.

Within the proposed framework, selection of the optimal fit is delegated to the user, who determines the most suitable approximation result by balancing the trade-off between RMSE reduction and the number of active degrees of freedom (i.e., control points). The framework is formulated on the basis of sum-of-squares error minimization for computational simplicity; however, alternative objective functions incorporating weighted combinations of approximation error and model complexity can be readily constructed to accommodate more specific requirements. An additional constraint in the form of a maximum number of allowable refinement levels is also admissible, and is inherently bounded by the data grid resolution $m_{0} \times m_{1}$ as defined in Eq. (13).

The THB-Spline representation, truncation mechanism and local least-squares refinement procedure used in this work are inherited from earlier studies on hierarchical spline fitting.^3,25,32 The novelty of the present work lies in coupling THB refinement with local parametric optimization used as a decision step prior to refinement, together with a complexity-oriented analysis and a topology-aware evaluation of the interaction between parameter correction and local hierarchical enrichment.

Although the framework is designed to maximize automation, several critical thresholds require user-defined input, most notably during the identification and definition of refinement subdomains. The incorporation of machine learning techniques to eliminate this dependency on expert judgment represents a promising direction for future research, with the potential to further enhance the degree of automation within the fitting procedure.

7 Conclusion

This work demonstrates that effective hierarchical spline reconstruction requires distinguishing geometric under-resolution from parametric distortion, a distinction that classical THB-Spline refinement does not make. Because standard refinement responds uniformly to all fitting error, it often introduces unnecessary subdivision, increasing the number of control points and the hierarchical depth. The proposed hybrid framework addresses this by first analysing the origin of the error and then applying an appropriate corrective action. Local parametric optimization is used to reduce mapping distortion, while refinement is invoked only when geometric limitations prevent further improvement. Through truncation, refinement remains confined to selected regions, so only the associated control points are recomputed. In practice, this leads to substantially fewer refinement operations, a lower number of degrees of freedom and more coherent spline representations, without changing the asymptotic complexity of THB refinement.

The numerical examples confirm that the method performs robustly on geometries with pronounced curvature variation, localized features and complex topology. The interaction between parameter optimization and localized refinement yields surfaces with improved local fidelity and more uniform control-point distribution, which is valuable in CAD and CAE workflows where model stability and interpretability are important. The framework can operate in a largely automatic mode, suitable for batch or pipeline integration, but also supports user-controlled refinement depth and subdomain selection for expert-driven reverse engineering and design tasks.

The present study is focused on validating the hybrid fitting strategy within the scope of parametric optimization and truncated hierarchical refinement, using a specific optimization and refinement configuration. A more detailed analysis of worst-case computational complexity, as well as systematic benchmarking across a broader range of industrial cases, is left for future work. Further extensions may include anisotropic refinement criteria, more advanced decision metrics for switching between optimization and refinement, and generalizations to volumetric or networked spline representations. Overall, the results show that integrating parametric optimization with localized hierarchical refinement provides an efficient, tunable and engineering-oriented approach to freeform surface reconstruction, achieving higher accuracy with lower model complexity than classical hierarchical fitting.

Footnotes

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Croatian Science Foundation under the project [HRZZ-IP-2022-10-8856].

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

ORCID iDs

Domagoj Samardžić

Andrijana Ćurković

Ivo Marinić-Kragić

Milan Ćurković

Appendix

References

Atri

Shojaee

. Meshfree truncated hierarchical refinement for isogeometric analysis. Comput Mech 2018; 62: 1583–1597.

Bracco

Lyche

Manni

, et al. Generalized spline spaces over T-meshes: dimension formula and locally refined generalized B-splines. Appl Math Comput 2016; 272: 187–198.

Bracco

Giannelli

Großmann

, et al. THB-Spline approximations for turbine blade design with local B-Spline approximations, 2022, 63–82.

Bressan

. Some properties of LR-splines. Comput Aided Geom Des 2013; 30: 778–794.

Buffa

Garau

. Refinable spaces and local approximation estimates for hierarchical splines. IMA J Numer Anal 2017; 37: 1125–1149.

Buffa

Cho

Sangalli

. Linear independence of the T-spline blending functions associated with some particular T-meshes. Comput Methods Appl Mech Eng 2010; 199: 1437–1445.

Carlson

. NURBS Surface fitting with Gauss-Newton. Luleå, Sweden: Luleå University of Technology, 2009.

Costa

Montemurro

Pailhès

. A general hybrid optimization strategy for curve fitting in the non-uniform rational basis spline framework. J Optim Theory Appl 2018; 176: 225–251.

Ćurković

Vučina

. Adaptive representation of large 3D point clouds for shape optimization. Cham: Springer, 2017, 547–553.

10.

Ćurković

Vučina

Ćurković

. Enhanced 3D parameterization for integrated shape synthesis by fitting parameter values to point sets. Integr Computer-Aided Eng 2017; 24: 241–260.

11.

Ćurković

Vučina

. Novel re-parameterization for shape optimization and comparison with knot-based gradient fitting method. Comput Methods Appl Mech Eng 2018a; 336: 304–332.

12.

Ćurković

Marinić-Kragić

Vučina

. A novel projection of open geometry into rectangular domain for 3D shape parameterization. Integr Comput Aided Eng 2018b; 25: 1–14.

13.

Dierckx

. An algorithm for smoothing, differentiation and integration of experimental data using spline functions. J Comput Appl Math 1975; 1: 165–184.

14.

Dierckx

. An algorithm for least-squares fitting of cubic spline surfaces to functions on a rectilinear mesh over a rectangle. J Comput Appl Math 1977; 3: 113–129.

15.

Dokken

Lyche

Pettersen

. Polynomial splines over locally refined box-partitions. Comput Aided Geom Des 2013; 30: 331–356.

16.

Elmidany

Elkeran

Galal

, et al. NURBS Surface reconstruction using rational B-spline neural networks. J Control Eng Technol 2011; 1: 34–38.

17.

Fan

Liu

, et al. A parametric and feature-based CAD dataset to support human-computer interaction for advanced 3D shape learning. Integr Computer-Aided Eng 2025; 32: 75–96.

18.

Farin

. Curves and surfaces for CAGD: A practical guide. San Francisco, CA, USA: Morgan Kaufmann, 2001.

19.

Floater

. Mean value coordinates. Comput Aided Geom Des 2003; 20: 19–27.

20.

Floater

Reimers

. Meshless parameterization and surface reconstruction. Comput Aided Geom Des 2001; 18: 77–92.

21.

Forsey

Bartels

. Hierarchical B-spline refinement. ACM SIGGRAPH Comp Graph 1988; 22: 205–212.

22.

Gálvez

Iglesias

. Efficient particle swarm optimization approach for data fitting with free knot B-splines. Computer-Aided Design 2011; 43: 1683–1692.

23.

Gálvez

Iglesias

. Particle swarm optimization for non-uniform rational B-spline surface reconstruction from clouds of 3D data points. Inf Sci (Ny) 2012; 192: 174–192.

24.

Gálvez

Iglesias

. Firefly algorithm for explicit B-spline curve fitting to data points. Math Probl Eng 2013; 2013: 1–12.

25.

Giannelli

Jüttler

Speleers

. THB-splines: the truncated basis for hierarchical splines. Comput Aided Geom Des 2012; 29: 485–498.

26.

Giannelli

Jüttler

Speleers

. Strongly stable bases for adaptively refined multilevel spline spaces. Adv Comput Math 2014; 40: 459–490.

27.

Giannelli

Imperatore

Mantzaflaris

, et al. Learning meshless parameterization with graph convolutional neural networks. In: Lecture notes in networks and systems 817. Singapore: Springer, 2024, pp.375–387.

28.

Giust

Jüttler

Mantzaflaris

. Local (T)HB-spline projectors via restricted hierarchical spline fitting. Comput Aided Geom Des 2020; 80: 101865.

29.

Han

Jin

Wang

, et al. A review of algorithms for filtering the 3D point cloud. Signal Process, Image Commun 2017; 57: 103–112.

30.

Huang

Cohen

. Affine-invariant B-spline moments for curve matching. IEEE Trans Image Process 1996; 5: 1473–1480.

31.

Kang

Chen

, et al. Knot calculation for spline fitting via sparse optimization. Computer-Aided Design 2015; 58: 179–188.

32.

Kiss

Giannelli

Zore

, et al. Adaptive CAD model (re-)construction with THB-splines. Graph Models 2014; 76: 273–288.

33.

Kiss

Giannelli

Jüttler

. Algorithms and data structures for truncated hierarchical B–splines. In: Mathematical methods for curves and surfaces: 8th international conference, Oslo. Norway, 2014, pp. 304–323. Springer Berlin Heidelberg.

34.

Kraft

. Hierarchical B-Splines. Citeseer, https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=189bbbe98584c6f1efcdc54529cb25a2f843970f (1994, accessed 1 February 2024).

35.

Kraft

. Adaptive and linearly independent multilevel B-Splines. In: Le Méhauté

Rabut

Schumaker

(eds) Surface fitting and multiresolution methods 1st ed. Nashville: Vanderbilt University Press, 1997, pp.209–218.

36.

Krassnig

Haselmann

Kremnitzer

, et al. Efficient surface defect detection in industrial screen printing with minimized labeling effort. Integr Computer-Aided Eng 2024; 32: 1–21.

37.

Laurent-Gengoux

Mekhilef

. Optimization of a NURBS representation. Computer-Aided Design 1993; 25: 699–710.

38.

T-H

Kim

D-J

Min

K-C

, et al. B-spline surface fitting using genetic algorithm. J Soc Naval Architects Korea 2009; 46: 87–95.

39.

Zheng

Sederberg

, et al. On linear independence of T-spline blending functions. Comput Aided Geom Des 2012; 29: 63–76.

40.

Liao

Liu

. Splines and subdivision. In: Subdivision surface modeling technology. Singapore: Springer Singapore, 2017, pp.15–63.

41.

Maillot

Yahiaz

Verroustz

. Interactive texture mapping. In Proceedings of the 20th annual conference on computer graphics and interactive techniques, SIGGRAPH 1993. Association for Computing Machinery, 1993, pp. 27–34.

42.

Naya-Varela

Faiña

Romero

, et al. Understanding joint range of motion development in robotic learning. Integr Computer-Aided Eng 2025; 32: 346–365.

43.

Nguyen

. 3D point cloud segmentation: A survey. In IEEE conference on robotics, automation and mechatronics, RAM - proceedings. IEEE Computer Society, 2013, pp. 225–230.

44.

Park

Lee

. B-spline curve fitting based on adaptive curve refinement using dominant points. Computer-Aided Design 2007; 39: 439–451.

45.

Park

Choi

Kim

, et al. Rapid simulation method for assessing seismic damage to building curtain walls on a regional scale using designed wind load capacity. Integr Computer-Aided Eng 2025; 32: 143–157.

46.

Piegl

Tiller

. The NURBS book. Berlin, Heidelberg: Springer Berlin Heidelberg, 1995.

47.

Remacle

Geuzaine

Compère

, et al. High-quality surface remeshing using harmonic maps. Int J Numer Methods Eng 2010; 83: 403–425.

48.

Rios

Jüttler

. LSPIA, (stochastic) gradient descent, and parameter correction. J Comput Appl Math 2022; 406: 113921.

49.

Saillot

Michel

Zidna

. B-spline curve approximation with transformer neural networks. Math Comput Simul 2024; 223: 275–287.

50.

Sederberg

Cardon

Finnigan

, et al. T-spline simplification and local refinement. ACM Trans Graph 2004; 23: 276–283.

51.

Speleers

. Hierarchical spline spaces: quasi-interpolants and local approximation estimates. Adv Comput Math 2017; 43: 235–255.

52.

Tian

Kong

Deming

, et al. An improved method for NURBS surface based on particle swarm optimization BP neural network. IEEE Access 2020; 8: 184656–184663.

53.

Ülker

Ulker

. NURBS CURVE FITTING USING ARTIFICIAL IMMUNE SYSTEM. Int J Innov Comp Inform Control ICIC Int c 2012; 8: 2875–2887.

54.

Weiss

Andor

Renner

, et al. Advanced surface fitting techniques. Comput Aided Geom Des 2002; 19: 19–42.

55.

Wen

Luo

Kang

. The deep neural network solver for B-spline approximation. Computer-Aided Design 2024; 169: 103668.

56.

Xian

Liu

Zhang

, et al. Gravity-Constrained simultaneous localization and mapping for suppressing map warping in complex large-scale environments. Integr Computer-Aided Eng 2025; 32: 229–243.

57.

Zhu

Iglesias

You

, et al. A review of 3D point clouds parameterization methods. In Lecture notes in computer science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 13352 LNCS. Springer, Cham, 2022, pp. 690–703.