doc update

README.md

@@ -1,126 +1,165 @@
 ## Spherical Triangular Tessellation (STT) Generator
 
+> **Note:** This document is automatically generated by Cursor AI.The function names,parameters,and usage instructions in the document may contain errors.Please refer to the actual function declarations in the header file(`.h`).In case of any inconsistencies,the source code shall prevail.
+
 ### Introduction
 
-The spherical triangular tessellation is a kind of partition of the spherical surface composed by only triangular cells. This program could generate the STT based on an icosahedron. The STT generated by the program could be refined around given points, lines, polygons and circles on the spherical surface. The exterior and interior outlines of the STT could also be customized.
+The spherical triangular tessellation (STT) is a method to partition a spherical surface into triangular cells. This program generates STT based on an icosahedron and provides various refinement options:
+
+1. Geometric refinement around:
+   - Points
+   - Lines
+   - Polygons
+   - Circles
+2. Topographic refinement based on elevation data
+3. Customizable exterior and interior boundaries
+
+### Technical Details
+
+#### Core Algorithm
+
+1. **Base Structure**
+   - Starts with an icosahedron (20 triangular faces)
+   - Uses quad-tree data structure for adaptive refinement
+   - Each face of the icosahedron becomes the root of a quad-tree
+
+2. **Refinement Process**
+   - Adaptive refinement based on various constraints
+   - Uses vector calculations for spherical geometry
+   - Supports multiple levels of refinement (controlled by depth parameter)
+
+3. **Constraint Types**
+   - Point constraints: Refines mesh around specific points
+   - Line constraints: Refines along great circle paths
+   - Polygon constraints: Refines within or around polygon boundaries
+   - Circle constraints: Refines within spherical caps
+   - Topographic constraints: Refines based on elevation data
+
+4. **Coordinate Systems**
+   - Supports multiple reference systems:
+     * WGS84 ellipsoid
+     * Spherical Earth
+     * Lunar sphere
+     * Custom ellipsoid
+     * Custom flattened sphere
+   - Configurable orientation of the base icosahedron
+
+### Algorithm Details
+
+1. **Initialization**
+   - Creates base icosahedron
+   - Sets up coordinate system and orientation
+   - Initializes quad-tree structure
+
+2. **Refinement Process**
+   - Reads constraint files
+   - For each triangle:
+     * Tests intersection with constraints
+     * Subdivides if needed based on depth and resolution
+     * Maintains neighbor relationships
+
+3. **Output Generation**
+   - Generates final mesh
+   - Computes vertex locations
+   - Creates neighbor lists
+   - Exports in various formats
+
+4. **Memory Management**
+   - Uses dynamic allocation for tree structures
+   - Automatically cleans up temporary data
+   - Efficient handling of large datasets
 
 ### Files and folders
 
