*Article* **An Empirical Study on the Effects of Temporal Trends in Spatial Patterns on Animated Choropleth Maps**

**Paweł Cybulski**

Department of Cartography and Geomatics, Faculty of Geographical and Geological Sciences, Adam Mickiewicz University Poznan, 61-712 Pozna ´n, Poland; p.cybulski@amu.edu.pl

**Abstract:** Animated cartographic visualization incorporates the concept of geomedia presented in this Special Issue. The presented study aims to examine the effectiveness of spatial pattern and temporal trend recognition on animated choropleth maps. In a controlled laboratory experiment with participants and eye tracking, fifteen animated maps were used to show a different spatial patterns and temporal trends. The participants' task was to correctly detect the patterns and trends on a choropleth map. The study results show that effective spatial pattern and temporal trend recognition on a choropleth map is related to participants' visual behavior. Visual attention clustered in the central part of the choropleth map supports effective spatio-temporal relationship recognition. The larger the area covered by the fixation cluster, the higher the probability of correct temporal trend and spatial pattern recognition. However, animated choropleth maps are more suitable for presenting temporal trends than spatial patterns. Understanding the difficulty in the correct recognition of spatio-temporal relationships might be a reason for implementing techniques that support effective visual searches such as highlighting, cartographic redundancy, or interactive tools. For end-users, the presented study reveals the necessity of the application of a specific visual strategy. Focusing on the central part of the map is the most effective strategy for the recognition of spatio-temporal relationships.

**Keywords:** geomedia; cartographic animation; choropleth map; temporal trend; spatial pattern; pattern recognition; eye tracking

### **1. Introduction**

Dynamic maps, so-called cartographic animations, are part of multimedia cartography [1–3]. The advancement of digital technology has brought significant developments. Existing animated maps have basic time navigation features such as a time sliders, using audio narration instead of a legend, or taking advantage of modern APIs [4,5] and JavaScript libraries [6,7]. Cartographic animation is consistent with the concept of geomedia, which means communicating about geographic space through multimedia related to spatio-temporal dimensions [8]. Owing to this, cartographic animation finds application not only in geovisualization [9], but also in augmented reality systems [10], movies [11], and computer games [12].

However, as Harrower [13,14] correctly noted, the problem with animated maps is their perception. Changes in space on maps require the user's memorization of spatial patterns, and not only comparison, observation, and interpretation. The problem of dynamic maps' perception relates to a multitude of changes such as disappearance, value changes, or movement [15]. Nevertheless, dynamic maps play an important role in visual analysis [16]. Hegarty et al. [17] and Lee et al. [18] gave evidence that people process animated maps through a series of static snapshots rather than a dynamic representation. Juergens [19] suggests that design issues such as map projections could be also a source of misleading interpretation of temporal geographic data.

Much scientific research has dealt with the animated map perception issue. Some research is concerned with the perception of smooth and abrupt changes [20–22]. Other

**Citation:** Cybulski, P. An Empirical Study on the Effects of Temporal Trends in Spatial Patterns on Animated Choropleth Maps. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, 273. https:// doi.org/10.3390/ijgi11050273

Academic Editor: Wolfgang Kainz

Received: 1 March 2022 Accepted: 15 April 2022 Published: 20 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

studies are concerned with comparing effectiveness against static maps [23,24]. Part of the research is related to searching for the most efficient solution in animated map design [25,26]. These studies focus on perception issues, e.g., change blindness [27–30].

Eye tracking technology allows one to record eye movements in real time. The main two components of this movement are saccades and fixations [31]. Saccades are rapid gaze movements from one location to another. The second component is fixation, which relates to a relatively stationary gaze on a single location. There is also a specific eye movement, i.e., smooth pursuit. This is slow eye tracking of a moving objects [32]. On dynamic maps, it concerns, among other things, the movement of point objects along the static route [33]. However, the dominant change in the choropleth map is the change in enumeration units [34]. Therefore, the user's visual strategy could rely on the specific sequence of fixations and saccades (e.g., fixation duration and saccadic amplitude). When users search for trends and patterns on a map, the visual strategy could be used to determine the effective performance. The connection between eye tracking measures and the cognitive processing of a map lies in the idea that fixating objects correlate with the cognitive processing of the fixated objects. This assumption is well established in cognitive studies such as Just and Carpenter [35] or Rayner [36] and also in recent cartographic research [37]. Thus, the following research questions arise: What is an effective visual strategy in pattern searching on a dynamic map? Can different spatial patterns be easier to detect? Are components of spatio-temporal relationships recognized with the same level of effectiveness?

Both spatial patterns and temporal trends are the subject of research in the context of the methodology adopted by Andrienko et al. [38]. In this regard, the main aim of this research was to assess the effectiveness of pattern and trend recognition on animated choropleth maps. Eye tracking could be used to help understand the differences between effective and ineffective recognition among study participants. Knowing that a specific spatial pattern or temporal trend is difficult to detect, end-users should be aware of a specific visual strategy that supports correct recognition. Additionally, we state the hypothesis that effective spatial pattern and temporal trend recognition on a choropleth map are related to participants' visual behavior. To verify this research assumption, a map-based empirical study with participants and eye tracking technology was conducted. Therefore, in Section 3, the adopted methodology is presented, in particular: participants, equipment, materials, procedure, and analysis.

### **2. Related Research**

Spatio-temporal data include data classifications according to changes occurring over time [39,40]. One can distinguish existential changes (disappearance/appearance), spatial property changes (visual variables), or thematic property changes (increase/decrease). Earlier, Peuquet [41] distinguished the space (where), time (when), and object (what) components of spatio-temporal data. Andrienko et al. [38] extended this issue. Identification and comparison tasks could serve as elementary searches or general searches since there is the possibility of distinguishing spatio-temporal map-based tasks. Such cartographic tasks do not contradict the phenomena taking place in geographical space. Moreover, many phenomena that change in space and time show regularities [42]. The correct identification of regularities is possible, among others, owing to an appropriate cartographic presentation. To estimate the correctness of the users' perception of these regularities, one needs to design appropriate map-based tasks [43,44].

Regularities in geographical space can be viewed as spatial patterns or temporal trends [45]. The spatial pattern is related to the distribution of a phenomenon in physical or administrative space [46]. In geographic research, one can distinguish several patterns: clustered, peripheral, grid, linear, radial, or dispersed patterns [47–50]. In animated maps, these patterns change over time. However, despite the changes, several trends occur. The most basic temporal trends are an increase in phenomena, a decrease in a phenomenon, or the absence of a trend [51,52]. Of course, a more detailed breakdown of temporal

trends considers increases and decreases, e.g., linear or logarithmic. In addition, it can also be cyclical.

Blok et al. [53] presented research on the visualization of the relationship between spatial patterns in time. They focused on synchronization, which is the main dynamic variable in cartographic animation. In their research, they designed an alternative visual exploration of spatio-temporal relationships. The results showed that one synchronous animation is sufficient to display juxtaposed data. However, it was found that, in some cases, the relationship between spatial patterns requires two alternative visualizations. Therefore, they concluded that empirical research is lacking in the determination of the effectiveness of the proposed tools.

Edsall et al. [54] provided tools for the visualization of spatial and temporal periodicity in geographic data. The proposed tools used both idealized and real-world environmental space–time data. Decomposing geographic time series involved using a three-dimensional Fourier transformation into spatio-temporal periodicities. Owing to these tools patterns and trends were easily observed. However, this study did not involve users' visual cognition to understand the difficulties of pattern and trend recognition.

Griffin et al. [55] compared the visual identification of space–time clusters on animated maps and static small-multiple maps. Their research proved that animated maps contributed to the correct identification of patterns, contrary to small-multiple maps. The animation pace supported different types of clusters. Griffin et al. (2006) found a threshold in the animation pace. They identified that if the animation pace is too slow, correct spatial pattern identification decreases among users. However, in their study, clustered patterns moved in time and space, and they did not test the identification of other types of spatial patterns.

Schiewe [56] conducted empirical research on the visual perception of spatial patterns in choropleth maps. He examined the intuitive ranking of color lightness, the neglect of small areas, and the attention to data classification when interpreting spatial patterns. The study participants detected extreme values, central tendencies, and homogeneities. The results confirmed that a legend improves the effective identification of extreme values. Additionally, the detection rate of extreme values in small areas was increased. However, most participants relied on the color lightness distribution and not on the data classification during pattern interpretation.

Karsznia et al. [57] studied effective spatial pattern recognition in choropleth maps with an optimal hexagonal enumeration unit size. They examined nine sizes of hexagonal units and showed that the enumeration unit size impacts users' spatial pattern recognition abilities. In conclusion, they argued that the largest and the smallest units' were not suitable for indicating spatial patterns. Their research delivered the optimal range of EUS (enumeration size units) without indicating sharp boundaries. However, their experimental procedure did not include changes in the spatial pattern over time.

Changes in the detection of global trends and local outliers on an animated choropleth map were studied by Traun et al. [58]. In a user testing study, they found that outlierpreserving value generalization in space supports local outlier identification. However, generalization did not support the participants' ability to correctly detect the spatial pattern movement. They conclude that spatial pattern recognition is related more to personal features (e.g., age) rather than specific map design.

Summarizing this section, it should be mentioned that spatial pattern cognition was well studied based on cartographic and psychological basics. The research gaps exist in the context of temporal trends in spatial patterns. There has been some indication of the effectiveness of patterns in trend recognition. However, the presented study focused more on trends in patterns, especially their effectiveness in terms of correct recognition and, additionally, users' visual behavior during a map-based task.

### **3. Methodology**

### *3.1. Participants*

Eighteen students of the Faculty of Geographic and Geological Sciences at Adam Mickiewicz University participated in the experiment. The participants described their characteristics: age, gender, and education. A total of 10 of them indicated male gender and 8 indicated female gender. The participants' age range was 21–22 years, and they were in the third year of their studies, taking geodetic and cartographic courses. They did not receive payment and participated voluntarily. They could resign from participation at any time in the study.

### *3.2. Equipment*

To display map stimuli, an 27" ASUS LCD with a 1980 × 1080 px screen resolution was used. A Gazepoint GP3 HD eye tracker with a sampling frequency of 150 Hz recorded the participants' eye movements. The visual angle accuracy was 0.5◦ , allowing up to 35 cm (horizontal) × 22 cm (vertical) movement. The average distance between the participant and the monitor was approximately 62 cm. The calibration procedure used a 9-point method using Gazepoint Control software (www.gazept.com, accessed 12 December 2021).

### *3.3. Materials*

Fifteen animated choropleth maps of the unemployment rate constituted the main cartographic material. Each choropleth had 63 enumeration units classified into 5 classes of a sequential single-hue greyscale according to ColorBrewer 2.0 (5-class Greys). Each of the animated maps lasted 30 seconds and showed changes over twelve months. The study participants viewed each map once.

In the first month of the animation (January), the maps legibly showed (using one color hue) one of the following patterns: radial, clustered, linear, grid, peripheral, or dispersed (Figure 1). In the following months, the enumeration units changed their hue randomly except for the pattern units. The pattern units did change their hue randomly, but changed in a specific direction. There were two possible directions of this temporal trend: an increase (from light to dark) or a decrease (from dark to light). This indicated the phenomenon intensity. The pattern units did not change their hue altogether, but consequently changed in one of the directions of the trend.

The only exception was the dispersed pattern. In this case, all of the enumeration units changed their hue randomly, even the pattern units. As a consequence, the pattern visible at the beginning of the animation was lost and did not follow any trend. Finally, in the last month of the animation (December), each pattern was again clearly visible with one color hue. Figure 2 present thes frames of the animated maps of linear patterns with an increasing trend and dispersed patterns with no trend.

### *3.4. Procedure*

The procedure consisted of presenting all fifteen maps in a random order to each participant. Before the presentation, each study participant was given an introduction to the purpose of the study. The experiment began with a 9-point calibration reaching a gaze sample score of 90% on average. Secondly, participants performed a set of familiarization tasks. Next, they watched fifteen animated maps. After each map, there were three questions regarding spatial pattern and temporal trend (Figure 3 presents an interactive response protocol that appeared after the animation ended):


5 of 16

**Figure 1.** Map stimuli used in the research. Spatial patterns: C: clustered, G: grid, L: linear, P: peripheral, R: radial, D: dispersed. Temporal trends: d: decrease, i: increase. **Figure 1.** Map stimuli used in the research. Spatial patterns: C: clustered, G: grid, L: linear, P: peripheral, R: radial, D: dispersed. Temporal trends: d: decrease, i: increase.

The only exception was the dispersed pattern. In this case, all of the enumeration units changed their hue randomly, even the pattern units. As a consequence, the pattern visible at the beginning of the animation was lost and did not follow any trend. Finally, in the last month of the animation (December), each pattern was again clearly visible with one color hue. Figure 2 present thes frames of the animated maps of linear patterns with

**Figure 2.** Selected frames of a linear and dispersed pattern. An increasing trend is in a linear pattern. The dispersed pattern has no trend. In D4, the spatial pattern is only visible at both ends of animation and did not appear in other months. Li—linear increasing, D4—dispersed pattern map number 4.

**Figure 3.** An interactive protocol for the participants' responses was presented after each animated map ended.

### *3.5. Analysis*

The analysis of participants' visual behavior included eye tracking metrics such as fixation location, fixation count, fixation duration, and saccadic magnitude. Fixation count describes how many times the participants gazed at a specific location. This might indicate the ease in understanding the spatio-temporal data. The more fixations on the timeline, the more focused attention it required. This is related to the fixation duration. Longer durations might indicate a longer duration of attention. Therefore, these issues are related to perception. Saccadic magnitude describes the length of rapid eye movement from one fixation to another. This metric might suggest the difficulty in map interpretation. Participants' shorter saccades correlated with a high number of long-duration fixations, which might result in ineffective spatio-temporal dependencies. The statistical analysis involved 6\*2 Friedman's test since the data had a non-parametrical distribution. The withinsubject factors were the spatial pattern type (six conditions) and the correct recognition of patterns and temporal trends (two conditions). Correct pattern recognition effectiveness was based on the number of correct spatial pattern indications after the animation was viewed. It was a binary response, where 0 was assumed when the participants incorrectly identified a spatial pattern in a given month. Additionally, analyzing the spatial distribution of fixations was based on fixation coordinates recorded with eye tracking and spatial statistics such as Average Nearest Neighbor and Kernel Density.

