Reconstructing Nerve Structures from Unorganized Points

Abstract

Realistic sensory feedback is paramount for amputees as it improves prosthetic limb control and boosts functionality, safety, and overall quality of life. This sensory restoration relies on the direct electrostimulation of residual peripheral nerves. Computational models are instrumental in simulating these neurostimulation effects, offering solutions to the complexities tied to extensive animal/human trials and costly materials. Central to these models is the detailed mapping of nerve geometry, necessitating the delineation of internal nerve structures, such as fascicles, across various cross-sections. In our modeling process, we faced the challenge of organizing an originally unstructured set of points into coherent contours. We introduced a parameter-free curve-reconstruction algorithm that combines valley-seeking clustering, an adaptive Kalman filter, and the nearest neighbor classification technique. While intuitively simple for humans, the task of reconstructing multiple open and/or closed lines with pronounced corners from a nonuniform point set is daunting for many algorithms. Additionally, the precise differentiation of adjacent curves, commonly encountered in realistic nerve models, remains a formidable challenge even for top-tier algorithms. Our proposed method adeptly navigates the complexities inherent to nerve structure reconstruction. While our algorithm is chiefly designed for closed curves, as dictated by nerve geometry, we believe it can be reconfigured with appropriate code adjustments to handle open curves. Beyond neuroprosthetics, our proposed model has the potential to be applied and spark innovations in biomedicine and a variety of other fields.

1. Introduction

There are a number of tasks we humans tend to solve effortlessly and almost intuitively, which machines find complicated and often require intricate algorithms for solving them. One such task is connecting seemingly randomized points into a meaningful shape. This type of 2D reconstruction is fundamental across a myriad of application areas: digital image processing, computer vision, reverse engineering, material quality control, 3D medical image reconstruction, and handwriting recognition, among others. In our case, the need for the development of such a technique for line reconstruction emerged during the in silico modeling of a human peripheral nerve.

A 2017 study showed that approximately 57.7-million people worldwide were living with limb amputation caused by a traumatic injury [1]. This study did not account for other common causes of extremity loss, such as diabetic neuropathy, cancer, infection, or burn injuries. Recent advances in the field of neurotechnology have enabled scientists to develop new methods that could help amputees [2]. Today, there are numerous prostheses available for both upper and lower limbs. However, there is significant room for improvement in restoring sensory feedback to amputees.

Nerve fibers are responsible for transmitting information between the Central Nervous System (CNS) and the peripheral effectors, and vice versa. These fibers bundle together to form nerve fascicles, which are essentially groups or bundles of nerve fibers within a nerve. Fascicles often exhibit intricate patterns of crossings and branchings. Recent studies have highlighted that direct electrostimulation of the remaining portions of a peripheral nerve can successfully generate physiologically relevant sensations, such as touch, pressure, vibration, and tingling [3,4,5,6,7]. Owing to the significant variability in peripheral nerve anatomy and its unique characteristics, it is unfeasible to standardize electrode positioning and stimulation protocols. Instead of resorting to clinical animal or human trials, which are generally complicated, risky, and bring a high cost, in silico modeling offers a promising alternative [3,8,9,10]. The combination of the Finite Element Method (FEM) with Neuron Modeling (hybrid electro-neural modeling) could simulate the effects of the electrical stimulation of residual peripheral nerves, i.e., the nerve fibers (axons) contained within, thus facilitating electrode positioning and stimulation paradigm optimization. Detailed and realistic nerve structure modeling would undoubtedly enhance the precision of these models. An in silico method for nerve contour reconstruction would certainly be an integral step. We introduced a technique specifically designed for this purpose, but one that could, with some minor adjustments, also find its use in a handful of other disciplines.

The persistent challenge of curve (re)construction from an unorganized set of points, coupled with the need for a reliable solution, has motivated numerous scientists to address this issue. However, this is not an easy task. Its complexity largely depends on the assignment, particularly the characteristics of the point set and the desired curve shape. A nonuniform point distribution, the high proximity of contours, as well as sharp corners they could contain are just some of the most-common problems such algorithms can encounter. However, few have successfully tackled these challenges. Existing techniques vary in their performance and in how they formulate and approach these tasks. More about line types each of them is able to reconstruct, their similarities, and their differences will be discussed in Section 2.

Section 3 presents the motivation for this study and gives a detailed insight into the problem that needed to be addressed during realistic nerve contour reconstruction. Section 4 describes the innovative way in which the authors approached the problem, which combines techniques such as valley-seeking clustering, the adaptive Kalman filter, and Nearest Neighbor (NN) classification.

Section 5 presents the reconstruction results generated in Matlab and the comparison with the representative state-of-the-art techniques, as well as outlines key considerations for future implementations of this algorithm. Finally, Section 6 sums up the advantages of the proposed algorithm and puts it into the context of potential future improvements and usages.

