stt/archived/QdTree_Icosa.h

#ifndef _QDTREE_ICOSA_H
#define _QDTREE_ICOSA_H
#include "dataStruct.h"
#include "progressBar_imp.h"

class QdTree_Icosa
{
public:
	QdTree_Icosa(){
		refradius = refRadius = defaultR;
		strcpy(referSystem,"NULL");
		strcpy(orientPara,"NULL");
		strcpy(pointFile,"NULL");
		strcpy(lineFile,"NULL");
		strcpy(polyFile,"NULL");
		strcpy(circleFile,"NULL");
		strcpy(outlineFile,"NULL");
		strcpy(holeFile,"NULL");
		strcpy(mshFile,"NULL");
		strcpy(sphFile,"NULL");
		strcpy(locFile,"NULL");
		strcpy(locOutFile,"NULL");
	}
	~QdTree_Icosa(){}

	void runtine(char*);
	int init_refer(char*);
	int set_orient(char*);
	void init_para(double,vertex); //初始化参数
	void create_branch(int,int,int,int,int,int,Qdtree_node**); //创建分枝
	void create_tree(int,int,int,Qdtree*);//创建树
	void delete_tree(Qdtree_node**);//清空整颗树
	void return_leaf(Qdtree_node**);//返回叶子
	void cut_outline(Qdtree_node**);//切割模型范围 为了保证后续处理中树形结构的完整性 我们只添加node的属性值来控制是否输出改节点
	void cut_holes(Qdtree_node**);
	void add_vertex(int,Qdtree_node**); //按照规则添加个额外顶点
	void close_surface(int,Qdtree_node**); //闭合模型表面
	void locating_triangle(vertex,_1iArray&,Qdtree_node**); //判断包含着球坐标位置的最小三角形
	//转到球坐标之后事情变得麻烦了一点点 策略 对于每一个三角平面建立一个新的坐标系 将待测试的点投影到新的坐标系然后判断 呵呵呵呵呵
	//int inTriangleDot_old(Qdtree_node*);//判断插入点是否在节点三角形内 直角坐标系默认使用x y坐标 球坐标默认使用坐标投影
	int inTriangleSpoint(vertex,Qdtree_node*); //判断一个单独点与节点三角形的关系
	int inTriangleDot(Qdtree_node*);//在球面上判断点和三角形关系的一种方法 直接使用矢量运算确定包含关系 更直接更简单
	//注意因为球面上的特殊关系 两个点之间的夹角不能大于等于180度 因为球面上总是沿着最短路径走 而且通常我们指的也是最短路径
	int inTriangleLine(Qdtree_node*);//判断插入线是否穿过节点三角形 使用的是球面下的方法 直接矢量计算
	int inTrianglePoly(Qdtree_node*);//判断多边形与三角形的关系
	int inTriangleCircle(Qdtree_node*);//判断圆与三角形的关系
	int outPolyOutline(Qdtree_node*);//判断多边形与三角形的关系 用于切割模型边界
	int inPolyOutline(Qdtree_node*);//判断多边形与三角形的关系 用于切割模型边界 挖洞
	int out_msh(char*,double,double);
	int out_sph(char*,double,double);
	int out_loc(char*);
	int get_extrapoint(char*); //读取额外的点
	int get_extraline(char*); //读取额外的线段
	int get_extrapoly(char*); //读取额外的多边形
	int get_extracircle(char*); //读取额外的圆
	int get_outline(char*); //读取模型范围文件
	int get_hole(char*); //读取模型范围文件
	int get_loc(char*); //读入一组位置 判断包含着一组位置的最小三角形
	//辅助函数
	void init_icosa(double,vertex); //初始化一个二十面体实例 需要给定一个默认半径值 二十面体顶点的经纬坐标 在init_para函数中调用
	void extract_vertex();//按return_*函数返回的三角形索引 提取顶点向量的子集 在outmsh函数中调用
	double mean_edgeDegree(Qdtree_node*,double);//计算节点三角形的平均边长相对于半径的角度 会使用点的x y z坐标
private:
	int tree_depth; //树深度
	int max_depth; //最大深度
	double refradius,refRadius;

	vertexArray vArray; //顶点向量与映射
	idMap mapId;
	idMap::iterator ivd;
	strMap mapStr;
	strMap::iterator ivm;
	outIdMap mapOutId;
	outIdMap::iterator ioim;

	Icosa ICOSA;
	Qdtree* FOREST[20]; //四叉树集合，每一个二十面体的三角面都确定为一个树根
	vertex orient; //二十面的顶点的朝向 由经纬度坐标给出 半径长度为defaultR 默认的二十面体顶点 第一个顶点 在地理北极点

	_1iArray out_iArray; //输出顶点向量索引
	triangleArray out_tArray; //输出三角形索引向量

	pointArray extrapArray; //额外顶点向量
	lineArray extralArray; //额外折线向量
	polygonArray extrapolyArray; //额外的多边形
	circleArray extracirArray; //额外的圆圈
	//模型范围多边形 我们使用一个多边形来规定需要保留的模型范围
	//如果outlineArray为空的话直接输出整个球面模型
	//而如果模型范围存在的话则删除位于模型范围之外的三角形节点
	//如果一个节点三角形的三个顶点均位于范围多边形之外 则删除之
	polygonArray outlineArray; //模型范围多边形 位于该多边形之外的三角形单元将被删去 从深到浅执行
	polygonArray holeArray; //挖洞多边形

	//输入的球坐标位置
	vertexArray sphLocations; //输入的球坐标位置
	_2iArray locTriangles; //输入球坐标位置所在的最小三角形的列表 一个位置可能被多个三角形包含 即在三角形的边 或 顶点上

	char referSystem[1024];
	char orientPara[1024];
	char pointFile[1024];
	char lineFile[1024];
	char polyFile[1024];
	char circleFile[1024];
	char outlineFile[1024];
	char holeFile[1024];
	char mshFile[1024];
	char sphFile[1024];
	char locFile[1024];
	char locOutFile[1024];
};