### **4. Results**

#### *4.1. Spatial Pattern and Temporal Trend Recognition Effectiveness* exception was the radial pattern. Its effectiveness achieved 28%. The effectiveness of iden-

**4. Results** 

Figure 4 presents the results of spatial pattern and temporal trend recognition effectiveness. The participants recognized the lack of a pattern the best (dispersed pattern); 73% of their answers were correct. There was a similar result in the case of recognizing the absence of a trend; 84% of the participants responded correctly. Correct identification at a random point (month) in the animation time was recorded for most of the studied patterns. The pattern recognition effectiveness ranged between 64% and 55%. The only exception was the radial pattern. Its effectiveness achieved 28%. The effectiveness of identifying temporal trends in spatial patterns ranged between 53% and 75%. tifying temporal trends in spatial patterns ranged between 53% and 75%. However, in some cases, the participants recognized the increasing trend more effectively than the decreasing trend. This was the case for the clustered pattern (recognition of 60% of the increasing trends vs. 40% of the decreasing trends), radial pattern (recognition of 60% of the increasing trends vs. 40% of the decreasing trends), and grid pattern (recognition of 58% of the increasing trends vs. 43% of the decreasing trends). In the case of a linear pattern and peripheral pattern, the trends' recognition rates were similar (recognition of 48% of the increasing trends vs. 52% of the decreasing trends).

*ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 7 of 16

such as Average Nearest Neighbor and Kernel Density.

*4.1. Spatial Pattern and Temporal Trend Recognition Effectiveness* 

perception. Saccadic magnitude describes the length of rapid eye movement from one fixation to another. This metric might suggest the difficulty in map interpretation. Participants' shorter saccades correlated with a high number of long-duration fixations, which might result in ineffective spatio-temporal dependencies. The statistical analysis involved 6\*2 Friedman's test since the data had a non-parametrical distribution. The within-subject factors were the spatial pattern type (six conditions) and the correct recognition of patterns and temporal trends (two conditions). Correct pattern recognition effectiveness was based on the number of correct spatial pattern indications after the animation was viewed. It was a binary response, where 0 was assumed when the participants incorrectly identified a spatial pattern in a given month. Additionally, analyzing the spatial distribution of fixations was based on fixation coordinates recorded with eye tracking and spatial statistics

Figure 4 presents the results of spatial pattern and temporal trend recognition effectiveness. The participants recognized the lack of a pattern the best (dispersed pattern); 73% of their answers were correct. There was a similar result in the case of recognizing the absence of a trend; 84% of the participants responded correctly. Correct identification at a random point (month) in the animation time was recorded for most of the studied patterns. The pattern recognition effectiveness ranged between 64% and 55%. The only

**Figure 4.** Results of spatial pattern and temporal trend recognition effectiveness calculated as a percentage of participants' correct responses.

However, in some cases, the participants recognized the increasing trend more effectively than the decreasing trend. This was the case for the clustered pattern (recognition of 60% of the increasing trends vs. 40% of the decreasing trends), radial pattern (recognition of 60% of the increasing trends vs. 40% of the decreasing trends), and grid pattern (recognition of 58% of the increasing trends vs. 43% of the decreasing trends). In the case of a linear pattern and peripheral pattern, the trends' recognition rates were similar (recognition of 48% of the increasing trends vs. 52% of the decreasing trends).

### *4.2. Participants' Visual Behavior*

Figure 5 presents Freidman's test results in the context of three eye tracking metrics describing the participants' visual behavior. Friedman's test detected statistically significant differences between all three eye tracking metrics. Fixation count varied substantially among the selected factors. The highest differences were among clustered patterns juxtaposed with correct pattern and trend recognition. The second large significant difference between correct and incorrect answers was in the peripheral and linear patterns. In both patterns, participants who answered correctly had more fixation on the map in contrast to participants who recognized them incorrectly. Additionally, the radial pattern had significant differences in the fixation count among participants who responded correctly or incorrectly. However, in the case of grid patterns, the participants who correctly recognized spatio-temporal dependencies had fewer fixations than those who responded incorrectly (χ <sup>2</sup> = 82.4; *p* < 0.000001). This shows that, in most cases, if participants had to fixate more

on the animated map, there was a higher probability of incorrect recognition of spatial patterns and temporal trends. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 9 of 16

**Figure 5.** 3 Spatial analysis of participants' visual behavior. **Figure 5.** 3 Spatial analysis of participants' visual behavior.

Based on Euclidean distance Average Nearest Neighbor was calculated for correct and incorrect recognition participants. The results are presented in Table 1, which contains all related values such as Observed Mean Distance (OMD), Expected Mean Distance (EMD), Nearest Neighbor Index (NNI), z-score (z), and *p*-value (*p*). Most fixations, both correct and incorrect recognition, showed statistically significant clustering visual behav-Fixation duration was differed significantly among participants' pattern and trend recognition in radial pattern conditions (χ <sup>2</sup> = 9.6; *p* < 0.00193). Interestingly, the radial pattern was the least effective recognized pattern by participants. Participants who incorrectly recognized spatio-temporal dependencies had longer fixations, while participants who responded correctly had a statistically significantly shorter fixation duration.

ior. The only exception was a dispersed pattern with incorrect recognition participants (random distribution of fixations). The calculated distance included the local coordinate system of the recorded screen. Therefore, the X and Y axes of the coordinate system were a percentage of the monitor height and width. **Table 1.** Average Nearest Neighbor Summary of the fixation distribution of participants' correct and incorrect recognition. **Average Nearest Neighbor Summary** Saccadic magnitude expressed in visual angle degrees was statistically significant among the selected within-subject factors (χ <sup>2</sup> = 3.28; *p* < 0.04). The saccadic magnitude was similar in the case of dispersed and peripheral patterns. The most crucial variation between participants' correct and incorrect recognition of spatio-temporal relations was visible in the case of the clustered, radial, and linear patterns. It is significant that these patterns had the lowest effectiveness in temporal trend detection. In the linear and radial pattern conditions, participants who incorrectly responded had longer saccades, while in the clustered pattern it was the opposite.

**Incorrect Recognition Correct Recognition OMD EMD NNI z** *p* **OMD EMD NNI z** *p* D1 0.0235 0.0299 0.785 −4.460 <0.0000 0,0099 0.0167 0,593 −16.922 <0.0000 D2 0.0148 0.0223 0.663 −7.327 <0.0000 0.0103 0.0159 0.652 −13.643 <0.0000 Inconsistency between the spatial patterns regarding correct and incorrect responses among the selected eye tracking metrics is visible. It should be noted that different spatial patterns require various visual behaviors for effective performance. The main factor crucial for determining effective visual behavior is the fixation count. Statistically significant

D3 0.0267 0.0285 0.936 −0.794 0.4272 0.0112 0.0168 0.667 −13.192 <0.0000 D4 0.0254 0.0292 0.871 −2.722 0.0065 0.0086 0.0136 0.637 −17.654 <0.0000 differences between the saccadic magnitude in the linear, clustered, and radial patterns might indicate that these are more complex patterns for visual analysis. Since the dispersed, peripheral, and grid patterns had similar saccadic magnitude differences, shorter saccades could be crucial for effective recognition of these patterns. This might be confirmed by the finding that these patterns had similar visual strategies according to the fixation count.

Based on Euclidean distance Average Nearest Neighbor was calculated for correct and incorrect recognition participants. The results are presented in Table 1, which contains all related values such as Observed Mean Distance (OMD), Expected Mean Distance (EMD), Nearest Neighbor Index (NNI), z-score (z), and *p*-value (*p*). Most fixations, both correct and incorrect recognition, showed statistically significant clustering visual behavior. The only exception was a dispersed pattern with incorrect recognition participants (random distribution of fixations). The calculated distance included the local coordinate system of the recorded screen. Therefore, the X and Y axes of the coordinate system were a percentage of the monitor height and width.

**Table 1.** Average Nearest Neighbor Summary of the fixation distribution of participants' correct and incorrect recognition.


This clustering of participant fixations was confirmed and located on a map stimulus with Kernel Density. In most cases, participants who correctly recognized the spatial patterns and temporal trends had different fixation concentration areas from those of participants with incorrect responses. As can be seen in Figure 6, correct recognition corresponds to a specific visual behavior. This consisted of the concentration of fixations near the spatial pattern. However, in particular, what distinguished correct recognition from incorrect recognition was that, if the recognitions were correct, the fixation clusters covered a larger area of the map. It is significant that the correct recognition fixation cluster covers a larger part of the spatial pattern than in the case of incorrect recognition. The only exception was the clustered pattern. For the decreasing temporal trends, the participants' correct recognition fixation cluster covers a larger area. However, it is slightly shifted from the pattern. Therefore, the cluster area of fixations is outside the pattern on the map. In the case of increasing temporal trends, the incorrect recognition fixation cluster covers a larger area than that of the correct recognition cluster. Therefore, there is a difference between the participants' visual behavior.

**Figure 6.** Fixation density clusters for correct and incorrect recognition of spatial patterns and temporal trends. **Figure 6.** Fixation density clusters for correct and incorrect recognition of spatial patterns and temporal trends.

When analyzing the lack of patterns (dispersed pattern), the results are not as clear as when there was a presence of a spatial pattern. Participants' fixation clusters occur in both incorrect and correct recognition. In the case of the D1, D2, and D4 maps, the cluster distribution is similar to the case of the occurrence of the spatial patterns and temporal trends. Correct recognition fixation clusters cover a larger area. The situation is different when analyzing the D3 and D5 maps. In both cases, the incorrect recognition fixation clusters cover a larger area than the correct recognition fixation clusters. Figure 7 presents the result of Kernel Density based on the participants' fixation and recognition. When analyzing the lack of patterns (dispersed pattern), the results are not as clear as when there was a presence of a spatial pattern. Participants' fixation clusters occur in both incorrect and correct recognition. In the case of the D1, D2, and D4 maps, the cluster distribution is similar to the case of the occurrence of the spatial patterns and temporal trends. Correct recognition fixation clusters cover a larger area. The situation is different when analyzing the D3 and D5 maps. In both cases, the incorrect recognition fixation clusters cover a larger area than the correct recognition fixation clusters. Figure 7 presents the result of Kernel Density based on the participants' fixation and recognition.

**Figure 7.** Fixation density clusters for correct and incorrect recognition of spatial patterns and temporal trends when analyzing dispersed patterns. **Figure 7.** Fixation density clusters for correct and incorrect recognition of spatial patterns and temporal trends when analyzing dispersed patterns.

#### **5. Discussion 5. Discussion**

The effectiveness, measured as the correct recognition of temporal trends and spatial patterns, differed among the analyzed stimuli. It appears that the study participants were better at recognizing temporal trends for most of the maps. This is especially beneficial for animated maps related to injury prevention. Cinnamon et al. [59] showed that participants could detect temporal trends relatively quickly, in comparison to static and interactive maps. Therefore, this study consists of empirical evidence supporting the previous hypothesis regarding the usefulness of animated products for the visualization of temporal trends [14,60]. Choropleth maps are more effective when the data have a peripheral spatial pattern or the lack of one (effectiveness above 64%). For the visualization of data with a radial pattern, one should search for other mapping methods due to the low choropleth map effectiveness (below 29%). The presented research had only one animation speed. However, according to Griffin et al. [55], different animation speeds may support the recognition of some patterns. Additionally, following Roth's [4] suggestions of including interactive tools could increase animated map effectiveness. Considering hexagonal enumeration units might also help in the recognition of less effective patterns (e.g., radial) [57]. Even some techniques of highlighting the most crucial information might be worthy of consideration [25]. The effectiveness, measured as the correct recognition of temporal trends and spatial patterns, differed among the analyzed stimuli. It appears that the study participants were better at recognizing temporal trends for most of the maps. This is especially beneficial for animated maps related to injury prevention. Cinnamon et al. [59] showed that participants could detect temporal trends relatively quickly, in comparison to static and interactive maps. Therefore, this study consists of empirical evidence supporting the previous hypothesis regarding the usefulness of animated products for the visualization of temporal trends [14,60]. Choropleth maps are more effective when the data have a peripheral spatial pattern or the lack of one (effectiveness above 64%). For the visualization of data with a radial pattern, one should search for other mapping methods due to the low choropleth map effectiveness (below 29%). The presented research had only one animation speed. However, according to Griffin et al. [55], different animation speeds may support the recognition of some patterns. Additionally, following Roth's [4] suggestions of including interactive tools could increase animated map effectiveness. Considering hexagonal enumeration units might also help in the recognition of less effective patterns (e.g., radial) [57]. Even some techniques of highlighting the most crucial information might be worthy of consideration [25].

The recognition of spatial patterns was found to be more difficult than the recognition of temporal trends. Based on the results, the most effective was the dispersed pattern. This proves previous assumptions that animated conditions support both the lack of pattern recognition and strong pattern recognition [55]. This is also true for the recognition The recognition of spatial patterns was found to be more difficult than the recognition of temporal trends. Based on the results, the most effective was the dispersed pattern. This proves previous assumptions that animated conditions support both the lack of pattern recognition and strong pattern recognition [55]. This is also true for the recognition of temporal trends, which did not appear while displaying dispersed pattern maps. It was found that most of the patterns analyzed in this study had similar effectiveness (from 55% to 64%). The only exception was the radial pattern, which had the lowest effectiveness (28%). This type of spatial pattern is related to road and railroad mapping [61].

