A Robust Skeletonization Method for High-Density Fringe Patterns in Holographic Interferometry Based on Parametric Modeling and Strip Integration

1. Introduction

Holographic interferometry is a cornerstone optical technique for full-field, non-contact measurement of displacements and deformations in experimental mechanics [1]. Its high sensitivity, enabling resolution of displacements on the order of the light wavelength, makes it invaluable for studying complex deformation phenomena. The accuracy of the method is intrinsically connected with the precision with which the resulting interference fringe patterns can be analyzed [2]. In principle, higher measurement resolution is achieved by analyzing interferograms with a high density of fringes. However, in modern practice, the outcome becomes critically dependent on the chosen digital image processing algorithm, as the manual detection of hundreds of fringes is unreliable [3].

The primary obstacle to automation is the presence of speckle noise – an unavoidable granular interference pattern generated by the random scattering of coherent light from surface micro-roughness, which is intrinsic to any physical hologram [4]. While post-processing filters (e.g., median or Gaussian) can reduce speckle visibility in a digital image, they require a delicate approach. Intensive smoothing can introduce systematic biases that shift the fringe locations, a particularly detrimental effect for high-density patterns where fine details are essential. Consequently, the development of robust, noise-resistant skeletonization algorithms, that can dispense with significant pre-filtering remains an actual problem in optical metrology.

A variety of general-purpose image processing techniques have been applied to fringe identification. For patterns with moderate noise, edge detection algorithms (e.g., Marr –- Hildreth [5], Canny [6] or Shen – Castan [7]) can be employed. However, for high-density patterns corrupted by strong speckle noise, these methods tend to produce fragmented, disconnected edge maps that require complex and often unreliable post-processing to assemble into complete contours.

Active contour models (snakes) and their variants (e.g., Gradient Vector Flow snakes [8,9]) offer a more flexible method by formulating the fringe identification process as a solution to an energy minimization problem [10]. This method has been successfully adapted for fringe analysis [11,12,13,14]; however, it has an inherent limitation when processing highly noisy, dense patterns, which lies in its high sensitivity to the snake’s initial position.

To overcome these limitations, we propose a novel skeletonization procedure that strategically integrates physical prior knowledge with a noise-robust numerical core. Our contribution is threefold:

• Physics-informed parametric modeling. We constrain the search for fringes to a purpose-built, finite-dimensional functional subspace ${\mathcal{A}}_{\mathit{\theta }}$, which is defined by specific parametric curves (e.g., trigonometric polynomials or splines). This subspace is constructed in such a way as to endow its elements with the required properties, such as smoothness, closure, and specific shape, guaranteeing that the identified fringes are physically plausible, connected curves.
• A robust strip-integration functional with local smoothing. For fringe localization we consider the specific functional in the form of a Gaussian-weighted integral of intensity over a narrow strip surrounding the candidate curve. This formulation allows the use of efficient gradient-based optimization techniques. To compute this functional, a continuously differentiable intensity field $\tilde{I}(x,y)$ is first obtained by local bicubic interpolation, described in Subsection 3.1. The resulting functional ${\mathcal{J}}^s\left(\mathit{\theta }\right)$ is smooth with respect to the curve parameters $\mathit{\theta }$, ensuring stable convergence.
• A recursive quasi-optimization algorithm. Exploiting the geometric similarity of adjacent fringes, the identification process proceeds recursively outward (or inward). The optimal parameters for an identified fringe seed the initial guess for the next, via a simple scaling transformation (quasi-optimization), followed by a full local refinement of all parameters. This strategy dramatically improves computational efficiency and reliability.

The performance of the proposed method is carefully validated. First, using synthetic interferograms with precisely controlled geometry and additive-multiplicative speckle noise, we demonstrate its acceptable accuracy and robustness compared to baseline methods. The error is estimated using two metrics: the Euclidean norm ($L^2$-norm) of the difference and its maximal absolute value ($L^{\infty }$-norm). Second, we apply the algorithm to a real, challenging interferogram from a bending experiment on a square plate, containing over 100 fringes. The algorithm successfully extracts the complete skeleton, enabling an accurate reconstruction of the displacement field $w(x,y)$. This result is shown to be in good agreement with an independent analytical solution, obtained in the paper [15], confirming the practical utility and accuracy of the entire procedure.

The remainder of this paper is organized as follows. Section 2 details the proposed skeletonization procedure. Section 3 describes key implementation aspects. Section 4 presents the quantitative validation on synthetic data. Section 5 showcases the application to real high-density interferograms. Finally, Section 6 discusses the results and outlines future work.

2. The Proposed Skeletonization Method

2.1. Physical Origins of Imperfections in Interferograms

To motivate the design choices of the proposed algorithm, it is essential to understand the physical nature of the imperfections present in a raw interferogram. These imperfections, which corrupt the ideal sinusoidal fringe pattern, arise from the optical setup imperfections and the nature of coherent light scattering. They can be categorized by their spatial scale, necessitating different corrective strategies in the digital processing procedure.

In a typical off-axis holographic setup (e.g., the Leith – Upatnieks scheme [2]), the recorded intensity results from the interference of a reference plane wave ${\mathit{E}}_r\left(\mathit{r}\right)$ and an object wave ${\mathit{E}}_o\left(\mathit{r}\right)$ reflected from the specimen surface. In complex notation, these waves are:

$$\begin{array}{cc}\hfill {\mathit{E}}_r\left(\mathit{r}\right)& =A_r{\mathit{e}}_r{e}^{i{\mathit{k}}_r·\mathit{r}},\hfill \\ \hfill {\mathit{E}}_o\left(\mathit{r}\right)& =A_o(x,y)\mathit{J}(x,y){\mathit{e}}_o{e}^{i{\mathit{k}}_o·\mathit{r}},\hfill \end{array}$$

where $A_r$ and $A_o(x,y)$ are the scalar amplitudes, ${\mathit{e}}_r$ and ${\mathit{e}}_o$ are the unit polarization vectors of the reference and initial object waves (directly from the laser source), ${\mathit{k}}_r$ and ${\mathit{k}}_o$ are the corresponding wave vectors, and $\mathit{J}(x,y)$ is the spatially varying Jones matrix accounting for polarization changes upon reflection from the deformed surface.

The interference of the wavefronts from the specimen’s reference and deformed states yields an intensity pattern of the form [4,16]:

$$I_{\mathrm{obs}}(x,y)\propto I_0(x,y)\left[1+V(x,y)cos\left(\Delta \varphi (x,y)\right)\right]R(x,y)+“\mathrm{noise}”.$$

(1)

Here, $\Delta \varphi (x,y)=(4\pi /\lambda )w(x,y)$ is the phase difference directly proportional to the deformable object out-of-plane displacement $w(x,y)$, $V(x,y)\propto |\left({\mathit{J}}_1{\mathit{e}}_o\right)·{\left({\mathit{J}}_2{\mathit{e}}_o\right)}^{*}|$ stands for the polarization-induced visibility term and $R(x,y)$ is the contrast variation, induced by coherent noise. A transition between consecutive dark (or bright) fringes corresponds to a displacement increment of $\lambda /2$.