void QdTree_Icosa::runtine(char* paras)
{
	if (!strcmp(paras,"NULL"))
	{
		cout << "set reference ellipsoid for output (input NULL to use defaults 1e+5/1e+5). use 0/0 to disable the function" << endl << "==> ";
		cin >> referSystem;
		cout << "set basic quad-tree depth. this is the depth all quad-tree nodes can achieve" << endl << "==> ";
		cin >> tree_depth;
		cout << "set maximal quad-tree depth. this is the depth no quad-tree node can surpass" << endl << "==> ";
		cin >> max_depth;
		cout << "set oriention of the basic icosahedron (lon/lat,input NULL to use defaults which are 0/90)" << endl << "==> ";
		cin >> orientPara;
		cout << "outline polygons locaiton file (input NULL if no)" << endl << "==> ";
		cin >> outlineFile;
		cout << "hole polygons locaiton file (input NULL if no)" << endl << "==> ";
		cin >> holeFile;
		cout << "extral points locaiton file (input NULL if no)" << endl << "==> ";
		cin >> pointFile;
		cout << "extral lines locaiton file (input NULL if no)" << endl << "==> ";
		cin >> lineFile;
		cout << "extral polygons locaiton file (input NULL if no)" << endl << "==> ";
		cin >> polyFile;
		cout << "extral circles locaiton file (input NULL if no)" << endl << "==> ";
		cin >> circleFile;
		cout << "extral spherical coordinates file (input NULL if no)" << endl << "==> ";
		cin >> locFile;
		cout << "set output Gmsh (.msh) file name (input NULL if no)" << endl << "==> ";
		cin >> mshFile;
		cout << "set output xyz (lon lat radius) file name (input NULL if no)" << endl << "==> ";
		cin >> sphFile;
		cout << "set output locations (lon lat radius triangles) file name (input NULL if no)" << endl << "==> ";
		cin >> locOutFile;
	}
	else
	{
		string temp_str;
		ifstream infile;
		if (open_infile(infile,paras)) return;
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#') continue;
			else
			{
				if (1==sscanf(temp_str.c_str(),"basic-depth=%d",&tree_depth)) cout << "basic tree depth:\t" << tree_depth << endl;
				else if (1==sscanf(temp_str.c_str(),"refer-system=%s",referSystem)) cout << "reference system:\t" << referSystem << endl;
				else if (1==sscanf(temp_str.c_str(),"max-depth=%d",&max_depth)) cout << "maximal tree depth:\t" << max_depth << endl;
				else if (1==sscanf(temp_str.c_str(),"orientation=%s",orientPara)) cout << "icosahedron orientation:\t" << orientPara << endl;
				else if (1==sscanf(temp_str.c_str(),"extra-points=%s",pointFile)) cout << "extral points:\t" << pointFile << endl;
				else if (1==sscanf(temp_str.c_str(),"extra-lines=%s",lineFile)) cout << "extral lines:\t" << lineFile << endl;
				else if (1==sscanf(temp_str.c_str(),"extra-polys=%s",polyFile)) cout << "extral polys:\t" << polyFile << endl;
				else if (1==sscanf(temp_str.c_str(),"extra-circles=%s",circleFile)) cout << "extral circles:\t" << circleFile << endl;
				else if (1==sscanf(temp_str.c_str(),"outline-polys=%s",outlineFile)) cout << "outline polys:\t" << outlineFile << endl;
				else if (1==sscanf(temp_str.c_str(),"hole-polys=%s",holeFile)) cout << "hole polys:\t" << holeFile << endl;
				else if (1==sscanf(temp_str.c_str(),"sph-locations=%s",locFile)) cout << "spherical locations:\t" << locFile << endl;
				else if (1==sscanf(temp_str.c_str(),"msh-save=%s",mshFile)) cout << "Gmsh save:\t" << mshFile << endl;
				else if (1==sscanf(temp_str.c_str(),"sph-save=%s",sphFile)) cout << "Xyz save:\t" << sphFile << endl;
				else if (1==sscanf(temp_str.c_str(),"loc-save=%s",locOutFile)) cout << "Xyz save:\t" << locOutFile << endl;
				else
				{
					cout << "Unknown parameter:\t" << temp_str << endl;
					return;
				}
			}
		}
	}

	if(init_refer(referSystem)) return;

	if(!get_outline(outlineFile))
		cout << "reading " << outlineFile << " ... done" << endl;
	if(!get_hole(holeFile))
		cout << "reading " << holeFile << " ... done" << endl;
	if(!get_extrapoint(pointFile))
		cout << "reading " << pointFile << " ... done" << endl;
	if(!get_extraline(lineFile))
		cout << "reading " << lineFile << " ... done" << endl;
	if (!get_extrapoly(polyFile))
		cout << "reading " << polyFile << " ... done" << endl;
	if (!get_extracircle(circleFile))
		cout << "reading " << circleFile << " ... done" << endl;
	if (!get_loc(locFile))
		cout << "reading " << locFile << " ... done" << endl;

	for (int i = 0; i < 20; i++)
	{
		FOREST[i] = new Qdtree;
		FOREST[i]->root = new Qdtree_node;
		FOREST[i]->root->tri = new triangle;
	}

	if (set_orient(orientPara)) return;
	init_para(defaultR,orient);

	ProgressBar *bar = new ProgressBar(20,"create model");
	for (int i = 0; i < 20; i++)
	{
		bar->Progressed(i);
		create_tree(ICOSA.tri[i].ids[0],ICOSA.tri[i].ids[1],ICOSA.tri[i].ids[2],FOREST[i]);
	}

	//如果没有额外的点 线 多边形 则跳过新增节点和封闭表面的操作
	if (strcmp(pointFile,"NULL") || strcmp(lineFile,"NULL") || strcmp(polyFile,"NULL") || strcmp(circleFile,"NULL"))
	{
		ProgressBar *bar2 = new ProgressBar(max_depth-tree_depth+1,"refine model");
		int curr_depth = tree_depth;
		while(curr_depth <= max_depth)
		{
			bar2->Progressed(curr_depth-tree_depth);
			for (int i = 0; i < 20; i++)
			{
				add_vertex(curr_depth,&(FOREST[i]->root));
			}
			for (int i = 0; i < 20; i++)
			{
				close_surface(curr_depth,&(FOREST[i]->root));
			}
			curr_depth++;
		}
	}

	//执行模型范围的剪切 在极点处有bug
	if (strcmp(outlineFile,"NULL")) //如果outlineFile存在 则执行剪裁流程
	{
		ProgressBar *bar3 = new ProgressBar(20,"cut model");
		for (int i = 0; i < 20; i++)
		{
			bar3->Progressed(i);
			cut_outline(&(FOREST[i]->root));
		}
	}

	//执行挖洞 在极点处有bug
	if (strcmp(holeFile,"NULL")) //如果outlineFile存在 则执行剪裁流程
	{
		ProgressBar *bar4 = new ProgressBar(20,"cut holes");
		for (int i = 0; i < 20; i++)
		{
			bar4->Progressed(i);
			cut_holes(&(FOREST[i]->root));
		}
	}

	for (int i = 0; i < 20; i++)
	{
		return_leaf(&(FOREST[i]->root));
	}

	//在确定节点三角形的输出id之后 我们再来确定输入位置所在的最小三角形
	if (strcmp(locFile,"NULL"))
	{
		//初始化locTriangles
		locTriangles.resize(sphLocations.size());
		ProgressBar *bar5 = new ProgressBar(sphLocations.size(),"locating");
		for (int i = 0; i < sphLocations.size(); i++)
		{
			bar5->Progressed(i);
			for (int j = 0; j < 20; j++)
			{
				locating_triangle(sphLocations[i],locTriangles[i],&(FOREST[j]->root));
			}
		}
	}

	if (!out_msh(mshFile,refradius,refRadius))
		cout << "file saved: " << mshFile << endl;
	if (!out_sph(sphFile,refradius,refRadius))
		cout << "file saved: " << sphFile << endl;
	if (!out_loc(locOutFile))
		cout << "file saved: " << locOutFile << endl;

	for (int i = 0; i < 20; i++)
	{
		delete_tree(&(FOREST[i]->root));
		if (FOREST[i] != NULL) delete FOREST[i];
	}
}

int QdTree_Icosa::init_refer(char* para)
{
	if (!strcmp(para,"NULL")) return 0;

	if (2 != sscanf(para,"%lf/%lf",&refradius,&refRadius))
	{
		cout << "reference system initialzation fail " << para << endl;
		return 1;
	}
	return 0;
}

int QdTree_Icosa::set_orient(char* para)
{
	if (!strcmp(para,"NULL"))
	{
		orient.posis.lon = 0;
		orient.posis.lat = 90;
		orient.posis.rad = defaultR;
		orient.set(orient.posis);
		return 0;
	}
	else
	{
		if (2!=sscanf(para,"%lf/%lf",&orient.posis.lon,&orient.posis.lat))
		{
			cout << "oriention initialzation fail " << para << endl;
			return 1;
		}
		else return 0;
	}
}

void QdTree_Icosa::init_para(double R,vertex ori)
{
	stringstream sstemp;
	string vert_id,mid_id;
	//初始化二十面体实例
	init_icosa(R,ori);
	for (int i = 0; i < 12; i++)
	{
		vArray.push_back(ICOSA.vert[i]); //添加二十面体顶点到顶点向量
		mapId[ICOSA.vert[i].id] = ICOSA.vert[i]; //保存顶点id映射
		//插入以坐标位置的字符串为键值的映射
		sstemp.str("");
		sstemp.clear();
		sstemp<<setprecision(16)<<ICOSA.vert[i].posic.x;
		sstemp>>vert_id;
		sstemp.str("");
		sstemp.clear();
		sstemp<<setprecision(16)<<ICOSA.vert[i].posic.y;
		sstemp>>mid_id;
		vert_id = vert_id + " " + mid_id;
		sstemp.str("");
		sstemp.clear();
		sstemp<<setprecision(16)<<ICOSA.vert[i].posic.z;
		sstemp>>mid_id;
		vert_id = vert_id + " " + mid_id;
		mapStr[vert_id] = ICOSA.vert[i];
	}
}

void QdTree_Icosa::create_branch(int upper_id,int order_id,int depth,int t_ids0,int t_ids1,int t_ids2,Qdtree_node** node)
{
	stringstream sstemp;
	string vert_id,mid_id;
	vertex localVert[6];
	Qdtree_node* currNode;

	*node = new Qdtree_node; //初始化一个新的四叉树节点
	currNode = *node;
	currNode->tri->ids[0] = t_ids0;//将上一节点的三角形顶点索引赋值给currNode内的triangle.ids,因此每一层节点实际上都保存了其本身的三角形顶点索引
	currNode->tri->ids[1] = t_ids1;
	currNode->tri->ids[2] = t_ids2;
	currNode->id = upper_id*10+order_id;//写入四叉树节点编号
	currNode->depth = depth;//记录四叉树深度

	//额外生长条件 满足其一即可生长 在局部加密模型的过程中 不同物理组的赋值顺序前后顺序为圈 多边形 线 点
	if ((depth < tree_depth //基本生长条件 所有节点都能达到的深度
	 || inTriangleCircle(currNode)
	 || inTrianglePoly(currNode)
	 || inTriangleLine(currNode)
	 || inTriangleDot(currNode))
	 && depth < max_depth) //最大深度限制 所有节点不能超过的深度
	{
		ivd = mapId.find(t_ids0);//利用map_ID映射找到四叉树节点的前三个点，这三个节点是上一层四叉树产生的，必然存在
		localVert[0] = ivd->second;
		ivd = mapId.find(t_ids1);
		localVert[1] = ivd->second;
		ivd = mapId.find(t_ids2);
		localVert[2] = ivd->second;

		for(int i = 0; i < 3; i++)//循环产生三个新的顶点，位于节点三角形的三条边的中点，给每一个新产生的节点赋索引值与坐标，并保存到顶点链表中
		{
			localVert[i+3].posic = middle_cpoint(localVert[i%3].posic,localVert[(i+1)%3].posic);//计算新顶点坐标,这里需要注意，如果顶点已经存在则需要将顶点索引重置，不增加顶点计数
			sstemp.str("");
			sstemp.clear();
			sstemp<<setprecision(16)<<localVert[i+3].posic.x;
			sstemp>>vert_id;
			sstemp.str("");
			sstemp.clear();
			sstemp<<setprecision(16)<<localVert[i+3].posic.y;
			sstemp>>mid_id;
			vert_id = vert_id + " " + mid_id;
			sstemp.str("");
			sstemp.clear();
			sstemp<<setprecision(16)<<localVert[i+3].posic.z;
			sstemp>>mid_id;
			vert_id = vert_id + " " + mid_id;
			ivm = mapStr.find(vert_id);
			if(ivm != mapStr.end())//利用map_vert查到当前顶点是否存在,这里需要注意，如果顶点已经存在则只需要将顶点索引置为已存在顶点的索引，不增加顶点计数
			{
				localVert[i+3].id = ivm->second.id;
			}
			else//若为新的顶点则将其增加到两个映射和一个链表中
			{
				localVert[i+3].set(vArray.size());//新的顶点索引等于顶点集的数量
				localVert[i+3].set(localVert[i+3].posic);
				vArray.push_back(localVert[i+3]);//将新产生的顶点保存到顶点链表中
				mapId[localVert[i+3].id] = localVert[i+3];//将新产生的顶点保存到顶点索引映射中
				mapStr[vert_id] = localVert[i+3];//将新产生的顶点保存到顶点位置映射中
			}
		}

		create_branch(currNode->id,1,depth+1,localVert[0].id,localVert[3].id,localVert[5].id,&(currNode->children[0]));
		create_branch(currNode->id,2,depth+1,localVert[1].id,localVert[4].id,localVert[3].id,&(currNode->children[1]));
		create_branch(currNode->id,3,depth+1,localVert[2].id,localVert[5].id,localVert[4].id,&(currNode->children[2]));
		create_branch(currNode->id,4,depth+1,localVert[3].id,localVert[4].id,localVert[5].id,&(currNode->children[3]));
	}
	return;
}

