initial upload

src/CMakeLists.txt

@@ -0,0 +1,8 @@
+aux_source_directory(. STT_SRC)
+
+add_executable(stt ${STT_SRC})
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+set_target_properties(stt PROPERTIES CXX_STANDARD 11)
+set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin)
+install(TARGETS stt RUNTIME DESTINATION sbin)

src/head_functions.cc

@@ -0,0 +1,65 @@
+#include "head_functions.h"
+
+// degree version of sin() and cos()
+inline double sind(double degree){
+	return sin(degree*Pi/180.0);
+}
+
+inline double cosd(double degree){
+	return cos(degree*Pi/180.0);
+}
+// calculate the semi-radius at a given latitude given the pole and equator radius of the reference system
+double EllipsoidRadius(double latitude,double pole_radius,double equator_radius){
+	return pole_radius*equator_radius/sqrt(pow(pole_radius,2)*pow(cosd(latitude),2)+pow(equator_radius,2)*pow(sind(latitude),2));
+}
+// Bilinear interpolation on the sphere
+/*
+ latitude
+   |
+   |
+value21-------value22
+   |             |
+   |             |
+   |      0      |
+   |             |
+   |             |
+value11-------value12----> longitude
+*/
+double SphBilinearInterpolation(
+	double latitude1,double latitude2,double longitude1,double longitude2,
+	double latitude0,double longitude0,
+	double value11,double value12,double value21,double value22){
+	double colatitude1 = 90.0 - latitude1;
+	double colatitude2 = 90.0 - latitude2;
+	double colatitude0 = 90.0 - latitude0;
+	double delta=(longitude2-longitude1)*(cosd(colatitude2)-cosd(colatitude1));
+	double A=(longitude1*(value12-value22)+longitude2*(value21-value11))/delta;
+	double B=(cosd(colatitude1)*(value21-value22)+cosd(colatitude2)*(value12-value11))/delta;
+	double C=(value11+value22-value21-value12)/delta;
+	double D=(longitude2*cosd(colatitude2)*value11-longitude2*cosd(colatitude1)*value21+longitude1*cosd(colatitude1)*value22-longitude1*cosd(colatitude2)*value12)/delta;
+	return A*cosd(colatitude0)+B*longitude0+C*longitude0*cosd(colatitude0)+D;
+}
+// Convert a string to stringstream
+stringstream Str2Ss(string input_string){
+	stringstream sstr;
+	sstr.str(""); sstr.clear(); sstr.str(input_string);
+	return sstr;
+}
+// Check the existence of a input file and return the running status
+int OpenInfile(ifstream &input_file,char* filename){
+	input_file.open(filename);
+	if (!input_file){
+		cerr << BOLDRED << "Error ==> " << RESET << "file not found: " << filename << endl;
+		return -1;
+	}
+	return 0;
+}
+// Check the existence of a output file and return the running status
+int OpenOutfile(ofstream &output_file,char* filename){
+	output_file.open(filename);
+	if (!output_file){
+		cerr << BOLDRED << "Error ==> " << RESET << "fail to create the file: " << filename << endl;
+		return -1;
+	}
+	return 0;
+}

src/head_functions.h

@@ -0,0 +1,71 @@
+#ifndef _HEAD_FUNCTIONS_H
+#define _HEAD_FUNCTIONS_H
+#include "iostream"
+#include "fstream"
+#include "sstream"
+#include "string.h"
+#include "cmath"
+#include "iomanip"
+#include "stdio.h"
+#include "stdlib.h"
+#include "unistd.h"
+#include "vector"
+#include "map"
+#include "algorithm"
+#include "ctime"
+
+using namespace std;
+
+// terminal controls
+#define BOLDRED "\033[1m\033[31m"
+#define BOLDGREEN "\033[1m\033[32m"
+#define BOLDBLUE "\033[1m\033[34m"
+#define UNDERLINE "\033[1m\033[4m"
+#define RESET "\033[0m"
+#define MOVEUP(x) printf("\033[%dA", (x))
+#define MOVEDOWN(x) printf("\033[%dB", (x))
+#define MOVELEFT(x) printf("\033[%dD", (x))
+#define MOVERIGHT(x) printf("\033[%dC", (x))
+#define MOVETO(y,x) printf("\033[%d;%dH", (y), (x))
+#define CLEARLINE "\033[K"
+#define CLEARALL "\033[2J"
+
+// define some mathematic constants
+#define DBL_MAX 1.0e+30
+#define DBL_MIN -1.0e+30
+#define ZERO 1.0e-20
+#define DefaultR 1e+5
+
+// define some physical constants
+// semi-radius (pole and equator) of the WGS84 reference system
+#define WGS84_r 6356752.314245179
+#define WGS84_R 6378137.0
+// mean radius of the Earth and Moon
+#define Earth_r 6371008.8
+#define Moon_r 1738000
+// Pi and the golden ratio
+#define Pi (4.0*atan(1.0))
+#define GoldenMean ((sqrt(5.0)+1)/2)
+// Universal gravitational constant
+#define G0 6.67408e-11
+
+// Macro functions
+#define MAX(a,b) (a>b?a:b)
+#define MIN(a,b) (a<b?a:b)
+#define SetToBox(low,high,input) (MAX(low,MIN(high,input)))
+
+// Define global functions
+// degree version of sin() and cos()
+inline double sind(double);
+inline double cosd(double);
+// calculate the semi-radius at a given latitude given the pole and equator radius of the reference system
+double EllipsoidRadius(double,double,double);
+// Bilinear interpolation on the spherical sphere
+double SphBilinearInterpolation(double,double,double,double,double,double,double,double,double,double);
+// Convert a string to stringstream
+stringstream Str2Ss(string);
+// Check the existence of a input file and return the running status
+int OpenInfile(ifstream&,char*);
+// Check the existence of a output file and return the running status
+int OpenOutfile(ofstream&,char*);
+#endif

src/main.cc

