TVM encompasses three stages: tensor representation, tensor communication, and tensor decomposition/encoding. ^{20 , 22} First, a tensor field was established by setting up a disk tensor for each LiDAR point, as proposed by You and Lin, ¹ and the geometric character of each LiDAR point was represented by a second symmetric tensor. Using tensor communication, each point collects the tensors of neighboring LiDAR points to form its final tensor. Consequently, each point has a strong geometric relationship with its neighboring points based on the final tensor, which is helpful for the extraction of ridge and valley lines in mountainous areas. After tensor communication, the final tensor of each point was decomposed into surface features, line features, and noise point features by tensor encoding. ²⁰ You and Lin ^{1 , 2} used a normal vector field to obtain the surface feature strengths using the TVM and applied the region-growing method to extract surface features from LiDAR point clouds for urban regions. For our purpose, the line feature strength (denoted as C here) is more suitable for extracting valley and river lines in mountainous regions than the other types of feature strength. This aspect is elaborated in the subsequent section.

C data can be obtained using the TVM. This section describes the characteristics of C data. Figure 1 shows the same mountainous region from various perspectives and data modes. In Figure 1(a), the ridge and valley regions appear in the 2.5D view, with an elevation range of approximately 1,400 meters. Figure 1(b) shows the elevation data from a 2D view. Figure 1(c) presents the C data from a 2D view.

To analyze the characteristics of the C data in relation to variations in elevation, this study established a section line. In Figure 2, the red line represents the elevation data and the blue line represents C data from the section lines shown in Figure 1(b) and (c). As depicted in Figure 2, C was higher in regions with large elevation gradients, and its magnitude was influenced by the extent of the elevation gradient. In areas with smaller elevation gradients, C was relatively stable. Hence, establishing a single threshold for identifying all points in the valley versus ridge regions is infeasible. In Figure 2, points near the ridge and valley regions display local maxima in the C data curves, confirming that C data can effectively assist in the identification of nodes of line features.

The eigenvectors for LiDAR point clouds can be derived using the TVM. Each point cloud possessed three associated eigenvectors, each with a distinct interpretation. For points located on a surface, the first eigenvector is the vector that is normal to the surface. Conversely, for points situated on a line, the third eigenvector is the tangent vector along the direction of the line, as shown by Medioni et al. ²⁰ and You and Lin. ^{1 , 2} Given that this study focuses on the extraction of ridge and valley lines, the tangent vector v for each point is considered. This vector is derived using the TVM and is crucial for eliminating extraneous point clouds.

Drawing inspiration from the region-growing method, this study established line growth criteria by choosing seed points and setting item thresholds. According to the next section, points located on the ridge or valley regions have higher C values than their surrounding points. Hence, points with local maxima C values serve as seeds for initiating line growth. Second, the v vectors of adjacent points located in the same valley or ridge line are similar. These two considerations guide the line-growth process, the specifics of which are explained in the following section.

This study establishes that such local maxima are directional in nature to provide a detailed explanation of how they are identified in 3D point cloud data. Vector v , derived from the TVM, is a critical factor. Points ( i = 1, 2, …, k) considered centers within a specified search radius are divided into right ( a = 1, …, m) and left ( b = 1, …, n) subsets based on v , which serves as a divide. A center point cloud is deemed to be a local maxima in terms of C when the average C values on the two sides are lower than the C value of the center point. This is formalized in the following equations: If ( ∑ a = 1 m ( C a ) m − C i ) < 0 and ( ∑ b = 1 n ( C b ) n − C i ) < 0 (1)

Subscript i indicates point i located on a line. C a represents the C value in point cloud a on the right-hand side of v ( a = 1, …, m). C b represents the C value in point cloud b on the left-hand side of v ( b = 1, …, n).

As depicted in Figure 3, points within the search radius are divided into right (R1, R2, R3, and R4) and left (L1 and L2) subsets based on v of P1. Point P1 is classified as the local maxima of C when it satisfies Equation (1). This criterion was employed to eliminate redundant points in the line growth process. Points satisfying Eq. (1) serve as candidate seeds for line growth.

When employing a large search radius, numerous point clouds not located on the line were included in the calculations. Conversely, when the search radius is small, many point clouds near the line segment are excluded, leading to incomplete line segments, particularly in the case of low-density point clouds. In this study, the radius used to identify point clouds with local maxima was consistent with the radius employed in TVM applications.