The terms $I_0(x,y)$, $V(x,y)$ and $R(x,y)$ describe large- to medium-scale corruptions:

• Non-uniform background intensity $I_0(x,y)$: Caused by uneven illumination, varying surface reflectivity, and polarization effects (${\mathit{J}}_1\ne {\mathit{J}}_2$). This leads to slow intensity variations across the image.
• Visibility term and contrast variation $V(x,y),R(x,y)$: The coherence noisy factor $R(x,y)$ (dependent on laser temporal/spectral properties and surface roughness) modulates the overall fringe contrast. The visibility term $V(x,y)$ induced by polarization and speckle decorrelation can further reduce contrast, even to zero.
• Geometric distortion: The non-coaxiality of reference and object beams in off-axis schemes, combined with lens imperfections, introduces a projective distortion between the object plane and the image sensor.

These low-frequency imperfections are corrected in the preprocessing stage (Section 3) via geometric unwarping and local intensity equalization.

The speckle-dependent terms in (1) present the most significant challenge. It is a fine-grained, high-frequency random pattern resulting from two physically unavoidable phenomena: (i) the finite spectral bandwidth of the laser source, and (ii) the random interference of waves scattered by surface micro-roughness. Unlike Gaussian additive noise, speckle is signal-dependent and exhibits characteristic correlation lengths (${\rho }_{\mathrm{speckle}}$). Its spatial frequency content often overlaps with that of the high-density fringes themselves, making linear filtering unsuitable as it erodes fringe edges and introduces systematic localization errors.

This multi-scale analysis dictates our algorithmic strategy: large-scale distortions are removed by explicit preprocessing, while the fine-scale speckle noise is addressed not by filtering, but by the main idea of the proposed method – the parametric modeling (Subsection 2.2) and the strip-integration functional (Subsection 2.4) – which provide robustness without distorting the signal.

2.2. Parametric Modeling and Problem Formulation

The goal of skeletonization is to extract a family of planar curves $\Gamma ={\left\{{\gamma }_i\right\}}_{i=1}^N$ that represent the level sets (fringes) of the phase map $\Delta \varphi (x,y)$ from a discretely sampled intensity image I. A fundamental premise of our approach is to treat fringes as continuous geometric entities (curves) rather than collections of disconnected pixels. Consequently, the search for an optimal curve $\gamma$ must be performed within a suitable functional space.

Let $\Omega \subset {\mathbb{R}}^2$ denote the rectangular image domain. We consider the space of admissible curves $\mathcal{V}$ to be a set of rectifiable, closed Jordan curves contained in $\Omega$. For any such curve $\gamma \in \mathcal{V}$, we define a functional as the normalized integral of the image intensity over its support:

$$\mathcal{J}\left(\gamma \right)={\displaystyle \frac{1}{\mu \left(\gamma \right)}}\underset{\gamma }{\int }Id\mu ,$$

where $\mu$ is an appropriate measure on $\gamma$. In the simplest case, where the support is the one-dimensional curve itself, $\mu$ is the arc-length measure ($d\mu =ds$), and the normalization factor $\mu \left(\gamma \right)$ is the total curve length $L\left(\gamma \right)$. This yields the line-integral functional:

$${\mathcal{J}}^l\left(\gamma \right)={\displaystyle \frac{1}{L\left(\gamma \right)}}\underset{\gamma }{\oint }Ids.$$

The choice of measure and its domain will be generalized in Subsection 2.4 to define more robust functionals. The core optimization problem is to find curves within a constrained set that extremize this functional. For an ideal, normalized interferogram where the low-frequency imperfections from (1) are absent ($I_0\equiv \mathrm{const}.$, $V\equiv 1,R\equiv 1$) and speckle noise is removed, the intensity is proportional to $1+cos(\Delta \varphi (x,y\left)\right)$. The bright fringes (intensity maxima) correspond to curves that maximize ${\mathcal{J}}^l$, while dark fringes (minima) correspond to its minimizers. In a real, noisy image, the functional ${\mathcal{J}}^l$ must be extremized despite the corruptions described in Subsection 2.1.

However, operating directly in the infinite-dimensional space $\mathcal{V}$ is intractable and ignores the crucial physical prior: the expected fringe family possesses a specific global structure. Instead of regularizing the problem within $\mathcal{V}$ (as active contours do), we constrain the solution a priori to a finite-dimensional subspace ${\mathcal{A}}_{\mathit{\theta }}\subset \mathcal{V}$. This subspace is defined by a parametric model:

$${\mathcal{A}}_{\mathit{\theta }}=\{\gamma \left(\mathit{\theta }\right)\mid \mathit{\theta }\in \Theta \subset {\mathbb{R}}^m\},$$

(2)

where $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_m)$ is a vector of parameters, and the mapping $\mathit{\theta }\mapsto \gamma \left(\mathit{\theta }\right)$ is continuously differentiable. The choice of the model ${\mathcal{A}}_{\mathit{\theta }}$ takes into account our a priori knowledge about the deformation. For instance, for the bending of a centrally loaded, clamped plate, the fringes are expected to be nested, closed, and quasi-concentric curves, which can be effectively modeled by a family of perturbed circles or ellipses (see Subsection 3.4).

By restricting the search to ${\mathcal{A}}_{\mathit{\theta }}$, the infinite-dimensional variational problem of extremizing ${\mathcal{J}}^l$ over $\mathcal{V}$ reduces to a finite-dimensional optimization problem over the parameter space $\Theta$:

$${\mathit{\theta }}^{*}=\underset{\mathit{\theta }\in \Theta }{arg\; extr}{\mathcal{J}}^l\left(\gamma \left(\mathit{\theta }\right)\right).$$

This is the core conceptual shift from other methods: the smoothness, connectivity, and global shape of the solution are built into the search space ${\mathcal{A}}_{\mathit{\theta }}$, rather than being enforced via penalty terms during optimization.

2.3. Discretization and Coordinate System

The image data I is provided as a discrete matrix $\mathit{I}[r,c]$ of size $r_{\max }\times c_{\max }$, where r and c are integer row and column indices. To connect the continuous formulation with the discrete data, we establish a fixed coordinate system. We assume a unit pixel pitch and place the origin of the Cartesian coordinates $(x,y)$ at the center of a chosen reference pixel with indices $(r_0,c_0)$. The correspondence between a pixel $P\left[\right[r,c\left]\right]$ and the coordinates of its center $(x_P,y_P)$ is then given by:

$$(x,y\in P\left[[r,c]\right])\iff \left(r=r_0-\lfloor y+1/2\rfloor ,c=c_0+\lfloor x+1/2\rfloor \right).$$

(3)

This mapping is illustrated in Figure 1.

2.4. Variants of the Functional