void QdTree_Icosa::create_tree(int t_ids0,int t_ids1,int t_ids2,Qdtree* pTree)
{
	if (tree_depth == 0)
	{
		pTree->root->id = 0;
		pTree->root->depth = 0;
		pTree->root->tri->ids[0] = t_ids0;//将上一节点的三角形顶点索引赋值给currNode内的triangle.ids,因此每一层节点实际上都保存了其本身的三角形顶点索引
		pTree->root->tri->ids[1] = t_ids1;
		pTree->root->tri->ids[2] = t_ids2;
		for(int i=0;i<4;i++)
		{
			pTree->root->children[i] = NULL;
		}
	}
	else
	{
		create_branch(0,0,0,t_ids0,t_ids1,t_ids2,&(pTree->root));//以根节点开始创建四叉树
	}
}

void QdTree_Icosa::delete_tree(Qdtree_node** pTree)
{
	Qdtree_node* currNode = *pTree;
	if (currNode == NULL)
	{
		return;
	}
	else
	{
		for (int i = 0; i < 4; i++)
		{
			delete_tree(&(currNode->children[i]));
		}
		delete currNode->tri; //构造函数初始化 一定存在
		delete currNode; //上面已经判断证否 一定存在
		currNode = NULL;
		return;
	}
}

void QdTree_Icosa::return_leaf(Qdtree_node** pTree)
{
	triangle temp_tri;
	Qdtree_node* currNode = *pTree;
	if (currNode->children[0]==NULL && 
		currNode->children[1]==NULL && 
		currNode->children[2]==NULL && 
		currNode->children[3]==NULL &&
		currNode->outOK==true)
	{
		for (int i = 0; i < 3; i++)
		{
			temp_tri.ids[i] = currNode->tri->ids[i];
			temp_tri.physic = currNode->tri->physic;
		}
		currNode->outId = out_tArray.size(); //为outId赋值
		out_tArray.push_back(temp_tri);
		return;
	}
	else
	{
		for (int i = 0; i < 4; i++)
		{
			if (currNode->children[i]!=NULL)
				return_leaf(&(currNode->children[i]));
		}
		return;
	}
}

void QdTree_Icosa::cut_outline(Qdtree_node** pTree)
{
	//切割范围多边形之外的四叉树节点 从深到浅执行
	Qdtree_node* currNode = *pTree;
	//当前节点是叶子节点 进行判断
	if (currNode->children[0]==NULL && 
		currNode->children[1]==NULL && 
		currNode->children[2]==NULL && 
		currNode->children[3]==NULL)
	{
		if (outPolyOutline(currNode)) //如果节点三角形在范围多边形之外 删除节点三角形 同时初始化指针
		{
			currNode->outOK = false;
			return;
		}
		else return;
	}
	else
	{
		for (int i = 0; i < 4; i++)
		{
			//顺序访问当前节点的四个子节点 先处理子节点的情况 当然前提是节点存在
			if (currNode->children[i] != NULL)
			{
				cut_outline(&(currNode->children[i]));
			}
			else continue;
		}
	}
	return;
}

void QdTree_Icosa::cut_holes(Qdtree_node** pTree)
{
	//切割范围多边形之外的四叉树节点 从深到浅执行
	Qdtree_node* currNode = *pTree;
	//当前节点是叶子节点 进行判断
	if (currNode->children[0]==NULL && 
		currNode->children[1]==NULL && 
		currNode->children[2]==NULL && 
		currNode->children[3]==NULL)
	{
		if (inPolyOutline(currNode)) //如果节点三角形在范围多边形之外 删除节点三角形 同时初始化指针
		{
			currNode->outOK = false;
			return;
		}
		else return;
	}
	else
	{
		for (int i = 0; i < 4; i++)
		{
			//顺序访问当前节点的四个子节点 先处理子节点的情况 当然前提是节点存在
			if (currNode->children[i] != NULL)
			{
				cut_holes(&(currNode->children[i]));
			}
			else continue;
		}
	}
	return;
}

void QdTree_Icosa::add_vertex(int vist_depth,Qdtree_node** pTree)
{
	cpoint middle_p;
	cpoint temp_p,temp_p2;
	vertex temp_v;
	stringstream sstemp;
	string vert_id,mid_id;
	Qdtree_node* currNode = *pTree;
	//访问深度为vist_depth的节点 节点深度不足 深入
	if (currNode->depth < vist_depth)
	{
		//访问存在的子节点
		for (int i = 0; i < 4; i++)
		{
			if (currNode->children[i] != NULL)
				add_vertex(vist_depth,&(currNode->children[i]));
			else continue;
		}
	}
	/*来到了目标节点 检查在当前节点三角形每条边的1/4和3/4位置上是否有已知顶点
	如有则检查在相对边的相应位置是否有已知的顶点 若无则添加*/
	else if (currNode->depth == vist_depth)
	{
		for (int i = 0; i < 3; i++)
		{
			middle_p = middle_cpoint(vArray.at(currNode->tri->ids[i%3]).posic,vArray.at(currNode->tri->ids[(i+1)%3]).posic);
			//检查1/4 3/4位置
			for (int n = 0; n < 2; n++)
			{
				temp_p = middle_cpoint(vArray.at(currNode->tri->ids[(i+n)%3]).posic,middle_p);
				sstemp.str("");
				sstemp.clear();
				sstemp<<setprecision(16)<<temp_p.x;
				sstemp>>vert_id;
				sstemp.str("");
				sstemp.clear();
				sstemp<<setprecision(16)<<temp_p.y;
				sstemp>>mid_id;
				vert_id = vert_id + " " + mid_id;
				sstemp.str("");
				sstemp.clear();
				sstemp<<setprecision(16)<<temp_p.z;
				sstemp>>mid_id;
				vert_id = vert_id + " " + mid_id;
				ivm = mapStr.find(vert_id);
				if(ivm!=mapStr.end())//利用map_vert查到1/4 3/4位置点存在 检查对边中点是否存在
				{
					temp_p2 = middle_cpoint(vArray.at(currNode->tri->ids[i%3]).posic,vArray.at(currNode->tri->ids[(i+2)%3]).posic);
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p2.x;
					sstemp>>vert_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p2.y;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p2.z;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					ivm = mapStr.find(vert_id);
					if(ivm==mapStr.end())//利用map_vert查到对边中点不存在 添加一个新的顶点到顶点向量
					{
						temp_v.set(vArray.size());//新的顶点索引等于顶点集的数量
						temp_v.set(temp_p2);
						vArray.push_back(temp_v);//将新产生的顶点保存到顶点链表中
						mapId[temp_v.id] = temp_v;//将新产生的顶点保存到顶点索引映射中
						mapStr[vert_id] = temp_v;//将新产生的顶点保存到顶点位置映射中
					}

					temp_p2 = middle_cpoint(vArray.at(currNode->tri->ids[(i+1)%3]).posic,vArray.at(currNode->tri->ids[(i+2)%3]).posic);
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p2.x;
					sstemp>>vert_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p2.y;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p2.z;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					ivm = mapStr.find(vert_id);
					if(ivm==mapStr.end())//利用map_vert查到对边中点不存在 添加一个新的顶点到顶点向量
					{
						temp_v.set(vArray.size());//新的顶点索引等于顶点集的数量
						temp_v.set(temp_p2);
						vArray.push_back(temp_v);//将新产生的顶点保存到顶点链表中
						mapId[temp_v.id] = temp_v;//将新产生的顶点保存到顶点索引映射中
						mapStr[vert_id] = temp_v;//将新产生的顶点保存到顶点位置映射中
					}
					break;
				}
			}
		}
	}
}

