# FDTD Complex Terrain Modeling Method Based on Papery Contour Map

Authors: Meng Chen; Bing Wei; Xinbo He
Journal: International Journal of Antennas and Propagation
 (2023)

---

## Abstract

This study aims to develop a finite difference time-domain (FDTD) modeling method for complex terrain using a paper contour map. First, the image preprocessing technology is employed to segment, denoise, and number the drawing area to arrive at independent figures and curve matrices. Subsequently, the support vector machine (SVM) classifier is implemented to identify the number and obtain the elevation value pertinent to the curve in the figure. Finally, the scattered points in the curve matrix are meshed by the Delaunay triangulation algorithm, and a three-dimensional triangular mesh model is established according to the recognition results. The FDTD calculation model is established by the method of intersection of Yee cell grid lines and triangle projection. The simulation results reveal that the recognition rate of the SVM classification model in this study is close to 100% for printed numerals. This approach is capable of quickly and accurately establishing the FDTD-based calculation model for complex terrain on the basis of the paper contour map.


---

## Body

## 1. Introduction

The establishment of a complex terrain calculation model is of substantial significance in the fields of radio wave propagation research and electromagnetic compatibility analysis. Currently, the common modeling scheme is based on the global geographic information system through the corresponding electronic information data and combined with various subdivision algorithms. This approach has a good effect in dealing with electronic maps; however, it is an undeniable fact that in the early days, a large number of paper contour maps of terrain are achieved through relatively primitive surveying and mapping methods. If the method of image recognition can be applied to combine contour maps with electromagnetic computational modeling, the modeling efficiency could be noticeably enhanced. Therefore, based on different image processing algorithms and support vector machine (SVM) classification models, the aforementioned problems are effectively unlocked by processing the paper contour maps.Among numerous numerical calculation methodologies, the finite difference time-domain (FDTD) approach is the most broadly utilized in the accurate analysis of radio wave propagation over a complex terrain. This method has natural parallelism and high computational accuracy [1–4]. The modeling of the FDTD algorithm chiefly includes the type value point method, analytical method, and projection intersection method [5–7]. The first two methods can only be exploited for simple object modeling, which is challenging to be applied to practical engineering, while the projection intersection method can be employed for various complex models.The study contour map is essentially composed of curves and figures, which should be extracted by artificial intelligence. The SVM classification-based model is based on the structural risk minimization principle and Vapnik–Chervonenkis theory. Compared with other classification-based models, this model is better at cracking the problems of a small number of samples, nonlinear separability, and high-dimensional pattern recognition. It also employs effective solutions to resolve the problem of “dimensionality disaster” and overfitting [8, 9]. This algorithm has been extensively utilized in signal feature recognition [10, 11], image detection [12–14], freight volume forecast [15], and other fields. In the current work, the SVM classifier is exploited for digital recognition to obtain the elevation value and image scatter matrix corresponding to the curve in the graph. On this basis, the Bowyer–Watson algorithm based on the Delaunay partition criterion is adopted to divide the image into triangles, and a three-dimensional (3D) triangle model of a complex terrain is established through elevation projection. On this basis, the FDTD calculation model is developed by the method of intersection of the Yee cell grid line and triangle projection.Herein, the FDTD calculation model of a complex terrain is developed on the basis of the paper contour map, and the electromagnetic coupling effect of the mountain model in the presence of high-altitude nuclear electromagnetic pulse (HEMP) incidence is calculated. The overall processing steps of a papery contour map are illustrated in Figure1. The simulation results reveal that the proposed SVM-based classification model has a recognition rate of nearly 100% for printed numerals. This approach could rapidly and precisely establish the FDTD calculation model of complex topographies on the basis of the paper contour map [16–18].Figure 1
The overall processing steps of the papery contour map.


## 2. Modeling Data Acquisition Based on a Contour Map

### 2.1. Image Preprocessing

The digital recognition process of the paper contour essentially includes image preprocessing and digital recognition. Image preprocessing is the foreshadowing work for digital recognition and subsequent mesh generation, which is implemented to deal with the problems of image noise, shadow, and curve scatter density.


#### 2.1.1. Connected Domain Segmentation

The connected domain refers to the set of pixels with the same value and adjacent positions in the pixel matrix, while connected domain segmentation is the process of marking and extracting each connected domain. During processing contour drawings, each curve and the corresponding numbers belong to different connected domains, as presented in Figure2. For the drawings with correct format and clear printing, the connected domain segmentation technology could be employed to distinguish the information such as numbers, symbols, and curves, as presented in Figure 3.Figure 2
A papery contour map.Figure 3
Demonstration of the connected domain segmentation technology: (a) curves and (b) numbers.
(a)(b)In addition to segmenting multiply connected domains in the image, it is also necessary to create an index between each curve and its corresponding number. Otherwise, the elevation associated with each curve cannot be achieved during the subsequent 3D modeling. Herein, the Euclidean distance between the matrices and the plane coordinate range of each connected domain is evaluated to assess whether an index should be created between two connected domains.


#### 2.1.2. Morphological Processing