2. Related Work

During the last several decades, numerous methods dedicated to line reconstruction from a set of unstructured points have been presented [11]. In general, these techniques can be classified based on their ability to handle lines, as well as sets of points, of different properties. In this paper, we highlight the most-prominent techniques, differentiated by their responses to the following questions:

• Is the point distribution uniform or nonuniform?
• Are the lines smooth, or do they contain sharp corners?
• Can the method handle both open and closed curves?

Additionally, while some algorithms can reconstruct only a single contour at a time, others are capable of detecting multiple lines. The majority of the techniques focus only on curves without self-intersection.

The earlier techniques, such as $\alpha$-shapes [12,13], $\beta$-skeletons [14], minimum spanning trees [15], and r-regular shapes [16], were limited to a uniform point distribution and could reconstruct only smooth lines, without sharp corners.

The following methods brought improvement in the aspect of working with a nonuniform set of points: Crust [17,18], Nearest Neighbor Crust (NN-Crust) [19], Conservative Crust [20], Gathan [21], GathanG [22], Discur [23], and Vicur [24], as well as the algorithms that are based on the Traveling Salesman Problem (TSP) [25,26].

The Crust algorithm, first introduced in [17], constructs a graph on a planar point set. It is designed for the reconstruction of smooth closed curves, under the condition that the points are dense enough in areas with higher curvature. Gold [18] introduced a simplified version of the original Crust algorithm. Similarly, Dey and Kumar [19] proposed their NN-Crust method, which introduces the nearest neighbor technique. The following method Conservative Crust [20] is suitable for the reconstruction of both open and closed lines.

However, the biggest difficulty for all the above-mentioned methods is the reconstruction of sharp angles. The problem occurs when two non-adjacent points are nearest neighbors, which leads to their incorrect connection (Figure 1). Giesen was the first to offer a solution to this problem [25], developing a reconstruction algorithm that could reconstruct open or closed non-smooth curves. This method, as well as the method proposed by Althaus and Mehlhorn [26] are based on the Traveling Salesman Problem (TSP). However, due to the way the TSP algorithm works, neither one of those two methods can reconstruct more than a single curve at a time, while the latter method can handle only closed curves.

Dey and Wenger proposed a method called Gathan [21], which can work with curves that contain sharp corners. The same authors then presented GathanG [22], which guarantees closedcurve reconstruction, with both smooth and corner parts, under specified conditions. Funke and Ramos [27] introduced a method able to reconstruct both open and/or closed single or multiple curves with or without sharp corners. However, its limited adoption is often attributed to specific sampling conditions and the lack of open-source code availability.

The algorithms Discur [23] and Vicur [24] bring a new perspective, inspired by the way human vision works: nearness (two nearest points are usually connected) and smoothness (the connected points tend to form a smooth curve). Both techniques enable reconstruction of single or multiple open and closed curvesand/or lines with sharp corners. Discur requires a very dense sampling near the sharp corner area, while Vicur brings an advancement in that aspect.

The method HNN-Crust, introduced in [28], deals with sampling conditions and improves them. Ohrhallinger et al. showed that their algorithm requires less-dense sampling than the prior methods of the same type (the Crust family of proximity-based reconstruction algorithms). The main deficiency of this method is it not being able to successfully reconstruct contours with sharp corners or two lines in close proximity.

3. Motivation of This Study

The method described in this paper was developed while working on a framework called ProprioStim [29], which focuses on planning biomimetic stimulation strategies to restore natural proprioceptive sensation, primarily for individuals with peripheral nerve damage, such as amputees. One of the framework’s key contributions is its departure from the traditional methods of generating in silico nerve models. The mentioned study demonstrated that the innovative approach of 3D modeling, which faithfully reconstructs peripheral nerves with consideration for their intricate complexity, branching, and internal nerve structures, yields significantly different results compared to the commonly used extruded models [8,9,30]. This approach substantially enhances the accuracy of the simulation.

During the development of ProprioStim, previously collected experimental data were utilized, including histological images of transverse sections of the sciatic nerves from two transfemoral amputees (as the one shown in Figure 2A) [4]. They were utilized for CAD modeling of the detailed fascicular nerve geometry (Figure 2B,C), which was then used to create an FEM model in the Comsol software (version 5.5). The spatial distribution of electrical potentials, obtained during the process, along with the nerve geometry were imported into Matlab for the purpose of calculating the recruitment of individual nerve fibers within the fascicle. Both nerves had similar cross-sectional dimensions (approximately 20 mm × 7.5 mm). The acquired histologies allowed for the modeling of nerve segments with a length of 20 mm along the longitudinal (z) axis. The nerve segments were discretized into 19 equidistant cross-sections for manipulation in Matlab, with a layer-to-layer spacing of 100 $\mathsf{\mu }$m.

