非创伤性股骨头坏死塌陷评估及保髋治疗进展与挑战

Abstract

非创伤性股骨头坏死（NONFH）是一种严重影响生活质量的致残性骨科疾病。人工全髋关节置换虽能有效治疗严重塌陷的终末期病变，但对于年轻、活跃患者而言，延缓或避免关节置换、保留自身关节功能的保髋治疗才是最终目标，这也是临床与基础研究的热点与难点。在此背景下，正确认识NONFH的临床特征与病理过程以及采用何种保髋治疗方法具有重要意义。本文立足于当前NONFH的研究现状，首先剖析了NONFH临床特征与病理过程的关联性，提出首次出现髋部疼痛预示着股骨头内已发生软骨下骨骨折或股骨头塌陷；基于“坏死-修复-塌陷-骨关节炎”的NONFH病理过程，发现了囊性变在NONFH的修复与塌陷中具有关键作用，靶向调控囊性变内的微环境来促进修复和避免塌陷是非手术保髋获得满意疗效的核心。其次总结了预测塌陷的影像学方法，基于股骨头前外侧负重区正常骨质保留量，提出股骨头前侧和外侧保留角概念，可以精确地评估塌陷风险以及塌陷是否会停止。最后分析了非手术保髋治疗方法、手术保髋治疗方法和新兴技术治疗方法目前的实际应用现状以及面临的挑战。提出NONFH的精准诊断（预测股骨头塌陷或塌陷是否会停止）对及时采取精确的保髋治疗至关重要，采取哪种保髋治疗取决于股骨头负重区是否会发生塌陷，而根据股骨头前侧和外侧保留角，既可以预测塌陷，也可以确定保髋治疗适应证。在适应证明确的前提下，运用现有治疗手段，联合多技术、多模式制定重建力学传导路径的个体化综合保髋治疗方案，是目前临床应用较成熟的方法。随着生物技术与再生医学的不断发展创新，以及通过高质量临床研究不断优化治疗路径，保髋治疗的成功率有望得到进一步提升。

Dupin cyclides are canal surfaces with circular lines of curvature and are widely used in geometric modelling applications such as Computer-Aided Design (CAD), surface fitting, and blending. In this research paper, we present a numerical approximation scheme for Dupin cyclide patches using rational bicubic Bézier surfaces. The boundary and interior control points are determined by enforcing C 1 continuity conditions, while the associated weights are optimized through curvature based variational fairness criterion. The complete Dupin cyclide is constructed by applying the suitable isometries to the constructed patches. The proposed approximation scheme interpolates the given cyclide data and ensures continuity across patch interfaces. The approximation quality is assessed using standard geometric error measures, including root mean square error and maximum relative error. Numerical experiments demonstrate that the maximum relative error attains the values of 0.6306 and 0.3982 for the two primary patches, confirming the computational efficiency of the proposed scheme.

Introduction

Cyclides are low-degree algebraic surfaces that are closed under offset, and their curvature lines are orthogonal circles [1]. Cyclides play a crucial role in the geometric modelling of free-form surfaces, architectural geometry, boundary representation in solid modelling, and surface blending (see references [1][2][3]). The cyclide surface serves as a fundamental component of both infrastructure development and industrial innovation. Through an analysis of cutting-edge designs and methods [4,5], recent developments contribute to the future of industry. Bo et al. [3] presented a numerical optimization technique to fit free-form shapes with a cyclide spline. Suitable values of corner points and orthogonal frames at these corner points were estimated by minimizing the distance between interpolating cyclide surfaces and the target surfaces. Hoxhaj et al. [6] characterized all Dupin cyclides passing through a section of bicircular quartic planar curve. All Dupin cyclides, admitting a given section, are classified as a surface of families that uniquely embeds into a Dupin cyclidic system. Garnier et al. [7] used Dupin cyclide for blending two quadrics of revolution. The authors used these blending techniques for geometric modelling of sea-horse and satellite antennas. Garnier et al. [8,9] developed the algorithm to subdivide Dupin cyclide (ring, horned, and spindle) using quadratic Bézier curves with mass points. Alcázar, Dahlb and Muntingh [10] characterized the symmetric canal surfaces using a rational spine curve and radius function. The symmetric blends between two canal surfaces were also discussed. Zube and Krasauskas [11] introduced Möbius invariant bilinear rational representation for parameterization of Dupin cyclide with quaternions weights. Many classical properties were represented by the quaternion formula. Garnier et al. [12] discussed the conversion of biquadratic rational Bézier surfaces into torus and double sphere (special cases of Dupin cyclide). Zhou and Straßer [13] derived a Non-Uniform Rational B-Splines (NURBS) representation of a cyclide. First, the cyclide was reparametrized to the domain [-1 , 1]. Then the blossom representation and control points at infinity were used to keep the representation compact and uniform.