For the segmented binarized pixel matrix, morphological processing should be applied for identification or meshing. The processing methodologies exploited in the present work include connected domain segmentation, erosion, dilation, and refinement. Erosion is to move the primitives along the boundary of the image, find the intersection between the moved area and the original image, and eliminate the intersection from the original image to obtain a new image. Dilation is to move the base element along the boundary of the image and obtain a new image from the union of the moving region and the original image. Graphical burrs can be removed through erosion and dilation, enhancing recognition accuracy and triangulation efficiency, as illustrated in Figure4.Figure 4
Demonstration of erosion and dilation: (a) original, (b) erosion, and (c) dilation.
(a)(b)(c)The thinning operation is utilized to extract the main “skeleton” of the image, which can convert the lines from multipixel width to single-pixel width. Using the thinning operation, the core frame of the contour line can be effectively extracted and the redundant pixels can be removed. In this way, in the subsequent mesh generation process, it can not only reduce the amount of calculation but also avoid generating too many invalid triangles, as illustrated in Figure5.Figure 5
Presentation of the thinning effect.


### 2.2. Digital Recognition

#### 2.2.1. SVM Classification Model

The SVM-based classification model can be implemented for analyzing both linear and nonlinear classification problems. For linear ones, the SVM is capable of maximizing the minimum distance between the sample points and the hyperplane on the basis of establishing the classification hyperplane.For the datasetX1,d1,X2,d2,⋯,Xi,di,⋯,XN,dN composed of N linearly separable samples, the i-th input sample Xi is expected to be output as di. If the hyperplane equation is set as WTX+b=0, where X represents the input vector, W is the weight vector, and b denotes the offset value, the classification problem could be transformed into finding the optimal solution of W and b. Calculate the distance from the sample point to the hyperplane, construct the Lagrange multiplier equation system, and combine the Karush–Kuhn–Tucker conditions to obtain the optimal weight vector:(1)W0=∑i=1Nλ0idiXi=∑Allsupportvectorsλ0idSXS,where λ0i represents the optimal solution of Qλ, XS denotes the input sample nearest to the optimal hyperplane, and dS is the expected output value associated with the support vector. Furthermore, if the optimal bias is given by b0=1−W0TXS, the optimal classification discriminant function for linear separable problems is provided as follows:(2)fX=sgn∑i=1Nλ0idiXiTX+b0,where sgnxrepresents a step function. When solving nonlinear problems, the input samples are mapped into a high-dimensional space to make them linearly separable in the high-dimensional space, and then, the optimal hyperplane is established in the high-dimensional space for classification based on the method exploited for solving linear separable problems [19, 20]. If ΦX=Φ1X,Φ2X,⋯,ΦMXT represents the nonlinear mapping of the input vector X from the original space to the M-dimensional feature space, the optimal classification discriminant function is expressed as follows:(3)fX=sgn∑i=1Nλ0idiΦTXiΦX+b0,where ΦTXiΦX denotes the inner product of the i-th input sample Xi and the input vector X mapped in the high-dimensional space. The operation of the vector inner product can be solved with the help of the kernel function. When a function satisfies the Mercer condition [21], the so-called kernel function, which is denoted by KXi,Xj, is exploited to receive input in the low-dimensional space and output in the high-dimensional space [22, 23]. The following linear kernel function is employed throughout all numerical studies in the present work: KXi,X=Xi⋅X.For the multiclassification problem of digital recognition, the error-correcting output codes are exploited to transform the multiclassification problem into multiple binary classification problems. Assuming that the number of binary classification problems isN, combining these N problems in pairs configures N (N−1)/2 classification tasks, in which one category in each classification task is taken as a positive example and the other as a negative example. N (N−1)/2 classifiers are attained through training. The samples are submitted to all classifiers at the same time to obtain N (N−1)/2 classification results, and then, the category with the most prediction is taken as the final result.


#### 2.2.2. Analysis of Digit Recognition Results

In this study, 800 printed digital images are chosen as training samples, and the image size is set as 106 × 106 pixels. The SVM-based classification model could be achieved by training these 800 images. The digitally connected domain matrix extracted from the contour map is input into the model, as presented in Table1. The input sample is demonstrated in Figure 6.Table 1
Digital recognition results of the contour map.
Input samplesSize (px)Samples numberRecognition numberRecognition rate (%)Test sample16 × 16100100100Image sample 1106 × 1062727100Image sample 2106 × 1061010100Figure 6
The input samples: (a) test sample and (b) image sample 1.
(a)(b)Table1 displays that the printed number recognition effect using the SVM-based classification model is accurate. Therefore, the suggested algorithm would be suitable for high-precision research work such as contour modeling.


### 2.3. Delaunay Triangulation Algorithm

According to the result of the connected domain extraction, the scatter set is triangulated. The Bowyer–Watson triangulation algorithm based on the Delaunay subdivision criterion is implemented to subdivide 1488 scatter points and finally arrive at 5908 Delaunay triangles. The triangles obtained by this algorithm not only meet the characteristics of minimum angle maximization and empty circle but also meet the established standard of the FDTD calculation model, as illustrated in Figure7.Figure 7
Display of Bowyer–Watson triangulation.


## 2.1. Image Preprocessing

The digital recognition process of the paper contour essentially includes image preprocessing and digital recognition. Image preprocessing is the foreshadowing work for digital recognition and subsequent mesh generation, which is implemented to deal with the problems of image noise, shadow, and curve scatter density.