@@ -0,0 +1,128 @@
+#include "stt_class.h"
+
+void disp_help(char* proname){
+	string exe_name = proname;
+	exe_name = " " + exe_name + 
+" -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]";
+	clog << proname << " - v1.3 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;
+	clog << "\t-r\tCoordinate reference system of the output STT, the default is 1e+5/1e+5." << endl;
+	clog << "\t-o\tOrientation of the top vertex of the base icosahedron, the default is 0/90." << endl;
+	clog << "\t-m\tOutput Gmsh(.msh) filename." << endl;
+	clog << "\t-v\tOutput vertex location(.txt) filename." << endl;
+	clog << "\t-t\tOutput triangle center location(.txt) filename." << endl;
+	clog << "\t-n\tOutput triangle neighbor(.txt) filename." << endl;
+	clog << "\t-p\tInput control point location(.txt) filename." << endl;
+	clog << "\t-l\tInput control line location(.txt) filename." << endl;
+	clog << "\t-g\tInput control polygon location(.txt) filename." << endl;
+	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-h\tShow help information." << endl;
+}
+
+int main(int argc, char* argv[]){
+	/*map of input options
+	0 -> tree depths
+	1 -> coordinate reference system
+	2 -> orientation of the base icosahedron
+	3 -> output filename for the constructed model (Gmsh .msh file)
+	4 -> output filename for vertex locations
+	5 -> output filename for triangle's center locations
+	6 -> output filename for triangle's neighboring index
+	7 -> input filename for point constraints
+	8 -> input filename for line constraints
+	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++){
+		strcpy(input_options[i],"NULL");
+	}
+
+	// show help information is no options is read
+	if (argc == 1){
+		disp_help(argv[0]);
+		return 0;
+	}
+
+	int curr, option_number;
+	while((curr = getopt(argc,argv,"hd:r:o:m:v:t:n:p:l:g:c:s:k:")) != -1){
+		// get option number
+		switch (curr){
+			case 'h': // show help information
+				disp_help(argv[0]); return 0;
+			case 'd':
+				option_number = 0; break;
+			case 'r':
+				option_number = 1; break;
+			case 'o':
+				option_number = 2; break;
+			case 'm':
+				option_number = 3; break;
+			case 'v':
+				option_number = 4; break;
+			case 't':
+				option_number = 5; break;
+			case 'n':
+				option_number = 6; break;
+			case 'p':
+				option_number = 7; break;
+			case 'l':
+				option_number = 8; break;
+			case 'g':
+				option_number = 9; break;
+			case 'c':
+				option_number =10; break;
+			case 's':
+				option_number =11; break;
+			case 'k':
+				option_number =12; break;
+			case '?': //处理未定义或错误参数
+				if (optopt == 'd' || optopt == 'r' || optopt == 'o' || optopt == 'm' || optopt == 'n'
+					|| optopt == 'v' || optopt == 't' || optopt == 'p' || optopt == 'l'
+					|| optopt == 'g' || optopt == 'c' || optopt == 's' || optopt == 'k'){
+					fprintf (stderr, "Option -%c requires an argument.\n", optopt);
+					return 0;
+				}
+				else if (isprint(optopt)){
+					fprintf (stderr, "Unknown option `-%c'.\n", optopt);
+					return 0;
+				}
+				else{
+					fprintf (stderr,"Unknown option character `\\x%x'.\n",optopt);
+					return 0;
+				} break;
+			default:
+				abort();
+		}
+
+		if (1!=sscanf(optarg,"%s",input_options[option_number])){
+			cerr<<BOLDRED<<"Error ==> "<<RESET<<"bad syntax: "<<optarg<<endl;
+			return 0;
+		}
+	}
+
+	SttGenerator instance;
+	// record commands
+	instance.set_command_record(argc,argv);
+	// generate stt model
+	instance.Routine(input_options);
+	return 0;
+}

src/progress_bar.cc

@@ -0,0 +1,114 @@
+//#ifdef _WINDOWS
+//#include <windows.h>
+//#else
+//#include <sys/ioctl.h>
+//#endif
+
+#include "progress_bar.h"
+
+ProgressBar::ProgressBar() {}
+
+ProgressBar::ProgressBar(unsigned long n_, const char* description_, std::ostream& out_){
+    
+    n = n_;
+    frequency_update = n_;
+    description = description_;
+    out = &out_;
+    
+    unit_bar = "\u2588";
+    unit_space = "-";
+    desc_width = std::strlen(description);  // character width of description field
+
+}
+
+void ProgressBar::SetFrequencyUpdate(unsigned long frequency_update_){
+    
+    if(frequency_update_ > n){
+        frequency_update = n;    // prevents crash if freq_updates_ > n_
+    }
+    else{
+        frequency_update = frequency_update_;
+    }
+}
+
+void ProgressBar::SetStyle(const char* unit_bar_, const char* unit_space_){
+    
+    unit_bar = unit_bar_;
+    unit_space = unit_space_;
+}
+
+int ProgressBar::GetConsoleWidth(){
+
+    int width;
+
+    #ifdef _WINDOWS
+        CONSOLE_SCREEN_BUFFER_INFO csbi;
+        GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
+        width = csbi.srWindow.Right - csbi.srWindow.Left;
+    #else
+        struct winsize win;
+        //注意！当我们使用pipe here-doc等通道获取程序参数时无法正确的获取窗口大小 此时我们将使用预定值
+        if (ioctl(0, TIOCGWINSZ, &win) != -1)
+            width = win.ws_col;
+        else width = 100;
+    #endif
+
+    return width;
+}
+
+int ProgressBar::GetBarLength(){
+
+    // get console width and according adjust the length of the progress bar
+
+    int bar_length = static_cast<int>((GetConsoleWidth() - desc_width - CHARACTER_WIDTH_PERCENTAGE) / 2.);
+
+    return bar_length;
+}
+
+void ProgressBar::ClearBarField(){
+
+    for(int i=0;i<GetConsoleWidth();++i){
+        *out << " ";
+    }
+    *out << "\r" << std::flush;
+}
+
+void ProgressBar::Progressed(unsigned long idx_)
+{
+    try{
+        if(idx_ > n) throw idx_;
+
+        // determines whether to update the progress bar from frequency_update
+        if ((idx_ != n-1) && ((idx_+1) % (n/frequency_update) != 0)) return;
+
+        // calculate the size of the progress bar
+        int bar_size = GetBarLength();
+    
+        // calculate percentage of progress
+        double progress_percent = idx_* TOTAL_PERCENTAGE/(n-1);
+
+        // calculate the percentage value of a unit bar 
+        double percent_per_unit_bar = TOTAL_PERCENTAGE/bar_size;
+
+        // display progress bar
+        *out << " " << description << " |";
+
+        for(int bar_length=0;bar_length<=bar_size-1;++bar_length){
+            if(bar_length*percent_per_unit_bar<progress_percent){
+                *out << unit_bar;
+            }
+            else{
+                *out << unit_space;
+            }
+        }
+
+        if(idx_ == n-1)
+            *out << "|" << std::setw(CHARACTER_WIDTH_PERCENTAGE + 1) << std::setprecision(1) << std::fixed << progress_percent << "%\r" << std::flush << std::endl;
+        else *out << "|" << std::setw(CHARACTER_WIDTH_PERCENTAGE + 1) << std::setprecision(1) << std::fixed << progress_percent << "%\r" << std::flush;
+    }
+    catch(unsigned long e){
+        ClearBarField();
+        std::cerr << "PROGRESS_BAR_EXCEPTION: _idx (" << e << ") went out of bounds, greater than n (" << n << ")." << std::endl << std::flush;
+    }
+
+}

src/progress_bar.h

@@ -0,0 +1,41 @@
+#ifndef _PROGRESS_BAR_
+#define _PROGRESS_BAR_
+#include <sys/ioctl.h>
+#include <iostream>
+#include <iomanip>
+#include <cstring>
+#include <thread>
+#include <chrono>
+
+#define TOTAL_PERCENTAGE 100.0
+#define CHARACTER_WIDTH_PERCENTAGE 4
+
+class ProgressBar
+{
+public: 
+
+    ProgressBar();
+    ProgressBar(unsigned long n_, const char *description_="", std::ostream& out_=std::cerr);
+
+    void SetFrequencyUpdate(unsigned long frequency_update_);
+    void SetStyle(const char* unit_bar_, const char* unit_space_);
+
+    void Progressed(unsigned long idx_);
+
+private:
+	
+    unsigned long n;
+    unsigned int desc_width;
+    unsigned long frequency_update;
+    std::ostream* out;
+		
+    const char *description;
+    const char *unit_bar;
+    const char *unit_space;
+		
+    void ClearBarField();
+    int GetConsoleWidth();
+    int GetBarLength();
+
+};
+#endif

src/struct_functions.cc

@@ -0,0 +1,157 @@
+#include "struct_functions.h"
+
+// mathematic functions regarding declared structures
+bool operator ==(Cpoint a,Cpoint b){
+	if(fabs(a.x-b.x)<ZERO&&fabs(a.y-b.y)<ZERO&&fabs(a.z-b.z)<ZERO){
+		return 1;
+	}
+	else return 0;
+}
+
+bool operator !=(Cpoint a,Cpoint b){
+	if(fabs(a.x-b.x)>=ZERO || fabs(a.y-b.y)>=ZERO || fabs(a.z-b.z)>=ZERO){
+		return 1;
+	}
+	else return 0;
+}
+
+Cpoint operator +(Cpoint a,Cpoint b){
+	Cpoint m;
+	m.x=a.x+b.x; m.y=a.y+b.y; m.z=a.z+b.z;
+	return m;
+}
+
+Cpoint operator -(Cpoint a,Cpoint b){
+	Cpoint m;
+	m.x=a.x-b.x; m.y=a.y-b.y; m.z=a.z-b.z;
+	return m;
+}
+
+Cpoint operator *(double sign,Cpoint b){
+	Cpoint m;
+	m.x=sign*b.x; m.y=sign*b.y; m.z=sign*b.z;
+	return m;
+}
+
+double DotProduct(Cpoint a,Cpoint b){
+	return a.x*b.x+a.y*b.y+a.z*b.z;
+}
+
+Cpoint CrossProduct(Cpoint a,Cpoint b){
+	Cpoint v;
+	v.x = a.y*b.z-a.z*b.y;
+	v.y = a.z*b.x-a.x*b.z;
+	v.z = a.x*b.y-a.y*b.x;
+	return v;
+}
+
+Cpoint CloudCenter(VertexArray input_verts){
+	Cpoint c;
+	c.x = 0; c.y = 0; c.z = 0;
+	if (!input_verts.empty()){
+		for (int i = 0; i < input_verts.size(); i++){
+			c.x += input_verts[i].posic.x;
+			c.y += input_verts[i].posic.y;
+			c.z += input_verts[i].posic.z;
+		}
+		c.x /= input_verts.size();
+		c.y /= input_verts.size();
+		c.z /= input_verts.size();
+	}
+	return c;
+}
+
+double ModuleLength(Cpoint v){
+	return sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
+}
+
+double ProjectAngle(Cpoint a,Cpoint b){
+	return acos((a.x*b.x+a.y*b.y+a.z*b.z)/(sqrt(a.x*a.x+a.y*a.y+a.z*a.z)*sqrt(b.x*b.x+b.y*b.y+b.z*b.z)))*180.0/Pi;
+}
+
+Cpoint Sphere2Cartesian(Spoint s){
+	Cpoint c;
+	c.x = s.rad*sin((0.5 - s.lat/180.0)*Pi)*cos((2.0 + s.lon/180.0)*Pi);
+	c.y = s.rad*sin((0.5 - s.lat/180.0)*Pi)*sin((2.0 + s.lon/180.0)*Pi);
+	c.z = s.rad*cos((0.5 - s.lat/180.0)*Pi);
+	return c;
+}
+
+Spoint Cartesian2Sphere(Cpoint c){
+	Spoint s;
+	s.rad = sqrt(pow(c.x,2)+pow(c.y,2)+pow(c.z,2));
+	//点距离原点极近 将点置于原点
+	if (fabs(s.rad)<ZERO){
+		s.lat = s.lon = s.rad = 0.0;
+	}
+	else{
+		s.lat = 90.0 - acos(c.z/s.rad)*180.0/Pi;
+		s.lon = atan2(c.y,c.x)*180.0/Pi;
+	}
+	return s;
+}
+
+Vertex RotateVertex(Vertex olda,Vertex newa,Vertex oldb){
+	Vertex newb;
+	Vertex temp_ref,temp_b;
+	double yz_angle = (newa.posis.lat - olda.posis.lat)*Pi/180.0;
+	//首先绕olda.lon即x轴旋转oldb
+	temp_b.posic.x = oldb.posic.x;
+	temp_b.posic.y = oldb.posic.y*cos(-1.0*yz_angle)+oldb.posic.z*sin(-1.0*yz_angle);
+	temp_b.posic.z = oldb.posic.z*cos(-1.0*yz_angle)-oldb.posic.y*sin(-1.0*yz_angle);
+	temp_b.posis = Cartesian2Sphere(temp_b.posic);
+	//计算绕x轴旋转后olda的位置 这是后一步旋转需要的参考值
+	temp_ref.posic.x = olda.posic.x;
+	temp_ref.posic.y = olda.posic.y*cos(-1.0*yz_angle)+olda.posic.z*sin(-1.0*yz_angle);
+	temp_ref.posic.z = olda.posic.z*cos(-1.0*yz_angle)-olda.posic.y*sin(-1.0*yz_angle);
+	temp_ref.posis = Cartesian2Sphere(temp_ref.posic);
+	//注意绕z轴旋转的经度参考位置为olda绕x轴旋转后的经度值
+	double xy_angle = (newa.posis.lon - temp_ref.posis.lon)*Pi/180.0;
+	//绕z轴旋转temp_b z值不变
+	newb.id = oldb.id;
+	newb.posic.x = temp_b.posic.x*cos(-1.0*xy_angle)+temp_b.posic.y*sin(-1.0*xy_angle);
+	newb.posic.y = temp_b.posic.y*cos(-1.0*xy_angle)-temp_b.posic.x*sin(-1.0*xy_angle);
+	newb.posic.z = temp_b.posic.z;
+	newb.posis = Cartesian2Sphere(newb.posic);
+	return newb;
+}
+
+Cpoint LineCrossPlane(Cpoint c,Cpoint normal,Cpoint p){
+	Cpoint m;
+	m.x = 0; m.y = 0; m.z = 0;
+	double t;
+	if (DotProduct(normal,p) != 0) //平面与矢量平行
+	{
+		t = DotProduct(normal,c)/DotProduct(normal,p);
+		m.x += p.x*t;
+		m.y += p.y*t;
+		m.z += p.z*t;
+	}
+	return m;
+}
+
+string GetStringIndex(Vertex input_vertex)
+{
+	stringstream sstemp;
+	string vert_id,mid_id;
+
+	sstemp.str(""); sstemp.clear();
+	sstemp << setprecision(16) << input_vertex.posic.x;
+	sstemp >> vert_id;
+	sstemp.str(""); sstemp.clear();
+	sstemp << setprecision(16) << input_vertex.posic.y;
+	sstemp >> mid_id;
+	vert_id = vert_id + " " + mid_id;
+	sstemp.str(""); sstemp.clear();
+	sstemp << setprecision(16) << input_vertex.posic.z;
+	sstemp >> mid_id;
+	vert_id = vert_id + " " + mid_id;
+	return vert_id;
+}
+
+int LocalIndex(int id, Triangle t)
+{
+	for (int i = 0; i < 3; i++)
+		if (id == t.ids[i])	return i;
+	return -1;
+}

src/struct_functions.h

@@ -0,0 +1,92 @@
+#ifndef _STRUCT_FUNCTIONS_H
+#define _STRUCT_FUNCTIONS_H
+#include "head_functions.h"
+
+// array structures
+typedef vector<int> IntArray1D; //1D int array
+typedef vector<vector<int> > IntArray2D; // 2D int array
+typedef vector<double> DoubleArray1D; // 1D double array
+// map structures
+typedef map<int,int> Int2IntMap; // int to int map
+// point structures in the Cartesian and spherical coordinates
+struct Cpoint{
+	double x = DBL_MAX, y = DBL_MAX, z = DBL_MAX;
+};
+typedef vector<Cpoint> CpointArray;
+// lon for longitude, lat for latitude and rad for radius
+struct Spoint{
+	double lon = DBL_MAX, lat = DBL_MAX, rad = DBL_MAX;
+};
+typedef vector<Spoint> SpointArray;
+// vertex structure
+struct Vertex{
+	int id = -1;
+	Cpoint posic; //position under the Cartesian coordinates
+	Spoint posis; //position under the sphere coordinates
+};
+typedef vector<Vertex> VertexArray;
+typedef map<int,Vertex> Int2VertexMap;
+typedef map<string,Vertex> String2VertexMap;
+// triangle structure
+struct Triangle{
+	int ids[3] = {-1,-1,-1}; // index s
+	int physic_group = 0;
+};
+typedef vector<Triangle> TriangleArray;
+// icosahedron structure
+struct Icosahedron{
+	Vertex vert[12]; // vert for vertex
+	Triangle tri[20]; // tir for triangle
+};
+// Quad-tree node structure
+struct QuadTreeNode{
+	int id = -1, depth = -1;
+	bool out_ok = true;
+	Triangle tri;
+	QuadTreeNode *children[4] = {nullptr,nullptr,nullptr,nullptr};
+	QuadTreeNode *neighbor[3] = {nullptr,nullptr,nullptr};
+};
+typedef vector<QuadTreeNode*> QuadTreeNodePointerArray;
+// quad-tree structure
+struct QuadTree{
+	QuadTreeNode *root = nullptr;
+};
+// control points, lines(polygons) and circles
+struct ControlPoint{
+	int id = -1, max_depth = -1, physic_group = 0;
+	double minimal_resolution = DBL_MAX;
+	Vertex vert; // vert for vertex
+};
+typedef vector<ControlPoint> ControlPointArray;
+
+struct ControlLine{
+	int id = -1, max_depth = -1, physic_group = 0;
+	double minimal_resolution = -1.0;
+	VertexArray vert; // verts for vertices
+};
+typedef vector<ControlLine> ControlLineArray;
+
+struct ControlCircle{
+	int id = -1, max_depth = -1, physic_group = 0;
+	double minimal_resolution = -1.0, circle_cap_degree = -1.0;
+	Vertex circle_center;
+};
+typedef vector<ControlCircle> ControlCircleArray;
+// mathematic functions regarding declared structures
+bool operator ==(Cpoint,Cpoint);
+bool operator !=(Cpoint,Cpoint);
+Cpoint operator +(Cpoint,Cpoint);
+Cpoint operator -(Cpoint,Cpoint);
+Cpoint operator *(double,Cpoint);
+double DotProduct(Cpoint,Cpoint);
+Cpoint CrossProduct(Cpoint,Cpoint);
+Cpoint CloudCenter(VertexArray);
+double ModuleLength(Cpoint);
+double ProjectAngle(Cpoint,Cpoint);
+Cpoint Sphere2Cartesian(Spoint);
+Spoint Cartesian2Sphere(Cpoint);
+Vertex RotateVertex(Vertex,Vertex,Vertex);
+Cpoint LineCrossPlane(Cpoint,Cpoint,Cpoint);
+string GetStringIndex(Vertex);
+int LocalIndex(int,Triangle);
+#endif

src/stt_class.h

@@ -0,0 +1,77 @@
+#ifndef _STT_CLASS_H
+#define _STT_CLASS_H
+#include "head_functions.h"
+#include "struct_functions.h"
+
+class SttGenerator
+{
+public:
+	SttGenerator()
+	{
+		for (int i = 0; i < 20; i++)
+		{
+			forest_[i] = nullptr;
+		}
+	}
+	~SttGenerator(){}
+	int set_command_record(int,char**);
+	int set_tree_depth(char*);
+	int set_pole_equator_radius(char*);
+	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 DeleteTree(QuadTreeNode**);//清空整颗树
+	void ReturnLeaf(QuadTreeNode**);//返回叶子
+	void ReturnDepth(QuadTreeNode**,int);
+	void SortNeighbor(QuadTreeNodePointerArray); // sort neighboring relationship for input triangles
+	void CutOutline(QuadTreeNode**);//切割模型范围 为了保证后续处理中树形结构的完整性 我们只添加node的属性值来控制是否输出该节点
+	void CutHole(QuadTreeNode**); //cut holes in the created STT
+	void CloseSurface(QuadTree**);
+	int ReturnBranchDepth(QuadTreeNode**); //返回当前枝桠的最大深度
+	int InTrianglePoint(QuadTreeNode*);//在球面上判断点和三角形关系的一种方法 直接使用矢量运算确定包含关系 更直接更简单
+	int InTriangleLine(QuadTreeNode*);//判断插入线是否穿过节点三角形 使用的是球面下的方法 直接矢量计算 注意因为球面上的特殊关系 两个点之间的夹角不能大于等于180度 因为球面上总是沿着最短路径走 而且通常我们指的也是最短路径
+	int InTrianglePolygon(QuadTreeNode*);//判断多边形与三角形的关系
+	int InTriangleCircle(QuadTreeNode*);//判断圆与三角形的关系
+	int OutPolyOutline(QuadTreeNode*);//判断多边形与三角形的关系 用于切割模型边界
+	int InPolyOutline(QuadTreeNode*);//判断多边形与三角形的关系 用于切割模型边界 挖洞
+	int OutputMshFile(char*,double,double);
+	int OutputVertexLocation(char*,double,double);
+	int OutputTriangleCenterLocation(char*,double,double);
+	int OutputNeighbor(char*);
+	int GetControlPoint(char*); //读取额外的点
+	int GetControlCircle(char*); //读取额外的圆
+	int GetControlLine(char*,ControlLineArray&); // Get control line arrays
+private:
+	// record input command line options for output records
+	string command_record_;
+	// minimal and maximal depths of quad-tree
+	int tree_depth_, max_depth_;
+	// pole and equator radius of the coordinate reference system
+	double pole_radius_, equator_radius_;
+	// orientation of the top vertex of the icosahedron
+	Vertex icosahedron_orient_;
+	// vertex array of the STT. This array defines the actual shape of the STT
+	VertexArray array_stt_vert_;
+	// map from the vertex's index to vertex. This map is used to find vertex by its index
+	Int2VertexMap map_id_vertex_;
+	Int2VertexMap::iterator ivd_;
+	// map from the vertex's position to vertex. This map is used to find vertex by its position
+	String2VertexMap map_str_vertex_;
+	String2VertexMap::iterator ivm_;
+	// base icosahedron used to construct the STT
+	Icosahedron base_icosahedron_;
+	// 20 quad-trees in which each of them represents the partition of one facet of the base icosahedron
+	QuadTree *forest_[20];
+	// 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)
+	ControlPointArray array_control_point_;
+	ControlCircleArray array_control_circle_;
+	ControlLineArray array_control_line_;
+	ControlLineArray array_control_polygon_;
+	ControlLineArray array_outline_polygon_;
+	ControlLineArray array_hole_polygon_;
+};
+#endif

src/stt_close_surface.cc

@@ -0,0 +1,221 @@
+#include "stt_class.h"
+#include "progress_bar.h"
+
+void SttGenerator::CloseSurface(QuadTree** forest){
+	int breaked,local_id,new_count;
+	int newids[4][3] = {{0,3,5},{1,4,3},{2,5,4},{3,4,5}};
+	Vertex temp_vert;
+	Vertex local_vert[6];
+	Triangle temp_tri,temp_tri2;
+	QuadTree* current_tree;
+
+	ProgressBar *bar = new ProgressBar(max_depth_,"Close surface");
+	//首先我们返回三角形列表再确定当前层的相邻关系
+	for (int d = 0; d < max_depth_; d++){
+		bar->Progressed(d);
+		if (!array_out_tri_pointer_.empty()) array_out_tri_pointer_.clear();
+		for (int t = 0; t < 20; t++){
+			current_tree = *(forest+t);
+			ReturnDepth(&(current_tree->root),d);
+		}
+		SortNeighbor(array_out_tri_pointer_);
+
+		//遍历所有三角形询问
+		for (int t = 0; t < array_out_tri_pointer_.size(); t++){
+			if (array_out_tri_pointer_[t]->children[0] != nullptr) continue; //如果节点三角形已经存在子节点则跳过
+			else{
+				//询问节点三角形的邻居
+				for (int n = 0; n < 3; n++){
+					//邻居必须存在
+					if (array_out_tri_pointer_[t]->neighbor[n] != nullptr){
+						//如果邻居最大深度大于等于三层当前节点的深度 则直接在所有邻居新建三个顶点 完成后跳出邻居循环
+						//这样做的主要原因是如果相邻节点深度差距太大的话无法保证在表面闭合的过程中三角形的质量 同时可以增加缝合带的宽度 有利于后续的建模使用
+						breaked = 0;
+						if (ReturnBranchDepth(&(array_out_tri_pointer_[t]->neighbor[n])) - array_out_tri_pointer_[t]->depth >= 3){
+							for (int q = 0; q < 3; q++){
+								if (array_out_tri_pointer_[t]->neighbor[(n+q)%3] != nullptr){
+									for (int l = 0; l < 3; l++)
+										temp_tri2.ids[l] = array_out_tri_pointer_[t]->neighbor[(n+q)%3]->tri.ids[l];
+									for (int h = 0; h < 3; h++){
+										//计算新顶点坐标
+										temp_vert.posic = 0.5*(array_stt_vert_[temp_tri2.ids[h]].posic + array_stt_vert_[temp_tri2.ids[(h+1)%3]].posic);
+										// map vertex to the reference sphere/ellipsoid
+										temp_vert.posis = Cartesian2Sphere(temp_vert.posic);
+										temp_vert.posis.rad = EllipsoidRadius(temp_vert.posis.lat, pole_radius_, equator_radius_);
+										temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+
+										ivm_ = map_str_vertex_.find(GetStringIndex(temp_vert));
+										//若为新的顶点则将其增加到两个映射和一个链表中
+										if(ivm_ == map_str_vertex_.end()){
+											temp_vert.id = array_stt_vert_.size();//新的顶点索引等于顶点集的数量
+											temp_vert.posis = Cartesian2Sphere(temp_vert.posic);
+											array_stt_vert_.push_back(temp_vert);//将新产生的顶点保存到顶点链表中
+											map_id_vertex_[temp_vert.id] = temp_vert;//将新产生的顶点保存到顶点索引映射中
+											map_str_vertex_[GetStringIndex(temp_vert)] = temp_vert;//将新产生的顶点保存到顶点位置映射中
+										}
+									}
+								}
+							}
+							breaked = 1;
+						}
+						if(breaked) break;
+						//如果邻居节点存在与该节点相邻的孙节点 则可能需要增加新的顶点 如果只存在子节点则不需要增加新的顶点 所以此时我们不考虑子节点的情况
+						//遍历邻居的所有子节点
+						for (int m = 0; m < 4; m++){
+							//要存在可能相邻的孙节点则必须存在子节点
+							if (array_out_tri_pointer_[t]->neighbor[n]->children[m] != nullptr){
+								//如果邻居的一个子节点存在子节点则需要判断这个子节点是否与当前的节点三角形相邻 相邻的条件是邻居的这个存疑的子节点应该与当前节点三角形共顶点
+								if (array_out_tri_pointer_[t]->neighbor[n]->children[m]->children[0] != nullptr){
+									//循环匹配当前节点三角形的三个顶点与存疑的子节点的三个顶点 如果共顶点则寻找当前节点三角形中唯一不与邻居三角形连接的顶点(外点)
+									//并在公共顶点与外顶点的中点尝试新建一个顶点 如果新的顶点不存在则添加 否则继续
+									for (int l = 0; l < 3; l++)
+										temp_tri.ids[l] = array_out_tri_pointer_[t]->neighbor[n]->children[m]->tri.ids[l];
+									for (int l = 0; l < 3; l++)
+										temp_tri2.ids[l] = array_out_tri_pointer_[t]->neighbor[n]->tri.ids[l];
+									breaked = 0;
+									for (int h = 0; h < 3; h++){
+										for (int f = 0; f < 3; f++){
+											//检验顶点索引 找到公共点 
+											if (array_out_tri_pointer_[t]->tri.ids[h] == temp_tri.ids[f]){
+												//找到外点
+												local_id = LocalIndex(temp_tri.ids[f],temp_tri2);
+												if (array_out_tri_pointer_[t]->tri.ids[(h+1)%3] != temp_tri2.ids[(local_id+1)%3] && 
+													array_out_tri_pointer_[t]->tri.ids[(h+1)%3] != temp_tri2.ids[(local_id+2)%3]){
+													//外点是array_out_tri_pointer_[t]->tri->ids[(h+1)%3]
+													temp_vert.posic = 0.5*(array_stt_vert_[array_out_tri_pointer_[t]->tri.ids[h]].posic
+													 + array_stt_vert_[array_out_tri_pointer_[t]->tri.ids[(h+1)%3]].posic);//计算新顶点坐标
+													// map vertex to the reference sphere/ellipsoid
+													temp_vert.posis = Cartesian2Sphere(temp_vert.posic);
+													temp_vert.posis.rad = EllipsoidRadius(temp_vert.posis.lat, pole_radius_, equator_radius_);
+													temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+
+													ivm_ = map_str_vertex_.find(GetStringIndex(temp_vert));
+													//若为新的顶点则将其增加到两个映射和一个链表中
+													if(ivm_ == map_str_vertex_.end()){
+														temp_vert.id = array_stt_vert_.size();//新的顶点索引等于顶点集的数量
+														temp_vert.posis = Cartesian2Sphere(temp_vert.posic);
+														array_stt_vert_.push_back(temp_vert);//将新产生的顶点保存到顶点链表中
+														map_id_vertex_[temp_vert.id] = temp_vert;//将新产生的顶点保存到顶点索引映射中
+														map_str_vertex_[GetStringIndex(temp_vert)] = temp_vert;//将新产生的顶点保存到顶点位置映射中
+													}
+												}
+												else
+												{
+													//外点是array_out_tri_pointer_[t]->tri->ids[(h+2)%3]
+													temp_vert.posic = 0.5*(array_stt_vert_[array_out_tri_pointer_[t]->tri.ids[h]].posic
+													 + array_stt_vert_[array_out_tri_pointer_[t]->tri.ids[(h+2)%3]].posic);//计算新顶点坐标
+													// map vertex to the reference sphere/ellipsoid
+													temp_vert.posis = Cartesian2Sphere(temp_vert.posic);
+													temp_vert.posis.rad = EllipsoidRadius(temp_vert.posis.lat, pole_radius_, equator_radius_);
+													temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+
+													ivm_ = map_str_vertex_.find(GetStringIndex(temp_vert));
+													//若为新的顶点则将其增加到两个映射和一个链表中
+													if(ivm_ == map_str_vertex_.end()){
+														temp_vert.id = array_stt_vert_.size();//新的顶点索引等于顶点集的数量
+														temp_vert.posis = Cartesian2Sphere(temp_vert.posic);
+														array_stt_vert_.push_back(temp_vert);//将新产生的顶点保存到顶点链表中
+														map_id_vertex_[temp_vert.id] = temp_vert;//将新产生的顶点保存到顶点索引映射中
+														map_str_vertex_[GetStringIndex(temp_vert)] = temp_vert;//将新产生的顶点保存到顶点位置映射中
+													}
+												}
+												breaked = 1; break; //跳出循环 并触发二次跳出
+											}
+										}
+										if(breaked) break;
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+		//第二次遍历所有三角形 按照未连接顶点个数补全下一层的三角形
+		for (int t = 0; t < array_out_tri_pointer_.size(); t++){
+			if (array_out_tri_pointer_[t]->children[0] != nullptr) continue; //如果节点三角形已经存在子节点则跳过
+			else{
+				for (int i = 0; i < 3; i++)
+					local_vert[i] = array_stt_vert_[array_out_tri_pointer_[t]->tri.ids[i]];
+				//查询当前节点三角形的各边是否有未连接的顶点并记录
+				new_count = 0;
+				for (int i = 0; i < 3; i++){
+					local_vert[i+3].posic = 0.5*(local_vert[i].posic + local_vert[(i+1)%3].posic);//计算新顶点坐标
+					// map vertex to the reference sphere/ellipsoid
+					local_vert[i+3].posis= Cartesian2Sphere(local_vert[i+3].posic);
+					local_vert[i+3].posis.rad = EllipsoidRadius(local_vert[i+3].posis.lat, pole_radius_, equator_radius_);
+					local_vert[i+3].posic = Sphere2Cartesian(local_vert[i+3].posis);
+
+					ivm_ = map_str_vertex_.find(GetStringIndex(local_vert[i+3]));
+					//若为新的顶点则将其增加到两个映射和一个链表中
+					if(ivm_!=map_str_vertex_.end()){
+						local_vert[i+3].id = ivm_->second.id;
+						new_count++;
+					}
+					else local_vert[i+3].id = -1;
+				}
+				//如果存在三个未连接的顶点则新建四个三角形
+				if (new_count == 3){
+					for (int i = 0; i < 4; i++){
+						array_out_tri_pointer_[t]->children[i] = new QuadTreeNode;
+						array_out_tri_pointer_[t]->children[i]->id = 10*(array_out_tri_pointer_[t]->id) + 1 + i;
+						array_out_tri_pointer_[t]->children[i]->depth = array_out_tri_pointer_[t]->depth + 1;
+						array_out_tri_pointer_[t]->children[i]->tri.ids[0] = local_vert[newids[i][0]].id;
+						array_out_tri_pointer_[t]->children[i]->tri.ids[1] = local_vert[newids[i][1]].id;
+						array_out_tri_pointer_[t]->children[i]->tri.ids[2] = local_vert[newids[i][2]].id;
+					}
+				}
+				else if (new_count == 2){
+					for (int i = 0; i < 3; i++){
+						if (local_vert[i+3].id == -1){
+							array_out_tri_pointer_[t]->children[0] = new QuadTreeNode;
+							array_out_tri_pointer_[t]->children[0]->id = 10*(array_out_tri_pointer_[t]->id) + 1;
+							array_out_tri_pointer_[t]->children[0]->depth = array_out_tri_pointer_[t]->depth + 1;
+							array_out_tri_pointer_[t]->children[0]->tri.ids[0] = local_vert[i].id;
+							array_out_tri_pointer_[t]->children[0]->tri.ids[1] = local_vert[(i+1)%3].id;
+							array_out_tri_pointer_[t]->children[0]->tri.ids[2] = local_vert[(i+1)%3+3].id;
+
+							array_out_tri_pointer_[t]->children[1] = new QuadTreeNode;
+							array_out_tri_pointer_[t]->children[1]->id = 10*(array_out_tri_pointer_[t]->id) + 2;
+							array_out_tri_pointer_[t]->children[1]->depth = array_out_tri_pointer_[t]->depth + 1;
+							array_out_tri_pointer_[t]->children[1]->tri.ids[0] = local_vert[i].id;
+							array_out_tri_pointer_[t]->children[1]->tri.ids[1] = local_vert[(i+1)%3+3].id;
+							array_out_tri_pointer_[t]->children[1]->tri.ids[2] = local_vert[(i+2)%3+3].id;
+
+							array_out_tri_pointer_[t]->children[2] = new QuadTreeNode;
+							array_out_tri_pointer_[t]->children[2]->id = 10*(array_out_tri_pointer_[t]->id) + 3;
+							array_out_tri_pointer_[t]->children[2]->depth = array_out_tri_pointer_[t]->depth + 1;
+							array_out_tri_pointer_[t]->children[2]->tri.ids[0] = local_vert[(i+2)%3].id;
+							array_out_tri_pointer_[t]->children[2]->tri.ids[1] = local_vert[(i+2)%3+3].id;
+							array_out_tri_pointer_[t]->children[2]->tri.ids[2] = local_vert[(i+1)%3+3].id;
+							break;
+						}
+					}
+				}
+				else if (new_count == 1){
+					for (int i = 0; i < 3; i++){
+						if (local_vert[i+3].id != -1){
+							array_out_tri_pointer_[t]->children[0] = new QuadTreeNode;
+							array_out_tri_pointer_[t]->children[0]->id = 10*(array_out_tri_pointer_[t]->id) + 1;
+							array_out_tri_pointer_[t]->children[0]->depth = array_out_tri_pointer_[t]->depth + 1;
+							array_out_tri_pointer_[t]->children[0]->tri.ids[0] = local_vert[i].id;
+							array_out_tri_pointer_[t]->children[0]->tri.ids[1] = local_vert[i+3].id;
+							array_out_tri_pointer_[t]->children[0]->tri.ids[2] = local_vert[(i+2)%3].id;
+
+							array_out_tri_pointer_[t]->children[1] = new QuadTreeNode;
+							array_out_tri_pointer_[t]->children[1]->id = 10*(array_out_tri_pointer_[t]->id) + 2;
+							array_out_tri_pointer_[t]->children[1]->depth = array_out_tri_pointer_[t]->depth + 1;
+							array_out_tri_pointer_[t]->children[1]->tri.ids[0] = local_vert[(i+2)%3].id;
+							array_out_tri_pointer_[t]->children[1]->tri.ids[1] = local_vert[i+3].id;
+							array_out_tri_pointer_[t]->children[1]->tri.ids[2] = local_vert[(i+1)%3].id;
+							break;
+						}
+					}
+				}
+			}
+		}
+	}
+
+	delete bar;
+	return;
+}

src/stt_create_branch.cc

@@ -0,0 +1,61 @@
+#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)
+{
+	Vertex local_vert[6];
+	QuadTreeNode* current_node;
+
+	*node = new QuadTreeNode; //初始化一个新的四叉树节点
+	current_node = *node;
+	current_node->tri.ids[0] = t_ids0;//将上一节点的三角形顶点索引赋值给current_node内的triangle.ids,因此每一层节点实际上都保存了其本身的三角形顶点索引
+	current_node->tri.ids[1] = t_ids1;
+	current_node->tri.ids[2] = t_ids2;
+	current_node->tri.physic_group = phy_group; //继承上层的物理组
+	current_node->id = upper_id*10+order_id;//写入四叉树节点编号
+	current_node->depth = depth;//记录四叉树深度
+
+	//额外生长条件 满足其一即可生长 在局部加密模型的过程中 不同物理组的赋值顺序前后顺序为圈 多边形 线 点
+	if ((depth < tree_depth_ //基本生长条件 所有节点都能达到的深度
+	 || InTriangleCircle(current_node)
+	 || InTrianglePolygon(current_node)
+	 || InTriangleLine(current_node)
+	 || InTrianglePoint(current_node))
+	 && depth < max_depth_) //最大深度限制 所有节点不能超过的深度
+	{
+		ivd_ = map_id_vertex_.find(t_ids0);//利用map_ID映射找到四叉树节点的前三个点，这三个节点是上一层四叉树产生的，必然存在
+		local_vert[0] = ivd_->second;
+		ivd_ = map_id_vertex_.find(t_ids1);
+		local_vert[1] = ivd_->second;
+		ivd_ = map_id_vertex_.find(t_ids2);
+		local_vert[2] = ivd_->second;
+
+		for(int i = 0; i < 3; i++)//循环产生三个新的顶点，位于节点三角形的三条边的中点，给每一个新产生的节点赋索引值与坐标，并保存到顶点链表中
+		{
+			local_vert[i+3].posic = 0.5*(local_vert[i%3].posic+local_vert[(i+1)%3].posic);//计算新顶点坐标,这里需要注意，如果顶点已经存在则需要将顶点索引重置，不增加顶点计数
+			// map vertex to the reference sphere/ellipsoid
+			local_vert[i+3].posis = Cartesian2Sphere(local_vert[i+3].posic);
+			local_vert[i+3].posis.rad = EllipsoidRadius(local_vert[i+3].posis.lat, pole_radius_, equator_radius_);
+			local_vert[i+3].posic = Sphere2Cartesian(local_vert[i+3].posis);
+
+			ivm_ = map_str_vertex_.find(GetStringIndex(local_vert[i+3]));
+			if(ivm_ != map_str_vertex_.end())//利用map_vert查到当前顶点是否存在,这里需要注意，如果顶点已经存在则只需要将顶点索引置为已存在顶点的索引，不增加顶点计数
+			{
+				local_vert[i+3].id = ivm_->second.id;
+			}
+			else//若为新的顶点则将其增加到两个映射和一个链表中
+			{
+				local_vert[i+3].id = array_stt_vert_.size();//新的顶点索引等于顶点集的数量
+				local_vert[i+3].posis = Cartesian2Sphere(local_vert[i+3].posic);
+				array_stt_vert_.push_back(local_vert[i+3]);//将新产生的顶点保存到顶点链表中
+				map_id_vertex_[local_vert[i+3].id] = local_vert[i+3];//将新产生的顶点保存到顶点索引映射中
+				map_str_vertex_[GetStringIndex(local_vert[i+3])] = local_vert[i+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]));
+		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]));
+	}
+	return;
+}

src/stt_create_tree.cc

@@ -0,0 +1,19 @@
+#include "stt_class.h"
+
+void SttGenerator::CreateTree(int tree_id,int t_ids0,int t_ids1,int t_ids2,QuadTree* p_tree){
+	if (max_depth_ == 0){
+		p_tree->root->id = 0;
+		p_tree->root->depth = 0;
+		p_tree->root->tri.ids[0] = t_ids0;//将上一节点的三角形顶点索引赋值给currNode内的triangle.ids,因此每一层节点实际上都保存了其本身的三角形顶点索引
+		p_tree->root->tri.ids[1] = t_ids1;
+		p_tree->root->tri.ids[2] = t_ids2;
+		for(int i=0;i<4;i++)
+		{
+			p_tree->root->children[i] = nullptr;
+		}
+	}
+	else
+	{
+		CreateBranch(0,tree_id,0,t_ids0,t_ids1,t_ids2,0,&(p_tree->root));//以根节点开始创建四叉树
+	}
+}

src/stt_cut_hole.cc

@@ -0,0 +1,26 @@
+#include "stt_class.h"
+
+void SttGenerator::CutHole(QuadTreeNode** p_tree){
+	//切割范围多边形之外的四叉树节点 从深到浅执行
+	QuadTreeNode* current_node = *p_tree;
+	//当前节点是叶子节点 进行判断
+	if (current_node->children[0]==nullptr && current_node->children[1]==nullptr && 
+		current_node->children[2]==nullptr && current_node->children[3]==nullptr){
+		//如果节点三角形在范围多边形之外 删除节点三角形 同时初始化指针
+		if (InPolyOutline(current_node)){
+			current_node->out_ok = false;
+			return;
+		}
+		else return;
+	}
+	else{
+		for (int i = 0; i < 4; i++){
+			//顺序访问当前节点的四个子节点 先处理子节点的情况 当然前提是节点存在
+			if (current_node->children[i] != nullptr){
+				CutHole(&(current_node->children[i]));
+			}
+			else continue;
+		}
+	}
+	return;
+}

src/stt_cut_outline.cc

@@ -0,0 +1,26 @@
+#include "stt_class.h"
+
+void SttGenerator::CutOutline(QuadTreeNode** p_tree){
+	//切割范围多边形之外的四叉树节点 从深到浅执行
+	QuadTreeNode* current_node = *p_tree;
+	//当前节点是叶子节点 进行判断
+	if (current_node->children[0]==nullptr && current_node->children[1]==nullptr && 
+		current_node->children[2]==nullptr && current_node->children[3]==nullptr){
+		//如果节点三角形在范围多边形之外 删除节点三角形 同时初始化指针
+		if (OutPolyOutline(current_node)){
+			current_node->out_ok = false;
+			return;
+		}
+		else return;
+	}
+	else{
+		for (int i = 0; i < 4; i++){
+			//顺序访问当前节点的四个子节点 先处理子节点的情况 当然前提是节点存在
+			if (current_node->children[i] != nullptr){
+				CutOutline(&(current_node->children[i]));
+			}
+			else continue;
+		}
+	}
+	return;
+}

src/stt_delete_tree.cc

@@ -0,0 +1,18 @@
+#include "stt_class.h"
+
+void SttGenerator::DeleteTree(QuadTreeNode **p_tree)
+{
+	QuadTreeNode *current_node = *p_tree;
+
+	if (current_node != nullptr)
+	{
+		for (int i = 0; i < 4; i++)
+		{
+			DeleteTree(&(current_node->children[i]));
+		}
+
+		delete current_node; current_node = nullptr;
+	}
+
+	return;
+}

src/stt_get_control_circle.cc

@@ -0,0 +1,36 @@
+#include "stt_class.h"
+
+int SttGenerator::GetControlCircle(char* filename)
+{
+	double node_eleva;
+	stringstream temp_ss;
+	string temp_str;
+	ControlCircle temp_circle;
+	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);
+				if (temp_ss >> temp_circle.circle_center.posis.lon >> temp_circle.circle_center.posis.lat >> temp_circle.circle_cap_degree
+				 >> temp_circle.max_depth >> temp_circle.minimal_resolution >> temp_circle.physic_group)
+				{
+					if (temp_circle.max_depth < 0) temp_circle.max_depth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
+					if (temp_circle.minimal_resolution < 0) temp_circle.minimal_resolution = -1.0; //这里直接给成-1
+					temp_circle.circle_center.posis.rad = DefaultR;
+					temp_circle.circle_center.id = array_control_circle_.size();
+					temp_circle.circle_center.posic = Sphere2Cartesian(temp_circle.circle_center.posis);
+					array_control_circle_.push_back(temp_circle);
+				}
+			}
+		}
+		infile.close();
+	}
+	return 0;
+}

src/stt_get_control_line.cc

@@ -0,0 +1,51 @@
+#include "stt_class.h"
+
+int SttGenerator::GetControlLine(char* filename,ControlLineArray& return_line_array)
+{
+	int count;
+	stringstream temp_ss;
+	string temp_str;
+	Vertex temp_vert;
+	ControlLine temp_line;
+	ifstream infile;
+
+	if (!strcmp(filename,"NULL")){
+		return 0;
+	}
+
+	//clear return array
+	if (!return_line_array.empty()){
+		for (int i = 0; i < return_line_array.size(); i++){
+			if (!return_line_array[i].vert.empty()) return_line_array[i].vert.clear();
+		}
+		return_line_array.clear();
+	}
+
+	if (OpenInfile(infile,filename)) return -1;
+	else{
+		while (getline(infile,temp_str)){
+			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
+			else{
+				if (!temp_line.vert.empty()) temp_line.vert.clear();
+				temp_ss = Str2Ss(temp_str);
+				temp_ss >> count >> temp_line.max_depth >> temp_line.minimal_resolution >> temp_line.physic_group;
+				if (temp_line.max_depth <= 0) temp_line.max_depth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
+				if (temp_line.minimal_resolution <= 0) temp_line.minimal_resolution = -1.0; //这里直接给成-1
+				for (int i = 0; i < count; i++){
+					getline(infile,temp_str);
+					temp_ss = Str2Ss(temp_str);
+					if (temp_ss >> temp_vert.posis.lon >> temp_vert.posis.lat){
+						temp_vert.posis.rad = DefaultR;
+						temp_vert.id = temp_line.vert.size();
+						temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+						temp_line.vert.push_back(temp_vert);
+					}
+				}
+				temp_line.id = return_line_array.size();
+				return_line_array.push_back(temp_line);
+			}
+		}
+		infile.close();
+	}
+	return 0;
+}

src/stt_get_control_point.cc

@@ -0,0 +1,39 @@
+#include "stt_class.h"
+
+int SttGenerator::GetControlPoint(char* filename)
+{
+	stringstream temp_ss;
+	string temp_str;
+	ControlPoint one_point;
+	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);
+				if (temp_ss >> one_point.vert.posis.lon 
+					>> one_point.vert.posis.lat 
+					>> one_point.max_depth 
+					>> one_point.minimal_resolution 
+					>> one_point.physic_group)
+				{
+					if (one_point.max_depth < 0) one_point.max_depth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
+					if (one_point.minimal_resolution < 0) one_point.minimal_resolution = -1.0; //这里直接给成-1
+					one_point.vert.posis.rad = DefaultR;
+					one_point.vert.id = array_control_point_.size();
+					one_point.vert.posic = Sphere2Cartesian(one_point.vert.posis);
+					array_control_point_.push_back(one_point);
+				}
+			}
+		}
+		infile.close();
+	}
+	return 0;
+}

src/stt_in_poly_outline.cc

@@ -0,0 +1,114 @@
+#include "stt_class.h"
+
+int SttGenerator::InPolyOutline(QuadTreeNode* node)
+{
+	//没有范围多边形 直接返回否
+	if (array_hole_polygon_.empty()){
+		return 0;
+	}
+	else{
+		int count, pnum;
+		Cpoint tri_nor;
+		Cpoint lineface_nor, edgeface_nor;
+		Cpoint intersect[2];
+		Cpoint angle_mid;
+		Cpoint polygon_mid;
+		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 < array_hole_polygon_.size(); i++){
+			pnum = array_hole_polygon_[i].vert.size();
+			for (int j = 0; j < array_hole_polygon_[i].vert.size(); j++){
+				//排除球形背面的点
+				if (DotProduct(tri_nor,array_hole_polygon_[i].vert[j].posic) > 0){
+					//多边形节点在当前节点三角形内
+					count = 0;
+					for (int k = 0; k < 3; k++){
+						cross_point = LineCrossPlane(array_stt_vert_[temp_tri.ids[k%3]].posic,tri_nor,array_hole_polygon_[i].vert[j].posic);
+						//依次判断前后两条边与待检测点的外法线是否同向 注意排除从球体背面穿射的情况 全为真则返回真
+						if (DotProduct(tri_nor,
+							CrossProduct(array_stt_vert_[temp_tri.ids[(k+1)%3]].posic-array_stt_vert_[temp_tri.ids[k%3]].posic,
+							 cross_point-array_stt_vert_[temp_tri.ids[k%3]].posic)) > 0)
+						{
+							count++;
+						}
+					}
+					//全部在左侧 多边形顶点至少有一个在节点三角形内 即节点三角形至少有一个顶点在多边形内 返回真
+					if (count == 3) return 1;
+				}
+				//多边形边与当前节点三角形相交
+				lineface_nor = CrossProduct(array_hole_polygon_[i].vert[j%pnum].posic,array_hole_polygon_[i].vert[(j+1)%pnum].posic);
+				angle_mid = 0.5*(array_hole_polygon_[i].vert[j%pnum].posic + array_hole_polygon_[i].vert[(j+1)%pnum].posic);
+				for (int n = 0; n < 3; n++){
+					edgeface_nor = CrossProduct(array_stt_vert_[temp_tri.ids[n%3]].posic,array_stt_vert_[temp_tri.ids[(n+1)%3]].posic);
+					//排除两个扇面在同一个平面的情况
+					if (fabs(DotProduct(lineface_nor,edgeface_nor))/(ModuleLength(lineface_nor)*ModuleLength(edgeface_nor)) != 1.0){
+						//两个扇面可能的交点矢量垂直于两个扇面的外法线矢量 正反两个矢量
+						intersect[0] = CrossProduct(lineface_nor,edgeface_nor);
+						intersect[1] = CrossProduct(edgeface_nor,lineface_nor);
+						for (int k = 0; k < 2; k++){
+							//交点矢量在两个线段端点矢量之间 注意端点先后顺序决定了大圆弧在球面上的范围 注意这里同样有从背面穿透的可能 因为我们不确定intersect中哪一个是我们想要的
+							//注意计算叉乘的时候 我们总是会走一个角度小于180的方向
+							//排除与angle_mid相反的半球上所有的三角形
+							if (DotProduct(CrossProduct(intersect[k],array_stt_vert_[temp_tri.ids[n%3]].posic),CrossProduct(intersect[k],array_stt_vert_[temp_tri.ids[(n+1)%3]].posic)) < 0
+							 && DotProduct(CrossProduct(intersect[k],array_hole_polygon_[i].vert[j%pnum].posic),CrossProduct(intersect[k],array_hole_polygon_[i].vert[(j+1)%pnum].posic)) < 0
+							 && DotProduct(angle_mid,tri_nor) > 0)
+							{
+								//多边形边与节点三角形相交 即节点三角形至少有一个顶点在多边形内 返回真
+								return 1;
+							}
+						}
+					}
+				}
+			}
+		}
+		//多边形的顶点和边与当前节点三角形不相交或者包含 判断三角形是否在多边形内
+		for (int i = 0; i < array_hole_polygon_.size(); i++){
+			pnum = array_hole_polygon_[i].vert.size();
+			polygon_mid = CloudCenter(array_hole_polygon_[i].vert);
+			//依次判断节点三角形的三条边与多边形边的交点个数
+			for (int k = 0; k < 3; k++){
+				count = 0;
+				//计算三角形边与球心的平面的法线矢量 只要任意一条边在多边形内 则三角形在多边形内
+				edgeface_nor = CrossProduct(array_stt_vert_[temp_tri.ids[(k)%3]].posic,array_stt_vert_[temp_tri.ids[(k+1)%3]].posic);
+				for (int j = 0; j < array_hole_polygon_[i].vert.size(); j++){
+					//多边形边与当前节点三角形相交
+					lineface_nor = CrossProduct(array_hole_polygon_[i].vert[j%pnum].posic,array_hole_polygon_[i].vert[(j+1)%pnum].posic);
+					angle_mid = 0.5*(array_hole_polygon_[i].vert[j%pnum].posic + array_hole_polygon_[i].vert[(j+1)%pnum].posic);
+					//排除两个扇面在同一个平面的情况
+					if (fabs(DotProduct(lineface_nor,edgeface_nor))/(ModuleLength(lineface_nor)*ModuleLength(edgeface_nor)) != 1.0){
+						//两个扇面可能的交点矢量垂直于两个扇面的外法线矢量 正反两个矢量
+						intersect[0] = CrossProduct(lineface_nor,edgeface_nor);
+						intersect[1] = CrossProduct(edgeface_nor,lineface_nor);
+						for (int n = 0; n < 2; n++){
+							/*注意 这里我们只判断交点是否在线段之间 或者一个点上 这里选择第一个点也可以选择第二点 但只能包含一个 不判断是不是在边之间
+							注意我们只统计从当前三角形顶点向左或者向右旋转180度中遇到的交点*/
+							//交点矢量在两个线段端点矢量之间 注意端点先后顺序决定了大圆弧在球面上的范围 注意这里同样有从背面穿透的可能 因为我们不确定intersect中哪一个是我们想要的
+							//注意计算叉乘的时候 我们总是会走一个角度小于180的方向
+							//排除与angle_mid相反的半球上所有的三角形
+							if (DotProduct(polygon_mid,tri_nor) > 0 //排除位于球背面的三角形
+								&& (DotProduct(CrossProduct(intersect[k],array_hole_polygon_[i].vert[j%pnum].posic),CrossProduct(intersect[k],array_hole_polygon_[i].vert[(j+1)%pnum].posic)) < 0
+								|| array_hole_polygon_[i].vert[j].posic == intersect[n]) //排除与多边形的边不相交的三角形边的延长线 这里包含了一个等于条件 即交点刚好在多边形的顶点上
+								&& DotProduct(angle_mid,intersect[n]) > 0 //排除位于球背面的多边形边与三角形边延长线的交点
+								&& DotProduct(edgeface_nor,CrossProduct(array_stt_vert_[temp_tri.ids[k%3]].posic,intersect[n])) > 0) //只取三角形边其中一则的延长线
+							{
+								count++;
+							}
+						}
+					}
+				}
+				//交点个数为奇数 边在多边形内 返回真
+				if (pow(-1,count) < 0) return 1;
+			}
+		}
+		//全不为真 返回假
+		return 0;
+	}
+}