The study results report statistically significant differences between correct (effective) and incorrect (ineffective) recognition of temporal trends and spatial patterns. These differences are related to participants' visual behavior. As Shiewe [56] showed, different choropleth mapping techniques could support the successful interpretation of cartographic data. However, this study used only one mapping method. Therefore, the spatial pattern classification is related to the color intensities, and a legend was rarely fixated upon. Differences in the participants' correct and incorrect recognition were visible in the saccadic magnitude in the case of the clustered, radial, and linear patterns. Two of these patterns were the only ones where the participants recognized more increasing trends than decreasing ones. In Cybulski's [33] study, different saccade lengths are related to variability in spatial abilities or other personal features. Additionally, Dong et al. [62] presented the role of previous experiences in differences between the spatial behavior of novices and experts. Finally, Trauns et al.'s [58] study confirmed the reported results regarding the role of map users' personal features, which may determine map effectiveness to a greater extent than the specific map design itself.

The spatial analysis of the visual behavior showed that the study participants fixated mostly on specific areas. Fixation clusters appeared among correct and incorrect recognitions. The fixation clusters were located in the central part of the map. The only exceptions were maps with a peripheral pattern and some dispersed patterns (D3), which generated clusters that were shifted to near the pattern. Previous studies by Dong et al. [63] and Cybulski and Krassanakis [64] have explained this type of visual behavior. To observe the whole map during animation, the participants mainly observed its central part. The participants' sporadically observed peripheral parts of the map. A larger area covered by longer duration fixations characterized the participants' correct recognition fixation clusters. This might be related to change blindness when a participant fails to notice a change on a dynamic map [63,65].

### **6. Conclusions**

The main aim of this research was to investigate the effectiveness of spatial pattern and temporal trend recognition assessment on a choropleth map. The presented study achieved this goal. The results show that spatial pattern recognition on a choropleth map varies, and that this depends on the pattern itself. The radial pattern was the least effective, and the peripheral pattern was the most effective. Despite this, the results were compared to those for maps with no pattern (dispersed pattern) and showed that recognizing the lack of a pattern is even more effective. Therefore, animated choropleth maps could support the visualization of specific spatial patterns. However, they are more suitable for presenting temporal trends. Temporal trends in spatial data had a recognition effectiveness of 64%.

The empirical results confirm the research hypothesis that effective spatial pattern and temporal trend recognition on a choropleth map is related to participants' visual behavior. Visual attention clustered in the central part of the choropleth map supports effective spatio-temporal relationship recognition. The larger the area covered by the fixation cluster, the higher the probability of correct temporal trend and spatial pattern recognition. This has specific implications for map makers and map users. For the visualization of temporal trends in spatial patterns using a choropleth map, one might consider a simple design since the interpretation of this type of data is more demanding than a static map. Understanding the difficulty in the correct recognition of spatio-temporal relationships might be a reason for implementing techniques that support effective visual searches such as highlighting, cartographic redundancy, or interactive tools. For end-users, the presented study reveals the necessity of the application of a specific visual strategy. Focusing on the central part of the map was found to be the most effective strategy for the recognition of spatio-temporal relationships.

The study limitations relate to the number of participants. Due to the COVID-19 pandemic, we experienced difficulties with user testing. This was related to legal restrictions and the necessity to maintain social distance and disinfect equipment and rooms. Therefore, it would be beneficial to increase the number of research participants in future studies.

In future studies, it would be interesting to compare the presented results with different spatial pattern and temporal trend mapping methods. This would help find the best cartographic solutions for spatio-temporal relationships, which are ineffective using a choropleth map. The presented study does not exhaust the topic of recognizing patterns and trends on maps. Users might experience more complex spatial patterns or more than one pattern on a map. Additionally, a temporal trend might initially be an increasing trend, and afterward become a decreasing trend. These more complicated spatio-temporal relationships require interdisciplinary studies, including cartography and psychology. Future cartographic visualization designers might benefit from these empirical results by taking into account the awareness that map users do not perceive space and time with the same level of effectiveness. The author believes that future temporal cartographic products need to be reconsidered in terms of temporal legend design to maintain effective communication with the user.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data that support the findings of this study are openly available in the "Harvard Dataverse" at https://doi.org/10.7910/DVN/NFJW6B.

**Acknowledgments:** The author would like to thank all anonymous reviewers for their fruitful comments and recommendations. The paper is the result of research on visualization methods carried out within statutory research in the Department of Cartography and Geomatics, Faculty of Geographical and Geological Sciences, Adam Mickiewicz University Poznan, Poland.

**Conflicts of Interest:** The author declares no conflict of interest.

### **References**


**An-Bo Li 1,2,3,\*, Hao Chen <sup>1</sup> , Xiao-Feng Du <sup>4</sup> , Guo-Kai Sun <sup>4</sup> and Xian-Yu Liu <sup>1</sup>**


**Abstract:** Most fabrication methods for three-dimensional (3D) geological symbols are limited to two types: directly increasing the dimensionality of a 2D geological symbol or performing appropriate modeling for an actual 3D geological situation. The former can express limited vertical information and only applies to the three-dimensional symbol-making of point mineral symbols, while the latter weakens the difference between 3D symbols and 3D geological models and has several disadvantages, such as high dependence on measured data, redundant 3D symbol information, and low efficiency when displayed in a 3D scene. Generating a 3D geological symbol is represented by the process of constructing a 3D geological model. This study proposes a parametric modeling method for 3D fold symbols according to the complexity and diversity of the fold structures. The method involves: (1) obtaining the location of each cross-section in the symbol model, based on the location parameters; (2) constructing the middle cross-section, based on morphological parameters and the Bezier curve; (3) performing affine transformation according to the morphology of the hinge zone and the middle section to generate the sections at both ends of the fold; (4) generating transition sections of the 3D symbol model, based on morphing interpolation; and (5) connecting the point sets of each transition section and stitching them to obtain a 3D fold-symbol model. Case studies for different typical fold structures show that this method can eliminate excessive dependence on geological survey data in the modeling process and realize efficient, intuitive, and abstract 3D symbol modeling of fold structures based on only a few parameters. This method also applies to the 3D geological symbol modeling of faults, joints, intrusions, and other geological structures and 3D geological modeling of typical geological structures with a relatively simple spatial morphology.

**Keywords:** fold structure; 3D geological symbols; parametric modeling; 3D geological modeling; 3D GIS; visual suggestiveness; symbols efficiency

### **1. Introduction**

In recent years, traditional geographic information system (GIS) technology and its applications have undergone revolutionary changes [1,2], and it has transformed from 2D GIS to 3D GIS, with a primary focus on 3D modeling, 3D scene display, and 3D spatial analysis [3]. In particular, 3D GIS has become one of the iconic contents of GIS technology today and for the future [4,5]. At present, Digital Earth systems, including Google Earth, NASA WorldWind and Cesium, have been widely embraced by geoscientists in various disciplines. These 3D GIS-related technologies are convenient tools for promoting the scientific research of different application scenarios and contribute to integrating geospatial information, improving visualization capabilities, exploring spatio-temporal changes, and communicating scientific results [6–8]. In this context, traditional 2D map symbols cannot easily meet the demand for the intuitive expression of complex geological objects in 3D

**Citation:** Li, A.-B.; Chen, H.; Du, X.-F.; Sun, G.-K.; Liu, X.-Y. Parametric Modeling Method for 3D Symbols of Fold Structures. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, 618. https://doi.org/ 10.3390/ijgi11120618

Academic Editors: Wolfgang Kainz, Beata Medynska-Gulij, David Forrest and Thomas P. Kersten

Received: 3 November 2022 Accepted: 12 December 2022 Published: 13 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

GIS. The research and application of 3D symbols have gradually become a hot topic in the current GIS field [9,10].

Similar to traditional 2D symbols, 3D symbols are also the language of map visualization, which aims to express the 3D spatial distribution and attribute information of geographical entities [11,12]. The essence of a 3D symbol is that it is a 3D model that directly reflects the spatial morphology of the corresponding geographical object. Therefore, the production of 3D symbols is the process of generating 3D models. Given the growing demand for 3D symbols in 3D GIS applications, such as virtual reality (VR)\augmented reality (AR)\mixed reality (MR), scholars have proposed various methods for creating 3D symbols. These methods can be roughly divided into dimension-raising methods for planar map symbols [13], 3D model construction methods [14–17], and parametric modeling methods [18,19]. In the following, a summary of the modeling methods for 3D geological symbols, from these three aspects, is provided.

Few studies have been conducted on 3D symbol-making methods based on the dimension-raising of planar map symbols. Zhu et al. [13] proposed an automatic dimensionraising scheme from 2D vector mineral point symbols to 3D symbols, based on the expression characteristics of the 2D vector point symbols. This method has a high degree of automation and is mainly applied to the 3D symbol fabrication of mineral point symbols. However, because of the failure to entirely use the essential geometric parameter information of related geological objects, this method is unsuitable for generating symbols of 3D geological structures.

Fabricating 3D symbols based on a 3D model construction method is the primary method of rapidly generating 3D symbols. Modeling methods can be further subdivided into manual 3D symbol construction methods, based on general 3D modeling software, and automatic 3D symbol construction methods, based on particular modeling software. Early research has primarily focused on manual modeling methods based on general 3D modeling software. Xu et al. [20] used computer aided design (CAD) and MDL models to build 3D solid model symbols for 3D urban scenes, and Gu et al. [21] used Sketchup to build a 3D solid symbol library containing point, line, face, and body symbols. With the continued maturation and development of 3D modeling technology, an automatic construction method based on special modeling software has become the primary method of performing 3D modeling in different fields. According to the source of modeling data, the methods of 3D geological modeling mainly include borehole-based modeling [14,22], section-based modeling [23–25], planar geological map-based modeling [26], and multisource data modeling [22,27–29].

Furthermore, based on the application of knowledge in the modeling process, 3D geological modeling methods can be divided into explicit modeling methods [16,30–32] and implicit modeling methods [17,33–35]. All of the 3D geological modeling methods described above can be used to rapidly construct 3D geological symbols. However, these methods require sufficient modeling data, and the model results are too elaborate. Directly applying the relevant 3D geological model to a 3D scene as a 3D symbol leads to low display efficiency, which contradicts the abstract principle of the symbol. Implicit modeling methods, which have developed rapidly in recent years, have dramatically reduced the dependence on modeling data by implicitly defining the geological interface as the isosurface of one or more scalar fields of 3D space [32,34,35]. However, owing to its refined modeling requirements, this method is unsuitable for creating many 3D geological symbols.

In the parametric modeling method, the characteristic parameters of natural objects are used to automatically control and construct 3D models based on parametric technology [18]. Parametric modeling was first invented by Rhino, which is a 3D draughting software that evolved from AutoCAD. The key advantage of parametric modeling is that, when setting up a 3D geometric model, the shape of the model's geometry can be changed as soon as the parameters, such as the dimensions or curvatures, are modified [19,36,37]. Parametric modeling has long been the primary method of machinery modeling in the CAD field. This method is suitable for structured entity modeling, owing to its low data dependency and

high modeling efficiency. In recent years, parametric modeling methods have been promoted and applied to 3D tree modeling [38,39], 3D human body modeling [40], 3D building (structure) modeling [36], and many other fields. However, owing to the heterogeneity and non-parametric characteristics of geological bodies, parametric modeling methods have not played a role in 3D geological modeling. Using a geological map sketch and symbol information, Amorim et al. [34] realized the rapid construction of 3D models of folds, faults, and other geological structures by quantitatively calculating the topological structure and stratigraphic contact relationship of geological maps. Through parametric analysis and 2D geological symbol application, this method shows a specific idea of parametric modeling and effectively reduces the dependence on modeling data. However, this method is not designed for 3D geological symbols; therefore, its relatively refined modeling requirements do not entirely eliminate the dependence on geological maps and fail to fully reflect the abstract and 3D characteristics of 3D geological symbols.

Map symbols have prominent parametric characteristics because of their high abstraction as either a 2D or multidimensional symbol [41]. The parametric modeling method should become the primary method of 3D symbol modeling in the future because it requires only a small amount of feature parameter data to support the construction of 3D models with abstract expression characteristics. Similarly, according to the characteristics and requirements of 3D geological symbols based on parametric modeling methods, studying targeted 3D geological symbol fabrication methods should be the primary trend of future 3D geological symbol modeling research and application. In addition, although geological bodies have various types and complex structures, they still have apparent spatial distribution characteristics, and their geometric shapes can be mathematically simulated [41]. Implicit modeling methods that have been developed in recent years, particularly geological structure 3D modeling methods based on geological map sketches and geological symbols proposed by Amorim et al. [34] have effectively verified the feasibility of this idea.

As one of the most common geological structures in the crust, folds are of great significance for studying oil and gas traps and other mineral resources [42]. Furthermore, their geometric shapes are complex, and their types are diverse. Therefore, this study intends to take fold structure as the study object and use the parametric modeling method to discuss the rapid construction method of three-dimensional geological structure symbols. The remainder of this paper is organized as follows: Section 2 introduces the methodology; Section 3 presents the experimental results; Section 4 presents a discussion; and Section 5 presents the conclusions and future work.

### **2. Methodology**

### *2.1. Modeling Parameters of the Fold Structures*

In general, fold structures not affected by strong denudation are wavy structures in shape [34] whose spatial morphology can be described by fold elements, such as limbs, hinges, and axial surfaces [42–44]. The limb refers to the rock strata on both sides of the core, the hinge is the connecting line of the maximum bending point on each cross-section of the fold, and the axial surface refers to the surfaces formed by connecting the hinge lines on each adjacent fold stratum (Figure 1).

The formation mechanism of the folds is closely related to their stress mode and deformation environment and the deformation behavior of the rock stratum [44]. The diversity of fold formation mechanisms determines the complexity of the fold structures and the diversity of the fold types. Different types of folds are mainly classified according to their position, axial occurrence, and section morphology (Table 1).