### 2.1.1. Connected Domain Segmentation

The connected domain refers to the set of pixels with the same value and adjacent positions in the pixel matrix, while connected domain segmentation is the process of marking and extracting each connected domain. During processing contour drawings, each curve and the corresponding numbers belong to different connected domains, as presented in Figure2. For the drawings with correct format and clear printing, the connected domain segmentation technology could be employed to distinguish the information such as numbers, symbols, and curves, as presented in Figure 3.Figure 2
A papery contour map.Figure 3
Demonstration of the connected domain segmentation technology: (a) curves and (b) numbers.
(a)(b)In addition to segmenting multiply connected domains in the image, it is also necessary to create an index between each curve and its corresponding number. Otherwise, the elevation associated with each curve cannot be achieved during the subsequent 3D modeling. Herein, the Euclidean distance between the matrices and the plane coordinate range of each connected domain is evaluated to assess whether an index should be created between two connected domains.


### 2.1.2. Morphological Processing

For the segmented binarized pixel matrix, morphological processing should be applied for identification or meshing. The processing methodologies exploited in the present work include connected domain segmentation, erosion, dilation, and refinement. Erosion is to move the primitives along the boundary of the image, find the intersection between the moved area and the original image, and eliminate the intersection from the original image to obtain a new image. Dilation is to move the base element along the boundary of the image and obtain a new image from the union of the moving region and the original image. Graphical burrs can be removed through erosion and dilation, enhancing recognition accuracy and triangulation efficiency, as illustrated in Figure4.Figure 4
Demonstration of erosion and dilation: (a) original, (b) erosion, and (c) dilation.
(a)(b)(c)The thinning operation is utilized to extract the main “skeleton” of the image, which can convert the lines from multipixel width to single-pixel width. Using the thinning operation, the core frame of the contour line can be effectively extracted and the redundant pixels can be removed. In this way, in the subsequent mesh generation process, it can not only reduce the amount of calculation but also avoid generating too many invalid triangles, as illustrated in Figure5.Figure 5
Presentation of the thinning effect.


## 2.1.1. Connected Domain Segmentation

The connected domain refers to the set of pixels with the same value and adjacent positions in the pixel matrix, while connected domain segmentation is the process of marking and extracting each connected domain. During processing contour drawings, each curve and the corresponding numbers belong to different connected domains, as presented in Figure2. For the drawings with correct format and clear printing, the connected domain segmentation technology could be employed to distinguish the information such as numbers, symbols, and curves, as presented in Figure 3.Figure 2
A papery contour map.Figure 3
Demonstration of the connected domain segmentation technology: (a) curves and (b) numbers.
(a)(b)In addition to segmenting multiply connected domains in the image, it is also necessary to create an index between each curve and its corresponding number. Otherwise, the elevation associated with each curve cannot be achieved during the subsequent 3D modeling. Herein, the Euclidean distance between the matrices and the plane coordinate range of each connected domain is evaluated to assess whether an index should be created between two connected domains.


## 2.1.2. Morphological Processing

For the segmented binarized pixel matrix, morphological processing should be applied for identification or meshing. The processing methodologies exploited in the present work include connected domain segmentation, erosion, dilation, and refinement. Erosion is to move the primitives along the boundary of the image, find the intersection between the moved area and the original image, and eliminate the intersection from the original image to obtain a new image. Dilation is to move the base element along the boundary of the image and obtain a new image from the union of the moving region and the original image. Graphical burrs can be removed through erosion and dilation, enhancing recognition accuracy and triangulation efficiency, as illustrated in Figure4.Figure 4
Demonstration of erosion and dilation: (a) original, (b) erosion, and (c) dilation.
(a)(b)(c)The thinning operation is utilized to extract the main “skeleton” of the image, which can convert the lines from multipixel width to single-pixel width. Using the thinning operation, the core frame of the contour line can be effectively extracted and the redundant pixels can be removed. In this way, in the subsequent mesh generation process, it can not only reduce the amount of calculation but also avoid generating too many invalid triangles, as illustrated in Figure5.Figure 5
Presentation of the thinning effect.


## 2.2. Digital Recognition

### 2.2.1. SVM Classification Model