src/stt_in_triangle_circle.cc

@@ -0,0 +1,41 @@
+#include "stt_class.h"
+
+int SttGenerator::InTriangleCircle(QuadTreeNode* node){
+	// If no constraint circle, return false
+	if (array_control_circle_.empty()){
+		return 0;
+	}
+	else{
+		int node_depth;
+		double node_resolution, center_degree;
+
+		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;
+
+		for (int i = 0; i < array_control_circle_.size(); i++){
+			for (int j = 0; j < 3; j++){
+				// calculate the geocentric angle between a vertex of the triangle and the center of the circle
+				center_degree = acos(DotProduct(array_control_circle_[i].circle_center.posic,array_stt_vert_[temp_tri.ids[j]].posic)
+					/(ModuleLength(array_control_circle_[i].circle_center.posic)*ModuleLength(array_stt_vert_[temp_tri.ids[j]].posic)))*180.0/Pi;
+
+				if (center_degree <= array_control_circle_[i].circle_cap_degree
+					&& array_control_circle_[i].max_depth >= node_depth && node_resolution >= array_control_circle_[i].minimal_resolution){
+					node->tri.physic_group = array_control_circle_[i].physic_group;
+					return 1;
+				}
+			}
+		}
+		return 0;
+	}
+}

src/stt_in_triangle_line.cc