void QdTree_Icosa::close_surface(int vist_depth,Qdtree_node** pTree)
{
	int pointNum;
	cpoint temp_p;
	stringstream sstemp;
	string vert_id,mid_id;
	vertex localVert[6];
	int newids[4][3] = {{0,3,5},{1,4,3},{2,5,4},{3,4,5}};
	Qdtree_node* currNode = *pTree;
	//访问深度为vist_depth的节点 节点深度不足 深入
	if (currNode->depth < vist_depth)
	{
		//访问存在的子节点
		for (int i = 0; i < 4; i++)
		{
			if (currNode->children[i] != NULL)
				close_surface(vist_depth,&(currNode->children[i]));
			else continue;
		}
	}
	/*来到了目标节点 首先检查在当前节点三角形第四个子节点是否存在 若存在则则前三个也存在 不用做操作*/
	else if (currNode->depth == vist_depth)
	{
		if (currNode->children[3] == NULL)
		{
			//循环确定有几个新顶点
			pointNum = 0;
			for (int j = 0; j < 3; j++)
			{
				temp_p = middle_cpoint(vArray.at(currNode->tri->ids[j%3]).posic,vArray.at(currNode->tri->ids[(j+1)%3]).posic);
				sstemp.str("");
				sstemp.clear();
				sstemp<<setprecision(16)<<temp_p.x;
				sstemp>>vert_id;
				sstemp.str("");
				sstemp.clear();
				sstemp<<setprecision(16)<<temp_p.y;
				sstemp>>mid_id;
				vert_id = vert_id + " " + mid_id;
				sstemp.str("");
				sstemp.clear();
				sstemp<<setprecision(16)<<temp_p.z;
				sstemp>>mid_id;
				vert_id = vert_id + " " + mid_id;
				ivm = mapStr.find(vert_id);
				if(ivm!=mapStr.end())//利用map_vert查到1/4 3/4位置点存在 检查对边中点是否存在
				{
					pointNum++;
				}
			}
			//按不同情况生成新的三角形
			if (pointNum == 1)
			{
				for (int j = 0; j < 3; j++)
				{
					temp_p = middle_cpoint(vArray.at(currNode->tri->ids[j%3]).posic,vArray.at(currNode->tri->ids[(j+1)%3]).posic);
					//计算点的位置ID
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p.x;
					sstemp>>vert_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p.y;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p.z;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					ivm = mapStr.find(vert_id);
					if(ivm!=mapStr.end())//利用map_vert查到当前顶点是否存在
					{
						//增加第一个三角形
						currNode->children[0] = new Qdtree_node;
						currNode->children[0]->id = 10*(currNode->id) + 1;
						currNode->children[0]->depth = currNode->depth + 1;
						currNode->children[0]->tri->ids[0] = currNode->tri->ids[j%3];
						currNode->children[0]->tri->ids[1] = ivm->second.id;
						currNode->children[0]->tri->ids[2] = currNode->tri->ids[(j+2)%3];
						//增加第二个三角形
						currNode->children[1] = new Qdtree_node;
						currNode->children[1]->id = 10*(currNode->id) + 2;
						currNode->children[1]->depth = currNode->depth + 1;
						currNode->children[1]->tri->ids[0] = currNode->tri->ids[(j+1)%3];
						currNode->children[1]->tri->ids[1] = currNode->tri->ids[(j+2)%3];
						currNode->children[1]->tri->ids[2] = ivm->second.id;
					}
				}
			}
			else if (pointNum == 2)
			{
				for (int i = 0; i < 3; i++)
				{
					temp_p = middle_cpoint(vArray.at(currNode->tri->ids[i%3]).posic,vArray.at(currNode->tri->ids[(i+1)%3]).posic);
					//计算点的位置ID
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p.x;
					sstemp>>vert_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p.y;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<temp_p.z;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					ivm = mapStr.find(vert_id);
					if(ivm==mapStr.end())//利用map_vert查到当前顶点是否存在
					{
						for(int k = 0; k < 2; k++)//循环产生三个新的顶点，位于节点三角形的三条边的中点，给每一个新产生的节点赋索引值与坐标，并保存到顶点链表中
						{
							localVert[k].posic = middle_cpoint(vArray.at(currNode->tri->ids[(i+k)%3]).posic,vArray.at(currNode->tri->ids[(i+2)%3]).posic);//计算新顶点坐标,这里需要注意，如果顶点已经存在则需要将顶点索引重置，不增加顶点计数
							sstemp.str("");
							sstemp.clear();
							sstemp<<setprecision(16)<<localVert[k].posic.x;
							sstemp>>vert_id;
							sstemp.str("");
							sstemp.clear();
							sstemp<<setprecision(16)<<localVert[k].posic.y;
							sstemp>>mid_id;
							vert_id = vert_id + " " + mid_id;
							sstemp.str("");
							sstemp.clear();
							sstemp<<setprecision(16)<<localVert[k].posic.z;
							sstemp>>mid_id;
							vert_id = vert_id + " " + mid_id;
							ivm = mapStr.find(vert_id);
							localVert[k].id = ivm->second.id;
						}
						currNode->children[0] = new Qdtree_node;
						currNode->children[0]->id = 10*(currNode->id) + 1;
						currNode->children[0]->depth = currNode->depth + 1;
						currNode->children[0]->tri->ids[0] = currNode->tri->ids[i%3];
						currNode->children[0]->tri->ids[1] = currNode->tri->ids[(i+1)%3];
						currNode->children[0]->tri->ids[2] = localVert[1].id;

						currNode->children[1] = new Qdtree_node;
						currNode->children[1]->id = 10*(currNode->id) + 2;
						currNode->children[1]->depth = currNode->depth + 1;
						currNode->children[1]->tri->ids[0] = currNode->tri->ids[i%3];
						currNode->children[1]->tri->ids[1] = localVert[1].id;
						currNode->children[1]->tri->ids[2] = localVert[0].id;

						currNode->children[2] = new Qdtree_node;
						currNode->children[2]->id = 10*(currNode->id) + 3;
						currNode->children[2]->depth = currNode->depth + 1;
						currNode->children[2]->tri->ids[0] = localVert[0].id;
						currNode->children[2]->tri->ids[1] = localVert[1].id;
						currNode->children[2]->tri->ids[2] = currNode->tri->ids[(i+2)%3];
					}
				}
			}
			else if (pointNum == 3)
			{
				for (int k = 0; k < 3; k++)
				{
					ivd = mapId.find(currNode->tri->ids[k]);//利用map_ID映射找到四叉树节点的前三个点，这三个节点是上一层四叉树产生的，必然存在
					localVert[k] = ivd->second;
				}
				for(int k = 0; k < 3; k++)//循环产生三个新的顶点，位于节点三角形的三条边的中点，给每一个新产生的节点赋索引值与坐标，并保存到顶点链表中
				{
					localVert[k+3].posic = middle_cpoint(localVert[k%3].posic,localVert[(k+1)%3].posic);//计算新顶点坐标,这里需要注意，如果顶点已经存在则需要将顶点索引重置，不增加顶点计数
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<localVert[k+3].posic.x;
					sstemp>>vert_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<localVert[k+3].posic.y;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					sstemp.str("");
					sstemp.clear();
					sstemp<<setprecision(16)<<localVert[k+3].posic.z;
					sstemp>>mid_id;
					vert_id = vert_id + " " + mid_id;
					ivm = mapStr.find(vert_id);
					localVert[k+3].id = ivm->second.id;
				}
				for (int i = 0; i < 4; i++)
				{
					currNode->children[i] = new Qdtree_node;
					currNode->children[i]->id = 10*(currNode->id) + 1 + i;
					currNode->children[i]->depth = currNode->depth + 1;
					currNode->children[i]->tri->ids[0] = localVert[newids[i][0]].id;
					currNode->children[i]->tri->ids[1] = localVert[newids[i][1]].id;
					currNode->children[i]->tri->ids[2] = localVert[newids[i][2]].id;
				}
			}
		}
	}
}

void QdTree_Icosa::locating_triangle(vertex sp,_1iArray& belongs,Qdtree_node** pTree)
{
	Qdtree_node* currNode = *pTree;
	//首先判断输入点是否在当前节点三角形内 如果不在则返回
	if (inTriangleSpoint(sp,currNode))
	{
		//如果节点三角形是叶子 则将其outId添加到belongs内 注意outOk应该也为真
		if (currNode->children[0]==NULL && 
		currNode->children[1]==NULL && 
		currNode->children[2]==NULL && 
		currNode->children[3]==NULL &&
		currNode->outOK==true)
		{
			belongs.push_back(currNode->outId);
			return;
		}
		else
		{
			for (int i = 0; i < 4; i++)
			{
				//顺序访问当前节点的四个子节点 先处理子节点的情况 当然前提是节点存在
				if (currNode->children[i] != NULL)
				{
					locating_triangle(sp,belongs,&(currNode->children[i]));
				}
				else continue;
			}
		}
	}
	else return;
}