A threshold of C value may be chosen to determine all candidate seeds whose C values are larger than that threshold, as in the conventional growing method. The result would lead to the extraction of the primary ridge and valley lines only, but might filter out other minor lines. Different from the conventional method, candidate seeds are identified by seeking local C maxima for the line-growing in this article. Using this line-growing method, more details of the ridge and valley lines can be obtained. Points with higher C values are more likely to be located on lines and are therefore prioritized for line growth. According to empirical observations, neighboring points on a line should have similar tangent vectors. The threshold angle between the two v vectors can then be defined.

To meet the criteria of the angle threshold, the five points with the largest C value are regarded as the center of the ball in the TV’s searching radius, and the tangent vectors of points that meet the local maxima C value within the range are the calculated angles between the center points. The average angle value of the five groups was used as the basic threshold (ω), and three times the average value (3ω) was used as the threshold value in this study. An excessively large threshold can introduce inaccuracies by leading to the inclusion of numerous point clouds that are not actually part of the line segment, skewing the final calculations. Conversely, a too-small threshold can lead to the omission of relevant clouds and fragmentation of lines; the procedure cuts off branch lines to multiple single lines.

To preserve the detailed line features, the proposed approach adopts a moving mask algorithm from computer vision and employs an overlapping mechanism. The size of the moving ball matches that of the search ball used in the TVM, and its moving distance equals the radius of the ball. Points that satisfy both the angle threshold and Equation (1) are considered, with their C values serving as weights because a higher C value implies a higher likelihood that the point is part of a line. The node of a line and its associated tangent vector were then calculated using the following equations: ( V i ) = ∑ ( C j ∙ V j ) ∑ C j when i = 1 … . n , j = 1 … . m (2) ( x y z ) i = ∑ ( C j ∙ ( x y z ) j ) ∑ C j when i = 1 … . n , j = 1 … . m (3)

i: node i of a line. j: point j within the search ball.

n the number of nodes in a line. m: number of points within the search ball.

( x y z ) i T and tangent vector ( V i ) . The subsequent node of the line was determined using the following equation: ( x y z ) i + 1 = ( x y z ) i + Radius ∙ ( V i | V i | ) (4)

In line with the techniques used in computer vision for line extraction, overlapping regions were maintained between the moving masks. In this study, approximately half of the search balls were designed to overlap adjacent balls.

P0 [ Figure 4(a)] was selected as the initial seed. Points that satisfy both Equation (1) and the angle threshold within searching ball O1 are evaluated using Equations (2) and (3), which yield new values for P1 and v1, respectively. Search ball O2 is generated through Equation (4), and the procedure is repeated to obtain P2, v2, and so on until no points within the search ball meet the criteria. At this juncture, P1, P2, and P3 become the nodes of the line, as illustrated in Figure 4(b), and new seeds are selected for the growth of additional lines until all the points satisfying Equation (1) are used.

Node-point coordinates were calculated by weighting them in accordance with their C values. Equation (4) describes the method used to determine the node coordinates. To ascertain the theoretical accuracy, error propagation calculations were conducted for each point i, as outlined in Equation (5). For simplicity, any correlations among x _i, y _i, and z _i were disregarded. σ x 2 = ∑ ( σ x i ( ∑ C i ) ) 2 σ y 2 = ∑ ( σ y i ( ∑ C i ) ) 2 σ z 2 = ∑ ( σ z i ( ∑ C i ) ) 2 (5)

Another quality indicator of a node point is interpolated accuracy using contour data obtained by Delauney triangulation. The formula is as follows: z contour i = ∑ ( d i ∙ d i ) n (6)

i: node i of a line. n the number of nodes in a line.

d i : elevation from point D to plane ABC in Figure 5.

Contour data were used to produce a Delauney triangulation network. The nodes of the line are interpolated to obtain the difference in the elevation data and calculate the interpolated accuracy for each line. Once the theoretical and interpolated accuracies of each node of a line segment have been determined, the error propagation technique can be used to obtain the error band of the segment ( Figure 6). Here, x and y are treated as independent uncorrelated random variables, as indicated by the following formula: σ x ′ 2 = σ x 2 . cos 2 θ + σ y 2 . sin 2 θ               σ y ′ 2 = σ x 2 . sin 2 θ + σ y 2 . cos 2 θ (7)

When ( x A , y A ) and ( x B , y B ) follow a normal distribution, the transformed coordinates ( x, y) and ( x’, y’) also have a normal distribution. Given a significance level of α, the confidence interval for the line on the y’-axis can be calculated as ( − Z α 2 . σ y ′   Z α 2 . σ y ′ ). Plotting and connecting all points within this confidence level yields Figure 6, which is similar to the GIS result. ²³

Ethical approval and consent were not required.