**Figure 1.** Diagram of fold structure characteristics. **Figure 1.** Diagram of fold structure characteristics.



series of cross‐cutting sections based on geological maps, and using the contour compar‐ ison algorithm [45] to build a 3D model of the fold structure, is the most feasible method. The two critical points of this method are to reasonably plan the location of the section Due to the lack of borehole data in the fold structure development area, generating a series of cross-cutting sections based on geological maps, and using the contour comparison algorithm [45] to build a 3D model of the fold structure, is the most feasible method. The two

Due to the lack of borehole data in the fold structure development area, generating a

critical points of this method are to reasonably plan the location of the section lines based on the location parameters and to efficiently generate a series of cross-sections based on the section morphology parameters. First, the location parameter of the cross-section (Table 2) was mainly used to determine the location and strike of the cross-sections. The relevant parameters can be directly obtained from geological maps or geological structure maps. Second, the morphological parameters of the cross-section (Table 2) represent the primary information that controls the cross-sectional morphology. The relevant parameters mainly include the axial dip angle, interlimb angle, limb occurrence, and stratigraphic information. These parameters control the local details of each stratum in the cross-sections of the fold structure. The geometric characteristics of the different fold types can be constrained and controlled by the different value ranges of the fold morphology parameters.


**Table 2.** Modeling parameters of the fold structure.

Based on the modeling parameters outlined above, this study presents a parametric 3D modeling method for fold symbols (Figure 2). The specific modeling steps primarily include: (1) the parametric generation of the cross-section lines based on the section location parameters; (2) the parametric generation of the middle cross-section based on the section morphology parameters; (3) the generation of the cross-cutting section at both ends based on the affine transformation; (4) the generation of the transition sections based on the morphing interpolation; and (5) the stitching and attribute assignment of the fold model.

**Figure 2.** Parametric modeling process of 3D fold symbols. **Figure 2.** Parametric modeling process of 3D fold symbols.

#### *2.2. Generation of Cross‐Section Lines 2.2. Generation of Cross-Section Lines*

The premise of generating a 3D fold‐symbol model based on a contour comparison algorithm is to generate cross‐sections by locating each cross‐section of the fold. The strat‐ igraphic information in the cross‐section of the fold structure should be comprehensive and reflect the typical characteristics of the fold structure. Therefore, the position of the cross‐sections should be determined according to the distribution of various strata of the fold structure and the direction of the hinge line (Figure 3a). This process makes the sym‐ bol model closer to the actual situation of the structure, which implies a higher degree of The premise of generating a 3D fold-symbol model based on a contour comparison algorithm is to generate cross-sections by locating each cross-section of the fold. The stratigraphic information in the cross-section of the fold structure should be comprehensive and reflect the typical characteristics of the fold structure. Therefore, the position of the cross-sections should be determined according to the distribution of various strata of the fold structure and the direction of the hinge line (Figure 3a). This process makes the symbol model closer to the actual situation of the structure, which implies a higher degree of model restoration. *ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 7 of 22

**Figure 3.** Determination of cross‐section locations: (**a**) extraction of location parameters based on a geological map and (**b**) location of sections based on location parameters. S1‐S3 are different stratum **Figure 3.** Determination of cross-section locations: (**a**) extraction of location parameters based on a geological map and (**b**) location of sections based on location parameters. S1-S3 are different stratum codes.

The construction of the middle cross‐section is essentially the process of generating

The stratum morphology of the fold structure is generally smooth and regular and shows specific function laws, and it can be described in terms of amplitude and wave‐ length [42]. Therefore, it is feasible to use a quadratic Bessel curve to simulate the stratum boundary of the fold structure. In addition, folds are usually formed by the hinge zone and its two connected limbs, and the stratum morphologies of the two limbs tend to be different. Therefore, dividing the boundary line of the fold stratum into two segments according to the hinge zone is appropriate, and the Bezier curve is used to express this

**Legend**

Referring stratum boundary Remaining stratum boundary

According to the morphological consistency of different stratum boundaries in the cross‐ section, and based on the premise of generating one stratum boundary based on the pa‐ rameters, this boundary can be used as the reference stratum boundary to deduce and generate other stratum boundaries. Therefore, the stratum boundary of the middle section can be divided into two types: the reference stratum boundary and the general stratum boundary (Figure 4). The boundary of the outermost stratum was generally selected as the reference stratum boundary. The generation of the stratum boundary mainly includes three parts: the generation of the reference stratum boundary; the deduction of the general

stratum boundary; and the trimming of the stratum boundary.

**Figure 4.** Diagrammatic sketch of stratum boundary.

2.3.1. Generation of Reference Stratum Boundary

*2.3. Construction of Middle Cross‐Section*

codes.

Furthermore, to break through the excessive dependence on geological survey data, this study applies the premise of extracting section location parameters based on geological maps (Figure 3a) and mainly determines the location of cross-sections according to the modeling parameters input by users (Figure 3b). The specific steps are as follows: (1) take the center point D of the middle section, according to the position parameters as the midpoint, and extend, perpendicularly, to the average hinge line strike (hinge line in Figure 3a) to obtain the positions of inflection points A and B of the middle section at d/2; and (2) calculate the affine transformation coefficient and generate the positions of the section lines at both ends, according to the distances (dist<sup>1</sup> and dist2) from the front and back sections to the middle section, the occurrence of the hinge line (including front and back strikes of hinge lines µ<sup>1</sup> and µ<sup>2</sup> and front and back dip angles of hinge lines δ<sup>1</sup> and δ2) and dip angle of the axial surface of the front and back sections γ<sup>1</sup> and γ2. **Figure 3.** Determination of cross‐section locations: (**a**) extraction of location parameters based on a geological map and (**b**) location of sections based on location parameters. S1‐S3 are different stratum codes.

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 7 of 22

#### *2.3. Construction of Middle Cross-Section 2.3. Construction of Middle Cross‐Section*

The construction of the middle cross-section is essentially the process of generating the boundary of each stratum based on the morphological parameters of the cross-section. According to the morphological consistency of different stratum boundaries in the crosssection, and based on the premise of generating one stratum boundary based on the parameters, this boundary can be used as the reference stratum boundary to deduce and generate other stratum boundaries. Therefore, the stratum boundary of the middle section can be divided into two types: the reference stratum boundary and the general stratum boundary (Figure 4). The boundary of the outermost stratum was generally selected as the reference stratum boundary. The generation of the stratum boundary mainly includes three parts: the generation of the reference stratum boundary; the deduction of the general stratum boundary; and the trimming of the stratum boundary. The construction of the middle cross‐section is essentially the process of generating the boundary of each stratum based on the morphological parameters of the cross‐section. According to the morphological consistency of different stratum boundaries in the cross‐ section, and based on the premise of generating one stratum boundary based on the pa‐ rameters, this boundary can be used as the reference stratum boundary to deduce and generate other stratum boundaries. Therefore, the stratum boundary of the middle section can be divided into two types: the reference stratum boundary and the general stratum boundary (Figure 4). The boundary of the outermost stratum was generally selected as the reference stratum boundary. The generation of the stratum boundary mainly includesthree parts: the generation of the reference stratum boundary; the deduction of the general stratum boundary; and the trimming of the stratum boundary.

**Figure 4.** Diagrammatic sketch of stratum boundary. **Figure 4.** Diagrammatic sketch of stratum boundary.

#### 2.3.1. Generation of Reference Stratum Boundary 2.3.1. Generation of Reference Stratum Boundary

The stratum morphology of the fold structure is generally smooth and regular and shows specific function laws, and it can be described in terms of amplitude and wave‐ length [42]. Therefore, it is feasible to use a quadratic Bessel curve to simulate the stratum boundary of the fold structure. In addition, folds are usually formed by the hinge zone and its two connected limbs, and the stratum morphologies of the two limbs tend to be different. Therefore, dividing the boundary line of the fold stratum into two segments according to the hinge zone is appropriate, and the Bezier curve is used to express this The stratum morphology of the fold structure is generally smooth and regular and shows specific function laws, and it can be described in terms of amplitude and wavelength [42]. Therefore, it is feasible to use a quadratic Bessel curve to simulate the stratum boundary of the fold structure. In addition, folds are usually formed by the hinge zone and its two connected limbs, and the stratum morphologies of the two limbs tend to be different. Therefore, dividing the boundary line of the fold stratum into two segments according to the hinge zone is appropriate, and the Bezier curve is used to express this division [46,47]. In addition, the fitting of the Bezier curve was realized through the endpoints and control points. Suppose that only one segment of the Bezier curve is used for the fitting; in this case, obtaining the Bezier control points according to the shape of the interlimb angle is difficult, which causes the morphology of the hinge zone to not be precisely controlled. Therefore, the generation of the reference stratum boundary can be transformed into the selection of the Bezier endpoints and control points of the two curves.

The generation of the reference stratum boundary mainly includes the following steps: (1) calculate the intersection point C of the tangent lines at the inflection points, according to the inflection points (A and B) of the middle section, and the dip angle of the two limbs (α and β) mentioned in the section morphological parameters; (2) find the intersection

point D between the angle bisector of angles ACB and AB, and CD is the axial surface line of the middle section; (3) cut point E proportionally on CD as the hinge point, based on the morphological parameter curve; (4) extend the vertical line of CD through point E and intersect AC and BC to points F and G, respectively; and (5) take points A and E and points E and B as the endpoints, and F and G as the Bezier control points of the left and right curve segments, respectively, and then conduct segmented Bezier interpolation to obtain the reference stratum boundary (Figure 5). tersection point D between the angle bisector of angles ACB and AB, and CD is the axial surface line of the middle section; (3) cut point E proportionally on CD as the hinge point, based on the morphological parameter curve; (4) extend the vertical line of CD through point E and intersect AC and BC to points F and G, respectively; and (5) take points A and E and points E and B as the endpoints, and F and G as the Bezier control points of the left and right curve segments, respectively, and then conduct segmented Bezier interpolation to obtain the reference stratum boundary (Figure 5).

division [46,47]. In addition, the fitting of the Bezier curve was realized through the end‐ points and control points. Suppose that only one segment of the Bezier curve is used for the fitting; in this case, obtaining the Bezier control points according to the shape of the interlimb angle is difficult, which causes the morphology of the hinge zone to not be pre‐ cisely controlled. Therefore, the generation of the reference stratum boundary can be transformed into the selection of the Bezier endpoints and control points of the two curves. The generation of the reference stratum boundary mainly includes the following steps: (1) calculate the intersection point C of the tangent lines at the inflection points, according to the inflection points (A and B) of the middle section, and the dip angle of the two limbs (α and β) mentioned in the section morphological parameters; (2) find the in‐

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 8 of 22

**Figure 5.** Diagrammatic sketch of reference stratum boundary generation. **Figure 5.** Diagrammatic sketch of reference stratum boundary generation.

2.3.2. Deduction of General Stratum Boundary 2.3.2. Deduction of General Stratum Boundary

According to the principle of stratigraphic superposition [48], the bottom surface of each stratum in the fold structure is generally morphologically consistent. Therefore, if the bottom deduction rules of the reference stratum are determined, all the general stra‐ tum boundaries can be obtained. However, the deduction rules for the general stratum boundaries differ for different types of fold structures. Taking the isopach fold and similar folds as examples, the deductive method of the general stratigraphic boundary is briefly According to the principle of stratigraphic superposition [48], the bottom surface of each stratum in the fold structure is generally morphologically consistent. Therefore, if the bottom deduction rules of the reference stratum are determined, all the general stratum boundaries can be obtained. However, the deduction rules for the general stratum boundaries differ for different types of fold structures. Taking the isopach fold and similar folds as examples, the deductive method of the general stratigraphic boundary is briefly described below.

described below. The isopach fold has typical characteristics and the "stratum thickness is equal eve‐ rywhere" [42]. Deducing the general stratum boundaries in such a structure can be real‐ ized by moving a certain distance along the normal direction, according to the stratum thickness based on the points in the reference stratum boundary (Figure 6a). The isopach fold has typical characteristics and the "stratum thickness is equal everywhere" [42]. Deducing the general stratum boundaries in such a structure can be realized by moving a certain distance along the normal direction, according to the stratum thickness based on the points in the reference stratum boundary (Figure 6a). *ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 9 of 22

**Figure 6.** Deduction of general stratum boundaries: (**a**) isopach fold and (**b**) similar fold. **Figure 6.** Deduction of general stratum boundaries: (**a**) isopach fold and (**b**) similar fold.

2.3.3. Trimming of Stratum Boundary As with 2D geological symbols, 3D fold symbols should also be designed to convey geological information as simply as possible and by focusing on the aesthetic expression of symbols. During stratum deduction, the strata of the isopach fold are inferred according to the normal direction of the reference stratum boundary, resulting in uneven stratum boundaries. Therefore, for anticlinal folds with low openness, the bottom of the entire The characteristic of a similar fold is that the curvature radius of each stratum boundary is similar, although a common curvature center is not observed [42]. In other words, the morphology of the stratum boundaries is the same, but their locations are different. Therefore, the general stratum boundaries can be obtained by translating the reference stratum boundary in the vertical direction, and the translation distance can be determined by the stratum thickness at the hinge zone (Figure 6b).

section must be trimmed based on the inflection points of the two limbs in the outermost stratum (Figure 7a,b). For the syncline, the topographic line at the section location was used as a trimming line to trim the top of each stratum to truly express the morphology of the stratum surface (Figure 7c,d). In addition, the core stratum should be calculated

tion according to the inflection points of the middle section; (2) eliminate the points out‐ side the trimming line in the stratum boundary of each stratum; (3) supplement the points on the trimming line in the stratum boundary of each stratum by linear interpolation; and (4) generate the core stratum boundary according to the area enclosed by the trimming

**Figure 7.** Trimming and optimization of the middle section, (**a**) bottom trimming line of the anti‐ cline, (**b**) supplement the core stratum of the anticline after trimming, (**c**) topographic line at the top