int QdTree_Icosa::inTriangleSpoint(vertex sp,Qdtree_node* node)
{
	int count;
	cpoint tri_nor;
	cpoint cross_point;
	cpoint edge_nor,edge_nor2;
	triangle temp_tri;

	for (int j = 0; j < 3; j++)
	{
		temp_tri.ids[j] = node->tri->ids[j];
	}

	//计算三角面元外法线矢量
	tri_nor = cross(vArray.at(temp_tri.ids[1]).posic-vArray.at(temp_tri.ids[0]).posic , vArray.at(temp_tri.ids[2]).posic-vArray.at(temp_tri.ids[0]).posic);
	if (dot(tri_nor,sp.posic) > 0)
	{
		//如果sp等于三角形顶点中的一个返回真
		for (int j = 0; j < 3; j++)
		{
			if (cpoint_angle(vArray.at(temp_tri.ids[j]).posic,sp.posic)<ZERO)
				return 1;
		}

		//我们需要判断点是否刚好落在三角形的一条边上
		for (int j = 0; j < 3; j++)
		{
			edge_nor = cross(sp.posic,vArray.at(temp_tri.ids[j]).posic);
			edge_nor2 = cross(vArray.at(temp_tri.ids[(j+1)%3]).posic,sp.posic);
			if (cpoint_angle(edge_nor,edge_nor2)<ZERO)
				return 1;
		}

		count = 0;
		for (int j = 0; j < 3; j++)
		{
			cross_point = lineOnPlane(vArray.at(temp_tri.ids[j]).posic,tri_nor,sp.posic);
			//依次判断前后两条边与待检测点的外法线是否同向 注意排除从球体背面穿射的情况 全为真则返回真
			if (dot(tri_nor,cross(vArray.at(temp_tri.ids[(j+1)%3]).posic-vArray.at(temp_tri.ids[j]).posic , cross_point-vArray.at(temp_tri.ids[j]).posic)) >= 0)
				count++;
		}
		if (count == 3) return 1;
	}
	return 0;
}

int QdTree_Icosa::inTriangleDot(Qdtree_node* node)
{
	if (extrapArray.empty()) //没有插入的点位 直接返回否
	{
		return 0;
	}
	else
	{
		int tempDepth;
		double tempDeg;
		int count;
		cpoint tri_nor;
		cpoint cross_point;
		triangle temp_tri;

		for (int j = 0; j < 3; j++)
		{
			temp_tri.ids[j] = node->tri->ids[j];
		}

		tempDepth = node->depth;
		tempDeg = 0;
		for (int i = 0; i < 3; i++)
		{
			tempDeg += acos(dot(vArray.at(temp_tri.ids[i]).posic,vArray.at(temp_tri.ids[(i+1)%3]).posic)
				/(length_cpoint(vArray.at(temp_tri.ids[i]).posic)*length_cpoint(vArray.at(temp_tri.ids[(i+1)%3]).posic)));
		}
		tempDeg = tempDeg*60/pi;

		//计算三角面元外法线矢量
		tri_nor = cross(vArray.at(temp_tri.ids[1]).posic-vArray.at(temp_tri.ids[0]).posic , vArray.at(temp_tri.ids[2]).posic-vArray.at(temp_tri.ids[0]).posic);
		for (int i = 0; i < extrapArray.size(); i++)
		{
			if (dot(tri_nor,extrapArray.at(i).vert.posic) > 0)
			{
				count = 0;
				for (int j = 0; j < 3; j++)
				{
					cross_point = lineOnPlane(vArray.at(temp_tri.ids[j]).posic,tri_nor,extrapArray.at(i).vert.posic);
					//依次判断前后两条边与待检测点的外法线是否同向 注意排除从球体背面穿射的情况 全为真则返回真
					if (dot(tri_nor,
						cross(vArray.at(temp_tri.ids[(j+1)%3]).posic-vArray.at(temp_tri.ids[j%3]).posic , cross_point-vArray.at(temp_tri.ids[j%3]).posic)) > 0)
					{
						count++;
					}
				}
				//全为真 当前节点深度小于允许的最大值 当前节点三角形角度大于允许的最小值 都为真则返回真
				if (count == 3)
				{
					node->tri->physic = extrapArray.at(i).physic;
					if (extrapArray.at(i).maxDepth >= tempDepth && tempDeg >= extrapArray.at(i).minDeg) return 1;
				}
			}
		}
		return 0;
	}
}

int QdTree_Icosa::inTriangleLine(Qdtree_node* node)
{
	if (extralArray.empty())
	{
		return 0;
	}
	else
	{
		int tempDepth;
		double tempDeg;
		int count;
		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];
		}

		tempDepth = node->depth;
		tempDeg = 0;
		for (int i = 0; i < 3; i++)
		{
			tempDeg += acos(dot(vArray.at(temp_tri.ids[i]).posic,vArray.at(temp_tri.ids[(i+1)%3]).posic)
				/(length_cpoint(vArray.at(temp_tri.ids[i]).posic)*length_cpoint(vArray.at(temp_tri.ids[(i+1)%3]).posic)));
		}
		tempDeg = tempDeg*60/pi;

		//计算三角面元外法线矢量
		tri_nor = cross(vArray.at(temp_tri.ids[1]).posic-vArray.at(temp_tri.ids[0]).posic , vArray.at(temp_tri.ids[2]).posic-vArray.at(temp_tri.ids[0]).posic);
		//遍历所有线段
		//判断线段端点是否在三角形内 如果是 返回真
		for (int i = 0; i < extralArray.size(); i++)
		{
			for (int j = 0; j < extralArray.at(i).vert.size(); j++)
			{
				//判断线段端点是否在三角形内 如果是 返回真
				if (dot(tri_nor,extralArray.at(i).vert.at(j).posic) > 0)
				{
					count = 0;
					for (int k = 0; k < 3; k++)
					{
						cross_point = lineOnPlane(vArray.at(temp_tri.ids[k%3]).posic,tri_nor,extralArray.at(i).vert.at(j).posic);
						//依次判断前后两条边与待检测点的外法线是否同向 注意排除从球体背面穿射的情况 全为真则返回真
						if (dot(tri_nor,
							cross(vArray.at(temp_tri.ids[(k+1)%3]).posic-vArray.at(temp_tri.ids[k%3]).posic , cross_point-vArray.at(temp_tri.ids[k%3]).posic)) > 0)
						{
							count++;
						}
					}
					//全部在左侧 返回真
					if (count == 3)
					{
						node->tri->physic = extralArray.at(i).physic;
						if (extralArray.at(i).maxDepth >= tempDepth && tempDeg >= extralArray.at(i).minDeg) return 1;
					}
				}
			}
		}
		//无端点在三角形内 判断球面大圆弧是否与三角形相交 三角形内任意一条边与大圆弧相交则两者相交
		//先计算线段球心扇面的法线矢量 因为点序不确定所以具体方向不确定 不过没关系 哈哈
		for (int i = 0; i < extralArray.size(); i++)
		{
			for (int j = 0; j < extralArray.at(i).vert.size()-1; j++)
			{
				lineface_nor = cross(extralArray.at(i).vert.at(j).posic,extralArray.at(i).vert.at(j+1).posic);
				angle_mid = middle_cpoint(extralArray.at(i).vert.at(j).posic,extralArray.at(i).vert.at(j+1).posic);
				for (int n = 0; n < 3; n++)
				{
					edgeface_nor = cross(vArray.at(temp_tri.ids[n%3]).posic,vArray.at(temp_tri.ids[(n+1)%3]).posic);
					edge_mid = middle_cpoint(vArray.at(temp_tri.ids[n%3]).posic,vArray.at(temp_tri.ids[(n+1)%3]).posic);
					//排除两个扇面在同一个平面的情况
					if (fabs(dot(lineface_nor,edgeface_nor))/(length_cpoint(lineface_nor)*length_cpoint(edgeface_nor)) != 1.0)
					{
						//两个扇面可能的交点矢量垂直于两个扇面的外法线矢量 正反两个矢量
						intersect[0] = cross(lineface_nor,edgeface_nor);
						intersect[1] = cross(edgeface_nor,lineface_nor);
						for (int k = 0; k < 2; k++)
						{
							//交点矢量在两个线段端点矢量之间 注意端点先后顺序决定了大圆弧在球面上的范围 注意这里同样有从背面穿透的可能 因为我们不确定intersect中哪一个是我们想要的
							//注意计算叉乘的时候 我们总是会走一个角度小于180的方向
							//排除与angle_mid相反的半球上所有的三角形
							if (dot(angle_mid,tri_nor) > 0
							 && dot(cross(intersect[k],vArray.at(temp_tri.ids[n%3]).posic),cross(intersect[k],vArray.at(temp_tri.ids[(n+1)%3]).posic)) < 0
							 && dot(cross(intersect[k],extralArray.at(i).vert.at(j).posic),cross(intersect[k],extralArray.at(i).vert.at(j+1).posic)) < 0
							 && dot(intersect[k],angle_mid) > 0
							 && dot(intersect[k],edge_mid) > 0)
							{
								node->tri->physic = extralArray.at(i).physic;
								if (extralArray.at(i).maxDepth >= tempDepth && tempDeg >= extralArray.at(i).minDeg) return 1;
							}
						}
					}
				}
			}
		}
		return 0;
	}
}