The SVM-based classification model can be implemented for analyzing both linear and nonlinear classification problems. For linear ones, the SVM is capable of maximizing the minimum distance between the sample points and the hyperplane on the basis of establishing the classification hyperplane.For the datasetX1,d1,X2,d2,⋯,Xi,di,⋯,XN,dN composed of N linearly separable samples, the i-th input sample Xi is expected to be output as di. If the hyperplane equation is set as WTX+b=0, where X represents the input vector, W is the weight vector, and b denotes the offset value, the classification problem could be transformed into finding the optimal solution of W and b. Calculate the distance from the sample point to the hyperplane, construct the Lagrange multiplier equation system, and combine the Karush–Kuhn–Tucker conditions to obtain the optimal weight vector:(1)W0=∑i=1Nλ0idiXi=∑Allsupportvectorsλ0idSXS,where λ0i represents the optimal solution of Qλ, XS denotes the input sample nearest to the optimal hyperplane, and dS is the expected output value associated with the support vector. Furthermore, if the optimal bias is given by b0=1−W0TXS, the optimal classification discriminant function for linear separable problems is provided as follows:(2)fX=sgn∑i=1Nλ0idiXiTX+b0,where sgnxrepresents a step function. When solving nonlinear problems, the input samples are mapped into a high-dimensional space to make them linearly separable in the high-dimensional space, and then, the optimal hyperplane is established in the high-dimensional space for classification based on the method exploited for solving linear separable problems [19, 20]. If ΦX=Φ1X,Φ2X,⋯,ΦMXT represents the nonlinear mapping of the input vector X from the original space to the M-dimensional feature space, the optimal classification discriminant function is expressed as follows:(3)fX=sgn∑i=1Nλ0idiΦTXiΦX+b0,where ΦTXiΦX denotes the inner product of the i-th input sample Xi and the input vector X mapped in the high-dimensional space. The operation of the vector inner product can be solved with the help of the kernel function. When a function satisfies the Mercer condition [21], the so-called kernel function, which is denoted by KXi,Xj, is exploited to receive input in the low-dimensional space and output in the high-dimensional space [22, 23]. The following linear kernel function is employed throughout all numerical studies in the present work: KXi,X=Xi⋅X.For the multiclassification problem of digital recognition, the error-correcting output codes are exploited to transform the multiclassification problem into multiple binary classification problems. Assuming that the number of binary classification problems isN, combining these N problems in pairs configures N (N−1)/2 classification tasks, in which one category in each classification task is taken as a positive example and the other as a negative example. N (N−1)/2 classifiers are attained through training. The samples are submitted to all classifiers at the same time to obtain N (N−1)/2 classification results, and then, the category with the most prediction is taken as the final result.


### 2.2.2. Analysis of Digit Recognition Results

In this study, 800 printed digital images are chosen as training samples, and the image size is set as 106 × 106 pixels. The SVM-based classification model could be achieved by training these 800 images. The digitally connected domain matrix extracted from the contour map is input into the model, as presented in Table1. The input sample is demonstrated in Figure 6.Table 1
Digital recognition results of the contour map.
Input samplesSize (px)Samples numberRecognition numberRecognition rate (%)Test sample16 × 16100100100Image sample 1106 × 1062727100Image sample 2106 × 1061010100Figure 6
The input samples: (a) test sample and (b) image sample 1.
(a)(b)Table1 displays that the printed number recognition effect using the SVM-based classification model is accurate. Therefore, the suggested algorithm would be suitable for high-precision research work such as contour modeling.


## 2.2.1. SVM Classification Model

The SVM-based classification model can be implemented for analyzing both linear and nonlinear classification problems. For linear ones, the SVM is capable of maximizing the minimum distance between the sample points and the hyperplane on the basis of establishing the classification hyperplane.For the datasetX1,d1,X2,d2,⋯,Xi,di,⋯,XN,dN composed of N linearly separable samples, the i-th input sample Xi is expected to be output as di. If the hyperplane equation is set as WTX+b=0, where X represents the input vector, W is the weight vector, and b denotes the offset value, the classification problem could be transformed into finding the optimal solution of W and b. Calculate the distance from the sample point to the hyperplane, construct the Lagrange multiplier equation system, and combine the Karush–Kuhn–Tucker conditions to obtain the optimal weight vector:(1)W0=∑i=1Nλ0idiXi=∑Allsupportvectorsλ0idSXS,where λ0i represents the optimal solution of Qλ, XS denotes the input sample nearest to the optimal hyperplane, and dS is the expected output value associated with the support vector. Furthermore, if the optimal bias is given by b0=1−W0TXS, the optimal classification discriminant function for linear separable problems is provided as follows:(2)fX=sgn∑i=1Nλ0idiXiTX+b0,where sgnxrepresents a step function. When solving nonlinear problems, the input samples are mapped into a high-dimensional space to make them linearly separable in the high-dimensional space, and then, the optimal hyperplane is established in the high-dimensional space for classification based on the method exploited for solving linear separable problems [19, 20]. If ΦX=Φ1X,Φ2X,⋯,ΦMXT represents the nonlinear mapping of the input vector X from the original space to the M-dimensional feature space, the optimal classification discriminant function is expressed as follows:(3)fX=sgn∑i=1Nλ0idiΦTXiΦX+b0,where ΦTXiΦX denotes the inner product of the i-th input sample Xi and the input vector X mapped in the high-dimensional space. The operation of the vector inner product can be solved with the help of the kernel function. When a function satisfies the Mercer condition [21], the so-called kernel function, which is denoted by KXi,Xj, is exploited to receive input in the low-dimensional space and output in the high-dimensional space [22, 23]. The following linear kernel function is employed throughout all numerical studies in the present work: KXi,X=Xi⋅X.For the multiclassification problem of digital recognition, the error-correcting output codes are exploited to transform the multiclassification problem into multiple binary classification problems. Assuming that the number of binary classification problems isN, combining these N problems in pairs configures N (N−1)/2 classification tasks, in which one category in each classification task is taken as a positive example and the other as a negative example. N (N−1)/2 classifiers are attained through training. The samples are submitted to all classifiers at the same time to obtain N (N−1)/2 classification results, and then, the category with the most prediction is taken as the final result.