The parametric modeling method of the cross‐section described above mainly ap‐ plies to types of folds where no inversion occurs. Thus, the distorted and deformed mor‐

The two limbs of the overturned fold inclined in the same direction, and one of the limbs was overturned. The stratum boundary generation method in Section 2.3.1 can be

of the syncline, and (**d**) supplement the core stratum of the syncline after trimming.

phology of overturned or recumbent folds must be optimized on this basis.

**Legend**

**Legend**

Clip line

Vertex

Vertex

Stratum core boundary Stratum boundary

Topographic line

Stratum core boundary Stratum boundary

line and the bottom stratum (Figure 7c,d).

2.3.4. Section Treatment for Special Fold Types

(1) Overturned fold

#### 2.3.3. Trimming of Stratum Boundary section must be trimmed based on the inflection points of the two limbs in the outermost

2.3.3. Trimming of Stratum Boundary

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 9 of 22

As with 2D geological symbols, 3D fold symbols should also be designed to convey geological information as simply as possible and by focusing on the aesthetic expression of symbols. During stratum deduction, the strata of the isopach fold are inferred according to the normal direction of the reference stratum boundary, resulting in uneven stratum boundaries. Therefore, for anticlinal folds with low openness, the bottom of the entire section must be trimmed based on the inflection points of the two limbs in the outermost stratum (Figure 7a,b). For the syncline, the topographic line at the section location was used as a trimming line to trim the top of each stratum to truly express the morphology of the stratum surface (Figure 7c,d). In addition, the core stratum should be calculated according to the principle of stratum superposition. stratum (Figure 7a,b). For the syncline, the topographic line at the section location was used as a trimming line to trim the top of each stratum to truly express the morphology of the stratum surface (Figure 7c,d). In addition, the core stratum should be calculated according to the principle of stratum superposition. The specific steps of this method are as follows: (1) calculate the trimming line equa‐ tion according to the inflection points of the middle section; (2) eliminate the points out‐ side the trimming line in the stratum boundary of each stratum; (3) supplement the points on the trimming line in the stratum boundary of each stratum by linear interpolation; and (4) generate the core stratum boundary according to the area enclosed by the trimming line and the bottom stratum (Figure 7c,d).

As with 2D geological symbols, 3D fold symbols should also be designed to convey geological information as simply as possible and by focusing on the aesthetic expression of symbols. During stratum deduction, the strata of the isopach fold are inferred according to the normal direction of the reference stratum boundary, resulting in uneven stratum boundaries. Therefore, for anticlinal folds with low openness, the bottom of the entire

**Figure 6.** Deduction of general stratum boundaries: (**a**) isopach fold and (**b**) similar fold.

**Figure 7.** Trimming and optimization of the middle section, (**a**) bottom trimming line of the anti‐ cline, (**b**) supplement the core stratum of the anticline after trimming, (**c**) topographic line at the top of the syncline, and (**d**) supplement the core stratum of the syncline after trimming. **Figure 7.** Trimming and optimization of the middle section, (**a**) bottom trimming line of the anticline, (**b**) supplement the core stratum of the anticline after trimming, (**c**) topographic line at the top of the syncline, and (**d**) supplement the core stratum of the syncline after trimming.

2.3.4. Section Treatment for Special Fold Types The parametric modeling method of the cross‐section described above mainly ap‐ plies to types of folds where no inversion occurs. Thus, the distorted and deformed mor‐ phology of overturned or recumbent folds must be optimized on this basis. (1) Overturned fold The two limbs of the overturned fold inclined in the same direction, and one of the The specific steps of this method are as follows: (1) calculate the trimming line equation according to the inflection points of the middle section; (2) eliminate the points outside the trimming line in the stratum boundary of each stratum; (3) supplement the points on the trimming line in the stratum boundary of each stratum by linear interpolation; and (4) generate the core stratum boundary according to the area enclosed by the trimming line and the bottom stratum (Figure 7c,d).

#### limbs was overturned. The stratum boundary generation method in Section 2.3.1 can be 2.3.4. Section Treatment for Special Fold Types

The parametric modeling method of the cross-section described above mainly applies to types of folds where no inversion occurs. Thus, the distorted and deformed morphology of overturned or recumbent folds must be optimized on this basis.

### (1) Overturned fold

The two limbs of the overturned fold inclined in the same direction, and one of the limbs was overturned. The stratum boundary generation method in Section 2.3.1 can be used to construct the stratum boundaries by treating the dip angle at the overturned limb as an obtuse angle (Figure 8a). However, owing to the compression and deformation between strata, most of the strata of the overturned folds will be distorted to varying degrees [49]. Therefore, when constructing the middle cross-section, the rotation angle of the section point set around vertex D can be calculated according to Formula (1). Distortion and deformation processing are then carried out to realize an approximate simulation of the overturned fold morphology.

$$\boldsymbol{\Theta} = \mathbf{a} \mathbf{m} \mathbf{a} \times (\mathbf{R} - \mathbf{r}) / \mathbf{R} \tag{1}$$

where amax is the maximum rotation angle, R is the maximum rotation radius, r is the distance from the current point to the rotation center D, and θ is the angle of rotation required for the current point.

used to construct the stratum boundaries by treating the dip angle at the overturned limb as an obtuse angle (Figure 8a). However, owing to the compression and deformation be‐ tween strata, most of the strata of the overturned folds will be distorted to varying degrees [49]. Therefore, when constructing the middle cross‐section, the rotation angle of the sec‐ tion point set around vertex D can be calculated according to Formula (1). Distortion and deformation processing are then carried out to realize an approximate simulation of the

used to construct the stratum boundaries by treating the dip angle at the overturned limb as an obtuse angle (Figure 8a). However, owing to the compression and deformation be‐ tween strata, most of the strata of the overturned folds will be distorted to varying degrees [49]. Therefore,when constructing the middle cross‐section, the rotation angle of the sec‐ tion point set around vertex D can be calculated according to Formula (1). Distortion and deformation processing are then carried out to realize an approximate simulation of the

where amax is the maximum rotation angle, R is the maximum rotation radius, r is the distance from the current point to the rotation center D, and θ is the angle of rotation

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 10 of 22

where amax is the maximum rotation angle, R is the maximum rotation radius, r is the distance from the current point to the rotation center D, and θ is the angle of rotation

θ = amax × (R − r)/R (1)

θ = amax × (R − r)/R (1)

**Figure 8.** Distortion treatment for overturned folds: (**a**) before distortion and (**b**) after distortion. **Figure 8.** Distortion treatment for overturned folds: (**a**) before distortion and (**b**) after distortion. **Figure 8.** Distortion treatment for overturned folds: (**a**) before distortion and (**b**) after distortion.

#### (2) Continuous fold belt (2) Continuous fold belt (2) Continuous fold belt

overturned fold morphology.

overturned fold morphology.

required for the current point.

required for the current point.

In general, typical geological structures with prominent characteristics are often ac‐ companied by large‐scale folding activities. One single symbol model of a fold cannot easily reflect the structural characteristics of a continuous fold belt. To enable fold symbols to effectively realize the abstract expression of geological characteristics under different geological scales, this study uses the periodic cycle method to repeat the middle section construction steps presented in Section 2.3 to generate the cross‐section of continuous fold belts (Figure 9). In general, typical geological structures with prominent characteristics are often accompanied by large-scale folding activities. One single symbol model of a fold cannot easily reflect the structural characteristics of a continuous fold belt. To enable fold symbols to effectively realize the abstract expression of geological characteristics under different geological scales, this study uses the periodic cycle method to repeat the middle section construction steps presented in Section 2.3 to generate the cross-section of continuous fold belts (Figure 9). In general, typical geological structures with prominent characteristics are often ac‐ companied by large‐scale folding activities. One single symbol model of a fold cannot easily reflect the structural characteristics of a continuous fold belt. To enable fold symbols to effectively realize the abstract expression of geological characteristics under different geological scales, this study uses the periodic cycle method to repeat the middle section construction steps presented in Section 2.3 to generate the cross‐section of continuous fold belts (Figure 9).

**Figure 9.** Periodic treatment for continuous fold belts. **Figure 9.** Periodic treatment for continuous fold belts.

#### *2.4. Generation of Cross‐Sections at Both Ends 2.4. Generation of Cross-Sections at Both Ends*

Based on the morphological consistency among cross‐sections of the fold structure, the cross‐sections at both ends can be obtained by affine transformation from the middle cross‐section. The affine transformation in this section includes translation, magnification, and rotation. The strike of the hinge affects the offset distance of the sections at both ends along the X‐axis, the dip angle of the hinge affects the magnification ratio of the sections at both ends, and the dip angle of the axial surfaces of the sections at both ends affects the *2.4. Generation of Cross‐Sections at Both Ends* Based on the morphological consistency among cross‐sections of the fold structure, the cross‐sections at both ends can be obtained by affine transformation from the middle cross‐section. The affine transformation in this section includes translation, magnification, and rotation. The strike of the hinge affects the offset distance of the sections at both ends along the X‐axis, the dip angle of the hinge affects the magnification ratio of the sections at both ends, and the dip angle of the axial surfaces of the sections at both ends affects the Based on the morphological consistency among cross-sections of the fold structure, the cross-sections at both ends can be obtained by affine transformation from the middle cross-section. The affine transformation in this section includes translation, magnification, and rotation. The strike of the hinge affects the offset distance of the sections at both ends along the X-axis, the dip angle of the hinge affects the magnification ratio of the sections at both ends, and the dip angle of the axial surfaces of the sections at both ends affects the rotation angle. The specific affine transformation parameters and corresponding calculation methods are listed in Table 3.



Note: Please refer to Figure 3b for the relevant parameters.

### *2.5. Generation of Transition Sections*

The cross-sections, including the front, middle, and back, were obtained. Theoretically, these cross-sections can be stitched directly to complete the rough construction of the entire

Offset distance along

Magnification scale ൌ ሺℎ tan ൈ ሻ/ℎ

<sup>X</sup> axis ൌ tan ൈ

symbol.

Rotation angle ൌ ′ െ

fold model. However, the spatial morphology of the hinge line is difficult to express using this method, which affects the accuracy and aesthetics of the model results to a large extent. Therefore, many transition sections need to be generated through morphing interpolation according to these cross-sections to realize the fine modeling of the fold symbol. Morphing is an interpolation technique that smoothly transitions an initial object into a target object. This technology has been widely used in animation production, virtual reality, image compression, and image reconstruction [50–52]. Three types of boundaries

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 11 of 22

the middle section.

the middle section.

Note: Please refer to Figure 3b for the relevant parameters.

tion methods are listed in Table 3.

**Table 3.** Affine transformation factor.

*2.5. Generation of Transition Sections*

**Factor Name Computing Formula Variable Description**

rotation angle. The specific affine transformation parameters and corresponding calcula‐

δ is the front or back dip angle of the hinge in the location

to the lower vertex D in the middle section (Figure 5).

µ is the front or back strike of the hinge in the location

The cross‐sections, including the front, middle, and back, were obtained. Theoreti‐

cally, these cross‐sections can be stitched directly to complete the rough construction of the entire fold model. However, the spatial morphology of the hinge line is difficult to express using this method, which affects the accuracy and aesthetics of the model results to a large extent. Therefore, many transition sections need to be generated through morph‐ ing interpolation according to these cross‐sections to realize the fine modeling of the fold

parameter, dist is the distance from the current section to the

parameters, and dist is the distance from the current section to

γ' is the front or back section dip angle of the axial surface in

location parameters, and γ is the dip angle of the axial surface in

middle section, and h is the height from the actual upper vertex E

Morphing is an interpolation technique that smoothly transitions an initial object into a target object. This technology has been widely used in animation production, virtual reality, image compression, and image reconstruction [50–52]. Three types of boundaries are necessary for morphing interpolation (Figure 10): (1) the starting boundary with direction (SRC); (2) the target boundary (DEST) with a direction that must be consistent with that of SRC; and (3) the first and last constraint boundaries with directions (FC and LC) [53]. are necessary for morphing interpolation (Figure 10): (1) the starting boundary with di‐ rection (SRC); (2) the target boundary (DEST) with a direction that must be consistent with that of SRC; and (3) the first and last constraint boundaries with directions (FC and LC) [53].

**Figure 10.** Principle of morphing interpolation. **Figure 10.** Principle of morphing interpolation.

2.5.1. Selection of Constraint Boundaries

2.5.1. Selection of Constraint Boundaries Considering that the monotony of a single fold element in a folded section will expe‐ rience one change at the hinge zone [54], this study considers the hinge line as the middle constraint and divides the strata of the fold sections into left and right parts. Morphing Considering that the monotony of a single fold element in a folded section will experience one change at the hinge zone [54], this study considers the hinge line as the middle constraint and divides the strata of the fold sections into left and right parts. Morphing interpolation with four boundaries intersected by the hinge line was performed for the two limbs. This method allows the hinge information to participate effectively in the morphing interpolation process, thus ensuring the accuracy of the overall interpolation result.

interpolation with four boundaries intersected by the hinge line was performed for the two limbs. This method allows the hinge information to participate effectively in the morphing interpolation process, thus ensuring the accuracy of the overall interpolation result. Among the four constraint boundaries of the morphing interpolation method, the Among the four constraint boundaries of the morphing interpolation method, the starting and ending stratum boundaries, SRC and DEST, can be obtained from the stratum boundary on the corresponding limb [46]. To make the section transition smoother, the constrained stratum boundaries, FC and LC, must be obtained by interpolation according to the corresponding characteristic points of each section. By connecting the points generated by the Lagrange interpolation results of A, E, and B (i.e., left inflection point, hinge point, and right inflection point, respectively), three smooth curves (AfrontAmidAback, EfrontEmidEback, and BfrontBmidBback) of each stratum surface can be obtained (Figure 11a).

starting and ending stratum boundaries, SRC and DEST, can be obtained from the stratum 2.5.2. Transition Section Generation Based on Morphing Interpolation