int QdTree_Icosa::inTrianglePoly(Qdtree_node* node)
{
	if (extrapolyArray.empty())
	{
		return 0;
	}
	else
	{
		int tempDepth;
		double tempDeg;
		int count;
		int 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];
		}

		tempDepth = node->depth;
		tempDeg = 0;
		for (int i = 0; i < 3; i++)
		{
			tempDeg += acos(dot(vArray.at(temp_tri.ids[i]).posic,vArray.at(temp_tri.ids[(i+1)%3]).posic)
				/(length_cpoint(vArray.at(temp_tri.ids[i]).posic)*length_cpoint(vArray.at(temp_tri.ids[(i+1)%3]).posic)));
		}
		tempDeg = tempDeg*60/pi;

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

int QdTree_Icosa::inTriangleCircle(Qdtree_node* node)
{
	if (extracirArray.empty()) //没有插入的圆圈 直接返回否
	{
		return 0;
	}
	else
	{
		int tempDepth;
		double tempDeg;
		double cenDeg;
		triangle temp_tri;

		for (int j = 0; j < 3; j++)
		{
			temp_tri.ids[j] = node->tri->ids[j];
		}

		tempDepth = node->depth;
		tempDeg = 0;
		for (int i = 0; i < 3; i++)
		{
			tempDeg += acos(dot(vArray.at(temp_tri.ids[i]).posic,vArray.at(temp_tri.ids[(i+1)%3]).posic)
				/(length_cpoint(vArray.at(temp_tri.ids[i]).posic)*length_cpoint(vArray.at(temp_tri.ids[(i+1)%3]).posic)));
		}
		tempDeg = tempDeg*60/pi;

		//计算节点三角形顶点与圆圈中心点的球心夹角
		for (int i = 0; i < extracirArray.size(); i++)
		{
			for (int j = 0; j < 3; j++)
			{
				cenDeg = acos(dot(extracirArray.at(i).cen.posic,vArray.at(temp_tri.ids[j]).posic)
					/(length_cpoint(extracirArray.at(i).cen.posic)*length_cpoint(vArray.at(temp_tri.ids[j]).posic)))*180/pi;
				if (cenDeg <= extracirArray.at(i).radDeg)
				{
					node->tri->physic = extracirArray.at(i).physic;
					if (extracirArray.at(i).maxDepth >= tempDepth && tempDeg >= extracirArray.at(i).minDeg) return 1;
				}
			}
		}
		return 0;
	}
}

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

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

int QdTree_Icosa::out_msh(char* filename,double refr,double refR)
{
	time_t now = time(0);
	char* dt = ctime(&now);
	vertex v;
	ofstream outfile;

	if (!strcmp(filename,"NULL"))
	{
		return -1;
	}

	if(open_outfile(outfile,filename)) return -1;

	extract_vertex();

	outfile << "$Comments" << endl << "this file is created by Qdtree_icosa.ex on " << dt;
	outfile << "basic quad-tree depth: " << tree_depth << endl;
	outfile << "maximal quad-tree depth: " << max_depth << endl;
	outfile << "outline: " << outlineFile << endl;
	outfile << "hole:" << holeFile << endl;
	outfile << "extral points: " << pointFile << endl;
	outfile << "extral lines: " << lineFile << endl;
	outfile << "extral polygons: " << polyFile << endl;
	outfile << "extral circles: " << circleFile << endl << "$EndComments" << endl;

	outfile<<"$MeshFormat"<<endl<<"2.2 0 8"<<endl<<"$EndMeshFormat"<<endl<<"$Nodes"<<endl<<out_iArray.size()<<endl;
	for (int i = 0; i < out_iArray.size(); i++)
	{
		if (refr == 0 && refR == 0)
		{
			outfile << i << " " << setprecision(16) << vArray.at(out_iArray.at(i)).posic.x << " " << vArray.at(out_iArray.at(i)).posic.y << " " << vArray.at(out_iArray.at(i)).posic.z << endl;
		}
		else if (refr != 0 && refR != 0)
		{
			v = vArray.at(out_iArray.at(i));
			v.posic = rescale_cpoint(vArray.at(out_iArray.at(i)).posic, REF_r(v.posis.lat,refr,refR));
			v.set(v.posic);
			outfile << i << " " << setprecision(16) << v.posic.x << " " << v.posic.y << " " << v.posic.z << endl;
		}
		else
		{
			cout << "error ==> reference system fail" << endl;
			return -1;
		}
	}
	outfile<<"$EndNodes"<<endl<<"$Elements"<<endl<<out_tArray.size()<<endl;
	for (int i = 0; i < out_tArray.size(); i++)
	{
		outfile<< i <<" 2 1 " << out_tArray.at(i).physic << " " << mapOutId[out_tArray.at(i).ids[0]] << " " << mapOutId[out_tArray.at(i).ids[1]] << " " << mapOutId[out_tArray.at(i).ids[2]] << endl;
	}
	outfile<<"$EndElements"<<endl;

	outfile.close();
	return 0;
}

int QdTree_Icosa::out_sph(char* filename,double refr,double refR)
{
	time_t now = time(0);
	char* dt = ctime(&now);
	vertex v;
	ofstream outfile;

	if (!strcmp(filename,"NULL"))
	{
		return -1;
	}

	if(open_outfile(outfile,filename)) return -1;

	extract_vertex();

	outfile << "# this file is created by Qdtree_icosa.ex on " << dt;
	outfile << "# basic quad-tree depth: " << tree_depth << endl;
	outfile << "# maximal quad-tree depth: " << max_depth << endl;
	outfile << "# outline: " << outlineFile << endl;
	outfile << "# hole: " << holeFile << endl;
	outfile << "# extral points: " << pointFile << endl;
	outfile << "# extral lines: " << lineFile << endl;
	outfile << "# extral polygons: " << polyFile << endl;
	outfile << "# extral circles: " << circleFile << endl;
	//outfile << "# lon(deg) lat(deg) rad(m)" << endl;
	outfile << "# node position lon(deg) lat(deg)" << endl;

	for (int i = 0; i < out_iArray.size(); i++)
	{
		if (refr == 0 && refR == 0) //not used
		{
			outfile << setprecision(16) << vArray.at(out_iArray.at(i)).posis.lon << " " << vArray.at(out_iArray.at(i)).posis.lat << " " << vArray.at(out_iArray.at(i)).posis.rad << endl;
		}
		else if (refr != 0 && refR != 0)
		{
			v = vArray.at(out_iArray.at(i));
			v.posic = rescale_cpoint(vArray.at(out_iArray.at(i)).posic, REF_r(v.posis.lat,refr,refR));
			v.set(v.posic);
			outfile << setprecision(16) << v.posis.lon << " " << v.posis.lat << " " << v.posis.rad << endl;
			//outfile << setprecision(16) << v.posis.lon << " " << v.posis.lat << endl;
		}
		else
		{
			cout << "error ==> reference system fail" << endl;
			return -1;
		}
	}

	outfile.close();
	return 0;
}

int QdTree_Icosa::out_loc(char* filename)
{
	time_t now = time(0);
	char* dt = ctime(&now);
	ofstream outfile;

	if (!strcmp(filename,"NULL"))
	{
		return -1;
	}

	if(open_outfile(outfile,filename)) return -1;

	outfile << "# this file is created by Qdtree_icosa.ex on " << dt;
	outfile << "# basic quad-tree depth: " << tree_depth << endl;
	outfile << "# maximal quad-tree depth: " << max_depth << endl;
	outfile << "# outline: " << outlineFile << endl;
	outfile << "# hole: " << holeFile << endl;
	outfile << "# extral points: " << pointFile << endl;
	outfile << "# extral lines: " << lineFile << endl;
	outfile << "# extral polygons: " << polyFile << endl;
	outfile << "# extral circles: " << circleFile << endl;
	outfile << "# node position lon(deg) lat(deg) Triangle-ids" << endl;

	for (int i = 0; i < sphLocations.size(); i++)
	{
		outfile << setprecision(16);
		outfile << sphLocations[i].posis.lon << " " << sphLocations[i].posis.lat << " ";
		if (locTriangles[i].empty())
		{
			outfile << "nan" << endl; //对于没有所属三角形的情况 输出-1
		}
		else
		{
			for (int j = 0; j < locTriangles[i].size(); j++)
			{
				outfile << locTriangles[i][j] << " ";
			}
			outfile << endl;
		}
	}
	outfile.close();
	return 0;
}

int QdTree_Icosa::get_extrapoint(char* filename)
{
	double node_eleva;
	stringstream temp_ss;
	string temp_str;
	point temp_v;
	ifstream infile;

	if (!strcmp(filename,"NULL"))
	{
		return 0;
	}

	if (open_infile(infile,filename)) return -1;
	else
	{
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
			else
			{
				temp_ss.str("");
				temp_ss.clear();
				temp_ss << temp_str;
				if (temp_ss >> temp_v.maxDepth >> temp_v.minDeg >> temp_v.physic >> temp_v.vert.posis.lon >> temp_v.vert.posis.lat >> node_eleva)
				{
					if (temp_v.maxDepth < 0) temp_v.maxDepth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
					if (temp_v.minDeg < 0) temp_v.minDeg = -1; //这里直接给成-1
					temp_v.vert.posis.rad = defaultR + node_eleva;
					temp_v.vert.set(extrapArray.size());
					temp_v.vert.set(temp_v.vert.posis);
					extrapArray.push_back(temp_v);
				}
			}
		}
		infile.close();
	}
	return 0;
}

