C++ MyMesh::point方法代码示例

void KnotsViewer::draw_curvature_error_domain()
{

    if(surfacedata==NULL)
        return;
    if ((surfacedata->get_domain()).empty())
        return;
    MyMesh polymesh = surfacedata->get_polymesh();
    if (polymesh.edges_empty()==true)
        return;
    MyMesh::ConstFaceIter        f_it(polymesh.faces_begin()),f_end(polymesh.faces_end());
    MyMesh::ConstFaceVertexIter  fv_it;
    //glLineWidth(0.5);
    for (; f_it!=f_end; ++f_it)
    {

        glColor3ubv(polymesh.color(*f_it).data());
        glBegin(GL_POLYGON);
        fv_it = polymesh.cfv_iter(*f_it);
        glVertex3dv(&polymesh.point(fv_it)[0]);
        ++fv_it;
        glVertex3dv(&polymesh.point(fv_it)[0]);
        ++fv_it;
        glVertex3dv(&polymesh.point(fv_it)[0]);
        ++fv_it;
        glVertex3dv(&polymesh.point(fv_it)[0]);
        glEnd();
    }


}

void BoundingBox::calculateAll(const MyMesh &mesh) {
    MyMesh::ConstVertexIter v_it=mesh.vertices_begin();
    minPoint.x = maxPoint.x = mesh.point(v_it)[0];
    minPoint.y = maxPoint.y = mesh.point(v_it)[1];
    minPoint.z = maxPoint.z = mesh.point(v_it)[2];
    ++v_it;

    for (; v_it!=mesh.vertices_end(); ++v_it){
        float x = mesh.point(v_it)[0];
        float y = mesh.point(v_it)[1];
        float z = mesh.point(v_it)[2];

        if (x < minPoint.x) minPoint.x = x;
        if (x > maxPoint.x) maxPoint.x = x;

        if (y < minPoint.y) minPoint.y = y;
        if (y > maxPoint.y) maxPoint.y = y;

        if (z < minPoint.z) minPoint.z = z;
        if (z > maxPoint.z) maxPoint.z = z;
    }

    calcDiagonal();
    calcCenterPoint();
}

void MyGL::setMesh(MyMesh mesh)
{
	
	petal_mesh = mesh;

	points.clear();


	//点
	for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); ++it){
		auto point = mesh.point(it.handle());
		coor t = { point.data()[0], point.data()[1], point.data()[2] };

		points.push_back(t);

	}

	//面
	face.clear();

	for (auto it = mesh.faces_begin(); it != mesh.faces_end(); ++it){
		OpenMesh::FaceHandle fh = it.handle();

		coor facet;//当前面索引
		int count = 0;//标记第几个点
		for (auto it1 = mesh.fv_begin(fh); it1 != mesh.fv_end(fh); ++it1){
			auto point = mesh.point(it1.handle());
			t = { point.data()[0], point.data()[1], point.data()[2] };

			//find
			vector<coor>::iterator itv;//记录该点对应Points中的位置
			itv = find_if(points.begin(), points.end(), [](coor const& obj){
				return (fabs(obj.x - t.x) < 0.01) && (fabs(obj.y - t.y) < 0.01) && (fabs(obj.z - t.z) < 0.01);
			});

			if (itv != points.end())//写入面facet
			{
				if (count == 0){
					facet.x = itv - points.begin();
				}
				else if (count == 1){
					facet.y = itv - points.begin();
				}
				else
					facet.z = itv - points.begin();
			}
			else
				break;
			count++;


		}
		face.push_back(facet);

	}

	return;
}

/*Add random Gaussian Noise to verteices*/
void AddNoise(double noise_standard_deviation,MyMesh &mesh)
{
	std::default_random_engine generator;
	std::normal_distribution<double> distribution(0.0,noise_standard_deviation); //Gaussian distribution: mean value = 0.0

	for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); ++it)
	{
		double Pt[3] = {};
		for (int d=0;d<3;d++)
		{
			Pt[d]=*(mesh.point(it).data()+d);
			double randn = distribution(generator);
			if ((randn>=-1.0)&&(randn<=1.0))//Gaussian distribution range [-1.0,1.0]							        
			{
				Pt[d]= Pt[d]*(1.0+randn);
				*(mesh.point(it).data()+d)=float(Pt[d]);
			}
		}
	}
	NOISE_CONTROL = false;
}