@@ -0,0 +1,85 @@
+#include "stt_class.h"
+
+int SttGenerator::InTriangleLine(QuadTreeNode* node){
+	// If no constraint line, return false
+	if (array_control_line_.empty()){
+		return 0;
+	}
+	else{
+		int count, node_depth;
+		double node_resolution;
+		Cpoint tri_nor;
+		Cpoint lineface_nor, edgeface_nor;
+		Cpoint intersect[2];
+		Cpoint angle_mid;
+		Cpoint edge_mid;
+		Cpoint cross_point;
+
+		Triangle temp_tri;
+		for (int i = 0; i < 3; i++){
+			temp_tri.ids[i] = node->tri.ids[i];
+		}
+
+		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;
+
+		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 < array_control_line_.size(); i++){
+			for (int j = 0; j < array_control_line_[i].vert.size(); j++){
+				if (DotProduct(tri_nor,array_control_line_[i].vert[j].posic) > 0){
+					count = 0;
+					for (int k = 0; k < 3; k++){
+						cross_point = LineCrossPlane(array_stt_vert_[temp_tri.ids[k]].posic, tri_nor, array_control_line_[i].vert[j].posic);
+						if (DotProduct(tri_nor,
+							CrossProduct(array_stt_vert_[temp_tri.ids[(k+1)%3]].posic - array_stt_vert_[temp_tri.ids[k]].posic,
+							 cross_point - array_stt_vert_[temp_tri.ids[k]].posic)) > 0){
+							count++;
+						}
+					}
+
+					if (count == 3 && array_control_line_[i].max_depth >= node_depth && node_resolution >= array_control_line_[i].minimal_resolution){
+						node->tri.physic_group = array_control_line_[i].physic_group;
+						return 1;
+					}
+				}
+			}
+		}
+
+		for (int i = 0; i < array_control_line_.size(); i++){
+			for (int j = 0; j < array_control_line_[i].vert.size()-1; j++){
+				lineface_nor = CrossProduct(array_control_line_[i].vert[j].posic, array_control_line_[i].vert[j+1].posic);
+				angle_mid = 0.5*(array_control_line_[i].vert[j].posic + array_control_line_[i].vert[j+1].posic);
+				for (int n = 0; n < 3; n++){
+					edgeface_nor = CrossProduct(array_stt_vert_[temp_tri.ids[n]].posic, array_stt_vert_[temp_tri.ids[(n+1)%3]].posic);
+					edge_mid = 0.5*(array_stt_vert_[temp_tri.ids[n]].posic + array_stt_vert_[temp_tri.ids[(n+1)%3]].posic);
+
+					if (fabs(DotProduct(lineface_nor,edgeface_nor))/(ModuleLength(lineface_nor) * ModuleLength(edgeface_nor)) != 1.0){
+						intersect[0] = CrossProduct(lineface_nor,edgeface_nor);
+						intersect[1] = CrossProduct(edgeface_nor,lineface_nor);
+						for (int k = 0; k < 2; k++){
+							if (DotProduct(angle_mid,tri_nor) > 0
+							 && DotProduct(CrossProduct(intersect[k],array_stt_vert_[temp_tri.ids[n]].posic),CrossProduct(intersect[k],array_stt_vert_[temp_tri.ids[(n+1)%3]].posic)) < 0
+							 && DotProduct(CrossProduct(intersect[k],array_control_line_[i].vert[j].posic),CrossProduct(intersect[k],array_control_line_[i].vert[j+1].posic)) < 0
+							 && DotProduct(intersect[k],angle_mid) > 0
+							 && DotProduct(intersect[k],edge_mid) > 0
+							 && array_control_line_[i].max_depth >= node_depth
+							 && node_resolution >= array_control_line_[i].minimal_resolution){
+								node->tri.physic_group = array_control_line_[i].physic_group;
+								return 1;
+							}
+						}
+					}
+				}
+			}
+		}
+		return 0;
+	}
+}