int QdTree_Icosa::get_extraline(char* filename) //点个数跟点列表
{
	int count;
	double node_eleva;
	stringstream temp_ss;
	string temp_str;
	vertex temp_v;
	line l;
	ifstream infile;

	if (!strcmp(filename,"NULL"))
	{
		return 0;
	}

	if (open_infile(infile,filename)) return -1;
	else
	{
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
			else
			{
				if (!l.vert.empty()) l.vert.clear();
				temp_ss.str("");
				temp_ss.clear();
				temp_ss << temp_str;
				temp_ss >> count >> l.maxDepth >> l.minDeg >> l.physic;
				if (l.maxDepth < 0) l.maxDepth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
				if (l.minDeg < 0) l.minDeg = -1; //这里直接给成-1
				for (int i = 0; i < count; i++)
				{
					getline(infile,temp_str);
					temp_ss.str("");
					temp_ss.clear();
					temp_ss << temp_str;
					if (temp_ss >> temp_v.posis.lon >> temp_v.posis.lat >> node_eleva)
					{
						temp_v.posis.rad = defaultR + node_eleva;
						temp_v.set(l.vert.size());
						temp_v.set(temp_v.posis);
						l.vert.push_back(temp_v);
					}
				}
				l.id = extralArray.size();
				extralArray.push_back(l);
			}
		}
		infile.close();
	}
	return 0;
}

int QdTree_Icosa::get_extrapoly(char* filename)
{
	int count;
	double node_eleva;
	stringstream temp_ss;
	string temp_str;
	vertex temp_v;
	polygon p;
	ifstream infile;

	if (!strcmp(filename,"NULL"))
	{
		return 0;
	}

	if (open_infile(infile,filename)) return -1;
	else
	{
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
			else
			{
				if (!p.vert.empty()) p.vert.clear();
				temp_ss.str("");
				temp_ss.clear();
				temp_ss << temp_str;
				temp_ss >> count >> p.maxDepth >> p.minDeg >> p.physic;
				if (p.maxDepth < 0) p.maxDepth = 100; //这里直接给一个很大的深度值 节点深度一定小于这个值
				if (p.minDeg < 0) p.minDeg = -1.0; //这里直接给成-1
				for (int i = 0; i < count; i++)
				{
					getline(infile,temp_str);
					temp_ss.str("");
					temp_ss.clear();
					temp_ss << temp_str;
					if (temp_ss >> temp_v.posis.lon >> temp_v.posis.lat >> node_eleva)
					{
						temp_v.posis.rad = defaultR + node_eleva;
						temp_v.set(p.vert.size());
						temp_v.set(temp_v.posis);
						p.vert.push_back(temp_v);
					}
				}
				p.id = extrapolyArray.size();
				extrapolyArray.push_back(p);
			}
		}
		infile.close();
	}
	return 0;
}

int QdTree_Icosa::get_extracircle(char* filename)
{
	double node_eleva;
	stringstream temp_ss;
	string temp_str;
	circle temp_v;
	ifstream infile;

	if (!strcmp(filename,"NULL"))
	{
		return 0;
	}

	if (open_infile(infile,filename)) return -1;
	else
	{
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
			else
			{
				temp_ss.str("");
				temp_ss.clear();
				temp_ss << temp_str;
				if (temp_ss >> temp_v.maxDepth >> temp_v.minDeg >> temp_v.physic >> temp_v.radDeg >> temp_v.cen.posis.lon >> temp_v.cen.posis.lat >> node_eleva)
				{
					if (temp_v.maxDepth < 0) temp_v.maxDepth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
					if (temp_v.minDeg < 0) temp_v.minDeg = -1; //这里直接给成-1
					temp_v.cen.posis.rad = defaultR + node_eleva;
					temp_v.cen.set(extracirArray.size());
					temp_v.cen.set(temp_v.cen.posis);
					extracirArray.push_back(temp_v);
				}
			}
		}
		infile.close();
	}
	return 0;
}

int QdTree_Icosa::get_outline(char* filename)
{
	int count;
	double node_eleva;
	stringstream temp_ss;
	string temp_str;
	vertex temp_v;
	polygon p;
	ifstream infile;

	if (!strcmp(filename,"NULL"))
	{
		return 0;
	}

	if (open_infile(infile,filename)) return -1;
	else
	{
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
			else
			{
				if (!p.vert.empty()) p.vert.clear();
				temp_ss.str("");
				temp_ss.clear();
				temp_ss << temp_str;
				//由于范围多边形对模型生成范围作出的是绝对性的约束 因此并不需要制定每个范围多边形的最大深度和最小尺寸等信息
				//这里暂时保留这两个参数 只是他们在后续的过程中并不起任何作用
				temp_ss >> count >> p.maxDepth >> p.minDeg >> p.physic; //后面三个没用
				if (p.maxDepth < 0) p.maxDepth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
				if (p.minDeg < 0) p.minDeg = -1; //这里直接给成-1
				for (int i = 0; i < count; i++)
				{
					getline(infile,temp_str);
					temp_ss.str("");
					temp_ss.clear();
					temp_ss << temp_str;
					if (temp_ss >> temp_v.posis.lon >> temp_v.posis.lat >> node_eleva)
					{
						temp_v.posis.rad = defaultR + node_eleva;
						temp_v.set(p.vert.size());
						temp_v.set(temp_v.posis);
						p.vert.push_back(temp_v);
					}
				}
				p.id = outlineArray.size();
				outlineArray.push_back(p);
			}
		}
		infile.close();
	}
	return 0;
}

int QdTree_Icosa::get_hole(char* filename)
{
	int count;
	double node_eleva;
	stringstream temp_ss;
	string temp_str;
	vertex temp_v;
	polygon p;
	ifstream infile;

	if (!strcmp(filename,"NULL"))
	{
		return 0;
	}

	if (open_infile(infile,filename)) return -1;
	else
	{
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
			else
			{
				if (!p.vert.empty()) p.vert.clear();
				temp_ss.str("");
				temp_ss.clear();
				temp_ss << temp_str;
				//由于范围多边形对模型生成范围作出的是绝对性的约束 因此并不需要制定每个范围多边形的最大深度和最小尺寸等信息
				//这里暂时保留这两个参数 只是他们在后续的过程中并不起任何作用
				temp_ss >> count >> p.maxDepth >> p.minDeg >> p.physic; //后面三个没用
				if (p.maxDepth < 0) p.maxDepth = 1e+3; //这里直接给一个很大的深度值 节点深度一定小于这个值
				if (p.minDeg < 0) p.minDeg = -1; //这里直接给成-1
				for (int i = 0; i < count; i++)
				{
					getline(infile,temp_str);
					temp_ss.str("");
					temp_ss.clear();
					temp_ss << temp_str;
					if (temp_ss >> temp_v.posis.lon >> temp_v.posis.lat >> node_eleva)
					{
						temp_v.posis.rad = defaultR + node_eleva;
						temp_v.set(p.vert.size());
						temp_v.set(temp_v.posis);
						p.vert.push_back(temp_v);
					}
				}
				p.id = holeArray.size();
				holeArray.push_back(p);
			}
		}
		infile.close();
	}
	return 0;
}

int QdTree_Icosa::get_loc(char* filename)
{
	stringstream temp_ss;
	string temp_str;
	vertex temp_vert;
	ifstream infile;

	if (!strcmp(filename,"NULL"))
	{
		return 0;
	}

	if (open_infile(infile,filename)) return -1;
	else
	{
		while (getline(infile,temp_str))
		{
			if (*(temp_str.begin()) == '#' || temp_str == "") continue;
			else
			{
				temp_ss.str("");
				temp_ss.clear();
				temp_ss << temp_str;
				temp_ss >> temp_vert.posis.lon >> temp_vert.posis.lat;
				temp_vert.posis.rad = defaultR;
				temp_vert.set(temp_vert.posis);
				sphLocations.push_back(temp_vert);
			}
		}
		infile.close();
	}
	return 0;
}