/*rotate mesh according to centre of mass*/
void RotateMesh(MyMesh &mesh)
{
	Point centroid(0.0,0.0,0.0);
	for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); ++it)
	{
		centroid.x() += mesh.point(it).data()[0];
		centroid.y() += mesh.point(it).data()[1];
		centroid.z() += mesh.point(it).data()[2];
	}
	centroid.x()/=float(mesh.n_vertices());
	centroid.y()/=float(mesh.n_vertices());
	centroid.z()/=float(mesh.n_vertices());

	double rotate_angle = M_PI/20.0;
	for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); ++it)
	{
		switch(ROTATE_CONTROL){
		case 1:	
			{
				double y = mesh.point(it).data()[1] - centroid.y();
				double z = mesh.point(it).data()[2] - centroid.z();

				*(mesh.point(it).data()+1) = y*cos(rotate_angle)-z*sin(rotate_angle) + centroid.y();
				*(mesh.point(it).data()+2) = y*sin(rotate_angle)+z*cos(rotate_angle) + centroid.z();
			}
			break;
		case 2:
			{
				double x = mesh.point(it).data()[0] - centroid.x();
				double z = mesh.point(it).data()[2] - centroid.z();

				*(mesh.point(it).data()+0) = x*cos(rotate_angle)-z*sin(rotate_angle) + centroid.x();
				*(mesh.point(it).data()+2) = x*sin(rotate_angle)+z*cos(rotate_angle) + centroid.z();
			} 
			break;
		case 3:
			{
				double x = mesh.point(it).data()[0] - centroid.x();
				double y = mesh.point(it).data()[1] - centroid.y();

				*(mesh.point(it).data()+0) = x*cos(rotate_angle)-y*sin(rotate_angle) + centroid.x();
				*(mesh.point(it).data()+1) = x*sin(rotate_angle)+y*cos(rotate_angle) + centroid.y();
			}
			break;
		}
	}

	ROTATE_CONTROL = 0;
}