Calculate any point *M* (*u*, *v*) inserted in the transition section according to the morphing formula (2), and then obtain the point set of the transition section.

$$\begin{cases} M(\mu, v) = H(\mu, v) + L(\mu, v) \\ H(\mu, v) = (1 - v)F(\mu) + vG(\mu) \\ L(\mu, v) = (1 - \mu)a + \mu b \\ a = P(v) - H(0, v) \\ b = Q(v) - H(1, v) \end{cases} \tag{2}$$

where *M* (*u*, *v*) represents the point to be inserted, and it is jointly controlled by the transverse path interpolation function *H* (*u*, *v*) and longitudinal constraint boundary function *L* (*u*, *v*). 0 ≤ *u* ≤ 1, 0 ≤ *v* ≤ 1. Moreover, as shown in Figure 10, *u* represents the close degree

between the point to be inserted and the FC or LC constraint lines and *v* represents the close degree between the point to be inserted and SRC or DEST lines. In addition, *F*(*u*) is the curve function of the SRC, *G*(*u*) is the curve function of the DEST, *P*(*v*) is a function of FC and *Q*(*v*) is a function of LC, *a* is the adjustment factor corresponding to the FC, and *b* is the adjustment factor corresponding to the LC. between the point to be inserted and the FC or LC constraint lines and *v* represents the close degree between the point to be inserted and SRC or DEST lines. In addition, *F(u)* is the curve function of the SRC, *G(u)* is the curve function of the DEST, *P(v)* is a function of FC and *Q(v)* is a function of LC, *a* is the adjustment factor corresponding to the FC, and *b* is the adjustment factor corresponding to the LC.

where *M (u, v)* represents the point to be inserted, and it is jointly controlled by the trans‐ verse path interpolation function *H (u, v)* and longitudinal constraint boundary function *L (u, v)*. 0 ≤ *u* ≤ 1, 0 ≤ *v* ≤ 1. Moreover, as shown in Figure 10, *u* represents the close degree

boundary on the corresponding limb [46]. To make the section transition smoother, the constrained stratum boundaries, FC and LC, must be obtained by interpolation according to the corresponding characteristic points of each section. By connecting the points gener‐ ated by the Lagrange interpolation results of A, E, and B (i.e., left inflection point, hinge point, and right inflection point, respectively), three smooth curves (AfrontAmidAback, Efron‐

Calculate any point *M (u, v)* inserted in the transition section according to the morph‐

(2)

tEmidEback, and BfrontBmidBback) of each stratum surface can be obtained (Figure 11a).

2.5.2. Transition Section Generation Based on Morphing Interpolation

ing formula (2), and then obtain the point set of the transition section.

ሺ, ሻ ൌ ሺ, ሻ ሺ, ሻ ሺ, ሻ ൌ ሺ1 െ ሻሺሻ ሺሻ ሺ, ሻ ൌ ሺ1 െ ሻ ൌ ሺሻ െ ሺ0, ሻ ൌ ሺሻ െ ሺ1, ሻ

⎩ ⎪ ⎨ ⎪ ⎧

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 12 of 22

**Figure 11.** Interpolation of transition sections based on morphing, (**a**) constraint boundaries for morphing interpolation, and (**b**) point set of transition sections. The point sets of two colors are the results of morphing interpolation, respectively. **Figure 11.** Interpolation of transition sections based on morphing, (**a**) constraint boundaries for morphing interpolation, and (**b**) point set of transition sections. The point sets of two colors are the results of morphing interpolation, respectively.

#### *2.6. Fold Model Construction and Attribute Assignment 2.6. Fold Model Construction and Attribute Assignment*

After the construction of all of the transition sections is completed, the boundary rep‐ resentation (Brep) surface model of the fold symbol can be built based on the contour comparison algorithm [55]. The difficulty of the contour algorithm lies in the treatment of branching, correspondence, meshing, and smoothing problems [56]. For the method used in this study, because the point sets between adjacent sections are all in one‐to‐one mode, the branching problem and corresponding problem are not observed. The smoothing problem can also be solved by generating transition sections. After the construction of all of the transition sections is completed, the boundary representation (Brep) surface model of the fold symbol can be built based on the contour comparison algorithm [55]. The difficulty of the contour algorithm lies in the treatment of branching, correspondence, meshing, and smoothing problems [56]. For the method used in this study, because the point sets between adjacent sections are all in one-to-one mode, the branching problem and corresponding problem are not observed. The smoothing problem can also be solved by generating transition sections. *ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 13 of 22

Considering the facet complexity and algorithm efficiency of the modeling result, the shortest diagonal method is selected [57] (Figure 12b). This is a locally optimal objective Considering the facet complexity and algorithm efficiency of the modeling result, the shortest diagonal method is selected [57] (Figure 12b). This is a locally optimal objective function method. The fold symbol model was finally obtained based on the contour comparison algorithm, as shown in Figure 12c. function method. The fold symbol model was finally obtained based on the contour com‐ parison algorithm, as shown in Figure 12c.

**Figure 12.** Fold surface model constructed based on the contour algorithm, (**a**) adjacent sections after morphing, (**b**) surface model constructed using the contour algorithm, and (**c**) fold symbol model result. Entities with different colors represent different strata. **Figure 12.** Fold surface model constructed based on the contour algorithm, (**a**) adjacent sections after morphing, (**b**) surface model constructed using the contour algorithm, and (**c**) fold symbol model result. Entities with different colors represent different strata.

After the model stitching is completed, attribute information, such as age, the litho‐ logical description of each stratum, and contact relationship between strata, can be stored in the corresponding JSON format model attribute file. The association query between the 3D model and the attribute information, 3D information annotation, and other functions can be effectively performed by binding the ID information of each object. After the model stitching is completed, attribute information, such as age, the lithological description of each stratum, and contact relationship between strata, can be stored in the corresponding JSON format model attribute file. The association query between the 3D model and the attribute information, 3D information annotation, and other functions can be effectively performed by binding the ID information of each object.

puter with the configuration of a 3.20 GHz Intel(R) Core (TM) i7‐8700 CPU, 16 GB RAM. The operation system is Windows 10 Professional. All algorithms are implemented using GDAL3.2 and Dotspatial 1.8 and compiled using Microsoft visual C# 2017 compiler. The number of control and transition sections in all cases is three and eight, respectively. In‐ formation for each fold is presented in Table 4. The parameters of the different cases are

**ID Name Fold Type Source** Dayue Mountain in Mount Lu Vertical horizontal anticline [58] Danaobo Mountain in Mount Lu Oblique horizontal anticline [59] Fangdou Mountain in East Sichuan Chevron fold [60] Scottish Highlands Recumbent fold [61]

Center point of middle section coordinate (1000, 0, 0) (1000, 0, 0) (1000, 0, 0) (1000, 0, 0)

Front and back dip angles of hinge line ° −10, −10 −5, −5 0, 0 0, 0 Front and back strikes of hinge line ° 0, 0 −10, −5 0, 0 0, 0 Axial dip angles of front and back sections ° 80, 80 70, 70 95, 95 20, 20

Dip angles of axial surface ° 80 70 95 20 Dip angles of inflection points ° 55, 35 40 120 150, 10 Morphological parameter 0.8 0.4 0.95 0.6 Collection of strata thickness M 50, 50, 50 50, 50, 50 30, 30, 30 40, 40, 40 Maximum rotation angle of distortion ° 0 0 0 30

middle section <sup>m</sup> 2000, <sup>2000</sup> 2000, <sup>2000</sup> 6000, <sup>6000</sup> 2000, <sup>2000</sup> Fold width m 2000 2000 2000 1000

**Mountain**

**Fangdou MOUNTAIN**

**Scottish Highlands**

**3. Case Studies**

summarized in Table 5.

Distance from front and back sections to

Location parameters

Morphological parameters

**Table 4.** Details of each experimental area.

**Type Name Unit Dayue Mountain Danaobo**

*3.1. Case 1: Anticline Modeling*

**Table 5.** Modeling parameters of each experimental area.

### **3. Case Studies**

This study selects four folds with different types for testing the performance of the 3D modeling of fold-structure symbols. The experiments were carried out using a computer with the configuration of a 3.20 GHz Intel(R) Core (TM) i7-8700 CPU, 16 GB RAM. The operation system is Windows 10 Professional. All algorithms are implemented using GDAL3.2 and Dotspatial 1.8 and compiled using Microsoft visual C# 2017 compiler. The number of control and transition sections in all cases is three and eight, respectively. Information for each fold is presented in Table 4. The parameters of the different cases are summarized in Table 5.

**Table 4.** Details of each experimental area.


### **Table 5.** Modeling parameters of each experimental area.


### *3.1. Case 1: Anticline Modeling*

Dayue Mountain is located southeast of Mount Lu in Jiangxi Province and represents the second highest peak in the Mount Lu region. The measured data from geological records of Mount Lu [58] indicate that the 3D symbol model of Dayue Mountain (Figure 13b), generated based on the above parameters, reflects the structural characteristics of Dayue Mountain, such as equal stratum thickness, vertical anticline, and rounded hinge zone. The number of vertices and triangles is 16,755 and 28,418, respectively.

The number of vertices and triangles is 16,755 and 28,418, respectively.

The number of vertices and triangles is 16,755 and 28,418, respectively.

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 14 of 22

**Figure 13.** 3D symbol modeling process of Dayue Mountain: (**a**) sections and hinge line of Dayue Mountain and (**b**) 3D symbol model of Dayue Mountain upright anticline. **Figure 13.** 3D symbol modeling process of Dayue Mountain: (**a**) sections and hinge line of Dayue Mountain and (**b**) 3D symbol model of Dayue Mountain upright anticline. **Figure 13.** 3D symbol modeling process of Dayue Mountain: (**a**) sections and hinge line of Dayue Mountain and (**b**) 3D symbol model of Dayue Mountain upright anticline.

Dayue Mountain is located southeast of Mount Lu in Jiangxi Province and represents the second highest peak in the Mount Lu region. The measured data from geological rec‐ ords of Mount Lu [58] indicate that the 3D symbol model of Dayue Mountain (Figure 13b), generated based on the above parameters, reflects the structural characteristics of Dayue Mountain, such as equal stratum thickness, vertical anticline, and rounded hinge zone.

Dayue Mountain is located southeast of Mount Lu in Jiangxi Province and represents the second highest peak in the Mount Lu region. The measured data from geological rec‐ ords of Mount Lu [58] indicate that the 3D symbol model of Dayue Mountain (Figure 13b), generated based on the above parameters, reflects the structural characteristics of Dayue Mountain, such as equal stratum thickness, vertical anticline, and rounded hinge zone.

#### *3.2. Case 2: Modeling of Syncline 3.2. Case 2: Modeling of Syncline*

*3.2. Case 2: Modeling of Syncline* Danaobo Mountain, as an oblique horizontal fold, is located in the Xifeng Cave and Danaobo mountain area at the northeast corner of the survey area (Mahui Ridge sheet (H‐ 50‐88‐D)). The axial trace is in the northeast direction, and the axis surface inclines to the southeast at a dip angle of approximately 70°. The hinge was nearly horizontal, the fold length was approximately 9 km, and the width was approximately 2.5 km. Among the strata involved, the core is the second and fourth member of the Hanyangfeng Formation and the two limbs are the first and second members of the Hanyangfeng Formation. The generated 3D symbol model of Danaobo Mountain is shown in Figure 14. The number of Danaobo Mountain, as an oblique horizontal fold, is located in the Xifeng Cave and Danaobo mountain area at the northeast corner of the survey area (Mahui Ridge sheet (H-50-88-D)). The axial trace is in the northeast direction, and the axis surface inclines to the southeast at a dip angle of approximately 70◦ . The hinge was nearly horizontal, the fold length was approximately 9 km, and the width was approximately 2.5 km. Among the strata involved, the core is the second and fourth member of the Hanyangfeng Formation and the two limbs are the first and second members of the Hanyangfeng Formation. The generated 3D symbol model of Danaobo Mountain is shown in Figure 14. The number of vertices and triangles is 14,539 and 29,062, respectively. Danaobo Mountain, as an oblique horizontal fold, is located in the Xifeng Cave and Danaobo mountain area at the northeast corner of the survey area (Mahui Ridge sheet (H‐ 50‐88‐D)). The axial trace is in the northeast direction, and the axis surface inclines to the southeast at a dip angle of approximately 70°. The hinge was nearly horizontal, the fold length was approximately 9 km, and the width was approximately 2.5 km. Among the strata involved, the core is the second and fourth member of the Hanyangfeng Formation and the two limbs are the first and second members of the Hanyangfeng Formation. The generated 3D symbol model of Danaobo Mountain is shown in Figure 14. The number of vertices and triangles is 14,539 and 29,062, respectively.

vertices and triangles is 14,539 and 29,062, respectively.

**Figure 14.** 3D symbol modeling process of the Danaobo mountain: (**a**) sections and hinge line of Danaobo mountain, (**b**) 3D symbol model of the Danaobo mountain upright syncline before trim‐ ming, and (**c**) 3D symbol model of the Danaobo mountain upright syncline after trimming. **Figure 14.** 3D symbol modeling process of the Danaobo mountain: (**a**) sections and hinge line of Danaobo mountain, (**b**) 3D symbol model of the Danaobo mountain upright syncline before trim‐ ming, and (**c**) 3D symbol model of the Danaobo mountain uprightsyncline after trimming. **Figure 14.** 3D symbol modeling process of the Danaobo mountain: (**a**) sections and hinge line of Danaobo mountain, (**b**) 3D symbol model of the Danaobo mountain upright syncline before trimming, and (**c**) 3D symbol model of the Danaobo mountain upright syncline after trimming.

#### *3.3. Case 3: Modeling of Chevron Anticline 3.3. Case 3: Modeling of Chevron Anticline*