The problem of constructing smooth, aesthetically pleasing and geometrically well behaved surfaces is closely related to broader research in surface approximation, surface fitting and fairness optimization. Huang et al. [14] aimed at the approximation of Weingarten surfaces using curvature based optimization technique. Zhu et al. [15] introduced a novel reconstruction of surfaces technique using partial differential equations and Coons patches. Tsuchie [16] proposed a curvature variation based fairness measure for curves for improving geometric evaluation and curve design by demonstrating automobile applications. Zafar and Hussain [17] discussed a novel method of construction of fair curves by controlling the curve length and curvature variation using rational cubic Said Ball curve. The shape parameters are obtained after optimizing the curvature variation and stretch energy as objective functional.

Despite the extensive theoretical study of Dupin cyclide and the development of various reparameterization, there are limited methods that translate cyclide into low degree rational Bézier patches with controllable shape parameters. The novel approximation technique presented in this work, is based on rational Bézier surfaces. It assures G 1 continuity, that is, calculates data values and derivatives at end points. Moreover, it has 8 free parameters which can improve the error of approximation as well as the geometric approximation of the surface. The geometric approximation strategy is adopted to compute the suitable value of the shape parameters, which has not been explicitly addressed in the existing literature. The shape design parameters are obtained by minimizing a fairness function that penalizes curvature variation to produce smooth surface. A suitable set of isometries is then applied to the approximated patches to generate a complete Dupin cyclide. The rational bicubic Bézier surface model is adopted in this work because it provides a low degree representation with sufficient degree of freedom through control points and weights. NURBS is the state of the art method for representation of surfaces. The NURBS and quaternion based approaches [11,13] do not address the patch-wise numerical approximation under interpolation and continuity constraints. In the approach [13], cyclides are re-parameterized by NURBS through blossom representation, which neither interpolates cyclide data nor provides information about the continuity of the resulting approximation. Furthermore, the numerical scheme presented in [13] cannot proceed without introducing additional knots. In contrast, the approximation scheme proposed in this paper not only interpolates the data but it guarantees continuity across the surface. In [11], the principal patches of Dupin cyclides were parameterized by a rational Bézier surface with quaternion weights. The proposed construction [11] was not affine invariant, but it was Möbius invariant. The representation produced an implicit equation, principal curvatures, and representation as canal surfaces. As compared to [11], the Dupin cyclide approximation scheme proposed in this research paper is affine invariant. The proposed scheme is direct and performs well for any specified number of data points. Moreover, introducing rational polynomial functions provides more degrees of freedom through the inclusion of the weights, thereby increasing its effectiveness in designing free-form surfaces.

The paper is organized as follows. Section 2 discusses preliminaries. In Section 3, the rational bicubic approximation of Dupin cyclide is elaborated. In Section 4, free parameters are optimized according to the selected fairness measure. Section 5 is dedicated to error analysis. Isometries are applied on the approximated patches to obtain designing a complete surface in Section 6. Concluding remarks are highlighted in Section 7.

Preliminaries

In this section, the terms to be used in the next sections are reviewed.

Dupin cyclides

Dupin cyclides are canal surfaces having focal conics as its directrix. These are the envelopes of the variable spheres on one of the pair of focal conics and touches the sphere whose center is on its other focal conic [18]. Cyclides are categorized corresponding to the nature of the focal conics. Cyclides consisting of an ellipse and a hyperbola as focal conics are named as elliptic cyclides (non-degenerate) and those having two parabolas as focal conics are termed as parabolic cyclides (degenerate form). In extreme cases, if the hyperbola degenerates to a line and ellipse to a circle then the cyclide is a torus. Parametric equations of Dupin cyclides (non-degenerate) are the following.