src/stt_in_triangle_point.cc

@@ -0,0 +1,50 @@
+#include "stt_class.h"
+
+int SttGenerator::InTrianglePoint(QuadTreeNode* node){
+	//没有插入的点位 直接返回否
+	if (array_control_point_.empty()){
+		return 0;
+	}
+	else{
+		int count, node_depth;
+		double node_resolution;
+		Cpoint tri_nor;
+		Cpoint cross_point;
+
+		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;
+
+		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 < array_control_point_.size(); i++){
+			if (DotProduct(tri_nor, array_control_point_[i].vert.posic) > 0){
+				count = 0;
+				for (int j = 0; j < 3; j++){
+					cross_point = LineCrossPlane(array_stt_vert_[temp_tri.ids[j]].posic, tri_nor, array_control_point_[i].vert.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 && array_control_point_[i].max_depth >= node_depth && node_resolution >= array_control_point_[i].minimal_resolution){
+					node->tri.physic_group = array_control_point_[i].physic_group;
+					return 1;
+				}
+			}
+		}
+		return 0;
+	}
+}

src/stt_in_triangle_polygon.cc

@@ -0,0 +1,159 @@
+#include "stt_class.h"
+
+int SttGenerator::InTrianglePolygon(QuadTreeNode* node){
+	// If no polygon, return false
+	if (array_control_polygon_.empty()){
+		return 0;
+	}
+	else{
+		int count, point_num, node_depth;
+		double node_resolution;
+		Cpoint tri_nor;
+		Cpoint lineface_nor, edgeface_nor;
+		Cpoint intersect[2];
+		Cpoint angle_mid;
+		Cpoint edge_mid;
+		Cpoint polygon_mid;
+		Cpoint cross_point;
+
+		// copy node triangle index to a local variable triangle
+		Triangle node_tri;
+		for (int j = 0; j < 3; j++){
+			node_tri.ids[j] = node->tri.ids[j];
+		}
+		// copy node depth to a local variable node_depth
+		node_depth = node->depth;
+		// calculate the spatial resolution of current node's triangle, which equals to the averaged geocentric angle of the triangle's edges.
+		// Note that using the minimal_resolution constraint may introduce uneven constrains at different parts of the spherical polygon. Use it 
+		// with caution.
+		node_resolution = 0;
+		for (int i = 0; i < 3; i++){
+			node_resolution += acos(DotProduct(array_stt_vert_[node_tri.ids[i]].posic,array_stt_vert_[node_tri.ids[(i+1)%3]].posic)
+				/(ModuleLength(array_stt_vert_[node_tri.ids[i]].posic)*ModuleLength(array_stt_vert_[node_tri.ids[(i+1)%3]].posic)));
+		}
+		node_resolution = node_resolution*60/Pi;
+		// calculate the outside normal vector of current node's triangle. The vector does not need to be normalized.
+		tri_nor = CrossProduct(array_stt_vert_[node_tri.ids[1]].posic-array_stt_vert_[node_tri.ids[0]].posic,
+		 array_stt_vert_[node_tri.ids[2]].posic-array_stt_vert_[node_tri.ids[0]].posic);
+		// To begin with, we deal with triangles that might be located on the edge of a polygon.
+		// loop all polygons
+		for (int i = 0; i < array_control_polygon_.size(); i++)
+		{
+			// vertex number of current polygon
+			point_num = array_control_polygon_[i].vert.size();
+			for (int j = 0; j < point_num; j++)
+			{
+				// Firstly, we determine if any vertex of a polygon locates inside of the node's triangle.
+				// The vertex must face the side as the same as the node's triangle. This condition rule out all vertexes that are located on the back of the sphere
+				if (DotProduct(tri_nor,array_control_polygon_[i].vert[j].posic) > 0)
+				{
+					count = 0;
+					for (int k = 0; k < 3; k++)
+					{
+						cross_point = LineCrossPlane(array_stt_vert_[node_tri.ids[k]].posic, tri_nor, array_control_polygon_[i].vert[j].posic);
+						// See if the intersection point is on the left side of the triangle's edges
+						if (DotProduct(tri_nor,
+							CrossProduct(array_stt_vert_[node_tri.ids[(k+1)%3]].posic-array_stt_vert_[node_tri.ids[k]].posic,
+							 cross_point-array_stt_vert_[node_tri.ids[k]].posic)) > 0){
+							count++;
+						}
+					}
+					// if the intersection point is on the left side of all three edges of the triangle, then the point is located inside the triangle
+					// Meanwhile, the node depth must be small than or equal to the maximal constraint depth of the polygon,
+					// or the node triangle's resolution must be bigger than or equal to the minimal constraint resolution of the polygon.
+					if (count == 3 && array_control_polygon_[i].max_depth >= node_depth && node_resolution >= array_control_polygon_[i].minimal_resolution)
+					{
+						// assign the polygon's physical group to the node's triangle as it is a part of the refined STT
+						node->tri.physic_group = array_control_polygon_[i].physic_group;
+						return 1;
+					}
+				}
+
+				// Secondly, we see if and edges of a polygon is intersected with the node's triangle
+				// calculate the outside normal vector of a triangle composed by two endpoints of a polygon's edge and the origin.
+				lineface_nor = CrossProduct(array_control_polygon_[i].vert[j].posic, array_control_polygon_[i].vert[(j+1)%point_num].posic);
+				// calculate the middle vector of two endpoints of a polygon's edge.
+				angle_mid = 0.5*(array_control_polygon_[i].vert[j].posic + array_control_polygon_[i].vert[(j+1)%point_num].posic);
+				// for edges of the node's triangle
+				for (int n = 0; n < 3; n++)
+				{
+					// calculate the outside normal vector of a triangle composed by two endpoints of a triangle's edge and the origin.
+					edgeface_nor = CrossProduct(array_stt_vert_[node_tri.ids[n]].posic,array_stt_vert_[node_tri.ids[(n+1)%3]].posic);
+					edge_mid = 0.5*(array_stt_vert_[node_tri.ids[n]].posic + array_stt_vert_[node_tri.ids[(n+1)%3]].posic);
+					// exclude the situation that the two edges are parallel to each other.
+					if (fabs(DotProduct(lineface_nor,edgeface_nor))/(ModuleLength(lineface_nor)*ModuleLength(edgeface_nor)) != 1.0)
+					{
+						// two intersection points of the edges are found using the cross product.
+						intersect[0] = CrossProduct(lineface_nor,edgeface_nor);
+						intersect[1] = CrossProduct(edgeface_nor,lineface_nor);
+						for (int k = 0; k < 2; k++)
+						{
+							// The two edges are intersected with each other if one intersection point is located between the polygon's edge and the triangle's edge.
+							// Note that the two edges should be facing the same hemisphere.
+							if (DotProduct(angle_mid, tri_nor) > 0
+							 && DotProduct(CrossProduct(intersect[k],array_stt_vert_[node_tri.ids[n]].posic),CrossProduct(intersect[k],array_stt_vert_[node_tri.ids[(n+1)%3]].posic)) < 0
+							 && DotProduct(CrossProduct(intersect[k],array_control_polygon_[i].vert[j].posic),CrossProduct(intersect[k],array_control_polygon_[i].vert[(j+1)%point_num].posic)) < 0
+							 && DotProduct(intersect[k], angle_mid) > 0
+							 && DotProduct(intersect[k], edge_mid) > 0
+							 && array_control_polygon_[i].max_depth >= node_depth
+							 && node_resolution >= array_control_polygon_[i].minimal_resolution)
+							{
+								node->tri.physic_group = array_control_polygon_[i].physic_group;
+								return 1;
+							}
+						}
+					}
+				}
+			}
+		}
+
+		// Now we found triangles that are inside the polygon.
+		for (int i = 0; i < array_control_polygon_.size(); i++)
+		{
+			point_num = array_control_polygon_[i].vert.size();
+			// find the center direction of the polygon
+			polygon_mid = CloudCenter(array_control_polygon_[i].vert);
+			// calculate number of intersection points for each edge of the triangle against the polygon.
+			// If any edge of the triangle is inside the polygon, then the triangle is inside the polygon.
+			for (int k = 0; k < 3; k++)
+			{
+				count = 0;
+				edgeface_nor = CrossProduct(array_stt_vert_[node_tri.ids[k]].posic,array_stt_vert_[node_tri.ids[(k+1)%3]].posic);
+				for (int j = 0; j < point_num; j++)
+				{
+					lineface_nor = CrossProduct(array_control_polygon_[i].vert[j].posic, array_control_polygon_[i].vert[(j+1)%point_num].posic);
+					angle_mid = 0.5*(array_control_polygon_[i].vert[j].posic + array_control_polygon_[i].vert[(j+1)%point_num].posic);
+					// exclude the situation that the two edges are parallel to each other.
+					if (fabs(DotProduct(lineface_nor,edgeface_nor))/(ModuleLength(lineface_nor)*ModuleLength(edgeface_nor)) != 1.0){
+						intersect[0] = CrossProduct(lineface_nor,edgeface_nor);
+						intersect[1] = CrossProduct(edgeface_nor,lineface_nor);
+						for (int n = 0; n < 2; n++)
+						{
+							// The two edges are intersected with each other if one intersection point is located between the polygon's edge and the triangle's edge.
+							// The two edges should be facing the same hemisphere.
+							// The intersection point is located between an edge of the polygon.
+							// Or the beginning point of an polygon's edge is the intersection point (you can choose to use the ending point too, but you should only use one of them).
+							// The intersection should be facing the same hemisphere as the polygon's edge.
+							// Only one side of the extended trajectory of the triangle's edge should be used.
+							if (DotProduct(polygon_mid,tri_nor) > 0
+								// Fixed an indexing error at the following line (change k to n) on 2021-10-13 by Yi Zhang
+								&& (DotProduct(CrossProduct(intersect[n],array_control_polygon_[i].vert[j].posic),CrossProduct(intersect[n],array_control_polygon_[i].vert[(j+1)%point_num].posic)) < 0
+								|| array_control_polygon_[i].vert[j].posic == intersect[n])
+								&& DotProduct(angle_mid,intersect[n]) > 0
+								&& DotProduct(edgeface_nor,CrossProduct(array_stt_vert_[node_tri.ids[k]].posic,intersect[n])) > 0){
+								count++;
+							}
+						}
+					}
+				}
+				// If the number of intersection is odd. The triangle is inside the polygon.
+				if (pow(-1,count) < 0 && array_control_polygon_[i].max_depth >= node_depth && node_resolution >= array_control_polygon_[i].minimal_resolution){
+					node->tri.physic_group = array_control_polygon_[i].physic_group;
+					return 1;
+				}
+			}
+		}
+		// All fails, the triangle is outside the polygon.
+		return 0;
+	}
+}