The motivation behind the method presented in this paper emerged as a result of the challenges encountered during the data transfer of nerve geometry from Comsol to Matlab. The vertices of the contours, which correspond to cross-sectional slices of nerve fascicles, initially became randomized upon import into Matlab, requiring a sorting process. Figure 3 depicts an example of a problem this paper focused on. The contour should be (re)constructed as precisely as possible from the nonuniformly distributed sample points (Figure 3A). Figure 3B illustrates its expected shape. If the points were connected according to the initial randomized order, the resulting contour would look as depicted in Figure 3C. Figure 3D shows an enlarged view of the reconstruction presented in Figure 3C, while four randomly chosen points are shown as references in both Figure 3C,D.

4. Materials and Methods

4.1. Initializing the Algorithm and the Clusters

The algorithm works recurrently and processes each cross-sectional level of each nerve fascicle individually. It starts by loading all the points from one level. As the points are initially unsorted, the goal is to organize and group them into clusters, i.e., connect them into lines. The initial point of the first cluster can be sampled from the input set randomly. Each following cluster from the same level is initialized by the first remaining, i.e., unsorted point from the input set.

The part of the algorithm that sorts points into corresponding clusters is based on the valley-seeking clustering method. This is a nonparametric technique that estimates the gradient of a density function [31,32,33]. Since a higher density of this function means more samples, our method’s objective is to find the regions with the highest density function and, thus, detect individual clusters.

The graph theoretic valley-seeking technique involves the following procedure: an area $G\left(X\right)$ around point X is selected (Figure 4):

$$G\left(X\right)=\{Y:d(Y,X)\le r\},$$

(1)

where r is the radius of the area $G\left(X\right)$. A measure of the distance $d(Y,X)$ is specified by:

$$d^2(Y,X)={(Y-X)}^T{A}^{-1}(Y-X),$$

(2)

where A is the metric to measure the distance. If the data about the distribution statistics are insufficient, the Euclidean metric $A=I$ can be used, in which case the area $G\left(X\right)$ will take the form of a circle. This assumption is applicable to the problem we tackled (Figure 4).

All the points that can be found inside the area $G\left(X\right)$ are taken into consideration as potential candidates for the point X’s successor. The local mean $M\left(X\right)$ can be calculated as:

$$M\left(X\right)=E\{(Y-X)\mid G\left(X\right)\}={\int }_{G\left(X\right)}(Y-X)\frac{f\left(Y\right)}{u_0}\mathrm{d}Y,$$

(3)

where $f\left(Y\right)$ is the a priori density function in the point Y and $u_0$ is the probability that Y belongs to the area $G\left(X\right)$ [32]. Under the assumption that the area $G\left(X\right)$ is small enough, the probability $u_0$ can be expressed as:

$$u_0={\int }_{G\left(X\right)}f\left(Y\right)\mathrm{d}Y\cong f\left(X\right)v,$$

(4)

where v would be the volume of the area $G\left(X\right)$. By implementing Equation (4) in Equation (3) and expanding $f\left(Y\right)$ from Equation (3) by a Taylor series, we can write the following:

$$\begin{array}{c}\hfill M\left(X\right)\cong {\int }_{G\left(X\right)}(Y-X)\frac{f\left(X\right)+{(Y-X)}^T▽f\left(X\right)}{f\left(X\right)v}\mathrm{d}Y=\frac{r^2}{n+2}A\frac{▽f\left(X\right)}{f\left(X\right)},\end{array}$$

(5)

where $n=dim\left\{X\right\}$ and $▽f\left(X\right)$ is the gradient of the density function in the area $G\left(X\right)$. From Equation (5), we can derive the normalized gradient:

$$\frac{▽f\left(X\right)}{f\left(X\right)}\cong \frac{n+2}{r^2}{A}^{-1}M\left(X\right),$$

(6)

where $M\left(X\right)$ can be estimated based on the sample of all the points contained in the area $G\left(X\right)$, as their mean value.

For finding the next point of the cluster, i.e., the successor of the point X, we can define a measure $S_{xy}$:

$$S_{xy}=∥\begin{array}{c}\frac{▽f\left(X\right)}{f\left(X\right)}\end{array}∥cos{\theta }_{xy},$$

(7)

where $S_{xy}$ is defined for each point Y from the area $G\left(X\right)$. $∥\begin{array}{c}\frac{▽f\left(X\right)}{f\left(X\right)}\end{array}∥$ is defined by Equation (6), while ${\theta }_{xy}$ is the angle between vectors $∥\begin{array}{c}▽f\left(X\right)\end{array}∥$ and $(Y-X)$. The point that maximizes the measure $S_{xy}$ is chosen as X’s successor. In practical terms, for all the analyzed points surrounding point X, the first portion of the expression in Equation (7), $∥\begin{array}{c}\frac{▽f\left(X\right)}{f\left(X\right)}\end{array}∥$, remains constant with each new point selection. This observation underscores that the measure $S_{xy}$ and the rationale behind defining the successor are solely influenced by the value of $cos{\theta }_{xy}$.