## 2.2.2. Analysis of Digit Recognition Results

In this study, 800 printed digital images are chosen as training samples, and the image size is set as 106 × 106 pixels. The SVM-based classification model could be achieved by training these 800 images. The digitally connected domain matrix extracted from the contour map is input into the model, as presented in Table1. The input sample is demonstrated in Figure 6.Table 1
Digital recognition results of the contour map.
Input samplesSize (px)Samples numberRecognition numberRecognition rate (%)Test sample16 × 16100100100Image sample 1106 × 1062727100Image sample 2106 × 1061010100Figure 6
The input samples: (a) test sample and (b) image sample 1.
(a)(b)Table1 displays that the printed number recognition effect using the SVM-based classification model is accurate. Therefore, the suggested algorithm would be suitable for high-precision research work such as contour modeling.


## 2.3. Delaunay Triangulation Algorithm

According to the result of the connected domain extraction, the scatter set is triangulated. The Bowyer–Watson triangulation algorithm based on the Delaunay subdivision criterion is implemented to subdivide 1488 scatter points and finally arrive at 5908 Delaunay triangles. The triangles obtained by this algorithm not only meet the characteristics of minimum angle maximization and empty circle but also meet the established standard of the FDTD calculation model, as illustrated in Figure7.Figure 7
Display of Bowyer–Watson triangulation.


## 3. Modeling Data Acquisition Based on a Contour Map

### 3.1. Convert Triangular Facets to the FDTD-Based Meshes

To construct the FDTD-based mesh model of the target from the data of the triangular face model, the key is to judge whether the FDTD-based cell is on the target. The projective intersection approach between the triangular face elements of the target surface and the FDTD-based mesh lines is exploited to examine the problem with the following steps:(i)
Determine the maximum and minimum values of the target shape on thex-axes and y-axes and then scan in this range.(ii)
Find out which triangular surface elements of the target surface intersect each line parallel to thez-axis and get the intersection point with the target surface. Hence, the maximum and minimum values of the external contour of the object in this direction are obtained.(iii)
The identical technique is applied to the projection intersection of lines parallel to thex-axes and y-axes.


#### 3.1.1. Solve the Intersection Point between the Line and the Plane Where the Triangular Plane Element Is Located

The coordinates of the three vertices of the triangle surface element are specified byAx1,y1,z1, Bx2,y2,z2, and Cx3,y3,z3 (Figure 8). In the rectangular coordinates system, the equation of the line l parallel to the z-axis could be stated as follows:(4)x=i·dx,y=j·dy.Figure 8
The intersection of the trigonometric elementABC and FDTD grid lines.The equation of the plane, in which any triangular surface element is located on the target surface, is given as follows:(5)y1z11y2z21y3z31x+z1x11z2x21z3x31y+x1y11x2y21x3y31z=x1y1z1x2y2z2x3y3z3.Let us define the following coefficients:(6)A=y1z11y2z21y3z31,B=z1x11z2x21z3x31y,C=x1y11x2y21x3y31,D=x1y1z1x2y2z2x3y3z3,then, equation (5) could be rewritten in a more compact form as follows:(7)Ax+By+Cz=D.Simultaneous application of equations (4) and (5) leads to(8)x=i·dx,y=j·dy,Ax+By+Cz=D.Thus,(9)z=1CD−Ai·dx−Bj·dy.The intersection of the line defined by equation (4), which is parallel to the z-axis, and the plane introduced by equation (7), where the triangular plane element is located, is then derived in the following form:(10)x=i·dx,y=j·dy,z=1CD−Ai·dx−Bj·dy.


#### 3.1.2. Determine Whether the Intersection Point Is Located in the Triangular Plane Element

The following is to further determine whether the intersection point is placed in the triangular plane. The specific idea is to determine whether the intersection point of the line provided by equation (10) and the coordinate plane xoy is located in the projection of the triangular plane on the coordinate plane xoy. Assume that the projection of a triangular plane element in the coordinate plane xoy is triangle ABC, and the projection of the intersection point between the line and the triangular plane element in the coordinate plane is denoted by P, as demonstrated in Figure 9. Figure 9(a) clearly represents the case where the intersection point is placed inside the triangular plane element; on the other hand, the case where the intersection point is located outside this specific zone is illustrated in Figure 9(b).Figure 9
The intersection point P is placed inside and outside of the projectionABC of the triangular plane element: (a) scenario 1 and (b) scenario 2.
(a)(b)The specific algorithm to assess whether the intersection point is located inside the triangular plane element is explained as follows. Consider the intersection point P in the planexoy, where the line and triangular plane element are located; then, the three vertices of the triangle could be constructed on the basis of three vectors PA,PB,andPC (Figure 9). Let us evaluate the cross product of the three vectors in turn as follows: PA×PB,PB×PC,andPC×PA; if the directions of the resulting vectors after their cross product are the same or one or two of them are equal to 0, then the intersection point of the line in equation (10) and the triangular plane element is located inside the triangular plane element; otherwise, the direction of the produced vectors after their cross product would be different, that is, one of them is different from the other two. The intersection point between the above line introduced by equation (10) and the triangular plane element is not inside the triangular plane element. According to the performed analysis, it can be assumed that the coordinates of the three vertices of the projection triangle of the triangular plane element would be Ax1,y1,0, Bx2,y2,0, and Cx3,y3,0, and the intersection point is Px,y,0; then, the three vectors PA,PB, and PC can be represented as follows:(11)PA=x1−xi^+y1−yj^,PB=x2−xi^+y2−yj^,PC=x3−xi^+y3−yj^,PA×PB=k^x1−xy2−y−y1−yx2−x,PB×PC=k^x2−xy3−y−y2−yx3−x,PC×PA=k^x3−xy1−y−y3−yx1−x,in which(12)t1=x1−xy2−y−y1−yx2−x,t2=x2−xy3−y−y2−yx3−x,t3=x3−xy1−y−y3−yx1−x.Ift1≥0,t2≥0,t3≥0 or t1≤0,t2≤0,t3≤0, then the intersection point is placed in the triangle element. The same method could be applied to the projection intersection of lines parallel to the x-axes and y-axes.