/*辅助函数*/
void QdTree_Icosa::init_icosa(double R,vertex ori)
{
	double constL,constZ;
	//计算二十面的十二个顶点位置 首先将顶点位置固定在地理北极点
	ICOSA.vert[0].id = 0;
	ICOSA.vert[0].posic.x = 0.0;
	ICOSA.vert[0].posic.y = 0.0;
	ICOSA.vert[0].posic.z = R;
	ICOSA.vert[0].set(ICOSA.vert[0].posic);
	ICOSA.vert[11].id = 11;
	ICOSA.vert[11].posic.x = 0.0;
	ICOSA.vert[11].posic.y = 0.0;
	ICOSA.vert[11].posic.z = -R;
	ICOSA.vert[11].set(ICOSA.vert[11].posic);
	
	constZ = R*(golden_mean*golden_mean - 1.0)/(golden_mean*golden_mean + 1.0);
	constL = R*sqrt(1.0 - pow((golden_mean*golden_mean - 1.0)/(golden_mean*golden_mean + 1.0),2));

	for(int i=1;i<=5;i++)
	{
		ICOSA.vert[i].id = i;
		ICOSA.vert[i].posic.x = cos(72.0*(i-1)*pi/180.0)*constL;
		ICOSA.vert[i].posic.y = sin(72.0*(i-1)*pi/180.0)*constL;
		ICOSA.vert[i].posic.z = constZ;
		ICOSA.vert[i].set(ICOSA.vert[i].posic);

		ICOSA.vert[i+5].id = i+5;
		ICOSA.vert[i+5].posic.x = cos((72.0*(i-1)+36.0)*pi/180.0)*constL;
		ICOSA.vert[i+5].posic.y = sin((72.0*(i-1)+36.0)*pi/180.0)*constL;
		ICOSA.vert[i+5].posic.z = -constZ;
		ICOSA.vert[i+5].set(ICOSA.vert[i+5].posic);
	}

	//给定二十面的面顶点索引，各个三角面顶点索引按逆时针排序
	ICOSA.tri[0].ids[0] = 0; ICOSA.tri[0].ids[1] = 1; ICOSA.tri[0].ids[2] = 2;
	ICOSA.tri[1].ids[0] = 0; ICOSA.tri[1].ids[1] = 2; ICOSA.tri[1].ids[2] = 3;
	ICOSA.tri[2].ids[0] = 0; ICOSA.tri[2].ids[1] = 3; ICOSA.tri[2].ids[2] = 4;
	ICOSA.tri[3].ids[0] = 0; ICOSA.tri[3].ids[1] = 4; ICOSA.tri[3].ids[2] = 5;
	ICOSA.tri[4].ids[0] = 0; ICOSA.tri[4].ids[1] = 5; ICOSA.tri[4].ids[2] = 1;
	ICOSA.tri[5].ids[0] = 1; ICOSA.tri[5].ids[1] = 6; ICOSA.tri[5].ids[2] = 2;
	ICOSA.tri[6].ids[0] = 2; ICOSA.tri[6].ids[1] = 6; ICOSA.tri[6].ids[2] = 7;
	ICOSA.tri[7].ids[0] = 2; ICOSA.tri[7].ids[1] = 7; ICOSA.tri[7].ids[2] = 3;
	ICOSA.tri[8].ids[0] = 3; ICOSA.tri[8].ids[1] = 7; ICOSA.tri[8].ids[2] = 8;
	ICOSA.tri[9].ids[0] = 3; ICOSA.tri[9].ids[1] = 8; ICOSA.tri[9].ids[2] = 4;
	ICOSA.tri[10].ids[0] = 4; ICOSA.tri[10].ids[1] = 8; ICOSA.tri[10].ids[2] = 9;
	ICOSA.tri[11].ids[0] = 4; ICOSA.tri[11].ids[1] = 9; ICOSA.tri[11].ids[2] = 5;
	ICOSA.tri[12].ids[0] = 5; ICOSA.tri[12].ids[1] = 9; ICOSA.tri[12].ids[2] = 10;
	ICOSA.tri[13].ids[0] = 5; ICOSA.tri[13].ids[1] = 10; ICOSA.tri[13].ids[2] = 1;
	ICOSA.tri[14].ids[0] = 1; ICOSA.tri[14].ids[1] = 10; ICOSA.tri[14].ids[2] = 6;
	ICOSA.tri[15].ids[0] = 6; ICOSA.tri[15].ids[1] = 11; ICOSA.tri[15].ids[2] = 7;
	ICOSA.tri[16].ids[0] = 7; ICOSA.tri[16].ids[1] = 11; ICOSA.tri[16].ids[2] = 8;
	ICOSA.tri[17].ids[0] = 8; ICOSA.tri[17].ids[1] = 11; ICOSA.tri[17].ids[2] = 9;
	ICOSA.tri[18].ids[0] = 9; ICOSA.tri[18].ids[1] = 11; ICOSA.tri[18].ids[2] = 10;
	ICOSA.tri[19].ids[0] = 10; ICOSA.tri[19].ids[1] = 11; ICOSA.tri[19].ids[2] = 6;

	//旋转二十面顶点的位置
	vertex ref_vert = ICOSA.vert[0]; //注意我们选取的参考点为z轴正方向
	for (int i = 0; i < 12; i++)
	{
		ICOSA.vert[i] = rotate_vertex(ref_vert,ori,ICOSA.vert[i]);
	}
	return;
}

void QdTree_Icosa::extract_vertex()
{
	vector<int>::iterator pos;
	for (int i = 0; i < out_tArray.size(); i++)
	{
		for (int j = 0; j < 3; j++)
		{
			out_iArray.push_back(out_tArray.at(i).ids[j]);
		}
	}
	sort(out_iArray.begin(),out_iArray.end()); //对顶点序列由小到大排序
	pos = unique(out_iArray.begin(), out_iArray.end()); //获取重复序列开始的位置
	out_iArray.erase(pos, out_iArray.end()); //删除重复顶点序列
	//初始化mapOutId;
	for (int i = 0; i < out_iArray.size(); i++)
	{
		mapOutId[out_iArray.at(i)] = i;
	}
}

double QdTree_Icosa::mean_edgeDegree(Qdtree_node* node,double refR)
{
	double l = 0;
	triangle temp_tri;
	for (int j = 0; j < 3; j++)
	{
		temp_tri.ids[j] = node->tri->ids[j];
	}
	for (int i = 0; i < 3; i++)
	{
		l += distance_cpoint(vArray.at(temp_tri.ids[i%3]).posic,vArray.at(temp_tri.ids[(i+1)%3]).posic);
	}
	return l*60.0/(pi*refR);
}

//判断点是否在三角形内 使用投影方法的函数 如要使用需要取消dataStruct.h中的相应部分
/*
int QdTree_Icosa::inTriangleDot_old(Qdtree_node* node)
{
	int count;
	planePara currPlane;
	point2d currTri[3];
	point2d currP2d;
	cpoint dir_v012;
	double xmin,xmax,ymin,ymax; //用于判断点是否在外接矩形外
	if (extrapArray.empty()) //没有插入的点位 直接返回否
	{
		return 0;
	}
	else
	{
		triangle temp_tri;
		for (int j = 0; j < 3; j++)
		{
			temp_tri.ids[j] = node->tri->ids[j];
		}
		//计算一样三角形的外法线矢量
		dir_v012 = cross(vArray.at(temp_tri.ids[1]).posic - vArray.at(temp_tri.ids[0]).posic,vArray.at(temp_tri.ids[2]).posic - vArray.at(temp_tri.ids[0]).posic);
		//计算三角面平面
		currPlane = get_planePara(vArray.at(temp_tri.ids[0]).posic,vArray.at(temp_tri.ids[1]).posic,vArray.at(temp_tri.ids[2]).posic);
		//计算三角形三个顶点在新坐标系中的坐标
		currTri[0].x = 0.0; currTri[0].y = 0.0; //新坐标系的原点
		currTri[1].x = distance_cpoint(vArray.at(temp_tri.ids[0]).posic,vArray.at(temp_tri.ids[1]).posic); currTri[1].y = 0.0; //新坐标系的x轴
		currTri[2] = newXY(vArray.at(temp_tri.ids[0]),vArray.at(temp_tri.ids[1]),vArray.at(temp_tri.ids[2]),vArray.at(temp_tri.ids[2]));
		//确定外接矩形
		xmin = MAX_DBL; ymin = MAX_DBL;
		xmax = MIN_BDL; ymax = MIN_BDL;
		for (int j = 0; j < 3; j++)
		{
			if (currTri[j].x < xmin) xmin = currTri[j].x;
			if (currTri[j].x > xmax) xmax = currTri[j].x;
			if (currTri[j].y < ymin) ymin = currTri[j].y;
			if (currTri[j].y > ymax) ymax = currTri[j].y;
		}

		for (int i = 0; i < extrapArray.size(); i++)
		{
			//首先要计算三角形外法线矢量与待检测点矢量的方向是否基本一致 需要排除反向穿越的情况
			if (dot(dir_v012,extrapArray.at(i).vert.posic) > 0)
			{
				//计算带监测点在新坐标系中的位置
				currP2d = newXY(vArray.at(temp_tri.ids[0]),vArray.at(temp_tri.ids[1]),vArray.at(temp_tri.ids[2]),dotOnPlane(currPlane,extrapArray.at(i).vert));
				//满足点在外接矩形内
				if (currP2d.x > xmin && currP2d.x < xmax && currP2d.y > ymin && currP2d.y < ymax)
				{
					//向右做射线 计算与线段的交点 待检验点横坐标如果小于交点横坐标 则表示射线穿过
					//如果穿过次数为奇数 则在多边形内
					count = 0;
					for (int j = 0; j < 3; j++)
					{
						//注意现在使用的条件不包含位于多边形的边和顶点上的点
						//前两组条件表示待检测点y坐标在线段内 单侧包含
						//后一组条件表示交点横坐标大于待检测点横坐标
						//注意这里线段两个端点的y值相同时分母为零的情况在前叙的if条件中已经规避了
						if (((currP2d.y < currTri[j%3].y && currP2d.y >= currTri[(j+1)%3].y) ||
							(currP2d.y > currTri[j%3].y && currP2d.y <= currTri[(j+1)%3].y)) &&
							currP2d.x < ((currP2d.y - currTri[j%3].y)*(currTri[(j+1)%3].x - currTri[j%3].x)/(currTri[(j+1)%3].y - currTri[j%3].y) + currTri[j%3].x))
						{
							count++;
						}
					}
					//至少有一个点在当前多边形内 直接返回真
					if (pow(-1,count) < 0) return 1;
				}
			}
		}
		//没有点在当前多边形内 返回假
		return 0;
	}
}
*/
#endif