The core of the skeletonization method is the optimization of a functional that measures how well a candidate curve aligns with a fringe. We trace its development from a naive discrete form to the final robust, continuous formulation.

2.4.1. Point-Sampling Along the Curve

The most direct approach is to rasterize the curve $\gamma \left(\mathit{\theta }\right)$ and sum the intensities of its pixels:

$${\mathcal{J}}_{\mathrm{point}}^l\left(\mathit{\theta }\right)={\displaystyle \frac{1}{N\left(\mathit{\theta }\right)}}\sum _{\mathcal{R}\left(\gamma \right(\mathit{\theta }\left)\right)}\mathbf{I}[r,c],$$

where $\mathcal{R}$ is the set of pixels that rasterize the curve, and $N\left(\mathit{\theta }\right)=\left|\mathcal{R}\right|$ is their number (an integer measure of curve length). This functional is discontinuous with respect to $\mathit{\theta }$, it means that a small parameter change can alter the set $\mathcal{R}$, causing a jump in the sum. This discontinuity precludes the use of gradient-based optimization.

2.4.2. Line-Integration with Interpolation

To obtain a smooth functional, we first construct a $C^1$-continuous intensity field $\tilde{I}(x,y)$ via bicubic interpolation (detailed in Subsection 3.1). Replacing the discrete sum with a normalized line integral yields:

$${\mathcal{J}}^l\left(\mathit{\theta }\right)={\displaystyle \frac{1}{L\left(\mathit{\theta }\right)}}\underset{\gamma \left(\mathit{\theta }\right)}{\oint }\tilde{I}(x,y)ds,$$

(4)

where $L\left(\mathit{\theta }\right)$ is the curve length. ${\mathcal{J}}^l\left(\mathit{\theta }\right)$ is now continuously differentiable, enabling gradient-based optimization. However, it remains sensitive to speckle noise, as it samples intensity only along an infinitesimally thin path.

2.4.3. Point-Sampling Over a Strip

To improve robustness, we consider a strip $\mathcal{S}(\gamma ,\delta )$ of width $2\delta$ around the curve. A discrete, robust functional can be defined by summing pixel intensities within this strip, weighted by their approximate Gaussian distance to the curve:

$${\mathcal{J}}_{\mathrm{point}}^s\left(\mathit{\theta }\right)={\displaystyle \frac{1}{N\left(\mathit{\theta }\right)}}\sum _{{\mathcal{R}}_{\mathcal{S}}}G\left({\displaystyle \frac{d_{rc}}{\sigma }}\right)\mathbf{I}[r,c],$$

where ${\mathcal{R}}_{\mathcal{S}}$ is the set of pixels whose centers lie within the strip $\mathcal{S}$, $d_{rc}$ is the minimal distance from the pixel center to the curve, and $G\left(u\right)=exp(-u^2/2)$. This choice of weight function is motivated by the statistical nature of speckle noise nature. The functional ${\mathcal{J}}_{\mathrm{point}}^s$ averages over many pixels, reducing the impact of speckles. However, it shares the discontinuity flaw of ${\mathcal{J}}_{\mathrm{point}}^l$ due to the discrete pixel set ${\mathcal{R}}_{\mathcal{S}}$.

2.4.4. Strip-Integration with Interpolation

The final, optimal form combines the robustness of strip averaging with the smoothness of continuous integration. We define the strip-integration functional as the normalized, weighted area integral over the continuous strip:

$${\mathcal{J}}^s\left(\mathit{\theta }\right)={\displaystyle \frac{{\displaystyle \underset{\mathcal{S}\left(\gamma \right(\mathit{\theta }),\delta )}{\int \int }G\left(\frac{d(\mathit{x},\gamma )}{\sigma }\right)\tilde{I}(x,y)dA}}{{\displaystyle \underset{\mathcal{S}\left(\gamma \right(\mathit{\theta }),\delta )}{\int \int }dA}}},$$

(5)

where, $d(\mathit{x},\gamma )$ is the minimal Euclidean distance between $\mathit{x}$ and the curve. The functional, ${\mathcal{J}}^s$, is both robust to noise (due to area averaging) and continuously differentiable in $\mathit{\theta }$ (due to the smoothness of $\tilde{I}$, G, and the parametric curve). It is the functional used in our final algorithm.

3. Numerical Implementation

The evaluation of ${\mathcal{J}}^l$ and ${\mathcal{J}}^s$ requires efficient numerical techniques for interpolation and integration.

3.1. Bicubic Interpolation Scheme

To obtain the $C^1$ intensity field $\tilde{I}(x,y)$ from the pixel matrix $\mathit{I}[r,c]$, we employ a compact bicubic formula. For a point $(x,y)$ in a grid cell with corners at pixel centers, the interpolation is:

$$\tilde{I}(x,y)=\left[\begin{array}{cccc}1& \Delta \left(x\right)& \Delta {\left(x\right)}^2& \Delta {\left(x\right)}^3\end{array}\right]{\mathit{M}}^{\mathsf{T}}\mathit{D}(i\left(y\right),j\left(x\right))\mathit{M}\left[\begin{array}{c}1\\ \Delta \left(y\right)\\ \Delta {\left(y\right)}^2\\ \Delta {\left(y\right)}^3\end{array}\right],$$

where $\Delta \left(x\right)=x+{\displaystyle \frac{1}{2}}-\lfloor x+{\displaystyle \frac{1}{2}}\rfloor$, $\Delta \left(y\right)=y+{\displaystyle \frac{1}{2}}-\lfloor y+{\displaystyle \frac{1}{2}}\rfloor$ are local coordinates within a cell of the interpolation grid. The integer-valued indices $i\left(y\right)$ and $j\left(x\right)$, which locate the pixel, are functions of the coordinates:

$$i\left(y\right)=r_0-⌊y+{\displaystyle \frac{1}{2}}⌋,j\left(x\right)=c_0+⌊x+{\displaystyle \frac{1}{2}}⌋.$$

(6)

The constant matrix $\mathit{M}$ and the matrix $\mathit{D}(i,j)$ (now with $i,j$ understood as the outputs of (6)) are:

$$\mathit{M}=\left(\begin{array}{cccc}1& 0& -3& 2\\ 0& 0& 3& -2\\ 0& 1& -2& 1\\ 0& 0& -1& 1\end{array}\right),\mathit{D}(i,j)=\left(\begin{array}{cccc}{f}_{ij}^{00}& f_{ij}^{01}& h_{ij}^{00}& h_{ij}^{01}\\ f_{ij}^{10}& f_{ij}^{11}& h_{ij}^{10}& h_{ij}^{11}\\ g_{ij}^{00}& g_{ij}^{01}& q_{ij}^{00}& q_{ij}^{01}\\ g_{ij}^{10}& g_{ij}^{11}& q_{ij}^{10}& q_{ij}^{11}\end{array}\right).$$