#### 3.1.3. Assign Values

The intersection points are marked and the electromagnetic parameters are assigned, and then, the grid profile data file for the FDTD-based calculation is generated. Note that in this step, the cell whose intersection is on the target region is marked as 1; otherwise, its mark is set as 0.


### 3.2. Establishment of the FDTD Calculation Model

The digital recognition result and the triangle set generated by triangulation are extracted, and the three-dimensional model of the mountain peak is constructed according to the digital-curve index relationship developed during the connected domain segmentation [24, 25]. A closed surface filling is performed for nonclosed surfaces. Finally, the triangular patches are shaded by the geometry shader of the OpenGL graphics library, the wireframe is filled, and the lighting is set to make the model close to the real effect (Figure 10).Figure 10
Three-dimensional model of the mountain: (a) grid model and (b) entity.
(a)(b)The projection intersection approach is exploited to transform the derived triangular patches into the required FDTD computational grid model, as demonstrated in Figure11.Figure 11
Yee grid model of the mountain.


### 3.3. Application Examples

The established mountain is used as a calculation model to analyze the field distribution of the mountain at various time steps and the waveforms of monitoring points at different heights. The HEMP is often exploited to examine the electromagnetic scattering characteristics adjacent to the ground and the coupling effect of cables [26, 27]. The medium parameters of the mountain are εr=4, σ=0.01S/m, and μr=1, and the FDTD discrete cell size is set as 10/10 m. The incident angle is specified by θ=180o, α=0o, and φ=0o. The field distribution of xoz and yoz surfaces of the mountain at 300-time steps is calculated under the HEMP incidence condition, as shown in Figure 12.Figure 12
Electric field distribution map of the mountain section: (a)yoz cross-section and (b) xoz cross-section.
(a)(b)Let us select two monitoring points at different heights to observe the time-domain waveform of the monitoring points. These two points with their coordinates are represented as follows: spot 1 (5, 5, 100) and spot 2 (5, 5, 80). The time-domain waveforms of the monitoring points are depicted in Figure13. The demonstrated results in Figure 13 indicate that the HEMP pulse reaches spot 1 first, and the peak field strength of the monitoring point attenuates 23.38% of the incident wave. After a certain time, the pulse reaches spot 2, and the peak field intensity at the monitoring point attenuates to 8.47% of the incident wave.Figure 13
Time-domain waveform of the electric field.


## 3.1. Convert Triangular Facets to the FDTD-Based Meshes

To construct the FDTD-based mesh model of the target from the data of the triangular face model, the key is to judge whether the FDTD-based cell is on the target. The projective intersection approach between the triangular face elements of the target surface and the FDTD-based mesh lines is exploited to examine the problem with the following steps:(i)
Determine the maximum and minimum values of the target shape on thex-axes and y-axes and then scan in this range.(ii)
Find out which triangular surface elements of the target surface intersect each line parallel to thez-axis and get the intersection point with the target surface. Hence, the maximum and minimum values of the external contour of the object in this direction are obtained.(iii)
The identical technique is applied to the projection intersection of lines parallel to thex-axes and y-axes.


### 3.1.1. Solve the Intersection Point between the Line and the Plane Where the Triangular Plane Element Is Located

The coordinates of the three vertices of the triangle surface element are specified byAx1,y1,z1, Bx2,y2,z2, and Cx3,y3,z3 (Figure 8). In the rectangular coordinates system, the equation of the line l parallel to the z-axis could be stated as follows:(4)x=i·dx,y=j·dy.Figure 8
The intersection of the trigonometric elementABC and FDTD grid lines.The equation of the plane, in which any triangular surface element is located on the target surface, is given as follows:(5)y1z11y2z21y3z31x+z1x11z2x21z3x31y+x1y11x2y21x3y31z=x1y1z1x2y2z2x3y3z3.Let us define the following coefficients:(6)A=y1z11y2z21y3z31,B=z1x11z2x21z3x31y,C=x1y11x2y21x3y31,D=x1y1z1x2y2z2x3y3z3,then, equation (5) could be rewritten in a more compact form as follows:(7)Ax+By+Cz=D.Simultaneous application of equations (4) and (5) leads to(8)x=i·dx,y=j·dy,Ax+By+Cz=D.Thus,(9)z=1CD−Ai·dx−Bj·dy.The intersection of the line defined by equation (4), which is parallel to the z-axis, and the plane introduced by equation (7), where the triangular plane element is located, is then derived in the following form:(10)x=i·dx,y=j·dy,z=1CD−Ai·dx−Bj·dy.