The constants a, b, c are defined by the focal conics. The constants a, b, c are related as c 2 = a 2b 2 because of the condition of focal conics i.e. one conic is passing through the foci of the other. m describes the position of the center of moving spheres to the radius of the sphere. Each category of cyclide is further subdivided into three types depending upon the parameters m, a and c. For ring cyclide c < m ≤ a, for horned cyclide 0 < m ≤ c and spindle cyclide is generated for m > a (Figure 1).

Quasi-Newton method

Quasi-Newton method [19] is a widely used algorithm for solving non-linear unconstrained optimization problems. It offers computational efficiency, particularly for the optimization problem involving large number of variables. In each iterative step, it updates the inverse of Hessian without solving it, making the quasi-Newton method less costly. For a problem of dimension n, the Quasi-Newton method converges in n-steps for the quadratic objective function. For the objective function of higher degree, the direction vector needs to be reinitialized in the direction of negative gradient generally after nth or (n + 1)th step. The process continues until the solution approaches the stopping criterion. As a gradient based technique, the Quasi-Newton method converges to a local minimum, not necessarily the global minimum.

Approximation of Dupin cyclide

This section introduces a numerical approximation framework for Dupin cyclide based on Rational Bicubic Bézier Surface (RBCBS) [12]. RBCBS is characterized by a set of control points, corresponding weights, and basic functions. The surface is defined as the following parametric expression:

Here (w i j ) 0≤i, j≤3 represent weights and (b i j ) 0≤i, j≤3 are the control points. P i (u) and P j (v) are cubic Bernstein polynomials. Rational bicubic form depends upon the control points and the associated weights. Thus, it has a total of 32 degrees of freedoms (control points and weights). A good approximation scheme requires an adequate number of such degrees of freedom. In the proposed scheme, we calculate the control points by applying continuity constraints while the weights are treated as optimization variables to further refine the surface. Cyclides are geometrically complex figures (continuously varying tangents and curvatures). Therefore, it is not possible to approximate whole cyclide by single RBCBS. Approximation is carried out by dividing the whole surface into small patches. An arbitrary patch of Dupin cyclide is presented here as

The rational bicubic Bézier representation of Dupin cyclide patch is obtained by imposing the following conditions.

(i) The four corner points of cyclide patch and rational bicubic Bézier surface are used to calculate the corresponding corner control points of rational bicubic Bézier surface.

(ii) Similarly, the remaining boundary control points are computed by equating the following end point tangent interpolation condition.

Here S θ and S φ represent the tangent vectors along the parameters θ and φ at the corner points of Dupin cyclide patch. P u (u, v) and P v (u, v) are the tangent vectors of Rational Bicubic Bézier Surface (RBCBS) along u and v directions.

(iii) Central control points of the Bézier patch are evaluated by using twist vectors.

Using conditions (3)-( 9), the following control points are calculated:

where f ′ i s are given in appendix.

An arbitrary patch of the Dupin cyclide

2π can be approximated using the calculated control points (10) of RBCBS. Substituting the x, y and z coordinates of the calculated control points (b i j ) 0≤i, j≤3 along with cubic Bernstein polynomials P i (u), P j (v) in equation ( 2), an RBCBS approximation of Dupin cyclide is obtained. A family of approximating surfaces is obtained by varying the weight parameters (w i j ) 0≤i, j≤3 . These free parameters provide additional degrees of freedom for controlling the shape of the approximated surface. How to choose the best RBCBS approximation for given Dupin cyclide patch is an open question. In this research paper, the weights (w i j ) 0≤i, j≤3 are computed using the surface fairness criterion in the following section. Two patches are approximated using this representation over the parametric domain

Curvature minimization