src/stt_initial_icosahedron.cc

@@ -0,0 +1,70 @@
+#include "stt_class.h"
+
+/**
+ * 初始化二十面实例ICOSA
+ * @param radius   二十面的顶点半径
+ * @param orient 二十面最高点顶点方位
+ */
+void SttGenerator::InitialIcosahedron(double radius,Vertex orient)
+{
+	double constL, constZ;
+	//计算二十面的十二个顶点位置 首先将顶点位置固定在地理北极点
+	base_icosahedron_.vert[0].id = 0;
+	base_icosahedron_.vert[0].posic.x = 0.0;
+	base_icosahedron_.vert[0].posic.y = 0.0;
+	base_icosahedron_.vert[0].posic.z = radius;
+	base_icosahedron_.vert[0].posis = Cartesian2Sphere(base_icosahedron_.vert[0].posic);
+	base_icosahedron_.vert[11].id = 11;
+	base_icosahedron_.vert[11].posic.x = 0.0;
+	base_icosahedron_.vert[11].posic.y = 0.0;
+	base_icosahedron_.vert[11].posic.z = -1.0*radius;
+	base_icosahedron_.vert[11].posis = Cartesian2Sphere(base_icosahedron_.vert[11].posic);
+
+	constZ = radius*(GoldenMean*GoldenMean - 1.0)/(GoldenMean*GoldenMean + 1.0);
+	constL = radius*sqrt(1.0 - pow((GoldenMean*GoldenMean - 1.0)/(GoldenMean*GoldenMean + 1.0),2));
+
+	for(int i = 1; i <= 5; i++)
+	{
+		base_icosahedron_.vert[i].id = i;
+		base_icosahedron_.vert[i].posic.x = cos(72.0*(i-1)*Pi/180.0)*constL;
+		base_icosahedron_.vert[i].posic.y = sin(72.0*(i-1)*Pi/180.0)*constL;
+		base_icosahedron_.vert[i].posic.z = constZ;
+		base_icosahedron_.vert[i].posis = Cartesian2Sphere(base_icosahedron_.vert[i].posic);
+
+		base_icosahedron_.vert[i+5].id = i+5;
+		base_icosahedron_.vert[i+5].posic.x = cos((72.0*(i-1)+36.0)*Pi/180.0)*constL;
+		base_icosahedron_.vert[i+5].posic.y = sin((72.0*(i-1)+36.0)*Pi/180.0)*constL;
+		base_icosahedron_.vert[i+5].posic.z = -constZ;
+		base_icosahedron_.vert[i+5].posis = Cartesian2Sphere(base_icosahedron_.vert[i+5].posic);
+	}
+
+	//给定二十面的面顶点索引，各个三角面顶点索引按逆时针排序
+	base_icosahedron_.tri[0].ids[0] = 0; base_icosahedron_.tri[0].ids[1] = 1; base_icosahedron_.tri[0].ids[2] = 2;
+	base_icosahedron_.tri[1].ids[0] = 0; base_icosahedron_.tri[1].ids[1] = 2; base_icosahedron_.tri[1].ids[2] = 3;
+	base_icosahedron_.tri[2].ids[0] = 0; base_icosahedron_.tri[2].ids[1] = 3; base_icosahedron_.tri[2].ids[2] = 4;
+	base_icosahedron_.tri[3].ids[0] = 0; base_icosahedron_.tri[3].ids[1] = 4; base_icosahedron_.tri[3].ids[2] = 5;
+	base_icosahedron_.tri[4].ids[0] = 0; base_icosahedron_.tri[4].ids[1] = 5; base_icosahedron_.tri[4].ids[2] = 1;
+	base_icosahedron_.tri[5].ids[0] = 1; base_icosahedron_.tri[5].ids[1] = 6; base_icosahedron_.tri[5].ids[2] = 2;
+	base_icosahedron_.tri[6].ids[0] = 2; base_icosahedron_.tri[6].ids[1] = 6; base_icosahedron_.tri[6].ids[2] = 7;
+	base_icosahedron_.tri[7].ids[0] = 2; base_icosahedron_.tri[7].ids[1] = 7; base_icosahedron_.tri[7].ids[2] = 3;
+	base_icosahedron_.tri[8].ids[0] = 3; base_icosahedron_.tri[8].ids[1] = 7; base_icosahedron_.tri[8].ids[2] = 8;
+	base_icosahedron_.tri[9].ids[0] = 3; base_icosahedron_.tri[9].ids[1] = 8; base_icosahedron_.tri[9].ids[2] = 4;
+	base_icosahedron_.tri[10].ids[0] = 4; base_icosahedron_.tri[10].ids[1] = 8; base_icosahedron_.tri[10].ids[2] = 9;
+	base_icosahedron_.tri[11].ids[0] = 4; base_icosahedron_.tri[11].ids[1] = 9; base_icosahedron_.tri[11].ids[2] = 5;
+	base_icosahedron_.tri[12].ids[0] = 5; base_icosahedron_.tri[12].ids[1] = 9; base_icosahedron_.tri[12].ids[2] = 10;
+	base_icosahedron_.tri[13].ids[0] = 5; base_icosahedron_.tri[13].ids[1] = 10; base_icosahedron_.tri[13].ids[2] = 1;
+	base_icosahedron_.tri[14].ids[0] = 1; base_icosahedron_.tri[14].ids[1] = 10; base_icosahedron_.tri[14].ids[2] = 6;
+	base_icosahedron_.tri[15].ids[0] = 6; base_icosahedron_.tri[15].ids[1] = 11; base_icosahedron_.tri[15].ids[2] = 7;
+	base_icosahedron_.tri[16].ids[0] = 7; base_icosahedron_.tri[16].ids[1] = 11; base_icosahedron_.tri[16].ids[2] = 8;
+	base_icosahedron_.tri[17].ids[0] = 8; base_icosahedron_.tri[17].ids[1] = 11; base_icosahedron_.tri[17].ids[2] = 9;
+	base_icosahedron_.tri[18].ids[0] = 9; base_icosahedron_.tri[18].ids[1] = 11; base_icosahedron_.tri[18].ids[2] = 10;
+	base_icosahedron_.tri[19].ids[0] = 10; base_icosahedron_.tri[19].ids[1] = 11; base_icosahedron_.tri[19].ids[2] = 6;
+
+	//旋转二十面顶点的位置
+	Vertex ref_vert = base_icosahedron_.vert[0]; //注意我们选取的参考点为z轴正方向
+	for (int i = 0; i < 12; i++)
+	{
+		base_icosahedron_.vert[i] = RotateVertex(ref_vert,orient,base_icosahedron_.vert[i]);
+	}
+	return;
+}

src/stt_out_poly_outline.cc