The proper choice of the area $G\left(X\right)$ is of great importance for the method’s initialization. The key parameter that affects the size of this area is r from (1). This parameter should be based on the points’ density. Based on their experience, the authors suggest that it should be estimated as:

$$r=k\times \frac{max(X_{max}-X_{min},Y_{max}-Y_{min})}{N/2},$$

(8)

where $(X_{max},X_{min})$ and $(Y_{max},Y_{min})$ are pairs of maximum and minimum values of all the points from the input set by the x-axis and y-axis, respectively, and N is the total number of points in the input set. The parameter k is calculated as:

$$k=\left\{\begin{array}{cc}4,\hfill & \mathrm{if}\frac{N}{averDist}<200\hfill \\ 5,\hfill & \mathrm{if}200\le \frac{N}{averDist}<240\hfill \\ 6,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(9)

where $averDist$ is defined as:

$$averDist=\frac{sum\left(minDist\right)}{N}$$

(10)

and where $minDist$ is an array, each element of which contains the distance from one point to its closest point in the set.

While the valley-seeking clustering method is highly applicable in many cases, it is not always able to successfully follow the dynamics of the trajectory. The problem it could encounter is a sudden change in the direction while forming the point order, as exemplified in Figure 5. To address this issue, we employed the classic valley-seeking method exclusively during the cluster initialization phase, specifically when determining the first two points of the cluster. This initial step defines the direction in which the cluster will be constructed, i.e., in which the next point will be searched for.

4.2. Trajectory Renewing

The objective of this part of the algorithm is forming an oriented graph that will, starting from the two points determined during the cluster initialization, result in a sequence of points. Once connected, the points from this sequence will form a line. The authors found that the way to accomplish this without the risk of diverging, as described in the previous section, is by combining the valley-seeking clustering method with an adaptive Kalman filter [34,35,36,37,38,39].

If we associate point-by-point contour reconstruction with target tracking, as in the example in Figure 6, the state vector ${\mathbf{x}}_{k+1}$ at time step $k+1$ can be expressed as:

$${\mathbf{x}}_{k+1}={\left[\begin{array}{c}{x}_{k+1}{\dot{x}}_{k+1}{y}_{k+1}{\dot{y}}_{k+1}\end{array}\right]}^T,$$

(11)

where $x_{k+1}$ and $y_{k+1}$ are the position coordinates, while ${\dot{x}}_{k+1}$ and ${\dot{y}}_{k+1}$ are the speeds along those two coordinates at time step $k+1$. The Kalman filter has the task of estimating the future states, i.e., the point positions, based on the available measurements.

If we adopt the constant velocity model, i.e., we assume that the velocity is constant between two consecutive sampling steps, we can write the following:

$${\mathbf{x}}_{k+1}=\mathbf{F}{\mathbf{x}}_k+{\mathbf{w}}_{k+1},$$

(12)

$${\mathbf{z}}_{k+1}=\mathbf{H}{\mathbf{x}}_{k+1}+{\mathbf{v}}_{k+1},$$

(13)

where $\mathbf{F}$ is the state transition matrix, ${\mathbf{z}}_{k+1}$ is the measurement at time step $k+1$, $\mathbf{H}$ is the measurement (i.e., observation) matrix, while ${\mathbf{w}}_{k+1}$ and ${\mathbf{v}}_{k+1}$ are the process noise and the measurement noise, respectively. The transition matrix is given by:

$$\mathbf{F}=\left[\begin{array}{cccc}1& △t& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& △t\\ 0& 0& 0& 1\end{array}\right]$$

(14)

where, for simplicity reasons, it was assumed that $△t=1$, and where the measurement matrix is:

$$\mathbf{H}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right].$$

(15)

With the nature of the particular points in mind, we can assume that the measurement noise ${\mathbf{v}}_{k+1}$ is negligible. On the other hand, the process noise ${\mathbf{w}}_{k+1}$ is very important, since its fluctuations indicate changes in trajectory, i.e., the presence of a maneuver.

Prediction can be defined as:

$${\widehat{\mathbf{x}}}_{k+1/k}=\mathbf{F}{\widehat{\mathbf{x}}}_{k/k},$$

(16)

$${\mathbf{P}}_{k+1/k}=\mathbf{F}{\mathbf{P}}_{k/k}{\mathbf{F}}^T+{\mathbf{Q}}_{k+1},$$

(17)

while updating can be expressed through the following set of equations:

$${\widehat{\mathbf{x}}}_{k+1/k+1}={\widehat{\mathbf{x}}}_{k+1/k}+{\mathbf{K}}_{k+1}({\mathbf{z}}_{k+1}-\mathbf{H}{\widehat{\mathbf{x}}}_{k+1/k}),$$

(18)

$${\mathbf{P}}_{k+1/k+1}=(\mathbf{I}-{\mathbf{K}}_{k+1}\mathbf{H}){\mathbf{P}}_{k+1/k},$$

(19)

$${\mathbf{K}}_{k+1}={\mathbf{P}}_{k+1/k}{\mathbf{H}}^T{(\mathbf{H}{\mathbf{P}}_{k+1/k}{\mathbf{H}}^T+\mathbf{R})}^{-1},$$

(20)

where ${\widehat{\mathbf{x}}}_{k+1/k}$ denotes the estimate for the following time step $k+1$, ${\mathbf{P}}_{k+1/k}$ and ${\mathbf{P}}_{k+1/k+1}$ denote the covariance matrices of the state estimations ${\widehat{\mathbf{x}}}_{k+1/k}$ and ${\widehat{\mathbf{x}}}_{k+1/k+1}$, respectively, ${\mathbf{Q}}_{k+1}$ denotes the covariance matrix of the process noise, $\mathbf{R}$ denotes the covariance matrix of the measurement noise, while ${\mathbf{K}}_{k+1}$ is the Kalman gain. Assuming that the individual components of the process noise and the measurement noise are uncorrelated, we can write that:

$$\mathbf{R}=\left[\begin{array}{cc}{\sigma }_x^2& 0\\ 0& {\sigma }_y^2\end{array}\right],$$

(21)

$${\mathbf{Q}}_{k+1}=\left[\begin{array}{cccc}\frac{q^3}{3}& \frac{q^2}{2}& 0& 0\\ \frac{q^2}{2}& q& 0& 0\\ 0& 0& \frac{q^3}{3}& \frac{q^2}{2}\\ 0& 0& \frac{q^2}{2}& q\end{array}\right],$$

(22)

where ${\sigma }_x^2$ and ${\sigma }_y^2$ are the variances of the x and y coordinate measurement noises, respectively, while the parameter q must be chosen in advance. The parameter q represents a deviation from the assumed constant velocity model. Since it is generally unknown, here, we adopted its value as $q=1$, while having in mind that the adaptation described in the following paragraphs will be able to indirectly model its unfixed value.

In order to detect the deviation from a linear trajectory, it is necessary to adapt the Kalman filter. We used the adaptation approach through the fading factor $S_{k+1}$ [36,37,38,39]. This implies changing the covariance matrix of the state estimation Equation (17):

$${\mathbf{P}}_{k+1/k}=S_{k+1}\mathbf{F}{\mathbf{P}}_{k/k}{\mathbf{F}}^T+{\mathbf{Q}}_{k+1}.$$

(23)

The fading factor $S_{k+1}$ represents the measure of model nonstationarity at time step $k+1$. It is based on the deviation of the theoretic covariance matrix (${\mathbf{M}}_{k+1}$) from the covariance matrix estimated on a sample (${\mathbf{N}}_{k+1}$) for the transformed prediction error sequence $\mathbf{H}\mathbf{F}\{{\mathbf{x}}_{k+1}-{\widehat{\mathbf{x}}}_{k+1/k}\}$:

$$S_{k+1}=max\{1,tr\left({\mathbf{N}}_{k+1}\right)/tr\left({\mathbf{M}}_{k+1}\right)\},$$

(24)

and where:

$${\mathbf{M}}_{k+1}=\mathbf{H}\mathbf{F}{\mathbf{P}}_{k/k}{\mathbf{F}}^T{\mathbf{H}}^T,$$

(25)

$${\mathbf{N}}_{k+1}={\mathbf{P}}_{V_{k+1}}-\mathbf{H}{\mathbf{Q}}_{k+1}{\mathbf{H}}^T-\mathbf{R},$$

(26)

where $P_{V_{k+1}}$ is the covariance matrix of the predicted residual vector. The predicted residual vector and its covariance matrix are given by:

$${\mathbf{V}}_{k+1}={\mathbf{z}}_{k+1}-\mathbf{H}{\mathbf{x}}_{k+1/k},$$

(27)

$${\mathbf{P}}_{V_{k+1}}=E\left(V_{k+1}{V}_{k+1}^T\right)\cong \frac{1}{k}\sum _{i=1}^k{\mathbf{V}}_{k+1}{\mathbf{V}}_{k+1}^T.$$

(28)

In an ideal case, when the trajectory is linear and the distances between the consecutive points are equal, ${\mathbf{M}}_k$ and ${\mathbf{N}}_k$ share the same value. In this case, the Kalman filter works in nominal mode ($S_k=1$). The adaptability of the Kalman filter is displayed through the situation when the measurements start to significantly differ from the predicted ones. This situation has a direct effect on the value of ${\mathbf{N}}_k$, which then as a result increases the value of the fading factor ($S_k>1$) and of the covariance matrix of the state estimation. As a result, the algorithm depends less on the previous and more on the new measurements.