Smoothness and fairness are central criteria in geometric modelling. Westgaard and Nowacki [20] provided detailed discussion of various second and higher order fairness measures. These measures based on different geometric properties of a surface, such as curvature or variation of curvature of the surface. In this study, Dupin cyclide is approximated using the curvature related fairness measure that depends upon the quadratic forms of partial derivatives of given surface P(u, v). The resulting approximation is obtained by minimizing the fairness criterion χ(P) = ∫∫ Ω (P It can be seen from the above calculation that each boundary control point depends on two weight parameters, the weight associated with the boundary control point and the weight corresponding to the adjacent corner of the rectangular patch. Therefore, without loss of generality, the four weights associated with the corner control points are set to one. The optimization problem (11) will provide the appropriate values of free parameters w ′ i j s for fair patch. Here P 2 uu and P 2 vv are the second order partial derivatives of RBCBS with the control points (b i j ) 1≤i, j≤3 and the weights (w i j ) 1≤i, j≤3 , calculated in Section 3. The integrals of P 2 uu and P 2 vv are evaluated by using Trapezoidal rule.

Here ∥.∥ is the Euclidean norm in R 3 . The computed values of M ′ i s and N ′ i s are provided in Appendix. The optimization problem (11) is solved in MATLAB using Quasi-Newton method discussed in Section 2. All the above discussion leads to the following algorithm.

Theorem 1 The Rational Bicubic Bézier Surface (RBCBS), defined in (1), approximate the Dupin cyclide if its control points are computed as and weights are computed by minimizing the surface fairness integral χ(P), defined in (11).

Algorithm 1. Details of the process shown in Figure 2 and steps are given below.

Step 1. First to approximate the Dupin cyclide by RBCBS over the rectangular patch

Step 2. Put values of θ 1 , φ 1 , θ 2 and φ 2 in (b i j ) 0≤i, j≤3 in set of equations (10), to calculate control points.

Step 3. Minimize χ(P) to obtain the fairest approximation of Dupin cyclide calculating most plausible values of weights. Algorithm 1 is applied to approximate the Dupin cyclide over two rectangular patches (see Figures 34). Table 1 shows the optimized values of the weights for each patch. The overall computational cost of the proposed approximation scheme is low and suitable for practical geometric modelling applications: as the evaluation of boundary and central control points involve only closed form algebraic expressions and only 12 weights are optimized using Quasi Newton method. As a result, the optimization converges quickly and each Dupin cyclide patch is generated within a fraction of a second in MATLAB. The efficiency of the proposed scheme is demonstrated by error analysis given in section 5.

Error analysis 5.1 Relative error

The efficiency of the proposed numerical scheme can be analyzed by evaluating relative radius error for each patch. Error function is defined as 2 is the radius function for the approximated patch. The reparameterization of θ , φ is considered corresponding to u, v for each patch to calculate the error. For the interval [0, π 2 ] × [0, π], the following reparameterization is used θ = u π 2 and φ = vπ for u, v ∈ [0, 1]. The maximum relative error for the approximated patch Ω 1 , is calculated to be 0.6306. Similarly, for the second patch defined over the interval [ π 2 , π] × [0, π], θ = (u + 1) π 2 and φ = vπ is used. The maximum relative error for second patch Ω 2 , is found to be 0.2982.

Root mean square error

The Euclidean distance is defined as D(u, v) = ∥S(θ , φ) -P(u, v)∥ over m × n sample points for the given domain Ω.

The root mean square error can be calculated as

, where r = m × n represents the total number of sample points.

Table 2 highlights the root mean square error and absolute maximum error of the proposed scheme for the two patches by taking 10 × 10 sampling points and using the reparameterization defined above. The table reports the magnitude of error before and after optimization of weights (whose optimal values are given in Table 1). The random values of weights before optimization are taken to be ω 01 = -1, ω 10 = 1.5, ω 02 = 1, ω 20 = 1, ω 11 = 2, ω 13 = 1, ω 12 = -1, ω 21 = 1.5, ω 22 = 2, ω 23 = 1, ω 31 = 1.5, ω 32 = -1. These metrics provide a measure of how closely rational bicubic Bézier surface approximates the Dupin cyclide across its domain. Figure 5 illustrates the RMSE surface map for the two patches Ω 1 and Ω 2 . The complete Dupin cyclide is approximated by applying the symmetry properties of Dupin cycle discussed in Section 6.

Construction of Dupin cyclide using isometries

An isometry g of R 3 is a mapping g: R 3 → R 3 such that it leaves the surface S (or curve C) invariant [7]. Identity mapping is the trivial isometry. Non-trivial isometries include axial symmetries (rotation about axes), mirror symmetries (reflection about a plane), translations and their compositions [8]. A remarkable property of Dupin cyclide is that it is symmetric about its planes containing focal conics. Here, to obtain complete Dupin surface, some of the nontrivial isometries are applied to two patches which are optimized by using the procedure defined in previous sections. Consider the first approximated cyclide patch, Ω 1 , with optimized weights over the interval, 6(a)). As Ω 1 is symmetric about xy-plane so by applying the isometric transformation g 1 (x, y, -z) the mirror image is obtained as shown in Figure 6(b). Isometry g 2 (x, -y, z) generates third patch of Dupin cyclide (Figure 6(c)). Composition of these two isometries g 3 (x, -y, -z) or half-turn (rotation by angle π) about the line intersecting the planes of symmetry generates the fourth patch (Figure 6(d)). Similarly, repeat the above-mentioned procedure for another cyclide patch, Ω 2 , which is approximated by RBCBs in Section 4 for the interval