The sixteen elements of $\mathit{D}$ are finite-difference approximations of the intensity and its derivatives at the four cell corners, computed from the $4\times 4$ neighborhood of the pixel $(i,j)$:

$$\begin{array}{cc}\hfill f_{ij}^{pk}& =\mathbf{I}[i-k,j+p],\hfill \\ \hfill g_{ij}^{pk}& ={\displaystyle \frac{1}{2}}\left(\mathbf{I}[i-k,j+p+1]-\mathbf{I}[i-k,j+p-1]\right),\hfill \\ \hfill h_{ij}^{pk}& ={\displaystyle \frac{1}{2}}\left(\mathbf{I}[i-k-1,j+p]-\mathbf{I}[i-k+1,j+p]\right),\hfill \\ \hfill q_{ij}^{pk}& ={\displaystyle \frac{1}{4}}(\mathbf{I}[i-k-1,j+p+1]-\mathbf{I}[i-k+1,j+p+1]\hfill \\ \hfill & -\mathbf{I}[i-k-1,j+p-1]+\mathbf{I}[i-k+1,j+p-1]),\hfill \end{array}$$

for $p,k\in \{0,1\}$. This scheme guarantees that $\tilde{I}(x,y)$ is a $C^1$ function across cell boundaries, providing the smoothness required for stable gradient computation. The example of interpolation is shown in Figure 2.

3.2. Numerical Quadrature

The line integral (4) is approximated by discretizing the curve into n segments. For a parametric curve $\left(x\right(\phi ),y(\phi \left)\right)$, $\phi \in [0,2\pi ]$, we use the trapezoidal rule:

$${\mathcal{J}}^l\left(\mathit{\theta }\right)\approx {\displaystyle \frac{{\displaystyle \sum _{k=1}^n}\left(\tilde{I}(x_k,y_k)+\tilde{I}(x_{k-1},y_{k-1})\right)\sqrt{{(x_k-x_{k-1})}^2+{(y_k-y_{k-1})}^2}}{2{\displaystyle \sum _{k=1}^n}\sqrt{{(x_k-x_{k-1})}^2+{(y_k-y_{k-1})}^2}}},$$

where $(x_k,y_k)=(x\left({\phi }_k\right)+1/2,y\left({\phi }_k\right)-1/2)$. The discretization step $\Delta \phi$ is chosen adaptively based on the curve’s local curvature to ensure at least one sample per intersected pixel:

$$\Delta ={({\dot{x}}^2+{\dot{y}}^2)}^{-1/2}.$$

In case of strip integration the denominator of (5) can be computed more convenient analytically:

$$\underset{\mathcal{S}\left(\gamma \right(\mathit{\theta }),\delta )}{\int \int }dA=\underset{\gamma }{\int }\underset{-\delta }{\overset{\delta }{\int }}Jd\xi d\phi ,$$

using the expansion for the Jacobian determinant, J, of the coordinate transformation from $(\phi ,\xi )$ (curve parameter and normal offset) to $(x,y)$:

$$J=\left|\begin{array}{cc}\dot{x}-\xi \left(\ddot{y}\Delta +\dot{y}\dot{\Delta }\right)& \dot{y}+\xi \left(\ddot{x}\Delta +\dot{x}\dot{\Delta }\right)\\ -\dot{y}\Delta & \dot{x}\Delta \end{array}\right|={\displaystyle \frac{1}{\Delta }}-\xi \left(\dot{x}\ddot{y}-\ddot{x}\dot{y}\right){\Delta }^2.$$

Substituting this relation into (5) brings it to a form similar to (4):

$${\displaystyle {\mathcal{J}}^s\left(\mathit{\theta }\right)=\frac{2\delta }{L\left(\mathit{\theta }\right)}\underset{\mathcal{S}\left(\gamma \right(\mathit{\theta }),\delta )}{\int \int }G\left(\frac{d(\mathit{x},\gamma )}{\sigma }\right)\tilde{I}(x,y)dA,}$$

This integral is approximated by summing contributions from triangular prisms formed by adjacent sample points:

$${\mathcal{J}}^s\left(\mathit{\theta }\right)\approx {\displaystyle \frac{2\delta }{L\left(\mathit{\theta }\right)}}\sum _{k=1}^{n-1}\sum _{p=-P}^{P-1}\left[V({\mathit{x}}_k^p,{\mathit{x}}_k^{p+1},{\mathit{x}}_{k+1}^p)+V({\mathit{x}}_k^{p+1},{\mathit{x}}_{k+1}^{p+1},{\mathit{x}}_{k+1}^p)\right],$$

where ${\mathit{x}}_k^p={\mathit{x}}_k+(p\delta /P){\mathit{n}}_k$ are the coordinates of a point along the normal ${\mathit{n}}_k$ at ${\mathit{x}}_k$, and $V(A,B,C)$ is the volume of a triangular prism with base $ABC$:

$$V(A,B,C)={\displaystyle \frac{G(d_A,\sigma )\tilde{I}\left(A\right)+G(d_B,\sigma )\tilde{I}\left(B\right)+G(d_C,\sigma )\tilde{I}\left(C\right)}{3}}·{\mathcal{S}}_{ABC}.$$

The first multiplier here is a height of the approximating prism, while the second one is $ABC$ triangle area, computed through coordinates of its corners:

$$S\left(A,B,C\right)={\displaystyle \frac{1}{2}}\left|\begin{array}{ccc}{A}_x& A_y& 1\\ B_x& B_y& 1\\ C_x& C_y& 1\end{array}\right|.$$

3.3. Discrete Differential Operators for Accurate Rasterization

The transition from continuous parametric curves to the discrete pixel grid is a critical step that, if done naively, can introduce systematic errors comparable to the inter-fringe spacing in high-density patterns. To achieve sub-pixel accuracy and topological correctness (4-connectivity), we introduce a unified approach based on integer-valued discrete differential operators. These operators analyze the local geometry of the curve directly in pixel-index space, enabling efficient and exact rasterization of both the curve and the surrounding strip.

3.3.1. First-Order Operator: Tangent Direction and Curve Rasterization

Given two consecutive pixel centers $P_{i-1}=(r_{i-1},c_{i-1})$ and $P_i=(r_i,c_i)$ along a discretized curve, we define the first-order discrete directional operator:

$$D_1(P_{i-1},P_i)=4+3{(\Delta r)}^3+{(\Delta c)}^3,\mathrm{where}\Delta r=r_i-r_{i-1},\Delta c=c_i-c_{i-1}.$$

This compact formula serves three purposes simultaneously:

• Direction defining: It maps the 9 possible neighbor vectors $(\Delta r,\Delta c)\in {\{-1,0,1\}}^2$ to unique integers in the range $[0,8]$. This mapping is shown in Figure 3.
• Neighborhood Filter: The sum of the cubic terms, $3{(\Delta r)}^3+{(\Delta c)}^3$ acts as a filter because for $(\Delta r,\Delta c)\in {\{-1,0,1\}}^2$, the identity $3{(\Delta r)}^3+{(\Delta c)}^3=3\Delta r+\Delta c$ holds, while if a parameter step is too large and causes a jump to a non-adjacent pixel, the sum makes $3{(\Delta r)}^3+{(\Delta c)}^3\ge 8)$ (within 2 pixels from $P_i$), flagging an invalid "skip".
• Adaptive Step Control: This flag triggers an adaptive adjustment of the parametric step $\Delta s$: $\Delta s\leftarrow \Delta s/k_{\mathrm{down}}$ if a skip is detected; $\Delta s\leftarrow \Delta s\times k_{\mathrm{up}}$ if $D_1=4$ (meaning the parameter advanced but the pixel did not change). Crucially, we use asymmetric coefficients ($k_{\mathrm{up}}=1.5$, $k_{\mathrm{down}}=1.4$) to prevent resonant cycling near stationary points, ensuring robust convergence.

The operator $D_1$ is the core of the curve rasterization algorithm. Its scheme is shown in the Appendix A (Algorithm 1). The algorithm marches along the continuous curve $\mathit{c}\left(s\right)$, adaptively adjusting $\Delta s$ to produce an ordered, 4-connected list of pixels $\mathcal{R}$. A key subtlety is the handling of diagonal moves: when the curve passes through a pixel corner, there are two equally valid 4-connected paths. In such cases, the algorithm use a special variable, which memorizes the choice made at the previous step, ensuring a smooth, consistent rasterization without jagged artifacts.

3.3.2. Second-Order Operator: Local Curvature and Strip Rasterization

Rasterizing a strip of width $2\delta$ around a curve requires understanding the local curvature to avoid gaps or overlaps in the coverage. Given three consecutive pixel centers $P_{i-1},P_i,P_{i+1}$ of the already rasterized curve, we define a second-order operator that classifies the local shape. Let $\Delta {\mathit{p}}_1=P_i-P_{i-1}$ and $\Delta {\mathit{p}}_2=P_{i+1}-P_i$. We define the signs of their components as integers: $s_x=\mathrm{sgn}(\Delta c_1)+2·\mathrm{sgn}(\Delta c_2)$ and $s_y=\mathrm{sgn}(\Delta r_1)+2·\mathrm{sgn}(\Delta r_2)$. The configuration is then defined as:

$$D_2=4s_y+s_x\in \{-15,-14,\dots ,14,15\}.$$

After shifting and clamping to the valid range of physically realizable 4-connected triples, we obtain the final formula used in the implementation:

$$D_2(P_{i-1},P_i,P_{i+1})=⌊{\displaystyle \frac{15-c_{i-1}-3c_i+4c_{i+1}-3r_{i-1}-9r_i+12r_{i+1}}{2}}⌋\in [0,15].$$

(7)

The operator $D_2$ evaluates to one of 16 integer values, each corresponding to a distinct local configuration of the three points (see Figure 4). Geometrically, it approximates the discrete curvature. For example:

• $D_2\in \{0,5,10,15\}$ indicates nearly linear motion (up, left, right, down), warranting a simple rectangular strip fill.
• $D_2\in \{1,2,4,7,8,11,13,14\}$ indicates a smooth turn, requiring a circular sector fill to cover the curved corner of the strip without gaps.
• $D_2\in \{3,6,9,12\}$ corresponds to inflection or "return" points, necessitating a combined fill (e.g., a half-circle plus strips).

Based on the value of $D_2$, the strip rasterization algorithm (Algorithm (2) in Appendix A) selects an optimal geometric primitive (line strip or circular arc) and fills the corresponding pixels within distance $\delta$ from the central pixel $P_i$. A distance mask ensures that if a pixel is covered by multiple primitives, the smallest distance to the central curve is retained, and the corresponding Gaussian weight $exp(-d^2/\left(2{\sigma }^2\right))$ is assigned. This approach guarantees a skip-free, accurate, and weight-aware rasterization of the continuous strip $\mathcal{S}\left(\gamma \right(\mathit{\theta }),\delta )$.

3.4. Parametric Models for Fringe Contours

The choice of the parametric subspace ${\mathcal{A}}_{\mathit{\theta }}$ is dictated by the expected physics of the deformation. For a broad class of problems involving axisymmetric or nearly axisymmetric bending of plates and membranes, the fringe contours are expected to be closed, nested, quasi-concentric curves. We introduce two related parametric families that efficiently described such shapes while offering sufficient flexibility to capture deviations from ideal geometry.

3.4.1. Perturbed Circle Model

The most natural model for quasi-circular fringes is a perturbed circle, where the constant radius of a circle is replaced by a periodic function of the polar angle $\phi$. In Cartesian coordinates with an origin at $(p_x,p_y)$, this is expressed as:

$$\begin{array}{cc}\hfill x\left(\phi \right)& =p_x+R\left(\phi \right)cos\phi ,\hfill \\ \hfill y\left(\phi \right)& =p_y+R\left(\phi \right)sin\phi ,\hfill \end{array}$$

(8)

where the radius function $R\left(\phi \right)$ is given by a trigonometric polynomial:

$$R\left(\phi \right)=r\left[1+\sum _{k=1}^n\left(a_{k}cos\left((k+1)\phi \right)+b_{k}sin\left((k+1)\phi \right)\right)\right].$$

Here, $r>0$ is a base radius, $(p_x,p_y)$ is the center, and ${\{a_k,b_k\}}_{k=1}^n$ are the perturbation coefficients. Setting all $a_k=b_k=0$ recovers a perfect circle of radius r centered at $(p_x,p_y)$. The shift in harmonic indices (using $(k+1)\phi$ instead of $k\phi$) ensures that the lowest-order perturbation already affects the shape, leaving the coefficients $a_0$ and $b_0$ implicitly absorbed into r and the center coordinates.

An equivalent, sometimes computationally convenient, representation expands (8) directly into a Fourier series:

$$\begin{array}{cc}\hfill x\left(\phi \right)& =A_0+\sum _{k=1}^{n+1}\left(A_{k}cos\left(k\phi \right)+B_{k}sin\left(k\phi \right)\right),\hfill \\ \hfill y\left(\phi \right)& =C_0+\sum _{k=1}^{n+1}\left(C_{k}cos\left(k\phi \right)+D_{k}sin\left(k\phi \right)\right),\hfill \end{array}$$

where the coefficients $\{A_k,B_k,C_k,D_k\}$ are linear combinations of $\{r,p_x,p_y,a_k,b_k\}$. It is simple to show that $A_0$, $C_0$, and the first harmonics are primarily governed by the “coarse” parameters $r,p_x,p_y$, while higher harmonics depend on the perturbation coefficients. This structure naturally suggests a two-stage optimization strategy:

• Quasi-optimization: Vary only the coarse parameters $(r,p_x,p_y)$ to quickly locate the approximate position and scale of the fringe.
• Full optimization: Refine all parameters $(r,p_x,p_y,\{a_k,b_k\})$ to capture fine shape details.