-1. **CMakeLists.txt** CMake project file;
-2. **\*.h and \*.cpp** Source files;
-3. **README.md** This file;
-4. **stt-example.png** Screen short of an example output;
-5. **archived** Old source files;
-6. **doc** Example files.
-
-#### Source file lists
-
-```shell
-head_functions.cc
-head_functions.h
-main.cc
-progress_bar.cc
-progress_bar.h
-struct_functions.cc
-struct_functions.h
-stt_class.h
-stt_close_surface.cc
-stt_create_branch.cc
-stt_create_tree.cc
-stt_cut_hole.cc
-stt_cut_outline.cc
-stt_delete_tree.cc
-stt_get_control_circle.cc
-stt_get_control_line.cc
-stt_get_control_point.cc
-stt_in_poly_outline.cc
-stt_in_triangle_circle.cc
-stt_in_triangle_line.cc
-stt_in_triangle_point.cc
-stt_in_triangle_polygon.cc
-stt_initial_icosahedron.cc
-stt_out_poly_outline.cc
-stt_output_msh_file.cc
-stt_output_neighbor.cc
-stt_output_triangle_center_location.cc
-stt_output_vertex_location.cc
-stt_return_branch_depth.cc
-stt_return_depth.cc
-stt_return_leaf.cc
-stt_routine.cc
-stt_set_command_record.cc
-stt_set_icosahedron_orient.cc
-stt_set_pole_equator_radius.cc
-stt_set_tree_depth.cc
-stt_sort_neighbor.cc
-```
+1. **CMakeLists.txt**: CMake project configuration file
+2. **src/**: Source code directory containing all implementation files
+3. **doc/**: Example files and test cases
+4. **README.md**: Documentation (this file)
+5. **archived/**: Legacy source files (for reference only)
 
 ### Installation
 
-This program is a toolkit that is developed and maintained by Dr. Yi Zhang (zhangyiss@icloud.com) , which could be compiled and installed using the [CMake](https://cmake.org) software. Follow the instructions bellow to compile this program:
+To compile and install using [CMake](https://cmake.org):
 
 ```shell
 mkdir build
 cd build
+cmake ..
+make
 make install
 ```
 
 ### Usage
 
 ```bash
-Usage: stt -d<minimal-depth>/<maximal-depth> [-r'WGS84'|'Earth'|'Moon'|<equator-radius>/<pole-radius>|<equator_radius>,<flat-rate>] [-o<orient-longitude>/<orient-latitude>] [-m<output-msh-filename>] [-v<output-vert-loc-filename>] [-t<output-tri-cen-filename>] [-n<output-tri-neg-filename>] [-p<control-point-filename>] [-l<control-line-filename>] [-g<control-poly-filename>] [-c<control-circle-filename>] [-s<outline-shape-filename>] [-k<hole-shape-filename>] [-h]
+Usage: stt -d<minimal-depth>/<maximal-depth> [options]
+
+Required:
+  -d<min>/<max>    Minimal and maximal depths of the quad-tree structure
+                   Example: -d3/7 for refinement from depth 3 to 7
+
+Optional:
+  -r<ref>          Coordinate reference system:
+                   - 'WGS84': WGS84 ellipsoid
+                   - 'Earth': Spherical Earth
+                   - 'Moon': Lunar sphere
+                   - <eq-rad>/<pole-rad>: Custom ellipsoid
+                   - <eq-rad>,<flat-rate>: Custom flattened sphere
+  -o<lon>/<lat>    Orientation of icosahedron top vertex
+                   Default: 0/90 (North Pole)
+  
+Output options:
+  -m<file>         Output Gmsh(.msh) mesh file
+  -v<file>         Output vertices' locations
+  -t<file>         Output triangle centers
+  -n<file>         Output triangle neighbors
+
+Refinement control:
+  -p<file>         Control points file
+  -l<file>         Control lines file
+  -g<file>         Control polygons file
+  -c<file>         Control circles file
+  -t<file>         Topography control file
+  -s<file>         Outline shape file
+  -k<file>         Hole shape file
+  -z<file>         Topography data file
+  
+Help:
+  -h               Show this help message
 ```
 
-#### Options
+### Input File Formats
 
-+ __-d__: Minimal and maximal depths of the quad-tree structures used to construct the STT.
-+ __-r__: Coordinate reference system of the output files.
-+ __-o__: Orientation of the top vertex of the base icosahedron.
-+ __-m__: Output filename of the Gmsh(.msh) file.
-+ __-v__: Output filename of the vertices' location.
-+ __-t__: Output filename of the triangles' center location.
-+ __-n__: Output filename of the triangles' neighbor.
-+ __-p__: Input filename of control points' location.
-+ __-l__: Input filename of control lines' location.
-+ __-g__: Input filename of control polygons' location.
-+ __-c__: Input filename of control circles' location.
-+ __-s__: Input filename of outline shapes' location.
-+ __-k__: Input filename of hole shapes' location.
-+ __-h__: show help information.
-
-#### Input File Formats
-
-##### Point format
-
-The format of the control points' location is a plain text file. Each line of the file has the information of one control point which contains the spherical coordinates, maximal quad-tree depth, minimal resolution and physical group of the point. The program takes both the tree depth and resolution to control the fineness of the refined STT. The refinement of the STT will stop which ever the two conditions has been reached. Note that any line that starts with '#' or any empty line will be skipped. An example file:
+#### Point Control Format
+Controls refinement around specific points. Each point is defined by its location and refinement parameters:
 
 ```bash
 # <longitude> <latitude> <maximal-depth> <minimal-resolution> <physical-group>
+# longitude, latitude: Coordinates in degrees
+# maximal-depth: Maximum refinement depth for this point
+# minimal-resolution: Minimum cell size in degrees
+# physical-group: Group identifier for the refined region
 -45 -45 5 1.0 7
 45 -45 5 1.0 7
-45 45 5 1.0 7
--45 45 5 1.0 7
 ```
 
-##### Circle format 
-
-The format of the control circles' location is a plain text file. Each line of the file has the information of one control circle which contains the spherical coordinates, spherical cap degree, maximal quad-tree depth, minimal resolution and physical group of the circle. The program takes both the tree depth and resolution to control the fineness of the refined STT. The refinement of the STT will stop which ever the two conditions has been reached. Note that any line that starts with '#' or any empty line will be skipped. An example file:
+#### Circle Control Format
+Controls refinement around spherical caps. Each circle is defined by its center, radius, and refinement parameters:
 
 ```bash
 # <longitude> <latitude> <spherical-cap-degree> <maximal-depth> <minimal-resolution> <physical-group>
+# spherical-cap-degree: Angular radius of the cap in degrees
 45 60 30 5 0.1 12
 -20 -45 20 6 0.1 13
 ```
 
-##### Line format
-
-The format of the control lines', polygons', outlines', and holes' location is a plain text file. Blocks separated by empty lines contain information of control units. For each block, the first line has the number of the spherical locations, maximal quad-tree depth, minimal resolution and physical group of the unit. Followed by spherical coordinates of the unit.The program takes both the tree depth and resolution to control the fineness of the refined STT. The refinement of the STT will stop which ever the two conditions has been reached. Note that any line that starts with '#' or any empty line will be skipped. An example file:
+#### Line/Polygon Control Format
+Controls refinement along lines or around polygons. Each shape is defined by a sequence of points and refinement parameters:
 
 ```bash
-# <number-of-points> <maximal-depth> <minimal-resolution> <physical-group>
-# <longitude> <latitude>
-# <longitude> <latitude>
-# ... ...
+# First line: <number-of-points> <maximal-depth> <minimal-resolution> <physical-group>
+# Following lines: <longitude> <latitude> of each point
+# Points are connected in sequence, last point connects to first for polygons
 4 6 0.1 5
 -10 10
 50 15
@@ -128,12 +167,40 @@ The format of the control lines', polygons', outlines', and holes' location is a
 -15 50
 ```
 
-### Examples
-
-An example of multi-resolution STT:
+#### Topography Control Format
+Controls refinement based on elevation data. The refinement is based on elevation variation within triangles:
 
 ```bash
+# First line: <maximum-STD> <maximal-depth> <minimal-resolution> <physical-group>
+# Following lines: <longitude> <latitude> <elevation(meters)>
+# maximum-STD: Maximum allowed standard deviation of elevation within a triangle
+200.0 10 -1 5
+-179.95 89.95 -4203.20
+-179.85 89.95 -4203.07
+-179.75 89.95 -4203.47
+```
+
+Note: maximum-STD represents the maximum standard deviation of elevation allowed in a triangle. A triangle will be refined if the elevation variation within it exceeds this threshold.
+
+### Examples
+
+1. Multi-resolution STT with geometric constraints:
+
+```bash
+# Creates a mesh with depth range 3-7, refined around lines, polygons and circles
 stt -d 3/7 -m example.msh -l doc/control_lines.txt -g doc/control_poly.txt -c doc/control_circle.txt
 ```
 
 ![stt-example](doc/stt-example.png)
+
+2. Topography-constrained STT:
+
+```bash
+# Creates a mesh with depth range 6-10, refined based on topography data
+stt -d 6/10 -m topo_example.msh -t doc/control_topo.txt
+```
+
+![topo-example](doc/topo_constraint.png)
+
+### Contact
+For any bug reports or you need some help. Please contact Dr. Yi Zhang ([yizhang-geo@zju.edu.cn](yizhang-geo@zju.edu.cn)).

doc/control_topo.txt

@@ -0,0 +1,10 @@
+200.0 10 -1 5
+-179.95 89.95 -4203.20800781
+-179.85 89.95 -4203.07958984
+-179.75 89.95 -4203.47851562
+-179.65 89.95 -4203.89404297
+-179.55 89.95 -4204.14257812
+-179.45 89.95 -4204.40771484
+-179.35 89.95 -4204.47558594
+-179.25 89.95 -4204.46972656
+-179.15 89.95 -4209.984375

src/main.cc

@@ -6,19 +6,27 @@ void disp_help(char* proname){
 " -d<minimal-depth>/<maximal-depth> \
 [-r'WGS84'|'Earth'|'Moon'|<equator-radius>/<pole-radius>|<equator_radius>,<flat-rate>] \
 [-o<orient-longitude>/<orient-latitude>] \
-[-m<output-msh-filename>] \
-[-v<output-vert-loc-filename>] \
-[-t<output-tri-cen-filename>] \
-[-n<output-tri-neg-filename>] \
-[-p<control-point-filename>] \
-[-l<control-line-filename>] \
-[-g<control-poly-filename>] \
-[-c<control-circle-filename>] \
-[-s<outline-shape-filename>] \
-[-k<hole-shape-filename>] \
+[-m<output-msh-file>] \
+[-v<output-vert-loc-file>] \
+[-t<output-tri-cen-file>] \
+[-n<output-tri-neg-file>] \
+[-p<control-point-file>] \
+[-l<control-line-file>] \
+[-g<control-poly-filen>] \
+[-c<control-circle-file>] \
+[-s<outline-shape-file>] \
+[-k<hole-shape-file>] \
+[-z<control-topography-file>] \
 [-h]";
-	clog << proname << " - v1.3 A generator of the Spherical Triangular Tessellation (STT)." << endl;
-	clog << "Usage: " << exe_name << endl;
+	clog << "       _    _   \n";
+	clog << "  ___ | |_ | |_ \n";
+	clog << " / __|| __|| __|\n";
+	clog << " \\__ \\| |_ | |_ \n";
+	clog << " |___/ \\__| \\__|\n";
+	clog << proname << " - v1.4 - A generator of the Spherical Triangular Tessellation (STT).\n\
+This program is distributed under a dual licensing scheme. It is free for academic use, but a commercial license is required for commercial use.\n\n";
+	clog << "Author: Dr. Yi Zhang (yizhang-geo@zju.edu.cn)\n\n";
+	clog << "Usage: " << exe_name << endl << endl;
 	clog << "Options:" << endl;
 	clog << "\t-d\tMinimal and maximal quad-tree depths of the output STT." << endl;
 	clog << "\t-r\tCoordinate reference system of the output STT, the default is 1e+5/1e+5." << endl;
@@ -33,6 +41,7 @@ void disp_help(char* proname){
 	clog << "\t-c\tInput control circle location(.txt) filename." << endl;
 	clog << "\t-s\tInput outline polygon location(.txt) filename." << endl;
 	clog << "\t-k\tInput hole polygon location(.txt) filename." << endl;
+	clog << "\t-z\tInput topography(.txt) filename." << endl;
 	clog << "\t-h\tShow help information." << endl;
 }
 
@@ -50,9 +59,10 @@ int main(int argc, char* argv[]){
 	9 -> input filename for polygon constraints
 	10-> input filename for circle constraints
 	11-> input filename for outline shape constraints
-	12-> input filename for hole shape constraints*/
-	char input_options[13][1024];
-	for (int i = 0; i < 13; i++){
+	12-> input filename for hole shape constraints
+	13-> input filename for topography constraints*/
+	char input_options[14][1024];
+	for (int i = 0; i < 14; i++){
 		strcpy(input_options[i],"NULL");
 	}
 
@@ -63,7 +73,7 @@ int main(int argc, char* argv[]){
 	}
 
 	int curr, option_number;
-	while((curr = getopt(argc,argv,"hd:r:o:m:v:t:n:p:l:g:c:s:k:")) != -1){
+	while((curr = getopt(argc,argv,"hd:r:o:m:v:t:n:p:l:g:c:s:k:z:")) != -1){
 		// get option number
 		switch (curr){
 			case 'h': // show help information
@@ -94,6 +104,8 @@ int main(int argc, char* argv[]){
 				option_number =11; break;
 			case 'k':
 				option_number =12; break;
+			case 'z':
+				option_number =13; break;
 			case '?': //处理未定义或错误参数
 				if (optopt == 'd' || optopt == 'r' || optopt == 'o' || optopt == 'm' || optopt == 'n'
 					|| optopt == 'v' || optopt == 't' || optopt == 'p' || optopt == 'l'

src/struct_functions.cc

@@ -155,3 +155,19 @@ int LocalIndex(int id, Triangle t)
 		if (id == t.ids[i])	return i;
 	return -1;
 }
+
+double data_std(const vector<double> &data)
+{
+	if (data.size() == 0) return 0.0;
+
+	double sum = 0.0;
+	for (int i = 0; i < data.size(); i++)
+		sum += data[i];
+	
+	double mean = sum / data.size();
+	double sq_sum = 0.0;
+	for (int i = 0; i < data.size(); i++)
+		sq_sum += (data[i] - mean)*(data[i] - mean);
+
+	return sqrt(sq_sum / data.size());
+}

src/struct_functions.h

@@ -59,6 +59,13 @@ struct ControlPoint{
 };
 typedef vector<ControlPoint> ControlPointArray;
 
+struct ControlTopo{
+	int id = -1, max_depth = -1, physic_group = 0;
+	double minimal_resolution = DBL_MAX;
+	VertexArray vert; // verts for topo
+	vector<double> topo;
+};
+
 struct ControlLine{
 	int id = -1, max_depth = -1, physic_group = 0;
 	double minimal_resolution = -1.0;
@@ -89,4 +96,5 @@ Vertex RotateVertex(Vertex,Vertex,Vertex);
 Cpoint LineCrossPlane(Cpoint,Cpoint,Cpoint);
 string GetStringIndex(Vertex);
 int LocalIndex(int,Triangle);
+double data_std(const vector<double> &data);
 #endif

src/stt_class.h

@@ -20,8 +20,8 @@ public:
 	int set_icosahedron_orient(char*);
 	int Routine(char [][1024]); // for a 2D array. you must specify enough dimensional information to make it unique
 	void InitialIcosahedron(double,Vertex); //初始化一个二十面体实例 需要给定一个默认半径值 二十面体顶点的经纬坐标 在init_para函数中调用
-	void CreateBranch(int,int,int,int,int,int,int,QuadTreeNode**); //创建分枝
-	void CreateTree(int,int,int,int,QuadTree*);//创建树
+	void CreateBranch(int,int,int,int,int,int,int,QuadTreeNode**,const ControlTopo&); //创建分枝
+	void CreateTree(int,int,int,int,QuadTree*,const ControlTopo&);//创建树
 	void DeleteTree(QuadTreeNode**);//清空整颗树
 	void ReturnLeaf(QuadTreeNode**);//返回叶子
 	void ReturnDepth(QuadTreeNode**,int);
@@ -34,6 +34,8 @@ public:
 	int InTriangleLine(QuadTreeNode*);//判断插入线是否穿过节点三角形 使用的是球面下的方法 直接矢量计算 注意因为球面上的特殊关系 两个点之间的夹角不能大于等于180度 因为球面上总是沿着最短路径走 而且通常我们指的也是最短路径
 	int InTrianglePolygon(QuadTreeNode*);//判断多边形与三角形的关系
 	int InTriangleCircle(QuadTreeNode*);//判断圆与三角形的关系
+	int InTriangleTopo(QuadTreeNode*, const ControlTopo& in_topo, double diff_threshold);//判断地形与三角形的关系
+	int InTriangleTopoSet(QuadTreeNode*, const ControlTopo& in_topo, ControlTopo& out_topo);
 	int OutPolyOutline(QuadTreeNode*);//判断多边形与三角形的关系 用于切割模型边界
 	int InPolyOutline(QuadTreeNode*);//判断多边形与三角形的关系 用于切割模型边界 挖洞
 	int OutputMshFile(char*,double,double);
@@ -43,6 +45,7 @@ public:
 	int GetControlPoint(char*); //读取额外的点
 	int GetControlCircle(char*); //读取额外的圆
 	int GetControlLine(char*,ControlLineArray&); // Get control line arrays
+	int GetControlTopography(char*); // Get control topography arrays
 private:
 	// record input command line options for output records
 	string command_record_;
@@ -67,11 +70,15 @@ private:
 	// pointer array of the extracted quad-tree nodes returned according to conditions
 	QuadTreeNodePointerArray array_out_tri_pointer_;
 	// external constraint information (point, line, polygons, circles, outline polygons and hole polygons)
+	// we add control topo for v1.4
+	ControlTopo control_topo_;
 	ControlPointArray array_control_point_;
 	ControlCircleArray array_control_circle_;
 	ControlLineArray array_control_line_;
 	ControlLineArray array_control_polygon_;
 	ControlLineArray array_outline_polygon_;
 	ControlLineArray array_hole_polygon_;
+	// threshold for topography constraint
+	double topo_max_diff_;
 };
 #endif

src/stt_create_branch.cc

@@ -1,6 +1,7 @@
 #include "stt_class.h"
 
-void SttGenerator::CreateBranch(int upper_id,int order_id,int depth,int t_ids0,int t_ids1,int t_ids2,int phy_group,QuadTreeNode** node)
+void SttGenerator::CreateBranch(int upper_id,int order_id,int depth,int t_ids0,int t_ids1,int t_ids2,
+	int phy_group,QuadTreeNode** node,const ControlTopo& in_topo)
 {
 	Vertex local_vert[6];
 	QuadTreeNode* current_node;
@@ -14,12 +15,16 @@ void SttGenerator::CreateBranch(int upper_id,int order_id,int depth,int t_ids0,i
 	current_node->id = upper_id*10+order_id;//写入四叉树节点编号
 	current_node->depth = depth;//记录四叉树深度
 
+	ControlTopo out_topo;
+	InTriangleTopoSet(current_node, in_topo, out_topo);
+
 	//额外生长条件 满足其一即可生长 在局部加密模型的过程中 不同物理组的赋值顺序前后顺序为圈 多边形 线 点
 	if ((depth < tree_depth_ //基本生长条件 所有节点都能达到的深度
 	 || InTriangleCircle(current_node)
 	 || InTrianglePolygon(current_node)
 	 || InTriangleLine(current_node)
-	 || InTrianglePoint(current_node))
+	 || InTrianglePoint(current_node)
+	 || InTriangleTopo(current_node, out_topo, topo_max_diff_))
 	 && depth < max_depth_) //最大深度限制 所有节点不能超过的深度
 	{
 		ivd_ = map_id_vertex_.find(t_ids0);//利用map_ID映射找到四叉树节点的前三个点，这三个节点是上一层四叉树产生的，必然存在
@@ -52,10 +57,13 @@ void SttGenerator::CreateBranch(int upper_id,int order_id,int depth,int t_ids0,i
 			}
 		}
 
-		CreateBranch(current_node->id,1,depth+1,local_vert[0].id,local_vert[3].id,local_vert[5].id,current_node->tri.physic_group,&(current_node->children[0]));
-		CreateBranch(current_node->id,2,depth+1,local_vert[1].id,local_vert[4].id,local_vert[3].id,current_node->tri.physic_group,&(current_node->children[1]));
-		CreateBranch(current_node->id,3,depth+1,local_vert[2].id,local_vert[5].id,local_vert[4].id,current_node->tri.physic_group,&(current_node->children[2]));
-		CreateBranch(current_node->id,4,depth+1,local_vert[3].id,local_vert[4].id,local_vert[5].id,current_node->tri.physic_group,&(current_node->children[3]));
+		CreateBranch(current_node->id,1,depth+1,local_vert[0].id,local_vert[3].id,local_vert[5].id,current_node->tri.physic_group,&(current_node->children[0]), out_topo);
+		CreateBranch(current_node->id,2,depth+1,local_vert[1].id,local_vert[4].id,local_vert[3].id,current_node->tri.physic_group,&(current_node->children[1]), out_topo);
+		CreateBranch(current_node->id,3,depth+1,local_vert[2].id,local_vert[5].id,local_vert[4].id,current_node->tri.physic_group,&(current_node->children[2]), out_topo);
+		CreateBranch(current_node->id,4,depth+1,local_vert[3].id,local_vert[4].id,local_vert[5].id,current_node->tri.physic_group,&(current_node->children[3]), out_topo);
 	}
+
+	out_topo.vert.clear();
+	out_topo.topo.clear();
 	return;
 }

src/stt_create_tree.cc

@@ -1,6 +1,7 @@
 #include "stt_class.h"
 
-void SttGenerator::CreateTree(int tree_id,int t_ids0,int t_ids1,int t_ids2,QuadTree* p_tree){
+void SttGenerator::CreateTree(int tree_id,int t_ids0,int t_ids1,int t_ids2,
+	QuadTree* p_tree,const ControlTopo& in_topo){
 	if (max_depth_ == 0){
 		p_tree->root->id = 0;
 		p_tree->root->depth = 0;
@@ -12,8 +13,7 @@ void SttGenerator::CreateTree(int tree_id,int t_ids0,int t_ids1,int t_ids2,QuadT
 			p_tree->root->children[i] = nullptr;
 		}
 	}
-	else
-	{
-		CreateBranch(0,tree_id,0,t_ids0,t_ids1,t_ids2,0,&(p_tree->root));//以根节点开始创建四叉树
+	else{
+		CreateBranch(0,tree_id,0,t_ids0,t_ids1,t_ids2,0,&(p_tree->root),in_topo);//以根节点开始创建四叉树
 	}
 }

src/stt_get_control_topo.cc

@@ -0,0 +1,45 @@
+#include "stt_class.h"
+
+int SttGenerator::GetControlTopography(char* filename)
+{
+    double one_topo;
+	stringstream temp_ss;
+	string temp_str;
+	Vertex temp_vert;
+	ifstream infile;
+
+	if (!strcmp(filename,"NULL")) return 0;
+
+	if (OpenInfile(infile,filename)) return -1;
+	else{
+        while (getline(infile,temp_str)){
+            if (*(temp_str.begin()) == '#' || temp_str == "") continue;
+            else{
+                temp_ss = Str2Ss(temp_str);
+				temp_ss >> topo_max_diff_ >> control_topo_.max_depth >> control_topo_.minimal_resolution >> control_topo_.physic_group;
+				if (control_topo_.max_depth <= 0) control_topo_.max_depth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
+				if (control_topo_.minimal_resolution <= 0) control_topo_.minimal_resolution = -1.0; //这里直接给成-1
+                break;
+            }
+        }
+
+		while (getline(infile,temp_str)){
+			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
+			else{
+                getline(infile,temp_str);
+                temp_ss = Str2Ss(temp_str);
+                if (temp_ss >> temp_vert.posis.lon >> temp_vert.posis.lat >> one_topo){
+                    temp_vert.posis.rad = DefaultR;
+                    temp_vert.id = control_topo_.vert.size();
+                    temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+                    control_topo_.vert.push_back(temp_vert);
+                    control_topo_.topo.push_back(one_topo);
+                }
+			}
+		}
+
+        control_topo_.id = 0;
+		infile.close();
+	}
+	return 0;
+}

src/stt_in_triangle_topo.cc

@@ -0,0 +1,33 @@
+#include "stt_class.h"
+
+int SttGenerator::InTriangleTopo(QuadTreeNode* node, const ControlTopo& in_topo, double diff_threshold){
+	//没有插入的地形 直接返回否
+	if (in_topo.vert.empty()){
+		return 0;
+	}
+	else{
+		int node_depth;
+		double node_resolution;
+        double max_topo = -1e+30, min_topo = 1e+30;
+
+		Triangle temp_tri;
+		for (int j = 0; j < 3; j++){
+			temp_tri.ids[j] = node->tri.ids[j];
+		}
+
+		node_depth = node->depth;
+
+		node_resolution = 0;
+		for (int i = 0; i < 3; i++){
+			node_resolution += acos(DotProduct(array_stt_vert_[temp_tri.ids[i]].posic,array_stt_vert_[temp_tri.ids[(i+1)%3]].posic)
+				/(ModuleLength(array_stt_vert_[temp_tri.ids[i]].posic)*ModuleLength(array_stt_vert_[temp_tri.ids[(i+1)%3]].posic)));
+		}
+		node_resolution = node_resolution*60/Pi;
+
+		// 将控制点的组别赋值给当前节点
+		node->tri.physic_group = in_topo.physic_group;
+
+        if (data_std(in_topo.topo) >= diff_threshold) return 1;
+		else return 0;
+	}
+}

src/stt_in_triangle_topo_set.cc

@@ -0,0 +1,52 @@
+#include "stt_class.h"
+
+int SttGenerator::InTriangleTopoSet(QuadTreeNode* node, const ControlTopo& in_topo, ControlTopo& out_topo){
+	//没有插入的地形 直接返回否
+	if (in_topo.vert.empty()){
+        out_topo.vert.clear();
+        out_topo.topo.clear();
+		return 0;
+	}
+	else{
+        out_topo.vert.clear();
+        out_topo.topo.clear();
+        out_topo.max_depth = in_topo.max_depth;
+        out_topo.minimal_resolution = in_topo.minimal_resolution;
+        out_topo.physic_group = in_topo.physic_group;
+
+		int count;
+		Cpoint tri_nor;
+		Cpoint cross_point;
+
+		Triangle temp_tri;
+		for (int j = 0; j < 3; j++){
+			temp_tri.ids[j] = node->tri.ids[j];
+		}
+
+		tri_nor = CrossProduct(array_stt_vert_[temp_tri.ids[1]].posic - array_stt_vert_[temp_tri.ids[0]].posic,
+		 array_stt_vert_[temp_tri.ids[2]].posic - array_stt_vert_[temp_tri.ids[0]].posic);
+
+        // 这一步可能会非常耗时
+		for (int i = 0; i < in_topo.vert.size(); i++){
+			// 控制点和三角形的法线方向相同（在同一个半球）
+			if (DotProduct(tri_nor, in_topo.vert[i].posic) > 0){
+				count = 0;
+				for (int j = 0; j < 3; j++){
+					cross_point = LineCrossPlane(array_stt_vert_[temp_tri.ids[j]].posic, tri_nor, in_topo.vert[i].posic);
+					if (DotProduct(tri_nor,
+						CrossProduct(array_stt_vert_[temp_tri.ids[(j+1)%3]].posic - array_stt_vert_[temp_tri.ids[j]].posic,
+						 cross_point - array_stt_vert_[temp_tri.ids[j]].posic)) > 0){
+						count++;
+					}
+				}
+
+                // 满足条件则 穿透点在三角形内
+				if (count == 3){
+                    out_topo.vert.push_back(in_topo.vert[i]);
+                    out_topo.topo.push_back(in_topo.topo[i]);
+				}
+			}
+		}
+		return 0;
+	}
+}

src/stt_routine.cc

@@ -12,6 +12,7 @@ int SttGenerator::Routine(char input_options[][1024]){
 	// get outline and hole polygons
 	if (GetControlPoint(input_options[7])) return -1;
 	if (GetControlCircle(input_options[10])) return -1;
+	if (GetControlTopography(input_options[13])) return -1;
 	if (GetControlLine(input_options[8],array_control_line_)) return -1;
 	if (GetControlLine(input_options[9],array_control_polygon_)) return -1;
 	if (GetControlLine(input_options[11],array_outline_polygon_)) return -1;
@@ -42,7 +43,7 @@ int SttGenerator::Routine(char input_options[][1024]){
 	for (int i = 0; i < 20; i++){
 		bar->Progressed(i);
 		// initialize the tree index starts from 50 to avoid possible repetition of vertex's index
-		CreateTree(i+50,base_icosahedron_.tri[i].ids[0],base_icosahedron_.tri[i].ids[1],base_icosahedron_.tri[i].ids[2], forest_[i]);
+		CreateTree(i+50,base_icosahedron_.tri[i].ids[0],base_icosahedron_.tri[i].ids[1],base_icosahedron_.tri[i].ids[2], forest_[i], control_topo_);
 	}
 
 	delete bar;