*3.3. Case 3: Modeling of Chevron Anticline* The Fangdou Mountain chevron anticline is located in the East Sichuan Fold Belt. The mountain is a narrow, high, and steep anticline in the Jura‐type folds of eastern Sichuan. The Fangdou Mountain structure is a composite anticline composed of fault‐bend folds and longitudinal‐bend folds controlled by detachment, and its surface section is an asym‐ metric fold with a sharp hinge zone, which gradually widens and slows to form an asym‐ metric box fold [62]. According to the 95‐23.5 survey line section data of Fangdou Moun‐ tain [60], the limb dip angles of the anticline are approximately 20°–25° and symmetrical and the axial surface is nearly vertical. Based on these parameters (Table 5), a 3D symbol model of Fangdou Mountain (Figure 15) is generated, which reflects the structural char‐ acteristics of its anticline, sharp edge, upright, and symmetrical characteristics. The num‐ The Fangdou Mountain chevron anticline is located in the East Sichuan Fold Belt. The mountain is a narrow, high, and steep anticline in the Jura‐type folds of eastern Sichuan. The Fangdou Mountain structure is a composite anticline composed of fault‐bend folds and longitudinal‐bend folds controlled by detachment, and its surface section is an asym‐ metric fold with a sharp hinge zone, which gradually widens and slows to form an asym‐ metric box fold [62]. According to the 95‐23.5 survey line section data of Fangdou Moun‐ tain [60], the limb dip angles of the anticline are approximately 20°–25° and symmetrical and the axial surface is nearly vertical. Based on these parameters (Table 5), a 3D symbol model of Fangdou Mountain (Figure 15) is generated, which reflects the structural char‐ acteristics of its anticline, sharp edge, upright, and symmetrical characteristics. The num‐ ber of vertices and triangles is 12,843 and 24,748, respectively. The Fangdou Mountain chevron anticline is located in the East Sichuan Fold Belt. The mountain is a narrow, high, and steep anticline in the Jura-type folds of eastern Sichuan. The Fangdou Mountain structure is a composite anticline composed of fault-bend folds and longitudinal-bend folds controlled by detachment, and its surface section is an asymmetric fold with a sharp hinge zone, which gradually widens and slows to form an asymmetric box fold [62]. According to the 95-23.5 survey line section data of Fangdou Mountain [60], the limb dip angles of the anticline are approximately 20◦–25◦ and symmetrical and the axial surface is nearly vertical. Based on these parameters (Table 5), a 3D symbol model of Fangdou Mountain (Figure15) is generated, which reflects the structural characteristics of its anticline, sharp edge, upright, and symmetrical characteristics. The number of vertices and triangles is 12,843 and 24,748, respectively.

ber of vertices and triangles is 12,843 and 24,748, respectively.

**Figure 15.** 3D symbol modeling process of Fangdou Mountain: (**a**) sections and hinge line of Fang‐ dou Mountain, and (**b**) 3D symbol model of the Fangdou Mountain chevron fold. **Figure 15.** 3D symbol modeling process of Fangdou Mountain: (**a**) sections and hinge line of Fangdou Mountain, and (**b**) 3D symbol model of the Fangdou Mountain chevron fold. **Figure 15.** 3D symbol modeling process of Fangdou Mountain: (**a**) sections and hinge line of Fang‐ dou Mountain, and (**b**) 3D symbol model of the Fangdou Mountain chevron fold.

#### *3.4. Case 4: Modeling of Recumbent Folds 3.4. Case 4: Modeling of Recumbent Folds*

*3.4. Case 4: Modeling of Recumbent Folds* A large number of overturned folds are observed in the Scottish Highlands [61]. This study selects a recumbent fold to demonstrate symbol modeling in the central south zone of the area. The approximate fold information can be obtained from the geological section map (Figure 16a): the dip angle of the left limb was approximately 150° and the right is at A large number of overturned folds are observed in the Scottish Highlands [61]. This study selects a recumbent fold to demonstrate symbol modeling in the central south zone of the area. The approximate fold information can be obtained from the geological section map (Figure 16a): the dip angle of the left limb was approximately 150◦ and the right is at approximately 10◦ . The generated 3D symbol model is shown in Figure 16b. The number of vertices and triangles is 7020 and 14,024, respectively. A large number of overturned folds are observed in the Scottish Highlands [61]. This study selects a recumbent fold to demonstrate symbol modeling in the central south zone of the area. The approximate fold information can be obtained from the geological section map (Figure 16a): the dip angle of the left limb was approximately 150° and the right is at approximately 10°. The generated 3D symbol model is shown in Figure 16b. The number of vertices and triangles is 7020 and 14,024, respectively.

**Figure 16.** Section map and modeling result of recumbent fold in the Scottish Highlands: (**a**) geolog‐ ical section map of the recumbent fold and (**b**) symbol model of the recumbent fold. Entities with different colors represent different strata. **Figure 16.** Section map and modeling result of recumbent fold in the Scottish Highlands: (**a**) geological section map of the recumbent fold and (**b**) symbol model of the recumbent fold. Entities with different colors represent different strata.

#### *3.5. Case 5: Other Types of Folds 3.5. Case 5: Other Types of Folds*

ure 17).

ure 17).

*3.5. Case 5: Other Types of Folds*