@@ -0,0 +1,112 @@
+#include "stt_class.h"
+
+int SttGenerator::OutPolyOutline(QuadTreeNode* node)
+{
+	//没有范围多边形 直接返回否
+	if (array_outline_polygon_.empty()){
+		return 0;
+	}
+	else{
+		int count, pnum;
+		Cpoint tri_nor;
+		Cpoint lineface_nor, edgeface_nor;
+		Cpoint intersect[2];
+		Cpoint angle_mid;
+		Cpoint polygon_mid;
+		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 < array_outline_polygon_.size(); i++){
+			pnum = array_outline_polygon_[i].vert.size();
+			for (int j = 0; j < array_outline_polygon_[i].vert.size(); j++){
+				//排除球形背面的点
+				if (DotProduct(tri_nor,array_outline_polygon_[i].vert[j].posic) > 0){
+					//多边形节点在当前节点三角形内
+					count = 0;
+					for (int k = 0; k < 3; k++){
+						cross_point = LineCrossPlane(array_stt_vert_[temp_tri.ids[k%3]].posic,tri_nor,array_outline_polygon_[i].vert[j].posic);
+						//依次判断前后两条边与待检测点的外法线是否同向 注意排除从球体背面穿射的情况 全为真则返回真
+						if (DotProduct(tri_nor,
+							CrossProduct(array_stt_vert_[temp_tri.ids[(k+1)%3]].posic-array_stt_vert_[temp_tri.ids[k%3]].posic,
+							 cross_point-array_stt_vert_[temp_tri.ids[k%3]].posic)) > 0){
+							count++;
+						}
+					}
+					//全部在左侧 多边形顶点至少有一个在节点三角形内 即节点三角形至少有一个顶点在多边形内 返回假
+					if (count == 3) return 0;
+				}
+				//多边形边与当前节点三角形相交
+				lineface_nor = CrossProduct(array_outline_polygon_[i].vert[j%pnum].posic,array_outline_polygon_[i].vert[(j+1)%pnum].posic);
+				angle_mid = 0.5*(array_outline_polygon_[i].vert[j%pnum].posic + array_outline_polygon_[i].vert[(j+1)%pnum].posic);
+				for (int n = 0; n < 3; n++){
+					edgeface_nor = CrossProduct(array_stt_vert_[temp_tri.ids[n%3]].posic,array_stt_vert_[temp_tri.ids[(n+1)%3]].posic);
+					//排除两个扇面在同一个平面的情况
+					if (fabs(DotProduct(lineface_nor,edgeface_nor))/(ModuleLength(lineface_nor)*ModuleLength(edgeface_nor)) != 1.0){
+						//两个扇面可能的交点矢量垂直于两个扇面的外法线矢量 正反两个矢量
+						intersect[0] = CrossProduct(lineface_nor,edgeface_nor);
+						intersect[1] = CrossProduct(edgeface_nor,lineface_nor);
+						for (int k = 0; k < 2; k++){
+							//交点矢量在两个线段端点矢量之间 注意端点先后顺序决定了大圆弧在球面上的范围 注意这里同样有从背面穿透的可能 因为我们不确定intersect中哪一个是我们想要的
+							//注意计算叉乘的时候 我们总是会走一个角度小于180的方向
+							//排除与angle_mid相反的半球上所有的三角形
+							if (DotProduct(CrossProduct(intersect[k],array_stt_vert_[temp_tri.ids[n%3]].posic),CrossProduct(intersect[k],array_stt_vert_[temp_tri.ids[(n+1)%3]].posic)) < 0
+							 && DotProduct(CrossProduct(intersect[k],array_outline_polygon_[i].vert[j%pnum].posic),CrossProduct(intersect[k],array_outline_polygon_[i].vert[(j+1)%pnum].posic)) < 0
+							 && DotProduct(angle_mid,tri_nor) > 0){
+								//多边形边与节点三角形相交 即节点三角形至少有一个顶点在多边形内 返回假
+								return 0;
+							}
+						}
+					}
+				}
+			}
+		}
+		//多边形的顶点和边与当前节点三角形不相交或者包含 判断三角形是否在多边形内
+		for (int i = 0; i < array_outline_polygon_.size(); i++){
+			pnum = array_outline_polygon_[i].vert.size();
+			polygon_mid = CloudCenter(array_outline_polygon_[i].vert);
+			//依次判断节点三角形的三条边与多边形边的交点个数
+			for (int k = 0; k < 3; k++){
+				count = 0;
+				//计算三角形边与球心的平面的法线矢量 只要任意一条边在多边形内 则三角形在多边形内
+				edgeface_nor = CrossProduct(array_stt_vert_[temp_tri.ids[(k)%3]].posic,array_stt_vert_[temp_tri.ids[(k+1)%3]].posic);
+				for (int j = 0; j < array_outline_polygon_[i].vert.size(); j++){
+					//多边形边与当前节点三角形相交
+					lineface_nor = CrossProduct(array_outline_polygon_[i].vert[j%pnum].posic,array_outline_polygon_[i].vert[(j+1)%pnum].posic);
+					angle_mid = 0.5*(array_outline_polygon_[i].vert[j%pnum].posic + array_outline_polygon_[i].vert[(j+1)%pnum].posic);
+					//排除两个扇面在同一个平面的情况
+					if (fabs(DotProduct(lineface_nor,edgeface_nor))/(ModuleLength(lineface_nor)*ModuleLength(edgeface_nor)) != 1.0){
+						//两个扇面可能的交点矢量垂直于两个扇面的外法线矢量 正反两个矢量
+						intersect[0] = CrossProduct(lineface_nor,edgeface_nor);
+						intersect[1] = CrossProduct(edgeface_nor,lineface_nor);
+						for (int n = 0; n < 2; n++){
+							/*注意 这里我们只判断交点是否在线段之间 或者一个点上 这里选择第一个点也可以选择第二点 但只能包含一个 不判断是不是在边之间
+							注意我们只统计从当前三角形顶点向左或者向右旋转180度中遇到的交点*/
+							//交点矢量在两个线段端点矢量之间 注意端点先后顺序决定了大圆弧在球面上的范围 注意这里同样有从背面穿透的可能 因为我们不确定intersect中哪一个是我们想要的
+							//注意计算叉乘的时候 我们总是会走一个角度小于180的方向
+							//排除与angle_mid相反的半球上所有的三角形
+							if (DotProduct(polygon_mid,tri_nor) > 0 //排除位于球背面的三角形
+								&& (DotProduct(CrossProduct(intersect[k],array_outline_polygon_[i].vert[j%pnum].posic),CrossProduct(intersect[k],array_outline_polygon_[i].vert[(j+1)%pnum].posic)) < 0
+								|| array_outline_polygon_[i].vert[j].posic == intersect[n]) //排除与多边形的边不相交的三角形边的延长线 这里包含了一个等于条件 即交点刚好在多边形的顶点上
+								&& DotProduct(angle_mid,intersect[n]) > 0 //排除位于球背面的多边形边与三角形边延长线的交点
+								&& DotProduct(edgeface_nor,CrossProduct(array_stt_vert_[temp_tri.ids[k%3]].posic,intersect[n])) > 0) //只取三角形边其中一则的延长线
+							{
+								count++;
+							}
+						}
+					}
+				}
+				//交点个数为奇数 边在多边形内 返回假
+				if (pow(-1,count) < 0) return 0;
+			}
+		}
+		//全不为假 返回真
+		return 1;
+	}
+}

src/stt_output_msh_file.cc

@@ -0,0 +1,58 @@
+#include "stt_class.h"
+
+int SttGenerator::OutputMshFile(char* filename,double pole_radius,double equator_radius)
+{
+	time_t now = time(0);
+	char* dt = ctime(&now);
+	IntArray1D array_vert_id;
+	Int2IntMap map_id_out_id;
+
+	if (!strcmp(filename,"NULL") || !strcmp(filename,""))
+		return -1;
+
+	ofstream outfile;
+	if(OpenOutfile(outfile,filename)) return -1;
+
+	vector<int>::iterator pos;
+	for (int i = 0; i < array_out_tri_pointer_.size(); i++)
+	{
+		for (int j = 0; j < 3; j++)
+		{
+			array_vert_id.push_back(array_out_tri_pointer_[i]->tri.ids[j]);
+		}
+	}
+	sort(array_vert_id.begin(), array_vert_id.end()); //对顶点序列由小到大排序
+	pos = unique(array_vert_id.begin(), array_vert_id.end()); //获取重复序列开始的位置
+	array_vert_id.erase(pos, array_vert_id.end()); //删除重复顶点序列
+	//初始化mapOutId;
+	for (int i = 0; i < array_vert_id.size(); i++)
+	{
+		map_id_out_id[array_vert_id[i]] = i;
+	}
+
+	Vertex temp_vert;
+	outfile << "$Comments" << endl << "This file is created by stt-generator.ex on " << dt;
+	outfile << "Commands: " << command_record_ << endl;
+	outfile << "$EndComments" << endl;
+	outfile << "$MeshFormat" << endl << "2.2 0 8" << endl << "$EndMeshFormat" << endl;
+	outfile << "$Nodes"<< endl << array_vert_id.size() << endl;
+	for (int i = 0; i < array_vert_id.size(); i++)
+	{
+		temp_vert = array_stt_vert_[array_vert_id[i]];
+		//temp_vert.posis.rad = EllipsoidRadius(temp_vert.posis.lat,pole_radius,equator_radius);
+		//temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+		outfile << i << " " << setprecision(16) << temp_vert.posic.x << " " << temp_vert.posic.y << " " << temp_vert.posic.z << endl;
+	}
+	outfile << "$EndNodes" << endl;
+	outfile << "$Elements" << endl << array_out_tri_pointer_.size() << endl;
+	for (int i = 0; i < array_out_tri_pointer_.size(); i++)
+	{
+		outfile<< i <<" 2 1 " << array_out_tri_pointer_[i]->tri.physic_group << " "
+		 << map_id_out_id[array_out_tri_pointer_[i]->tri.ids[0]] << " "
+		  << map_id_out_id[array_out_tri_pointer_[i]->tri.ids[1]] << " "
+		   << map_id_out_id[array_out_tri_pointer_[i]->tri.ids[2]] << endl;
+	}
+	outfile<<"$EndElements"<<endl;
+	outfile.close();
+	return 0;
+}

src/stt_output_neighbor.cc

@@ -0,0 +1,34 @@
+#include "stt_class.h"
+
+int SttGenerator::OutputNeighbor(char* filename)
+{
+	time_t now = time(0);
+	char* dt = ctime(&now);
+	ofstream outfile;
+	Int2IntMap map_out_triangle_id_;
+
+	if (!strcmp(filename,"NULL")) return 0;
+	if (OpenOutfile(outfile,filename)) return -1;
+
+	for (int i = 0; i < array_out_tri_pointer_.size(); i++){
+		map_out_triangle_id_[array_out_tri_pointer_[i]->id] = i;
+	}
+
+	outfile << "# This file is created by stt-generator.ex on " << dt;
+	outfile << "# Commands: " << command_record_ << endl;
+	outfile << "# Triangle number: "<< array_out_tri_pointer_.size() << endl;
+	outfile << "# Invalid number: -1" << endl;
+	outfile << "# triangle_id neighbor_id1 neighbor_id2 neighbor_id3" << endl;
+	for (int i = 0; i < array_out_tri_pointer_.size(); i++){
+		outfile << i;
+		for (int j = 0; j < 3; j++){
+			if (array_out_tri_pointer_[i]->neighbor[j] != nullptr){
+				outfile << " " << map_out_triangle_id_[array_out_tri_pointer_[i]->neighbor[j]->id];
+			}
+			else outfile << " -1";
+		}
+		outfile << endl;
+	}
+	outfile.close();
+	return 0;
+}

src/stt_output_triangle_center_location.cc

@@ -0,0 +1,32 @@
+#include "stt_class.h"
+
+int SttGenerator::OutputTriangleCenterLocation(char* filename,double pole_radius,double equator_radius)
+{
+	time_t now = time(0);
+	char* dt = ctime(&now);
+
+	if (!strcmp(filename,"NULL") || !strcmp(filename,""))
+		return -1;
+
+	ofstream outfile;
+	if(OpenOutfile(outfile,filename)) return -1;
+
+	Vertex temp_vert;
+	outfile << "# This file is created by stt-generator.ex on " << dt;
+	outfile << "# Commands: " << command_record_ << endl;
+	outfile << "# Vertex number: "<< array_out_tri_pointer_.size() << endl;
+	outfile << "# x y z longitude latitude radius (meter)" << endl;
+	for (int i = 0; i < array_out_tri_pointer_.size(); i++)
+	{
+		temp_vert.posic = 1.0/3.0*(array_stt_vert_[array_out_tri_pointer_[i]->tri.ids[0]].posic 
+			+ array_stt_vert_[array_out_tri_pointer_[i]->tri.ids[1]].posic 
+			+ array_stt_vert_[array_out_tri_pointer_[i]->tri.ids[2]].posic);
+		temp_vert.posis = Cartesian2Sphere(temp_vert.posic);
+		temp_vert.posis.rad = EllipsoidRadius(temp_vert.posis.lat,pole_radius,equator_radius);
+		temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+		outfile << setprecision(16) << temp_vert.posic.x << " " << temp_vert.posic.y << " " << temp_vert.posic.z 
+		<< " " << temp_vert.posis.lon << " " << temp_vert.posis.lat << " " << temp_vert.posis.rad << endl;
+	}
+	outfile.close();
+	return 0;
+}

src/stt_output_vertex_location.cc

@@ -0,0 +1,42 @@
+#include "stt_class.h"
+
+int SttGenerator::OutputVertexLocation(char* filename,double pole_radius,double equator_radius)
+{
+	time_t now = time(0);
+	char* dt = ctime(&now);
+	IntArray1D array_vert_id;
+
+	if (!strcmp(filename,"NULL") || !strcmp(filename,""))
+		return -1;
+
+	ofstream outfile;
+	if(OpenOutfile(outfile,filename)) return -1;
+
+	vector<int>::iterator pos;
+	for (int i = 0; i < array_out_tri_pointer_.size(); i++)
+	{
+		for (int j = 0; j < 3; j++)
+		{
+			array_vert_id.push_back(array_out_tri_pointer_[i]->tri.ids[j]);
+		}
+	}
+	sort(array_vert_id.begin(), array_vert_id.end()); //对顶点序列由小到大排序
+	pos = unique(array_vert_id.begin(), array_vert_id.end()); //获取重复序列开始的位置
+	array_vert_id.erase(pos, array_vert_id.end()); //删除重复顶点序列
+
+	Vertex temp_vert;
+	outfile << "# This file is created by stt-generator.ex on " << dt;
+	outfile << "# Commands: " << command_record_ << endl;
+	outfile << "# Vertex number: "<< array_vert_id.size() << endl;
+	outfile << "# x y z longitude latitude radius (meter)" << endl;
+	for (int i = 0; i < array_vert_id.size(); i++)
+	{
+		temp_vert = array_stt_vert_[array_vert_id[i]];
+		//temp_vert.posis.rad = EllipsoidRadius(temp_vert.posis.lat,pole_radius,equator_radius);
+		//temp_vert.posic = Sphere2Cartesian(temp_vert.posis);
+		outfile << setprecision(16) << temp_vert.posic.x << " " << temp_vert.posic.y << " " << temp_vert.posic.z 
+		<< " " << temp_vert.posis.lon << " " << temp_vert.posis.lat << " " << temp_vert.posis.rad << endl;
+	}
+	outfile.close();
+	return 0;
+}

src/stt_return_branch_depth.cc

@@ -0,0 +1,19 @@
+#include "stt_class.h"
+
+int SttGenerator::ReturnBranchDepth(QuadTreeNode** p_tree){
+	int node_depth, max_depth = 0;
+	QuadTreeNode* current_node = *p_tree;
+	if (current_node->children[0]==nullptr && current_node->children[1]==nullptr && 
+		current_node->children[2]==nullptr && current_node->children[3]==nullptr){
+		return current_node->depth;
+	}
+	else{
+		for (int i = 0; i < 4; i++){
+			if (current_node->children[i] != nullptr){
+				node_depth = ReturnBranchDepth(&(current_node->children[i]));
+				if(node_depth > max_depth) max_depth = node_depth;
+			}
+		}
+		return max_depth;
+	}
+}

src/stt_return_depth.cc