After each state vector estimation ${\widehat{\mathbf{x}}}_{k+1/k}$, an area $G\left(X\right)$ is formed around each newly found point, i.e., graph node, according to Relations (1), (2), and (8). For each unassigned point from $G\left(X\right)$, a measure of accordance between the points and the trajectory is set. It is calculated in the manner of (7):

$$S_{xy}=cos{\alpha }_{xy},$$

(29)

where ${\alpha }_{xy}$ is the angle between the vector $(Y-X)$ and the speed vector ${\dot{x}}_{k+1}{\mathbf{i}}_x+{\dot{y}}_{k+1}{\mathbf{i}}_y$, where ${\mathbf{i}}_x$ and ${\mathbf{i}}_y$ are orthogonal unit vectors along the x- and y-axis, respectively. The point that maximizes the measure $S_{xy}$ is chosen as the successor of the point X.

4.3. Adding Omitted Points

The previously described procedure for arranging points into a cluster implies using the angle between the vector ($Y-X$) and the vector $∥\begin{array}{c}▽f\left(X\right)\end{array}∥$ during the initialization or the speed vector during the trajectory-renewing process. However, this leaves the possibility of missing some points that should be placed in the oriented graph between the nodes X and Y. The solution to this problem can be found in the implementation of the Nearest Neighbor (NN) algorithm after each step of the previously described procedure. In order to detect the omitted points and organize them correctly, we define a sector of a circle with its center in the point X:

$$\Xi \left(X\right)=\{Z:∥\begin{array}{c}X-Z\end{array}∥\le r_1,\angle (Y-X,Z-X)<\alpha \},$$

(30)

where $r_1=∥\begin{array}{c}X-Y\end{array}∥$ is the radius and the angle $\alpha$ is one half of the sector’s central angle (Figure 7). The angle $\alpha$ must be small enough so that the points that should belong to the other clusters are not encapsulated and added to the current cluster/graph. On the other hand, $\alpha$ must also be large enough to prevent any of the points that should belong to the current cluster being left outside. All the detected points should be placed between the points X and Y and sorted by ascending distance from the point X.

4.4. Finalizing the Individual Clusters and Forming the Final Contours

The parameter $S_{xy}$ from the formula (29), or the formula (7), is used as a criterion for finalizing the current cluster/graph. $S_{xy}$ with a positive value indicates that the successor Y of the point X was found. A negative value of $S_{xy}$ suggests that the found node does not have its successor, i.e., that the algorithm has reached the root of the graph. The situation where $S_{xy}$ is zero represents a borderline case. In terms of the probability density function, this case would correspond to a “flat” local maximum, where more than one potential maximum value can be found.

At the end of each iteration, which involves the analysis and incorporation of point(s) into the current cluster, the algorithm checks for any unsorted points remaining in the area $G\left(Y\right)$ surrounding the cluster’s last point (Y). If such points are detected, the previously described procedure is repeated, focusing on the point Y (as illustrated in Figure 8, and the accompanying pseudocode, which can be found in Appendix A). If the control does not return/find any potential candidates in the area $G\left(Y\right)$, the current cluster is finalized. New clusters are formed as long as there are still some unsorted points left.

Once all the clusters are formed, the algorithm checks which ones should be connected, in order to get the final nerve contours. For that purpose, around each cluster’s terminal node X and its predecessor $X_p$, an area in the form of a circle sector is made (as shown in Figure 9):

$$\Xi \left(X\right)=\{Z:∥\begin{array}{c}X-Z\end{array}∥\le r_1,\angle (X-X_p,Z-X)<\theta \},$$

(31)

where $r_1$ and $\theta$ are empirically derived values of the sector’s radius and of one half of the sector’s central angle, respectively. When some other cluster’s terminal node Y belongs to $\Xi \left(X\right)$, that cluster is considered to be an extension/continuation of the current cluster.

Remarks: Two evaluations used before the final contours are formed will be reported in the following paragraphs.

The initial evaluation indirectly assesses the dimensions of the selected area $G\left(X\right)$. If, in post-cluster formation, there is a pronounced number of extremely small clusters (consisting of three points or fewer, determined empirically), as in example in Figure 10, the algorithm deduces that the chosen $G\left(X\right)$ area was undersized. As a result, the entire process restarts with an expanded radius for $G\left(X\right)$.

In the second evaluation, clusters containing fewer than n points, where n was empirically set to 1 or 2 points, are considered for elimination. These smaller clusters, often located close to larger ones, typically do not significantly influence the final reconstruction, but can pose challenges during the final contour-formation stage.