The set of isometries is again applied to obtain the remaining patches (symmetries of Ω 2 ) of Dupin cyclide as illustrated in Figure 7. Combination of two approximated patches Ω 1 , Ω 2 and their symmetries approximate the complete Dupin cyclide surface (Figure 8).

Conclusion

This paper presents a rational bicubic Bézier surface approach to approximate Dupin cyclide patches defined over the domains

For each patch, the boundary and interior control points are calculated using C 1 continuity conditions. Optimization of the Bézier weights through a curvature based variational criterion provides additional flexibility and smoothness to each patch. Individual cyclide patches are constructed over the prescribed domain and suitable isometries are applied to obtain complete Dupin surface. Quantitative error analysis, including RMSE and absolute error metrics, confirmed the accuracy and robustness of the proposed approximation scheme.

The method provides a practical computational framework for applications in Computer-Aided Design (CAD), blending, architecture geometry and solid modelling, where smooth, low degree rational representation of cyclides are required [1][2][3]. However, the approach may have reduced flexibility for degenerate cyclide surfaces and the effectiveness of the fairness optimization depends on the choice of a suitable optimization strategy and appropriate initial values. Future work may extend the proposed framework to rational biquartic surface to increase modelling capacity and integrating adaptive patch refinement strategies for complex surface regions.

Appendix

The f ′ i s (for equation (10)) are given below:

) ,

) ,

) ,

) ,

) ,

) , f 13 = (g 13 (θ 1 , φ 1 ) , g 14 (θ 1 , φ 1 ) , g 15 (θ 1 , φ 1 )) , f 14 = (g 13 (θ 2 , φ 1 ) , g 14 (θ 2 , φ 1 ) , g 15 (θ 2 , φ 1 )) , f 15 = (g 13 (θ 1 , φ 2 ) , g 14 (θ 1 , φ 2 ) , g 15 (θ 1 , φ 2 )) , f 16 = (g 13 (θ 2 , φ 2 ) , g 14 (θ 2 , φ 2 ) , g 15 (θ 2 , φ 2 )) , g 1 (θ 1 , φ 1 ) = m (c -a cos cos θ 1 cos φ 1 ) + b 2 cos cos θ 1 , g 2 (θ 1 , φ 1 ) = b sin sin θ 1 (am cos cos φ 1 ) ,

) ,

The computed values of M ′ i S and )

,

Volume 7 Issue 2|2026| 1967Contemporary Mathematics

Contemporary Mathematics

Volume 7 Issue 2|2026| 1971Contemporary Mathematics

Volume 7 Issue 2|2026| 1975Contemporary Mathematics

Volume 7 Issue 2|2026| 1985Contemporary Mathematics

Volume 7 Issue 2|2026| 1985Contemporary Mathematics