### 3.1.2. Determine Whether the Intersection Point Is Located in the Triangular Plane Element

The following is to further determine whether the intersection point is placed in the triangular plane. The specific idea is to determine whether the intersection point of the line provided by equation (10) and the coordinate plane xoy is located in the projection of the triangular plane on the coordinate plane xoy. Assume that the projection of a triangular plane element in the coordinate plane xoy is triangle ABC, and the projection of the intersection point between the line and the triangular plane element in the coordinate plane is denoted by P, as demonstrated in Figure 9. Figure 9(a) clearly represents the case where the intersection point is placed inside the triangular plane element; on the other hand, the case where the intersection point is located outside this specific zone is illustrated in Figure 9(b).Figure 9
The intersection point P is placed inside and outside of the projectionABC of the triangular plane element: (a) scenario 1 and (b) scenario 2.
(a)(b)The specific algorithm to assess whether the intersection point is located inside the triangular plane element is explained as follows. Consider the intersection point P in the planexoy, where the line and triangular plane element are located; then, the three vertices of the triangle could be constructed on the basis of three vectors PA,PB,andPC (Figure 9). Let us evaluate the cross product of the three vectors in turn as follows: PA×PB,PB×PC,andPC×PA; if the directions of the resulting vectors after their cross product are the same or one or two of them are equal to 0, then the intersection point of the line in equation (10) and the triangular plane element is located inside the triangular plane element; otherwise, the direction of the produced vectors after their cross product would be different, that is, one of them is different from the other two. The intersection point between the above line introduced by equation (10) and the triangular plane element is not inside the triangular plane element. According to the performed analysis, it can be assumed that the coordinates of the three vertices of the projection triangle of the triangular plane element would be Ax1,y1,0, Bx2,y2,0, and Cx3,y3,0, and the intersection point is Px,y,0; then, the three vectors PA,PB, and PC can be represented as follows:(11)PA=x1−xi^+y1−yj^,PB=x2−xi^+y2−yj^,PC=x3−xi^+y3−yj^,PA×PB=k^x1−xy2−y−y1−yx2−x,PB×PC=k^x2−xy3−y−y2−yx3−x,PC×PA=k^x3−xy1−y−y3−yx1−x,in which(12)t1=x1−xy2−y−y1−yx2−x,t2=x2−xy3−y−y2−yx3−x,t3=x3−xy1−y−y3−yx1−x.Ift1≥0,t2≥0,t3≥0 or t1≤0,t2≤0,t3≤0, then the intersection point is placed in the triangle element. The same method could be applied to the projection intersection of lines parallel to the x-axes and y-axes.


### 3.1.3. Assign Values

The intersection points are marked and the electromagnetic parameters are assigned, and then, the grid profile data file for the FDTD-based calculation is generated. Note that in this step, the cell whose intersection is on the target region is marked as 1; otherwise, its mark is set as 0.


## 3.1.1. Solve the Intersection Point between the Line and the Plane Where the Triangular Plane Element Is Located

The coordinates of the three vertices of the triangle surface element are specified byAx1,y1,z1, Bx2,y2,z2, and Cx3,y3,z3 (Figure 8). In the rectangular coordinates system, the equation of the line l parallel to the z-axis could be stated as follows:(4)x=i·dx,y=j·dy.Figure 8
The intersection of the trigonometric elementABC and FDTD grid lines.The equation of the plane, in which any triangular surface element is located on the target surface, is given as follows:(5)y1z11y2z21y3z31x+z1x11z2x21z3x31y+x1y11x2y21x3y31z=x1y1z1x2y2z2x3y3z3.Let us define the following coefficients:(6)A=y1z11y2z21y3z31,B=z1x11z2x21z3x31y,C=x1y11x2y21x3y31,D=x1y1z1x2y2z2x3y3z3,then, equation (5) could be rewritten in a more compact form as follows:(7)Ax+By+Cz=D.Simultaneous application of equations (4) and (5) leads to(8)x=i·dx,y=j·dy,Ax+By+Cz=D.Thus,(9)z=1CD−Ai·dx−Bj·dy.The intersection of the line defined by equation (4), which is parallel to the z-axis, and the plane introduced by equation (7), where the triangular plane element is located, is then derived in the following form:(10)x=i·dx,y=j·dy,z=1CD−Ai·dx−Bj·dy.


## 3.1.2. Determine Whether the Intersection Point Is Located in the Triangular Plane Element