5. Results and Discussion

The combination of methods and procedures reported in the previous section enables the reconstruction of one (Figure 9) or several contours (Figure 11).

A complete reconstruction of one complex (i.e., that has at least one branching along its length) nerve fascicle is shown in Figure 12. Fascicle shapes, dimensions, and branching locations differ across nerves and individuals. The example from Figure 12 shows a fascicle branching from a single section (in the lower part) into two sections (in the upper part of the figure). Indeed, nerve structures such as this one can present the biggest challenge to the algorithm, due to the complexity of nerve cross-sections at the branching points. This intricacy can lead to semi-connected contours with pronounced angles or contours in near proximity to one another, etc.

Figure 13 displays a segment (due to the clarity of the representation) of the reconstructed fascicles from one of the two subjects.

In the context of comparing our algorithm with existing state-of-the-art methods for contour reconstruction from 2D point sets, we primarily focused on analyzing two algorithms: HNN-Crust [28] and Vicur [24]. Our selection of these two was informed by a thorough review paper [11], where they particularly distinguished themselves from other methods in an experimental study. We also selected these two algorithms because, as highlighted in Section 2, they represent cutting-edge techniques within their respective algorithm categories. HNN-Crust is the latest iteration in the Crust algorithm lineage, introducing new enhancements to this widely used algorithm category. On the other hand, Vicur is an evolved version of the acclaimed Discur algorithm. Given that it was developed by authors previously involved with Discur, we opted to analyze this advanced iteration, Vicur.

While HNN-Crust and Vicur adeptly handle a broad range of reconstruction scenarios, they encountered challenges with specific nerve contour cases presented in this study. Particularly in situations requiring the differentiation of two closely situated contours, our newly introduced algorithm demonstrated a higher success rate than Vicur (Figure 14) and HNN-Crust (Figure 15).

A comparative analysis between our proposed method, HNN-Crust, and the Vicur algorithm is presented in

Table 1. Statistical comparison between methods *.

Algorithm	RMSE	% of Omitted Points	No. of Gaps per Fascicle
our method	0.020	1.94 (%)	0
HNN-Crust	0.043	0.4 (%)	9.5
Vicur	0.019	0 (%)	0.2

* Data for branching fascicles.

. Given that we observed minimal variation among the three methods except in instances of fascicle branching, we opted to specifically analyze and present statistics for fascicles exhibiting branching throughout their length. The results represent an average taken from all the cross-sections of such fascicles. This indicates that, while most of the analyzed contours are straightforward, a few within each fascicle are intricate, characterized by sharp corners and notably closely spaced contours. The first column displays the Root-Mean-Squared Error (RMSE) values of each of the methods in comparison to the expert curve. The second column indicates the percentage of points that were omitted during reconstruction, relative to the total number of points in the initial set. The third column presents the count of unintended gaps or interruptions in the reconstructed curve.

Table 1 reveals that, when considering the RMSE, Vicur and our method produce closely aligned results. However, in the context of missed points, Vicur excels over both our strategy and the HNN-Crust approach, seamlessly retaining every point from the input set. This is, in fact, the main liability of our method. Conversely, in terms of faithfully reconstructing continuous curves, both HNN-Crust and Vicur demonstrate inconsistencies, especially in the scenarios illustrated in Figure 14 and Figure 15.

Emphasizing the importance of the precise reconstruction of neighboring contours, particularly fascicle branches,

Table 2. Statistical details concerning two realistic nerve models used in this study.

Subject	Total No. of Fascicles	No. of Branchings	No. of Complex Fascicles (%)
Subject 1	49	20	≈33%
Subject 2	47	21	≈28%

showcases statistical data from realistic nerve models derived from [29] relevant to this study. The data reveal the count of fascicles recognized in the available nerve sections, as well as the overall number of fascicle branchings occurring along their paths. Notably, about a third of all detected fascicles in the examined nerve sections exhibited at least one branching, with some presenting multiple branches. It us noteworthy that a heightened resolution in the count of cross-sectional layers amplifies the chances of encountering such scenarios, as depicted in Figure 14 and Figure 15.

While the presented algorithm effectively performs reconstructions for its intended and similar applications, it is important to acknowledge its potential limitations:

• Given that our algorithm adopts the valley-seeking clustering method, the primary limitation of our approach lies in its inherent characteristic of omitting certain points deemed non-essential for reconstruction. Hence, it is not particularly tailored for reconstructions reliant on a small number of points. In such cases, there is an increased risk of overlooking or disregarding points that the algorithm determines as not aligning with the intended “line” (see Figure 16). However, this shortcoming diminishes with reconstructions based on a larger point set.
• Considering the sensitivity of the results to the selection of specific parameter values (k, defined in Equation (9), $\alpha$, from Equation (30), or $\theta$, from Equation (31)), it is important to reevaluate and adjust these parameter values when dealing with a significantly different dataset. Further details regarding the values of these parameters will be provided in the subsequent paragraphs.