/*Get normal vectors from OpenMesh library
-> bilateral filter denoising of vertices*/
void BilateralDenoise(MyMesh &mesh)
{
	/*ANN kd-tree find nearest point*/
	ANNpointArray	meshArray;				// mesh points array
	ANNpoint		Pt;						// point
	ANNkd_tree*		kdTree;					// search structure

	int PtNum = mesh.n_vertices();
	meshArray = annAllocPts(PtNum, dim);
	Pt = annAllocPt(dim);

	//assign mesh points to ANN array
	for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); ++it)
	{   
		//Pt get the space of data array
		double getPt[3] = {};

		//Pt get the coordinates of mesh point
		int index = it->idx();
		Pt = meshArray[index];
		for(int d = 0;d < dim; d++)
		{
			getPt[d] = *(mesh.point(it).data()+d);
			Pt[d] = getPt[d];
		}
		//assign Pt coordinates to data array
		meshArray[index] = Pt;
	}

	//build kd-tree
	kdTree = new ANNkd_tree(	// build search structure
		meshArray,				// the data points
		PtNum,					// number of points
		dim);					// dimension of space

	/*Request vertex normal*/
	// request vertex normals, so the mesh reader can use normal information
	// if available
	mesh.request_vertex_normals();

	OpenMesh::IO::Options opt;

	// If the file did not provide vertex normals, then calculate them
	if ( !opt.check( OpenMesh::IO::Options::VertexNormal ) )
	{
		// we need face normals to update the vertex normals
		mesh.request_face_normals();
		// let the mesh update the normals
		mesh.update_normals();

		// dispose the face normals, as we don't need them anymore
		mesh.release_face_normals();
	}

	/*Bilateral filtering*/
	//initial parameters for denoising
	double radius  = 0.0;
	double sigma_c = 0.05; //depend on the distance of the points
	double sigma_s = 0.0;
	int neighbour_num = 12;

	//get certain neighbour and denoise
	for (auto it = mesh.vertices_begin(); it != mesh.vertices_end(); ++it)
	{   
		//Pt get the coordinates from mesh array
		int index = it->idx();
		Pt = meshArray[index];
		std::vector<MyMesh::Point> neighbours;
		vector<double> distVector;

		//Find K(neighbours for normal calculation) nearest neighbours, 
		//return neighbours,distVector,radius
		FindNeighbours(neighbour_num,kdTree,meshArray,Pt,
			neighbours,distVector,radius);

		//get height from current vertex to tangent plane
		//update neighbour number
		neighbour_num = neighbours.size();
		vector<double> height;
		double mean_height = 0.0;
		for(int i=0;i<neighbour_num;i++)
		{
			MyMesh::Point  current_neighbour;
			current_neighbour = neighbours.at(i);

			//get normal vector from normal mesh
			double normal_vector[3]={};
			normal_vector[0]= *mesh.normal(it).data();
			normal_vector[1]= *(mesh.normal(it).data()+1);
			normal_vector[2]= *(mesh.normal(it).data()+2);

			//calculate height
			height.push_back((Pt[0]-current_neighbour.data()[0])*normal_vector[0]+(Pt[1]-current_neighbour.data()[1])*normal_vector[1]+(Pt[2]-current_neighbour.data()[2])*normal_vector[2]);
			mean_height += height.at(i);
		}

		//Calculate standard deviation(sigma_s)
		mean_height /= neighbour_num;
//.........这里部分代码省略.........

/*rasterize the model to 2R*2R*2R voxel grid*/
void RasterizeMesh(MyMesh &mesh,vector<Point> &grid_points)
{
	//reserve space 
	grid_points.reserve(10000);

	//create grid 
	int grid_size = 2*RADIUS*2*RADIUS*2*RADIUS;

	bool *grid;
	grid = new bool [grid_size];
	memset(grid,false,grid_size*sizeof(bool));

	for(MyMesh::FaceIter f_it=mesh.faces_begin();f_it!=mesh.faces_end();++f_it)
	{
		vector<Point> face_points;

		for(MyMesh::FaceVertexIter v_it=mesh.fv_iter(f_it);v_it;++v_it)
		{
			Point temp_face_point;
			//scale back from range [0,1] to [0,2R]
			temp_face_point.x() = *(mesh.point(v_it).data()+0)*2*RADIUS;
			temp_face_point.y() = *(mesh.point(v_it).data()+1)*2*RADIUS;
			temp_face_point.z() = *(mesh.point(v_it).data()+2)*2*RADIUS;
			face_points.push_back(temp_face_point);

			//get grid coordinate
			int grid_coordinate	= int(temp_face_point.x()*2*RADIUS*2*RADIUS + temp_face_point.y()*2*RADIUS + temp_face_point.z());

			//discard points out of range
			if(grid_coordinate>(2*RADIUS*2*RADIUS*2*RADIUS-1) || grid_coordinate<0)  continue; 

			//if not registered
			if(*(grid+grid_coordinate)!=true)
			{
				*(grid+grid_coordinate) = true;
				grid_points.push_back(temp_face_point);
			}
		}

		//fill the edges for a triangle face
		vector<Point> inter12,inter23,inter13;
		int n_points12 = fillGridLine(face_points.at(0),face_points.at(1),inter12,grid,grid_points);
		int n_points23 = fillGridLine(face_points.at(1),face_points.at(2),inter23,grid,grid_points);
		int n_points13 = fillGridLine(face_points.at(0),face_points.at(2),inter13,grid,grid_points);

		//fill in the face
		if(n_points12>=n_points13 && n_points13>0)
		{
			for(unsigned int n=0;n<n_points13;n++)
			{
				vector<Point> temp_inter;
				fillGridLine(inter12.at(n),inter13.at(n),temp_inter,grid,grid_points);
			}
			for(unsigned int n=n_points13;n<n_points12;n++)
			{
				vector<Point> temp_inter;
				fillGridLine(inter12.at(n),inter13.at(inter13.size()-1),temp_inter,grid,grid_points);
			}
		} 
		else if(n_points13>=n_points12 && n_points12>0)
		{
			for(unsigned int n=0;n<n_points12;n++)
			{
				vector<Point> temp_inter;
				fillGridLine(inter12.at(n),inter13.at(n),temp_inter,grid,grid_points);
			}
			for(unsigned int n=n_points12;n<n_points13;n++)
			{
				vector<Point> temp_inter;
				fillGridLine(inter12.at(inter12.size()-1),inter13.at(n),temp_inter,grid,grid_points);
			}
		}
	}

	//get centroid again
	Point centroid_after(0.0,0.0,0.0);
	for(unsigned int p_it = 0;p_it<grid_points.size();p_it++)
	{
		centroid_after.x()+=grid_points.at(p_it).x();
		centroid_after.y()+=grid_points.at(p_it).y();
		centroid_after.z()+=grid_points.at(p_it).z();

	}
	centroid_after.x()/=double(grid_points.size());
	centroid_after.y()/=double(grid_points.size());
	centroid_after.z()/=double(grid_points.size());

	//move grid to the origin
	for(unsigned int p_it = 0;p_it<grid_points.size();p_it++)
	{
		grid_points.at(p_it).x() -= centroid_after.x();
		grid_points.at(p_it).y() -= centroid_after.y();
		grid_points.at(p_it).z() -= centroid_after.z();
	}

	delete [] grid;

	RASTERIZE_CONTROL = FALSE;
}

/*normalize the model */
void NormalizeMesh(MyMesh &mesh)
{
	double x_max,y_max,z_max,x_min,y_min,z_min;
	FindMaxMin(mesh,x_max,y_max,z_max,x_min,y_min,z_min);

	double distance_x = x_max - x_min;
	double distance_y = y_max - y_min;
	double distance_z = z_max - z_min;
	double max_distance = distance_x;

	if (distance_y > max_distance) max_distance = distance_y;
	if (distance_z > max_distance) max_distance = distance_z;

	/*normalize*/
	//initial centroid
	Point centroid(0.0,0.0,0.0);
	for (MyMesh::VertexIter v_it = mesh.vertices_begin();v_it!=mesh.vertices_end(); ++v_it)
	{
		//move to positive, normalize to 1 
		double x_normalize = (mesh.point(v_it).data()[0] - x_min)/max_distance;
		double y_normalize = (mesh.point(v_it).data()[1] - y_min)/max_distance;
		double z_normalize = (mesh.point(v_it).data()[2] - z_min)/max_distance;

		*(mesh.point(v_it).data()+0) = x_normalize;
		*(mesh.point(v_it).data()+1) = y_normalize;
		*(mesh.point(v_it).data()+2) = z_normalize;

		//get centroid
		centroid.x() += double(x_normalize);
		centroid.y() += double(y_normalize);
		centroid.z() += double(z_normalize);
	}
	centroid.x()/=double(mesh.n_vertices());
	centroid.y()/=double(mesh.n_vertices());
	centroid.z()/=double(mesh.n_vertices());

	//get distance
	float mean_dist = 0.0;
	for (MyMesh::VertexIter v_it = mesh.vertices_begin();v_it!=mesh.vertices_end(); ++v_it)
	{
		//get distance between vertex and origin
		double current_dist = sqrt(pow(*(mesh.point(v_it).data()+0)-centroid.x(),2.0) 
			+ pow(*(mesh.point(v_it).data()+1)-centroid.y(),2.0) + pow(*(mesh.point(v_it).data()+2)-centroid.z(),2.0));
		mean_dist += current_dist;
	}
	mean_dist/=double(mesh.n_vertices());

	//scale and make the average distance to center of mass is R/2
	float scale_ratio = (double(RADIUS)/2.0)/mean_dist;
	for (MyMesh::VertexIter v_it = mesh.vertices_begin();v_it!=mesh.vertices_end(); ++v_it)
	{
		*(mesh.point(v_it).data()+0) *= scale_ratio;
		*(mesh.point(v_it).data()+1) *= scale_ratio;
		*(mesh.point(v_it).data()+2) *= scale_ratio;

		//scale back to range [0,1] for better viewing 
		*(mesh.point(v_it).data()+0) /= double(2*RADIUS);
		*(mesh.point(v_it).data()+1) /= double(2*RADIUS);
		*(mesh.point(v_it).data()+2) /= double(2*RADIUS);
	}

	NORMALIZE_CONTROL = FALSE;
}

void SolveLinearSystem(MyMesh* mesh,std::vector<int> controlPoint_indexes){
	if(mesh==NULL)
	{
		return ;
	}
	//Ax=B
	LinearSystemLib::GeneralSparseMatrix GA;//Ax=B 的 A
	double** B_X = new double*[1];//Ax=B 的 B 針對x座標解
	double** B_Y = new double*[1];//Ax=B 的 B 針對y座標解
	double** B_Z = new double*[1];//Ax=B 的 B 針對z座標解
	double** solution_X = 0;//Ax=B x的解 針對x座標解
	double** solution_Y = 0;//Ax=B x的解 針對y座標解
	double** solution_Z = 0;//Ax=B x的解 針對z座標解
	//
	int n_vertices = mesh->n_vertices();//Mesh's total vertex
	int n_control_point = controlPoint_indexes.size();//total control point
	//Set matrix size
	GA.Create(n_vertices+n_control_point, n_vertices);
	B_X[0] = new double[n_vertices+n_control_point];
	B_Y[0] = new double[n_vertices+n_control_point];
	B_Z[0] = new double[n_vertices+n_control_point];
	//填上半部矩陣 點對點
	for(MyMesh::VIter v_it = mesh->vertices_begin(); v_it!=mesh->vertices_end(); ++v_it)
	{
		MyMesh::VHandle vh = v_it.handle();
		int row,col;//要放入的r,c
		int vv_total = 0;//點對點的總量
		double element;//要放入矩陣的元素

		row = v_it.handle().idx();//Set row index
		GA.SetElement(row,row,1);//if row == column , set 1 to (row,column)

		//先獲得點對點總量
		for(MyMesh::VVIter vv_it = mesh->vv_iter(v_it.handle());vv_it;++vv_it)
		{
			vv_total++;
		}
		//要放入矩陣的元素 為 -1/vv_total
		element = double(-1) / double(vv_total);
		//放入矩陣的元素
		for(MyMesh::VVIter vv_it = mesh->vv_iter(v_it.handle());vv_it;++vv_it)
		{
			col = vv_it.handle().idx();//Set column index
			GA.SetElement(row,col,element);//Set element
			//printf("row : %d , col : %d , element : %f\n",row,col,element);
		}

		//Set B's element 上半部為0
		B_X[0][row] = 0.0;
		B_Y[0][row] = 0.0;
		B_Z[0][row] = 0.0;
	}
	//Control point
	for(int i=0;i<controlPoint_indexes.size();i++)
	{
		//填下半部矩陣 control point
		GA.SetElement(n_vertices+i,controlPoint_indexes[i],1);
		//Set control point to B's 下半部
		MyMesh::VHandle vh = mesh->vertex_handle(controlPoint_indexes[i]);
		B_X[0][ n_vertices+i ] = mesh->point(vh)[0];
		B_Y[0][ n_vertices+i ] = mesh->point(vh)[1];
		B_Z[0][ n_vertices+i ] = mesh->point(vh)[2];
	}
	//Solve linear system
	LinearSystemLib::SparseLinearSystem sls_X(new LinearSystemLib::StableSparseMatrix(GA), B_X, 1);
	LinearSystemLib::SparseLinearSystem sls_Y(new LinearSystemLib::StableSparseMatrix(GA), B_Y, 1);
	LinearSystemLib::SparseLinearSystem sls_Z(new LinearSystemLib::StableSparseMatrix(GA), B_Z, 1);
	try
	{
		bool result;
		result = LinearSystemLib::LeastSquareSparseLSSolver::GetInstance()->Solve( &sls_X, solution_X );
		result = LinearSystemLib::LeastSquareSparseLSSolver::GetInstance()->Solve( &sls_Y, solution_Y );
		result = LinearSystemLib::LeastSquareSparseLSSolver::GetInstance()->Solve( &sls_Z, solution_Z );
	} 
	catch( std::exception e )
	{
		printf("Error : %s",e.what());
		//cerr << e.what() << endl;
		return ;
	}
	//Set point to mesh
	for(int i=0;i<n_vertices;i++)
	{
		MyMesh::VHandle vh = mesh->vertex_handle(i);
		mesh->point(vh)[0] =  solution_X[0][i];
		mesh->point(vh)[1] =  solution_Y[0][i];
		mesh->point(vh)[2] =  solution_Z[0][i];
	}

	printf("Solve success");
}