**Figure 16.** Section map and modeling result of recumbent fold in the Scottish Highlands: (**a**) geolog‐ ical section map of the recumbent fold and (**b**) symbol model of the recumbent fold. Entities with different colors represent different strata. In addition to the fold structures in the specific areas mentioned above, this method can also be applied to a more general fold symbol modeling based on the inclusion of the appropriate modeling parameters. When the morphological parameter cur is set to 0.8, 0.5, and 0.2, the symbol models of the chevron fold, circular fold, and box fold can be generated. When the interlimb angle θ is set to 30°, 50°, 100°, and 160°, the symbol models of a tight fold, closed fold, broad fold, and gentle fold can be generated, respectively (Fig‐ In addition to the fold structures in the specific areas mentioned above, this method can also be applied to a more general fold symbol modeling based on the inclusion of the appropriate modeling parameters. When the morphological parameter cur is set to 0.8, 0.5, and 0.2, the symbol models of the chevron fold, circular fold, and box fold can be generated. When the interlimb angle θ is set to 30◦ , 50◦ , 100◦ , and 160◦ , the symbol models of a tight fold, closed fold, broad fold, and gentle fold can be generated, respectively (Figure 17).

In addition to the fold structures in the specific areas mentioned above, this method

can also be applied to a more general fold symbol modeling based on the inclusion of the appropriate modeling parameters. When the morphological parameter cur is set to 0.8,

generated. When the interlimb angle θ is set to 30°, 50°, 100°, and 160°, the symbol models of a tight fold, closed fold, broad fold, and gentle fold can be generated, respectively (Fig‐

**Figure 17.** Other types of fold symbols obtained by setting the cur and θ parameters to different values, (**a**) cur = 0.8, chevron fold, (**b**) cur = 0.5, circular fold, (**c**) cur = 0.2, box fold, (**d**) θ = 30°, tight fold, (**e**) θ = 50°, closed fold, (**f**) θ = 100°, broad fold, and (**g**) θ = 160°, gentle fold. For a description of parameter cur and θ, please refer to Table 2. Entities with different colors represent different strata. **Figure 17.** Other types of fold symbols obtained by setting the cur and θ parameters to different values, (**a**) cur = 0.8, chevron fold, (**b**) cur = 0.5, circular fold, (**c**) cur = 0.2, box fold, (**d**) θ = 30◦ , tight fold, (**e**) θ = 50◦ , closed fold, (**f**) θ = 100◦ , broad fold, and (**g**) θ = 160◦ , gentle fold. For a description of parameter cur and θ, please refer to Table 2. Entities with different colors represent different strata.

### **4. Discussion**

#### **4. Discussion** *4.1. Analysis of the Influencing Factors*

*4.1. Analysis of the Influencing Factors* The use of a greater number of control sections increases the accuracy of the gener‐ ated model and ensures that the model is able to better restore the morphological charac‐ teristics of the fold structure. For fold structures with small‐scale and prominent charac‐ teristics, fewer control sections are required to reflect the hinge morphology. In contrast, for fold types with large‐scale and relatively complex hinge morphologies, more control sections are required. The typical fold structure selected in this experiment has a small The use of a greater number of control sections increases the accuracy of the generated model and ensures that the model is able to better restore the morphological characteristics of the fold structure. For fold structures with small-scale and prominent characteristics, fewer control sections are required to reflect the hinge morphology. In contrast, for fold types with large-scale and relatively complex hinge morphologies, more control sections are required. The typical fold structure selected in this experiment has a small scale; therefore, only three sections in the middle and at both ends were generated.

scale; therefore, only three sections in the middle and at both ends were generated. The use of a greater number of transition sections leads to the generation of a smoother model and requires a greater amount of model data. However, the running ef‐ ficiency tends to be lower when a model with a large amount of data displays a large scene. For the fold type with a relatively smooth hinge zone, fewer transition sections are required to effectively reflect the 3D morphology of the fold structure. Otherwise, more The use of a greater number of transition sections leads to the generation of a smoother model and requires a greater amount of model data. However, the running efficiency tends to be lower when a model with a large amount of data displays a large scene. For the fold type with a relatively smooth hinge zone, fewer transition sections are required to effectively reflect the 3D morphology of the fold structure. Otherwise, more transition sections would be required to effectively reflect the 3D morphology of the fold structure.

transition sections would be required to effectively reflect the 3D morphology of the fold structure. The different selection methods of the location parameters of the cross‐sections sig‐ nificantly influence the model morphology. The location parameters of the cross‐sections can be measured based on the outer rectangle of the fold distribution range or the outer The different selection methods of the location parameters of the cross-sections significantly influence the model morphology. The location parameters of the cross-sections can be measured based on the outer rectangle of the fold distribution range or the outer rectangle of the core formation. These methods can be applied as long as they can effectively reflect the morphological characteristics of the fold structures.

rectangle of the core formation. These methods can be applied as long as they can effec‐ tively reflect the morphological characteristics of the fold structures. The quality of the middle cross‐section is a critical factor for determining the model‐ ing quality of the fold symbol. The parametric modeling method is effective for the middle cross‐section with a relatively simple stratum boundary morphology. However, the effect of parametric modeling is poor for a more complex middle cross‐section. In this case, a The quality of the middle cross-section is a critical factor for determining the modeling quality of the fold symbol. The parametric modeling method is effective for the middle cross-section with a relatively simple stratum boundary morphology. However, the effect of parametric modeling is poor for a more complex middle cross-section. In this case, a better modeling quality can be obtained by replacing the cross-sections generated by the parameters with the measured cross-section.

#### better modeling quality can be obtained by replacing the cross‐sections generated by the *4.2. Factors That Influence the Symbol Display Effect*

parameters with the measured cross‐section. *4.2. Factors That Influence the Symbol Display Effect* The spatial information of the constructed models is stored in an OBJ format file, while the attribute information is in JSON. The constructed 3D fold symbols can be placed on a 3D Earth display platform (such as Cesium or Google Earth). When importing a model into Cesium, the process mainly includes: (1) converting the OBJ format model file to FBX format file, which contains information such as material and lighting; (2) The spatial information of the constructed models is stored in an OBJ format file, while the attribute information is in JSON. The constructed 3D fold symbols can be placed on a 3D Earth display platform (such as Cesium or Google Earth). When importing a model into Cesium, the process mainly includes: (1) converting the OBJ format model file to FBX format file, which contains information such as material and lighting; (2) converting the FBX format model files to GLB format file, which is the 3D model data format supported by the Cesium platform; (3) directly loading the glb file into the Cesium scene.

scene.

Consistent with the expression requirements of 2D symbols [63], when the constructed 3D folding symbols are placed on the 3D earth display platform, it is also necessary to set the appropriate symbol size to optimize the display effect. When the scale of the fold structure is small and the internal characteristics of the fold elements are apparent, the effect is better when the ratio of the symbol to the actual fold is 1:1. When the scale of the fold structure is large, the internal characteristics of the fold elements will not be sufficiently prominent, and the effect of the proportional magnification will be better. In addition, in 2D electronic maps, the size of related symbols usually needs to be adaptively adjusted at different scales. This rule also applies to symbol-size settings in 3D scenes. For example, when displayed on the Cesium platform, the size of the corresponding model can be changed for adaptive scaling and adjustment to ensure that the relative size of the symbol model remains unchanged at different scales (Figure 18). Consistent with the expression requirements of 2D symbols [63], when the con‐ structed 3D folding symbols are placed on the 3D earth display platform, it is also neces‐ sary to set the appropriate symbol size to optimize the display effect. When the scale of the fold structure is small and the internal characteristics of the fold elements are appar‐ ent, the effect is better when the ratio of the symbol to the actual fold is 1:1. When the scale of the fold structure is large, the internal characteristics of the fold elements will not be sufficiently prominent, and the effect of the proportional magnification will be better. In addition, in 2D electronic maps, the size of related symbols usually needs to be adaptively adjusted at different scales. This rule also applies to symbol‐size settings in 3D scenes. For example, when displayed on the Cesium platform, the size of the corresponding model can be changed for adaptive scaling and adjustment to ensure that the relative size of the symbol model remains unchanged at different scales (Figure 18).

converting the FBX format model files to GLB format file, which is the 3D model data format supported by the Cesium platform; (3) directly loading the glb file into the Cesium

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 17 of 22

**Figure 18.** Display effect of fold symbols with different scaling ratios: (**a**) scaling down and (**b**) equal ratio. **Figure 18.** Display effect of fold symbols with different scaling ratios: (**a**) scaling down and (**b**) equal ratio.

When the symbol is smaller than the actual geological structure, it should be placed directly above the center point of the actual geological structure. When the scale is 1:1, the symbol should be placed in the actual position of the structure and the height should be consistent with the actual structure. In addition, when symbols are placed, the strike of When the symbol is smaller than the actual geological structure, it should be placed directly above the center point of the actual geological structure. When the scale is 1:1, the symbol should be placed in the actual position of the structure and the height should be consistent with the actual structure. In addition, when symbols are placed, the strike of the symbols must be consistent with that of the actual structure.

the symbols must be consistent with that of the actual structure. Different proportions may occur between the actual length along the hinge direction and the width of the fold structure. According to the ratio of longitudinal length to trans‐ verse width, folds can be divided into linear, long‐axis, short‐axis, or equiaxed folds (dome structure and tectonic basin) [64]. However, if the symbol model is constructed according to the actual proportion of the geological structures, then the structural infor‐ mation conveyed by the fold elements is weakened and the information about the hinge zone is highlighted. Moreover, the typical geological characteristics of some ample linear Different proportions may occur between the actual length along the hinge direction and the width of the fold structure. According to the ratio of longitudinal length to transverse width, folds can be divided into linear, long-axis, short-axis, or equiaxed folds (dome structure and tectonic basin) [64]. However, if the symbol model is constructed according to the actual proportion of the geological structures, then the structural information conveyedby the fold elements is weakened and the information about the hinge zone is highlighted. Moreover, the typical geological characteristics of some ample linear folds are likely to be reflected in the spatial information of the hinge; thus, comprehensive restoration of the entire fold is also required. Therefore, intercepting a part of the fold in the hinge direction or reflecting the 3D morphology of the entire fold depends on the characteristics of the fold and preferences of the user (Figure 19).

characteristics of the fold and preferences of the user (Figure 19).

folds are likely to be reflected in the spatial information of the hinge; thus, comprehensive

**Figure 19.** Examples of fold symbol models with different hinge lengths: (**a**) local morphology and (**b**) whole morphology. **Figure 19.** Examples of fold symbol models with different hinge lengths: (**a**) local morphology and (**b**) whole morphology. **Figure 19.** Examples of fold symbol models with different hinge lengths: (**a**) local morphology and (**b**) whole morphology.

folds are likely to be reflected in the spatial information of the hinge; thus, comprehensive restoration of the entire fold is also required. Therefore, intercepting a part of the fold in the hinge direction or reflecting the 3D morphology of the entire fold depends on the

Stratigraphic legends usually use color to indicate age and texture to indicate lithol‐ ogy. Therefore, by using an appropriate legend to decorate the corresponding strata, the fold symbols can better express and transmit the age and lithology information. Stratigraphic legends usually use color to indicate age and texture to indicate lithology. Therefore, by using an appropriate legend to decorate the corresponding strata, the fold symbols can better express and transmit the age and lithology information. Stratigraphic legends usually use color to indicate age and texture to indicate lithol‐ ogy. Therefore, by using an appropriate legend to decorate the corresponding strata, the fold symbols can better express and transmit the age and lithology information.

#### *4.3. Analysis of Applicability 4.3. Analysis of Applicability 4.3. Analysis of Applicability*

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 18 of 22

characteristics of the fold and preferences of the user (Figure 19).

The premise of parametric modeling is to obtain relevant modeling parameters effec‐ tively. In this method, the modeling parameters of the fold structure can be directly ob‐ tained from the regional geological survey report or measured from corresponding maps, such as geological maps, structural maps, or section maps. Under exceptional circum‐ stances, the parameters required for fold modeling can also be empirically assigned ac‐ cording to the fold type and structural knowledge. The premise of parametric modeling is to obtain relevant modeling parameters effectively. In this method, the modeling parameters of the fold structure can be directly obtained from the regional geological survey report or measured from corresponding maps, such as geological maps, structural maps, or section maps. Under exceptional circumstances, the parameters required for fold modeling can also be empirically assigned according to the fold type and structural knowledge. The premise of parametric modeling is to obtain relevant modeling parameters effec‐ tively. In this method, the modeling parameters of the fold structure can be directly ob‐ tained from the regional geological survey report or measured from corresponding maps, such as geological maps, structural maps, or section maps. Under exceptional circum‐ stances, the parameters required for fold modeling can also be empirically assigned ac‐ cording to the fold type and structural knowledge.

This study's 3D parametric modeling method of fold symbol construction can realize symbol model construction when the geological data are insufficient. Although the sym‐ bol has an abstract requirement for the model structure, less dependence on geological data also means that the model results will be simplified in terms of morphology and large discrepancies with the actual structure may occur. In the case of complex sections, if the measured cross‐sections are available, they can be used to replace the cross‐sections gen‐ erated by the parameters, so that more accurate modeling of the fold symbol can be carried out [23–25]. For example, by selecting the cross‐section (Figure 20a) of the Danaobo moun‐ tain syncline as the middle section, a more detailed symbol model of the Danaobo moun‐ tain syncline can be built (Figure 20b). This experiment shows that the symbol model of the fold structure based on the measured cross‐section can restore the actual morphology of the current fold structure. This method also has some shortcomings, such as relying on geological section data and being unable to directly reflect the original structural mor‐ phology before being subjected to external forces. In specific applications, the appropriate method for the cross‐section can be selected according to the geological data and require‐ ments to achieve the best data utilization and model expression effect. This study's 3D parametric modeling method of fold symbol construction can realize symbol model construction when the geological data are insufficient. Although the symbol has an abstract requirement for the model structure, less dependence on geological data also means that the model results will be simplified in terms of morphology and large discrepancies with the actual structure may occur. In the case of complex sections, if the measured cross-sections are available, they can be used to replace the cross-sections generated by the parameters, so that more accurate modeling of the fold symbol can be carried out [23–25]. For example, by selecting the cross-section (Figure 20a) of the Danaobo mountain syncline as the middle section, a more detailed symbol model of the Danaobo mountain syncline can be built (Figure 20b). This experiment shows that the symbol model of the fold structure based on the measured cross-section can restore the actual morphology of the current fold structure. This method also has some shortcomings, such as relying on geological section data and being unable to directly reflect the original structural morphology before being subjected to external forces. In specific applications, the appropriate method for the cross-section can be selected according to the geological data and requirements to achieve the best data utilization and model expression effect. This study's 3D parametric modeling method of fold symbol construction can realize symbol model construction when the geological data are insufficient. Although the sym‐ bol has an abstract requirement for the model structure, less dependence on geological data also means that the model results will be simplified in terms of morphology and large discrepancies with the actual structure may occur. In the case of complex sections, if the measured cross‐sections are available, they can be used to replace the cross‐sections gen‐ erated by the parameters, so that more accurate modeling of the fold symbol can be carried out [23–25]. For example, by selecting the cross‐section (Figure 20a) of the Danaobo moun‐ tain syncline as the middle section, a more detailed symbol model of the Danaobo moun‐ tain syncline can be built (Figure 20b). This experiment shows that the symbol model of the fold structure based on the measured cross‐section can restore the actual morphology of the current fold structure. This method also has some shortcomings, such as relying on geological section data and being unable to directly reflect the original structural mor‐ phology before being subjected to external forces. In specific applications, the appropriate method for the cross‐section can be selected according to the geological data and require‐ ments to achieve the best data utilization and model expression effect.

**Figure 20.** Fold structure modeling based on measured cross-section: (**a**) section of the Danaobo mountain syncline at Hanyang Peak of Mount Lu and (**b**) symbol model of the Danaobo mountain syncline.

A box fold refers to a fold in which the strata are all overturned, and the hinge zone is flat and wide, similar to a box [64]. The box fold is a component of the world-famous Jura-type folds [65], which were developed in East Sichuan in China as well as in the syncline.

Alps in Europe. In contrast to general folds, box folds typically have two hinges and a pair of conjugate axial surfaces. The parametric modeling method of fold symbols proposed in this paper cannot easily and directly realize the symbol modeling of box folds; therefore, a particular modeling method must be designed for box folds. According to the morphological characteristics of "two hinges and a pair of conjugate axial surfaces" of box folds, the middle cross-section of a box fold can be divided into left and right parts for modeling, and the topological transition at the middle junction should be guaranteed (Figure 21). The construction method of the stratum boundaries in the middle section of the box fold is as follows: (1) obtain the positions of inflection points A and B and the corresponding extension lines along the dip angles, according to the position of the inflection points and the occurrence of both limbs input by the user; (2) calculate the distance h from A to the extension line of B, take the midpoint K of AB, and make CK⊥AB so that the length of segment CK is equal to h; and (3) take point E on CK, according to the morphological parameter cur entered by the user, make FG//AB so that the dip extension lines of the two inflection points intersect with points F and G. After that, Bezier interpolation is carried out according to points A, F, E, E, G, and B to obtain the stratum boundary of the box fold; namely, the arc AEB. Jura‐type folds [65], which were developed in East Sichuan in China as well as in the Alps in Europe. In contrast to general folds, box folds typically have two hinges and a pair of conjugate axial surfaces. The parametric modeling method of fold symbols proposed in this paper cannot easily and directly realize the symbol modeling of box folds; therefore, a particular modeling method must be designed for box folds. According to the morpho‐ logical characteristics of "two hinges and a pair of conjugate axial surfaces" of box folds, the middle cross‐section of a box fold can be divided into left and right parts for modeling, and the topological transition at the middle junction should be guaranteed (Figure 21). The construction method of the stratum boundaries in the middle section of the box fold is as follows: (1) obtain the positions of inflection points A and B and the corresponding extension lines along the dip angles, according to the position of the inflection points and the occurrence of both limbs input by the user; (2) calculate the distance h from A to the extension line of B, take the midpoint K of AB, and make CK⊥AB so that the length of segment CK is equal to h; and (3) take point E on CK, according to the morphological parameter cur entered by the user, make FG//AB so that the dip extension lines of the two inflection points intersect with points F and G. After that, Bezier interpolation is carried out according to points A, F, E, E, G, and B to obtain the stratum boundary of the box fold;

**Figure 20.** Fold structure modeling based on measured cross‐section: (**a**) section of the Danaobo mountain syncline at Hanyang Peak of Mount Lu and (**b**) symbol model of the Danaobo mountain

A box fold refers to a fold in which the strata are all overturned, and the hinge zone is flat and wide, similar to a box [64]. The box fold is a component of the world‐famous

*ISPRS Int. J. Geo‐Inf.* **2022**, *11*, 618 19 of 22

**Figure 21.** Construction process for the stratum boundaries of box folds. **Figure 21.** Construction process for the stratum boundaries of box folds.

A lack of apparent regularity is observed among stratum boundaries for more com‐ plex disharmonious folds. Rather than being deduced from a specific stratum boundary, the boundary of each stratum must be generated using the method described in Section 2.3.1. Therefore, the morphological parameters of each stratum of the disharmonious fold (such as the occurrence of inflection points, dip angles of the axial surfaces, and interlimb angles of each stratum) must be given in detail. A lack of apparent regularity is observed among stratum boundaries for more complex disharmonious folds. Rather than being deduced from a specific stratum boundary, the boundary of each stratum must be generated using the method described in Section 2.3.1. Therefore, the morphological parameters of each stratum of the disharmonious fold (such as the occurrence of inflection points, dip angles of the axial surfaces, and interlimb angles of each stratum) must be given in detail.

#### **5. Conclusions 5. Conclusions**

namely, the arc AEB.

This study aims to develop a novel parametric modeling method for the 3D symbols of fold structures using cross‐sections generated by fold parameters and a contour com‐ parison algorithm. The cross‐section generation method based on parameters and the Bezier curve effectively supports the parametric modeling of cross‐sections of various fold structures. Transition section generation based on morphing satisfies the smooth require‐ ment of the symbol model. The symbol‐modeling method can control the model accuracy and data size by setting an appropriate number and quality of cross‐sections, which may be measured or generated by fold parameters. This study aims to develop a novel parametric modeling method for the 3D symbols of fold structures using cross-sections generated by fold parameters and a contour comparison algorithm. The cross-section generation method based on parameters and the Bezier curve effectively supports the parametric modeling of cross-sections of various fold structures. Transition section generation based on morphing satisfies the smooth requirement of the symbol model. The symbol-modeling method can control the model accuracy and data size by setting an appropriate number and quality of cross-sections, which may be measured or generated by fold parameters.

In this study, different types of fold structures were selected as experimental objects for modeling experiments. The experimental results show that the generated 3D fold sym‐ bol models can abstractly express the morphologies of fold structures developed in the experimental area and can express most fold structure subdivision types. The parametric In this study, different types of fold structures were selected as experimental objects for modeling experiments. The experimental results show that the generated 3D fold symbol models can abstractly express the morphologies of fold structures developed in the experimental area and can express most fold structure subdivision types. The parametric modeling method for the 3D symbols of fold structures follows the traditional symbol creation rules [66] and cartographic rules [67]. The modeling method in this study is simple and efficient, avoids excessive dependence on geological survey data, and fits well with the abstract character of the geological symbol. Furthermore, the controllable characteristics of model accuracy and data size better meet the needs for efficient display in a wide range of 3D scenes.

In view of the inherent non-parametric geometric features of geological structures [41], this method breaks through the application limitation of parametric modeling methods in the field of 3D geological modeling. The parametric modeling method used for 3D foldsymbol modeling in this study can be applied for 3D symbol modeling of other geological structure types, such as faults, joints, and intrusions. It also applies for the 3D modeling of typical geological structures with relatively simple spatial morphology. The 3D geological symbols generated by this method have excellent research significance and application value for the 3D symbol expression of geological structures in digital Earth, digital cities, AR\VR\MR, geography teaching, science popularization, and other applications and research areas. In addition, the method in this paper has laid a sound data foundation for future research on the construction method of 3D spatio-temporal geological symbols (4D geological symbols).

**Author Contributions:** An-Bo Li conceived the original idea and is the primary author of the manuscript. Hao Chen developed the main modules of the prototype system and is the key contributor to the methodology. Xian-Yu Liu and Xiao-Feng Du developed partial modules of the prototype system. Guo-Kai Sun processed the relative data. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Key R&D Program of China (No. 2022YFB3904101, 2022YFB3904104) and by the National Natural Science Foundation of China (Project No. 41971068, 41771431).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

