version update

2025-01-15 22:42:06 +08:00 · 2025-01-15 22:42:06 +08:00 · b0ab1fe089
commit b0ab1fe089
parent d9ae10d0e4
14 changed files with 307 additions and 114 deletions
--- a/README.md
+++ b/README.md
@ -2,109 +2,86 @@

 ### 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

 ### 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:
+This program is developed and maintained by Dr. Yi Zhang (yizhang-geo@zju.edu.cn). 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
+
+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
+  
+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:

 ```bash
 # <longitude> <latitude> <maximal-depth> <minimal-resolution> <physical-group>
 -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:

 ```bash
 # <longitude> <latitude> <spherical-cap-degree> <maximal-depth> <minimal-resolution> <physical-group>
@ -112,15 +89,12 @@ The format of the control circles' location is a plain text file. Each line of t
 -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:

 ```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
 4 6 0.1 5
 -10 10
 50 15
@ -128,12 +102,32 @@ The format of the control lines', polygons', outlines', and holes' location is a
 -15 50
 ```

+#### Topography Control Format
+Controls refinement based on elevation data:
+
+```bash
+# First line: <maximum-STD> <maximal-depth> <minimal-resolution> <physical-group>
+# Following lines: <longitude> <latitude> <elevation(meters)>
+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.
+
 ### Examples

-An example of multi-resolution STT:
+1. Multi-resolution STT with geometric constraints:

 ```bash
 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
+stt -d 6/10 -m topo_example.msh -t doc/control_topo.txt
+```
+![topo-example](doc/topo_constraint.png)
--- a/STT操作说明书.docx
+++ b/STT操作说明书.docx
--- a/doc/control_topo.txt
+++ b/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
--- a/doc/topo_constraint.png
+++ b/doc/topo_constraint.png
--- a/src/main.cc
+++ b/src/main.cc
@ -6,18 +6,19 @@ 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 << proname << " - v1.4 A generator of the Spherical Triangular Tessellation (STT)." << endl;
 	clog << "Usage: " << exe_name << endl;
 	clog << "Options:" << endl;
 	clog << "\t-d\tMinimal and maximal quad-tree depths of the output STT." << endl;
@ -33,6 +34,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 +52,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 +66,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 +97,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'
--- a/src/struct_functions.cc
+++ b/src/struct_functions.cc
@ -154,4 +154,20 @@ int LocalIndex(int id, Triangle t)
 	for (int i = 0; i < 3; i++)
 		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());
 }
--- a/src/struct_functions.h
+++ b/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
--- a/src/stt_class.h
+++ b/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
--- a/src/stt_create_branch.cc
+++ b/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;
 }
--- a/src/stt_create_tree.cc
+++ b/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);//以根节点开始创建四叉树
 	}
 }
--- a/src/stt_get_control_topo.cc
+++ b/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;
+}
--- a/src/stt_in_triangle_topo.cc
+++ b/src/stt_in_triangle_topo.cc
@ -0,0 +1,47 @@
+#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;
+/*
+		for (int i = 0; i < in_topo.vert.size(); i++){
+			// 1. 允许的最大深度大于当前节点的深度
+			// 2. 允许的最小分辨率小于当前节点的分辨率
+			if (in_topo.max_depth >= node_depth && node_resolution >= in_topo.minimal_resolution){
+				// 统计符合条件的地形的最大值与最小值
+				if (in_topo.topo[i] > max_topo) max_topo = in_topo.topo[i];
+				if (in_topo.topo[i] < max_topo) min_topo = in_topo.topo[i];
+			}
+		}
+
+		if (std::abs(max_topo - min_topo) >= diff_threshold) return 1;
+		return 0;
+*/
+	}
+}
--- a/src/stt_in_triangle_topo_set.cc
+++ b/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;
+	}
+}
--- a/src/stt_routine.cc
+++ b/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;