@@ -0,0 +1,16 @@
+#include "stt_class.h"
+
+void SttGenerator::ReturnDepth(QuadTreeNode** p_tree,int back_depth){
+	QuadTreeNode* current_node = *p_tree;
+	if (current_node->depth == back_depth && back_depth <= max_depth_){
+		array_out_tri_pointer_.push_back(current_node);
+		return;
+	}
+	else{
+		for (int i = 0; i < 4; i++){
+			if (current_node->children[i] != nullptr)
+				ReturnDepth(&(current_node->children[i]),back_depth);
+		}
+		return;
+	}
+}

src/stt_return_leaf.cc

@@ -0,0 +1,24 @@
+#include "stt_class.h"
+
+void SttGenerator::ReturnLeaf(QuadTreeNode** p_tree){
+	QuadTreeNode* current_node = *p_tree;
+	if (current_node->children[0]==nullptr && current_node->children[1]==nullptr && 
+		current_node->children[2]==nullptr && current_node->children[3]==nullptr && 
+		current_node->out_ok==true){
+		// bug fix 这里将输出的节点的邻居重置为nullptr 否则在输出局部stt的邻居列表时将产生一个bug 即边缘位置的三角形会默认与0号三角形相邻
+		// 这是因为在之前的闭合三角面过程中我们已全部对邻居列表赋值，所以所有节点三角形的邻居列表全不为空，因此一定会输出三个邻居，但对局部
+		// stt而言，因为在输出环节的邻居排序中并不能找到对应的三角形索引，所以会输出一个默认值，即0。
+		current_node->neighbor[0] = nullptr;
+		current_node->neighbor[1] = nullptr;
+		current_node->neighbor[2] = nullptr;
+		array_out_tri_pointer_.push_back(current_node);
+		return;
+	}
+	else{
+		for (int i = 0; i < 4; i++){
+			if (current_node->children[i]!=nullptr)
+				ReturnLeaf(&(current_node->children[i]));
+		}
+		return;
+	}
+}

src/stt_routine.cc

@@ -0,0 +1,101 @@
+#include "stt_class.h"
+#include "progress_bar.h"
+
+int SttGenerator::Routine(char input_options[][1024]){
+	// set values of tree depths, terminate the program if failed
+	if (set_tree_depth(input_options[0])) return -1;
+	// set values of pole_radius_ and equator_radius_, terminate the program if failed
+	if (set_pole_equator_radius(input_options[1])) return -1;
+	// set orientation of the base icosahedron, terminate the program if failed
+	if (set_icosahedron_orient(input_options[2])) return -1;
+	// get extra-control information for control points, lines, polygons and circles. Terminate the program if failed
+	// get outline and hole polygons
+	if (GetControlPoint(input_options[7])) return -1;
+	if (GetControlCircle(input_options[10])) 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;
+	if (GetControlLine(input_options[12],array_hole_polygon_)) return -1;
+
+	// initial spaces for tree root
+	for (int i = 0; i < 20; i++){
+		forest_[i] = new QuadTree;
+		forest_[i]->root = new QuadTreeNode;
+	}
+	// initial the base icosahedron
+	// the radius of the base icosahedron is set to DefaultR
+	InitialIcosahedron(DefaultR,icosahedron_orient_);
+	// map vertex to the reference sphere/ellipsoid
+	for (int i = 0; i < 12; i++){
+		base_icosahedron_.vert[i].posis = Cartesian2Sphere(base_icosahedron_.vert[i].posic);
+		base_icosahedron_.vert[i].posis.rad = EllipsoidRadius(base_icosahedron_.vert[i].posis.lat, pole_radius_, equator_radius_);
+		base_icosahedron_.vert[i].posic = Sphere2Cartesian(base_icosahedron_.vert[i].posis);
+	}
+	// add vertices of the base icosahedron to array_stt_vert_, map_id_vertex_ and map_str_vertex_
+	for (int i = 0; i < 12; i++){
+		array_stt_vert_.push_back(base_icosahedron_.vert[i]);
+		map_id_vertex_[base_icosahedron_.vert[i].id] = base_icosahedron_.vert[i];
+		map_str_vertex_[GetStringIndex(base_icosahedron_.vert[i])] = base_icosahedron_.vert[i];
+	}
+
+	ProgressBar *bar = new ProgressBar(20,"Initialize STT");
+	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]);
+	}
+
+	delete bar;
+
+	// close surface after construction
+	CloseSurface(&forest_[0]);
+
+	// if outline polygon exists
+	if (!array_outline_polygon_.empty()){
+		ProgressBar *bar3 = new ProgressBar(20,"Cut outline");
+		for (int i = 0; i < 20; i++){
+			bar3->Progressed(i);
+			CutOutline(&(forest_[i]->root));
+		}
+		delete bar3;
+	}
+
+	// if hole polygon exists
+	if (!array_hole_polygon_.empty())
+	{
+		ProgressBar *bar4 = new ProgressBar(20,"Cut holes");
+		for (int i = 0; i < 20; i++){
+			bar4->Progressed(i);
+			CutHole(&(forest_[i]->root));
+		}
+		delete bar4;
+	}
+
+	// return leafs and prepare for outputs
+	if (!array_out_tri_pointer_.empty()) array_out_tri_pointer_.clear();
+	for (int i = 0; i < 20; i++)
+		ReturnLeaf(&(forest_[i]->root));
+
+	if (!OutputMshFile(input_options[3],pole_radius_,equator_radius_))
+		clog << "file saved: " << input_options[3] << endl;
+	if (!OutputVertexLocation(input_options[4],pole_radius_,equator_radius_))
+		clog << "file saved: " << input_options[4] << endl;
+	if (!OutputTriangleCenterLocation(input_options[5],pole_radius_,equator_radius_))
+		clog << "file saved: " << input_options[5] << endl;
+	if (strcmp(input_options[6],"NULL")){
+		SortNeighbor(array_out_tri_pointer_);
+		if (!OutputNeighbor(input_options[6]))
+			clog << "file saved: " << input_options[6] << endl;
+	}
+
+	for (int i = 0; i < 20; i++)
+	{
+		DeleteTree(&(forest_[i]->root));
+		if (forest_[i] != nullptr)
+		{
+			delete forest_[i];
+			forest_[i] = nullptr;
+		}
+	}
+	return 0;
+}

src/stt_set_command_record.cc

@@ -0,0 +1,12 @@
+#include "stt_class.h"
+
+int SttGenerator::set_command_record(int argv_num,char** argvs)
+{
+	command_record_ = argvs[0];
+	for (int i = 1; i < argv_num; i++)
+	{
+		command_record_ += " ";
+		command_record_ += argvs[i];
+	}
+	return 0;
+}

src/stt_set_icosahedron_orient.cc

@@ -0,0 +1,19 @@
+#include "stt_class.h"
+
+int SttGenerator::set_icosahedron_orient(char* input_parameter){
+	// if input_parameter is NULL or empty, set icosahedron_orient_ to the positive direction of z-axis
+	if (!strcmp(input_parameter,"NULL") || !strcmp(input_parameter,"")){
+		icosahedron_orient_.posis.lon = 0.0; icosahedron_orient_.posis.lat = 90.0;
+	}
+	// try to read the input parameter as icosahedron_orient_.posis.lon/icosahedron_orient_.posis.lat along with some boundary checks
+	else if (2!=sscanf(input_parameter,"%lf/%lf",&icosahedron_orient_.posis.lon,&icosahedron_orient_.posis.lat) 
+		|| icosahedron_orient_.posis.lon < -180.0 || icosahedron_orient_.posis.lon > 180.0
+		|| icosahedron_orient_.posis.lat < -90.0 || icosahedron_orient_.posis.lat > 90.0){
+		cerr << BOLDRED << "Error ==> " << RESET << "fail to initialize the orient of the base icosahedron: " << input_parameter << endl;
+		return -1;
+	}
+
+	icosahedron_orient_.posis.rad = DefaultR;
+	icosahedron_orient_.posic = Sphere2Cartesian(icosahedron_orient_.posis);
+	return 0;
+}

src/stt_set_pole_equator_radius.cc

@@ -0,0 +1,37 @@
+#include "stt_class.h"
+
+int SttGenerator::set_pole_equator_radius(char* input_parameter){
+	double ratio;
+	// if input_parameter is NULL or empty, set pole_radius_ and equator_radius_ to DefaultR
+	if (!strcmp(input_parameter,"NULL") || !strcmp(input_parameter,"")){
+		pole_radius_ = equator_radius_ = DefaultR;
+	}
+	// use predefined values
+	else if (!strcmp(input_parameter,"WGS84")){
+		pole_radius_ = WGS84_r; equator_radius_ = WGS84_R;
+	}
+	else if (!strcmp(input_parameter,"Earth")){
+		pole_radius_ = equator_radius_ = Earth_r;
+	}
+	else if (!strcmp(input_parameter,"Moon")){
+		pole_radius_ = equator_radius_ = Moon_r;
+	}
+	// first try to read the input parameter as equator_radius_/pole_radius_
+	// note that equator_radius_ and pole_radius_ must be bigger than zero
+	else if (2 != sscanf(input_parameter,"%lf/%lf",&equator_radius_,&pole_radius_) || pole_radius_ <= 0.0 || equator_radius_ <= 0.0){
+		// then try to read it as equator_radius_,ratio in which ratio = pole_radius_/equator_radius_
+		// note the ratio must be bigger than zero but not necessarily smaller than or equal to one
+		// However, for reality account, we set the limit of the ratio as 2.0
+		if (2 == sscanf(input_parameter,"%lf,%lf",&equator_radius_,&ratio) && equator_radius_ > 0.0 && ratio > 0.0 && ratio <= 2.0){
+			pole_radius_ = ratio * equator_radius_;
+		}
+		// all attempts fail, return -1
+		else{
+			cerr << BOLDRED << "Error ==> " << RESET << "fail to initialize the coordinate reference system: " << input_parameter << endl;
+			return -1;
+		}
+	}
+	// test output
+	//clog << equator_radius_ << " " << pole_radius_ << endl;
+	return 0;
+}

src/stt_set_tree_depth.cc

@@ -0,0 +1,15 @@
+#include "stt_class.h"
+
+int SttGenerator::set_tree_depth(char* input_parameter){
+	// if input_parameter is NULL or empty, set tree_depth_ and max_depth_ to zero
+	if (!strcmp(input_parameter,"NULL") || !strcmp(input_parameter,"")){
+		tree_depth_ = max_depth_ = 0;
+	}
+	// try to read the input parameter as tree_depth_/max_depth_ along with some boundary checks
+	else if (2!=sscanf(input_parameter,"%d/%d",&tree_depth_,&max_depth_) 
+		|| tree_depth_ > max_depth_ || max_depth_ < 0){
+		cerr << BOLDRED << "Error ==> " << RESET << "fail to initialize the minimal and maximal quad-tree depths: " << input_parameter << endl;
+		return -1;
+	}
+	return 0;
+}

src/stt_sort_neighbor.cc

@@ -0,0 +1,78 @@
+#include "stt_class.h"
+#include "progress_bar.h"
+
+void SttGenerator::SortNeighbor(QuadTreeNodePointerArray input_pointers)
+{
+	int local_id1,local_id2;
+	IntArray1D vert_index; //当前层的顶点索引列表
+	vector<int>::iterator pos; //整型向量的迭代器
+	IntArray2D vert_neighbor; //对应于当前层顶点索引的相邻三角形索引列表
+	IntArray2D neighbor_index; //三角形的相邻索引
+	Int2IntMap map_vert_id;
+
+	//确定当前层的所有顶点索引列表
+	if (!vert_index.empty()) vert_index.clear();
+	for (int i = 0; i < input_pointers.size(); i++){
+		for (int t = 0; t < 3; t++){
+			vert_index.push_back(input_pointers[i]->tri.ids[t]);
+		}
+	}
+	//去除所有重复点
+	sort(vert_index.begin(),vert_index.end()); //对顶点序列由小到大排序
+	pos = unique(vert_index.begin(),vert_index.end()); //获取重复序列开始的位置
+	vert_index.erase(pos,vert_index.end()); //删除重复点
+	//初始化map_vert_id
+	if (!map_vert_id.empty()) map_vert_id.clear();
+	for (int i = 0; i < vert_index.size(); i++){
+		map_vert_id[vert_index[i]] = i;
+	}
+	//确定与顶点相连的三角形索引列表
+	if (!vert_neighbor.empty()){
+		for (int i = 0; i < vert_neighbor.size(); i++){
+			if (!vert_neighbor[i].empty()) vert_neighbor[i].clear();
+		}
+	}
+	vert_neighbor.resize(vert_index.size());
+	for (int i = 0; i < input_pointers.size(); i++){
+		for (int j = 0; j < 3; j++){
+			vert_neighbor[map_vert_id[input_pointers[i]->tri.ids[j]]].push_back(i);
+		}
+	}
+
+	if (!neighbor_index.empty()){
+		for (int i = 0; i < neighbor_index.size(); i++){
+			if (!neighbor_index[i].empty()) neighbor_index[i].clear();
+		}
+	}
+	neighbor_index.resize(input_pointers.size());
+	for (int i = 0; i < neighbor_index.size(); i++)
+		neighbor_index[i].reserve(3);
+
+	//确定相邻的三角形
+	for (int i = 0; i < vert_neighbor.size(); i++){
+		//顺序取其中两个三角形 size为1时循环直接跳过
+		for (int m = 0; m < vert_neighbor[i].size()-1; m++){
+			local_id1 = LocalIndex(vert_index[i], input_pointers[vert_neighbor[i][m]]->tri);
+			for (int n = m+1; n < vert_neighbor[i].size(); n++){
+				local_id2 = LocalIndex(vert_index[i], input_pointers[vert_neighbor[i][n]]->tri);
+				//使用单边相邻进行判定 避免重复添加
+				if (input_pointers[vert_neighbor[i][m]]->tri.ids[(local_id1+1)%3] == input_pointers[vert_neighbor[i][n]]->tri.ids[(local_id2+2)%3]){
+					//相互添加
+					neighbor_index[vert_neighbor[i][m]].push_back(vert_neighbor[i][n]);
+					neighbor_index[vert_neighbor[i][n]].push_back(vert_neighbor[i][m]);
+				}
+			}
+		}
+	}
+
+
+	//最后为对应的节点添加neighbor
+	for (int i = 0; i < neighbor_index.size(); i++)
+	{
+		for (int j = 0; j < neighbor_index[i].size(); j++)
+		{
+			input_pointers[i]->neighbor[j] = input_pointers[neighbor_index[i][j]];
+		}
+	}
+	return;
+}