The following is to further determine whether the intersection point is placed in the triangular plane. The specific idea is to determine whether the intersection point of the line provided by equation (10) and the coordinate plane xoy is located in the projection of the triangular plane on the coordinate plane xoy. Assume that the projection of a triangular plane element in the coordinate plane xoy is triangle ABC, and the projection of the intersection point between the line and the triangular plane element in the coordinate plane is denoted by P, as demonstrated in Figure 9. Figure 9(a) clearly represents the case where the intersection point is placed inside the triangular plane element; on the other hand, the case where the intersection point is located outside this specific zone is illustrated in Figure 9(b).Figure 9
The intersection point P is placed inside and outside of the projectionABC of the triangular plane element: (a) scenario 1 and (b) scenario 2.
(a)(b)The specific algorithm to assess whether the intersection point is located inside the triangular plane element is explained as follows. Consider the intersection point P in the planexoy, where the line and triangular plane element are located; then, the three vertices of the triangle could be constructed on the basis of three vectors PA,PB,andPC (Figure 9). Let us evaluate the cross product of the three vectors in turn as follows: PA×PB,PB×PC,andPC×PA; if the directions of the resulting vectors after their cross product are the same or one or two of them are equal to 0, then the intersection point of the line in equation (10) and the triangular plane element is located inside the triangular plane element; otherwise, the direction of the produced vectors after their cross product would be different, that is, one of them is different from the other two. The intersection point between the above line introduced by equation (10) and the triangular plane element is not inside the triangular plane element. According to the performed analysis, it can be assumed that the coordinates of the three vertices of the projection triangle of the triangular plane element would be Ax1,y1,0, Bx2,y2,0, and Cx3,y3,0, and the intersection point is Px,y,0; then, the three vectors PA,PB, and PC can be represented as follows:(11)PA=x1−xi^+y1−yj^,PB=x2−xi^+y2−yj^,PC=x3−xi^+y3−yj^,PA×PB=k^x1−xy2−y−y1−yx2−x,PB×PC=k^x2−xy3−y−y2−yx3−x,PC×PA=k^x3−xy1−y−y3−yx1−x,in which(12)t1=x1−xy2−y−y1−yx2−x,t2=x2−xy3−y−y2−yx3−x,t3=x3−xy1−y−y3−yx1−x.Ift1≥0,t2≥0,t3≥0 or t1≤0,t2≤0,t3≤0, then the intersection point is placed in the triangle element. The same method could be applied to the projection intersection of lines parallel to the x-axes and y-axes.


## 3.1.3. Assign Values

The intersection points are marked and the electromagnetic parameters are assigned, and then, the grid profile data file for the FDTD-based calculation is generated. Note that in this step, the cell whose intersection is on the target region is marked as 1; otherwise, its mark is set as 0.


## 3.2. Establishment of the FDTD Calculation Model

The digital recognition result and the triangle set generated by triangulation are extracted, and the three-dimensional model of the mountain peak is constructed according to the digital-curve index relationship developed during the connected domain segmentation [24, 25]. A closed surface filling is performed for nonclosed surfaces. Finally, the triangular patches are shaded by the geometry shader of the OpenGL graphics library, the wireframe is filled, and the lighting is set to make the model close to the real effect (Figure 10).Figure 10
Three-dimensional model of the mountain: (a) grid model and (b) entity.
(a)(b)The projection intersection approach is exploited to transform the derived triangular patches into the required FDTD computational grid model, as demonstrated in Figure11.Figure 11
Yee grid model of the mountain.


## 3.3. Application Examples

The established mountain is used as a calculation model to analyze the field distribution of the mountain at various time steps and the waveforms of monitoring points at different heights. The HEMP is often exploited to examine the electromagnetic scattering characteristics adjacent to the ground and the coupling effect of cables [26, 27]. The medium parameters of the mountain are εr=4, σ=0.01S/m, and μr=1, and the FDTD discrete cell size is set as 10/10 m. The incident angle is specified by θ=180o, α=0o, and φ=0o. The field distribution of xoz and yoz surfaces of the mountain at 300-time steps is calculated under the HEMP incidence condition, as shown in Figure 12.Figure 12
Electric field distribution map of the mountain section: (a)yoz cross-section and (b) xoz cross-section.
(a)(b)Let us select two monitoring points at different heights to observe the time-domain waveform of the monitoring points. These two points with their coordinates are represented as follows: spot 1 (5, 5, 100) and spot 2 (5, 5, 80). The time-domain waveforms of the monitoring points are depicted in Figure13. The demonstrated results in Figure 13 indicate that the HEMP pulse reaches spot 1 first, and the peak field strength of the monitoring point attenuates 23.38% of the incident wave. After a certain time, the pulse reaches spot 2, and the peak field intensity at the monitoring point attenuates to 8.47% of the incident wave.Figure 13
Time-domain waveform of the electric field.


## 4. Conclusion

This study deals with an implementation scheme from the papery contour map to the three-dimensional triangular patch model and finally to the FDTD calculation model. The scheme combines image processing, SVM pattern recognition, and the Bowyer–Watson triangulation algorithm. Each step of the algorithm in the current work is realized by the self-controllable code compiled by the author’s team. The method proposed in this study reduces the dependence on foreign commercial software, and the established model exhibits a good visualization effect and could be seamlessly connected with the subsequent calculation program. It is anticipated that the obtained results from this work will provide a solid foundation for developing the FDTD electromagnetic simulation software and exploiting it to deal with the interaction between electromagnetic waves and complex terrains.



---

*Source: 6986225-2023-01-30.xml*