The total number of parameters is $m=3+2n$. In our experiments with plate bending, $n=15$ (giving 33 parameters) proved sufficient to accurately represent fringes even under significant nonlinear deformation.

3.4.2. Perturbed Ellipse Model

For specimens with pronounced anisotropy or non-axisymmetric boundary conditions, the fringe contours may exhibit a preferred elliptical orientation. The model can be generalized to a perturbed ellipse by introducing semi-axes a and b ($a\ge b>0$) and an orientation angle $\psi$:

$$\begin{array}{cc}\hfill x\left(\phi \right)& =p_x+\rho (\phi ;a,b,\psi )\left[1+\sum _{k=1}^n\left(a_{k}cos\left((k+1)(\phi -\psi )\right)+b_{k}sin\left((k+1)(\phi -\psi )\right)\right)\right]cos\phi ,\hfill \\ \hfill y\left(\phi \right)& =p_y+\rho (\phi ;a,b,\psi )\left[1+\sum _{k=1}^n\left(a_{k}cos\left((k+1)(\phi -\psi )\right)+b_{k}sin\left((k+1)(\phi -\psi )\right)\right)\right]sin\phi ,\hfill \end{array}$$

where

$$\rho (\phi ;a,b,\psi )={\displaystyle \frac{ab}{\sqrt{(b^2-a^2){cos}^2(\phi -\psi )+a^2}}}.$$

When $a=b=r$, this relation reduces to $\rho =r$, and the model reverts to the perturbed circle (8). The perturbation is now applied relative to the elliptical base shape and is oriented with it (via the phase shift $\phi -\psi$). The coarse parameter set for this model is $(a,b,\psi ,p_x,p_y)$, totaling $5+2n$ parameters.

3.4.3. Parameter Constraints and Initialization

To ensure physical plausibility and improve optimization convergence, constraints are applied: $r,a,b>0$; perturbation coefficients are bounded ($|a_k|,|b_k|<0.5$ to prevent self-intersection); and for the ellipse, $a\ge b$. The initial guess for the first (innermost) fringe is obtained via a coarse search: the image center is estimated from intensity function cross-sections ($\phi =\mathrm{const}$), then r (and $a,b$ later) is varied until a local extremum of ${\mathcal{J}}^s$ is found. For subsequent fringes, the parameters of the previous fringe are scaled outward by a factor slightly greater than one to provide a good starting point for the quasi-optimization step.

These parametric models, combined with the strip-integration functional and the rasterization operators, constitute the proposed skeletonization framework. However, before proceeding to the experimental validation of its performance in Section 4, it is necessary to prepare the intensity matrix to reduce the influence of speckle noise.

3.5. Local Intensity Equalization Algorithm

Raw interferograms often exhibit slow spatial variations in background intensity and contrast due to uneven illumination and polarization effects (see Subsection 2.1). To normalize these variations without distorting the fringe positions, we employ a local intensity equalization procedure based on quantile statistics within sliding windows.

3.5.1. Local Quantile Computation

Before processing the intensity matrix it is necessary to remove areas of the image without fringes, that contain no fringes, retaining only the useful region. Furthermore geometric distortion correction must be applied to compensate for the misalignment between the camera lens axis and the object beam. The example of this procedure, which we will call the geometric correction, is shown in Figure 5.

Let $\mathbf{I}[r,c]$ denote the geometrically corrected intensity matrix of size $r_{max}\times c_{max}$. For each pixel $(r_q,c_q)$ considered as a potential window center, we define a local block $M_{r_q,c_q,q}$ of size $q\times q$:

$$M_{r_q,c_q,q}[i,j]=\mathbf{I}\left[R(r_q+i-\lfloor q/2\rfloor ),C(c_q+j-\lfloor q/2\rfloor )\right],i,j=1,\dots ,q,$$

where $R\left(k\right)$ and $C\left(k\right)$ are reflection-padding functions that handle boundary pixels:

$$R\left(k\right)=\left\{\begin{array}{cc}k,\hfill & 1\le k\le r_{max}\hfill \\ 1-k,\hfill & k<1\hfill \\ 2r_{max}-k,\hfill & k>r_{max}\hfill \end{array}\right.,C\left(k\right)=\left\{\begin{array}{cc}k,\hfill & 1\le k\le c_{max}\hfill \\ 1-k,\hfill & k<1\hfill \\ 2c_{max}-k,\hfill & k>c_{max}\hfill \end{array}\right..$$

This padding ensures that windows near image boundaries remain fully populated without introducing artificial discontinuities.

For each block $M_{r_q,c_q,q}$, we compute two statistics: the lower quantile $a[r_q,c_q]$ and the upper quantile $b[r_q,c_q]$ at probability levels $\nu$ and $1-\nu$, respectively (typically $\nu =0.05$). These quantiles estimate the local minimum and maximum intensity while being resistant to outlier pixels caused by speckle noise:

$$a[r_q,c_q]=Q_{\nu }\left(M_{r_q,c_q,q}\right),b[r_q,c_q]=Q_{1-\nu }\left(M_{r_q,c_q,q}\right).$$

3.5.2. Interpolation and Normalization

The quantiles $a[r,c]$ and $b[r,c]$ are computed only on a coarse grid with spacing $d=\lfloor q/2\rfloor$ to reduce computational cost. To obtain values at every pixel $(i,j)$, we use linear interpolation:

$$\widehat{a}\left(i,j\right)=\mathrm{Int}\left(\begin{array}{cccc}\mathrm{a}\left[1,1\right]& \mathrm{a}\left[1,\mathrm{d}+1\right]& \dots & \mathrm{a}\left[1,{\mathrm{c}}_{\max }\right]\\ \mathrm{a}\left[\mathrm{d}+1,1\right]& \mathrm{a}\left[\mathrm{d}+1,\mathrm{d}+1\right]& \dots & \mathrm{a}\left[\mathrm{d}+1,{\mathrm{c}}_{\max }\right]\\ \dots & \dots & \dots & \dots \\ \mathrm{a}\left[{\mathrm{r}}_{\max },1\right]& \mathrm{a}\left[{\mathrm{r}}_{\max },\mathrm{d}+1\right]& \dots & \mathrm{a}\left[{\mathrm{r}}_{\max },{\mathrm{c}}_{\max }\right]\end{array}\right),$$

$$\widehat{b}\left(i,j\right)=\mathrm{Int}\left(\begin{array}{cccc}\mathrm{b}\left[1,1\right]& \mathrm{b}\left[1,\mathrm{d}+1\right]& \dots & \mathrm{b}\left[1,{\mathrm{c}}_{\max }\right]\\ \mathrm{b}\left[\mathrm{d}+1,1\right]& \mathrm{b}\left[\mathrm{d}+1,\mathrm{d}+1\right]& \dots & \mathrm{b}\left[\mathrm{d}+1,{\mathrm{c}}_{\max }\right]\\ \dots & \dots & \dots & \dots \\ \mathrm{b}\left[{\mathrm{r}}_{\max },1\right]& \mathrm{b}\left[{\mathrm{r}}_{\max },\mathrm{d}+1\right]& \dots & \mathrm{b}\left[{\mathrm{r}}_{\max },{\mathrm{c}}_{\max }\right]\end{array}\right).$$

The normalized intensity at each pixel is then obtained by linear rescaling:

$$\tilde{\mathbf{I}}[i,j]={\displaystyle \frac{\mathbf{I}[i,j]-\widehat{a}(i,j)}{\widehat{b}(i,j)-\widehat{a}(i,j)}}.$$

This operation maps the local intensity range $[\widehat{a}(i,j),\widehat{b}(i,j)]$ approximately to $[0,1]$, effectively compensating for low-frequency inhomogeneities while preserving the high-frequency fringe structure.

3.5.3. Parameter Selection

The window size q is chosen to be several times larger than the expected inter-fringe distance but smaller than characteristic scales of illumination non-uniformity. In our experiments with square plates, $q=31$ pixels (covering about 3–4 fringes) worked well. The quantile level $\nu =0.05$ provides a good compromise between rejecting speckle outliers and retaining true fringe extrema. The equalization is applied once after geometric correction and, if needed, repeated after optional Fourier filtering to compensate for global intensity shifts introduced by the filter.

3.6. Complete Skeletonization Algorithm

Integrating all components described in the previous subsections, we present the complete procedure for automated fringe skeletonization. The algorithm proceeds recursively from the innermost to the outermost fringe, taking into account the quasi-similarity property. The main steps are illustrated in Figure 6. The algorithm iteratively identifies fringes until one of the three termination conditions is met:

• The last computed intensity integral is equal to zero;
• The last identified fringe lies outside the image boundaries;
• The required (or preset) number of fringes has been identified;

3.7. Synthetic Interferogram Generation

To quantitatively validate the proposed algorithm under controlled conditions, we developed a virtual interferogram generator that produces images with known pattern geometry and physically realistic noise characteristics. The generator taking into account two key aspects: (1) the typical fringe pattern of a bent square plate, and (2) the corruption mechanisms described in Subsection 2.1.

3.7.1. Ideal Fringe Geometry via Morphing

The underlying displacement field is modeled by a family of closed curves that morph smoothly from a circle at the center to a square at the boundary, parameterized by a morphing parameter $\alpha \in [0,1]$. In the first octant ($\phi \in [0,\pi /4]$), the curve consists of a straight segment and a curved segment that ensures $C^1$ continuity:

$$\begin{array}{cccc}\hfill r_{\mathrm{s}}(\phi ;\alpha )& =\alpha sec\phi ,\hfill & \hfill \phi & \in \left(0,\alpha \pi /4\right],\hfill \\ \hfill r_{\mathrm{c}}(\phi ;\alpha )& =\alpha \left[{\displaystyle \frac{2}{\pi (1-\alpha )}}\left({\displaystyle \frac{\alpha \pi }{4}}-\phi \right)\left({\displaystyle \frac{\alpha \pi }{4}}+\phi -{\displaystyle \frac{\pi }{2}}\right)tan{\displaystyle \frac{\alpha \pi }{4}}+1\right]sec{\displaystyle \frac{\alpha \pi }{4}},\hfill & \hfill \phi & \in \left(\alpha \pi /4,\pi /4\right].\hfill \end{array}$$

The full curve over all angles is obtained by symmetric replication. Equivalently, it can be expressed compactly as:

$$r(\phi ;\alpha )=\left\{\begin{array}{cc}\alpha r_1\left(\phi mod{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi \alpha }{4}}\right),\hfill & 0<\phi mod{\displaystyle \frac{\pi }{2}}\le {\displaystyle \frac{\pi }{4}},\hfill \\ \alpha r_2\left(\phi mod{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{4}}(2-\alpha )\right),\hfill & {\displaystyle \frac{\pi }{4}}<\phi mod{\displaystyle \frac{\pi }{2}}\le {\displaystyle \frac{\pi }{2}},\hfill \end{array}\right.$$

where

$$\begin{array}{cc}\hfill & r_1(\phi ,\alpha )=\left\{\begin{array}{cc}\left[{\displaystyle \frac{\alpha -\phi }{\pi -4\alpha }}(2\alpha +2\phi -\pi )tan\alpha +1\right]sec\alpha ,\hfill & \phi >\alpha \hfill \\ sec\phi \hfill & \phi \le \alpha \hfill \end{array}\right.,\hfill \\ \hfill & r_2(\phi ,\alpha )=\left\{\begin{array}{cc}csc\phi ,\hfill & \phi >\alpha \hfill \\ \left[{\displaystyle \frac{\alpha -\phi }{\pi -4\alpha }}(\pi -2\alpha -2\phi )cot\alpha +1\right]csc\alpha ,\hfill & \phi \le \alpha \hfill \end{array}\right..\hfill \end{array}$$

The ideal, noise-free intensity pattern is then defined as a periodic function of $\alpha$:

$$I_{\mathrm{ideal}}(x,y)={sin}^2\left(k\pi \alpha (x,y)\right),$$

where k controls the fringe density. The mapping $(x,y)\mapsto \alpha$ is obtained by numerically solving the transcendental equation:

$$r(\phi ,\alpha )=\sqrt{x^2+y^2},$$

for $\alpha$ at each grid point.

3.7.2. Physically Motivated Noise Model

The ideal intensity is corrupted by two noise sources derived from the physical model in Subsection 2.1, which presents the general formula for corrupted intensity (1); however the terms describes speckle noisy, were not specified. Accurately defining of these terms is a complicated problem, which was considered in a huge amount of works (e. g. [16,17,18]). For this paper, it is enough to mark that speckle noise consists of three parts, which related with reflected surface factor, depended on its roughness, spatial and temporal coherence:

$$N_{\mathrm{temp}}=1+{\left(2\Delta L\Delta \nu /c\right)}^2,N_{\mathrm{spat}}=1+{\left(D/{\rho }_{\mathrm{speckle}}\right)}^2,$$

where c is the speed of light in vacuum, $\Delta \nu$ is the laser spectral linewidth, $\Delta L$ the path length difference, D is the diameter of the illuminated object area and ${\rho }_{\mathrm{speckle}}$ is the characteristic speckle size:

$${\rho }_{\mathrm{speckle}}\approx max\left(\lambda z/D,L_{c}z/D\right).$$

Here, $\lambda$ is the laser wavelength, z is the distance from object to observation plane and $L_c$ is the reflected surface roughness parameter. The resulting speckle noise can be modeled as the product of different coherence terms and the surface factor, which is obtained empirically.

In practice, for rapid commutation and parameter studies, it is convenient to employ a simplified combined noise model, constructed as stochastic perturbations of phase difference and amplitude:

$$\tilde{I}(x,y)=I_0(x,y)+k\nu (x,y)+2\sqrt{k\nu (x,y)I_0(x,y)}cos\left(\delta \pi \right),$$

where k controls the overall noise level and $\nu (x,y)$, $\delta (x,y)$ are independent stochastic processes with correlation:

$$\langle f\left(\mathit{x}\right)f\left({\mathit{x}}^{\prime }\right)\rangle ={\left[2J_1\left(\pi |\mathit{x}-{\mathit{x}}^{\prime }|/{\rho }_{\mathrm{speckle}}\right)/\left(\pi |\mathit{x}-{\mathit{x}}^{\prime }|/{\rho }_{\mathrm{speckle}}\right)\right]}^2,$$

where $J_1(·)$ is the Bessel function. This model approximates the statistical properties of speckle while remaining computationally efficient.

3.7.3. Error Metrics

For each synthetic image, the exact curves $\left\{{\gamma }_i^{\mathrm{true}}\right\}$ are known analytically. This allows us to define rigorous error metrics for any skeletonization result $\left\{{\gamma }_i^{\mathrm{res}}\right\}$:

$$\begin{array}{cccc}\hfill & \mathbf{Euclidean} \mathbf{metric}:\hfill & \hfill & {\epsilon }_{\mathrm{Euc}}=\sqrt{{\displaystyle \frac{1}{L^{\mathrm{true}}}}\underset{0}{\overset{2\pi }{\int }}\left[{\left(x^{\mathrm{true}}\left(\phi \right)-x^{\mathrm{res}}\left(\phi \right)\right)}^2+{\left(y^{\mathrm{true}}\left(\phi \right)-y^{\mathrm{res}}\left(\phi \right)\right)}^2\right]d\phi },\hfill \\ \hfill & \mathbf{Max} \mathbf{error}:\hfill & \hfill & {\epsilon }_{max}={\displaystyle \frac{\underset{\phi \in \left[0,2\pi \right)}{max}\sqrt{{\left(x^{\mathrm{true}}\left(\phi \right)-x^{\mathrm{res}}\left(\phi \right)\right)}^2+{\left(y^{\mathrm{true}}\left(\phi \right)-y^{\mathrm{res}}\left(\phi \right)\right)}^2}}{\underset{\phi \in \left[0,2\pi \right)}{max}\sqrt{{\left(x^{\mathrm{true}}\left(\phi \right)\right)}^2+{\left(y^{\mathrm{true}}\left(\phi \right)\right)}^2}}}.\hfill \end{array}$$

These metrics are used in Section 4 to quantitatively compare the performance of different methods.

3.8. Implementation Details.

The described skeletonization procedure has been implemented in a custom C++ program. This program is capable of performing all steps of the proposed algorithm (see Figure 6):

• Geometric correction;
• Filtering in the frequency domain with different kernels;
• Intensity equalization, as described in Subsection 3.5;
• Localization of the image center;
• Recursive identification of fringes, which includes:

  (a)

  Initializing of parametric curve (Subsection 3.4.1 and Section 3.4.2);

  (b)

  Computing of the intensity integral (Subsection 2.4) using numerical quadratures (Subsection 3.2);

  (c)

  Defining of the approximation curve’s parameters by one of the optimization methods:

  • Advanced Coordinate Descent with Nesterov’s Acceleration for fastest convergence (typically 20-30 iterations per fringe) (e. g. [19]).
  • Nelder – Mead simplex method when derivative information is not readily available (robust but slightly slower).
  • Conjugate Gradient method for highest precision in final refinement (used only when the highest accuracy is required, despite longer computation time).

An example of the preprocessing of a real interferogram, performed using the developed program, is shown in Figure 7.

Apart from that, the program contains a set of standard tools for working with images and includes a custom generator of synthetic fringe patterns, a described in Subsection 3.7. The generator accepts the following control parameters:

• Fringe density k (number of fringes across the field);
• Speckle noise, dependent on surface roughness parameter and the laser spectral linewidth;
• Noise models the influence of non-uniform illumination;

For each generated image, the exact curves $\left\{{\gamma }_i^{\mathrm{true}}(\alpha =i/k)\right\}$ are known analytically. This allows for the quantification of the algorithm errors.

The main computational cost of the algorithms associated with the strip-integrating procedure. Each evaluation of ${\mathcal{J}}^s\left(\mathit{\theta }\right)$ requires $O\left(M\right)$ operations, where M is the number of quadrature points in the strip discretization. With N fringes and an average of K optimization iterations per fringe, the total cost is $O\left(KNM\right)$. In practice, for a $1000\times 1000$ image with 100 fringes, the algorithm completes in 2 – 3 minutes on a standard desktop CPU, making it suitable for practical laboratory use.

4. Validation on Synthetic Interferograms

We compared the efficiency of fringe identification with different variants of the intensity integral, using the metrics defined in Subsection 3.7.3. Figure 8 demonstrates that slightly noisy images, both methods have similar performance, whereas for highly noisy patterns the stripe-integration method has a significant advantage over the line-integration. Moreover, as it is shown in Figure 9, the stripe-integration method is capable of identifying fringes in images with noise levels as high as 95%, where it is hardly possible to differ the fringes even visually.

5. Application to Real High-Density Interferograms

To validate the proposed algorithms, it was tested on a real hologram, obtained from a bending experiment on a square plate. The plate (copper, side length 6 cm, thickness 184 $\mu$m) was subjected to a uniform transverse load in a holographic interferometry setup (Leith – Upatnieks off-axis scheme with a solid-state laser, $\lambda =532$ nm). Interferograms were recorded for several load increments. The most challenging case, with a deflection generating over 100 fringes (shown in Figure 10i)), was selected for analysis.

The set of identified fringes (red and blue lines in Figure 10ii) and Figure 10iii)) was compared with the results (Figure 10iv)), obtained from theoretical modeling, presented in [15]. As shown in Figure 10, the identified fringes are in good agreement with isolines of deformed surface, demonstrating the acceptable accuracy of the proposed procedure.

6. Conclusion

This paper presents a novel, robust method for skeletonizing fringe patterns in holographic interferometry. By constraining the solution to a physics-informed parametric subspace and employing a strip-integration functional, the method achieves high accuracy and stability in the presence of strong speckle noise, where conventional techniques may fail. Quantitative validation on synthetic data demonstrated a significant reduction in error compared to baseline methods. The practical utility was confirmed by successfully processing a real interferogram with over 100 fringes, the results of which closely matched those of an independent theoretical model.

The main limitation of the algorithm is the need to select an appropriate parametric model ${\mathcal{A}}_{\theta }$ beforehand. Thus, the method is most effective when the fringe family topology is known to be a set of nested, simply connected curves. It may be less effective for patterns featuring fringe bifurcations, interruptions, or multiple disconnected families without significant modifications. In the future, the concept of quasi-similarity and recursive search may be modified to a more powerful paradigm for handling complex topologies.