The parameter k is defined in Equation (9). It is worth noting that, for future parameter selection, a value of 5 was suitable in roughly 93% of cases examined in this study. However, values of 4 and 6 reflect special cases and should be considered when setting up this method with initial values.

The parameter $\alpha$, as introduced in Equation (30), plays a vital role in defining the area where omitted points will be sought. Achieving an optimal value for $\alpha$ must be balanced due to the reasons and logic explained in Section 4.3. In our experimentation, we explored a range of values from 2${}^{\circ }$ to 90${}^{\circ }$, as these values align with the intended purpose; values below 2${}^{\circ }$ result in a region that is too narrow, while values exceeding 90${}^{\circ }$ yield regions that are excessively wide. It is noteworthy that, for nearly 97% of cases, any value within this range yields satisfactory reconstructions. However, in special cases, such as the one illustrated in the lower row of Figure 15, an angle between 5${}^{\circ }$ and 25${}^{\circ }$ proves necessary. In the context of nerve contour reconstruction, we determined that the optimal value for $\alpha$ is 8${}^{\circ }$.

The parameter $\theta$, as outlined in Equation (31), plays a pivotal role in consolidating clusters to formulate the final contour. When deliberating on the suitable value for this parameter, there was no predefined rationale regarding its permissible range. Unlike $\alpha$, which intuitively had a “forward” direction, $\theta$ possessed greater flexibility in its potential values. Consequently, we conducted experiments involving both the lower and upper limits of the angles that permit successful reconstruction. We systematically recorded the values that resulted in successful reconstructions, and these findings are summarized in

Table 3. Statistical analysis of minimum and maximum angle ranges for successful reconstruction with parameter $\theta$.

	Minimum Angles for Successful Reconstruction	Maximum Angles for Successful Reconstruction
Mean value	17.34${}^{\circ }$	252.5${}^{\circ }$
Standard deviation	13.55${}^{\circ }$	115.97${}^{\circ }$
Minimum	5${}^{\circ }$	110${}^{\circ }$
Maximum	45${}^{\circ }$	360${}^{\circ }$

, which presents statistical information regarding the range, i.e., the minimum and maximum angle values. As a practical rule of thumb, we identified the optimal value for $\theta$ to be 45${}^{\circ }$.

6. Conclusions

In this paper, we present an algorithm designed to reconstruct realistic nerve cross-sectional contours from unsorted points. Our approach addresses the complex geometry of human nerves’ internal structures, characterized by frequent fascicle crossings and branching, presenting the algorithm with challenges such as handling sharp corners and closely positioned individual contours that require precise differentiation.

This innovative method has significant potential contributions to the field of neuroprostheses, particularly in simulating natural sensations in amputated limbs. We recognize the ongoing advancements needed in this area. Additionally, this study built upon our previous work in modeling nerve geometry, as detailed in [29].

Our newly developed method offers several key advantages when compared to existing approaches for curve reconstruction from unarranged points:

• Handling nonuniform point sets: Our method excels in working with nonuniform sets of points and is capable of reconstructing one or more contours, even when they are not necessarily smooth.
• Reconstructing closely positioned contours: An essential advantage of our algorithm is its proficiency in reconstructing closely positioned contours, a task where many existing methods face challenges.
• Automation: Once the parameters are set for the input values and remain relatively consistent, our algorithm operates without the need for further manual input from the user.
• Dynamics and fascicle model: Unlike many existing methods that primarily focus on geometric considerations, our algorithm stands out by incorporating dynamics and implicitly integrating a fascicle model through an adaptive Kalman filter.
• Innovative combination: To our knowledge, there is currently no algorithm for curve reconstruction from randomized points that combines valley-seeking clustering with a Kalman filter and NN classification method.

A comparative analysis with state-of-the-art methods demonstrates that our approach excels in cases involving intricate internal nerve structures, making it particularly suitable for such scenarios. While it may not be ideal for cases with a small number of input points, it outperforms other methods in situations similar to those presented here.

In this paper, we primarily focused on reconstructing closed contours due to the nature of the contours examined. However, our method has the potential to reconstruct open lines with minor algorithm modifications. In future work, we plan on expanding the algorithm’s capabilities to handle open curves and, potentially, make the parameters more robust to the input values.

Beyond its primary purpose, this innovative method has broader potential applications and can serve as inspiration or a foundational technique in various fields. These encompass 2D and 3D image processing, target tracking (see Figure 6), and spatial analysis of linearly arranged points across diverse domains, including biomedicine. Its impact extends not only to the field of simulation planning in neuroprosthetics, but also to virtual surgeries and many other specialized domains.