C++ Vertex_iterator类代码示例

void alignMesh( Polyhedron & poly )
{
    int num = poly.size_of_vertices();

    // initialization here is very important!!
    Point3 ave( 0.0, 0.0, 0.0 );
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	ave = ave + ( vi->point() - CGAL::ORIGIN );
    }
    ave = CGAL::ORIGIN + ( ave - CGAL::ORIGIN )/( double )num;

    unsigned int dim	= 3;
    double * data	= new double [ num*dim ];

    int nPoints = 0;
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	data[ nPoints * dim + 0 ] = vi->point().x() - ave.x();
	data[ nPoints * dim + 1 ] = vi->point().y() - ave.y();
	data[ nPoints * dim + 2 ] = vi->point().z() - ave.z();
	nPoints++;
    }

    assert( nPoints == ( int )num );

    /*****************************************
      analyze the engenstructure of X^T X
     *****************************************/
    /* define the matrix X */
    gsl_matrix_view X = gsl_matrix_view_array( data, num, dim );

    /* memory reallocation */
    // proj = ( double * )realloc( proj, sizeof( double ) * PDIM * num );

    /* calculate the covariance matrix B */
    gsl_matrix * B              = gsl_matrix_alloc( dim, dim );
    gsl_blas_dgemm( CblasTrans, CblasNoTrans, 1.0,
                    &(X.matrix),
                    &(X.matrix), 0.0,
                    B );
       
    /* divided by the number of samples */
    gsl_matrix_scale( B, 1.0/(double)num );

    gsl_vector * eVal       = gsl_vector_alloc( dim );
    gsl_matrix * eVec       = gsl_matrix_alloc( dim, dim );
    gsl_eigen_symmv_workspace * w
                                = gsl_eigen_symmv_alloc( dim );

    // eigenanalysis of the matrix B
    gsl_eigen_symmv( B, eVal, eVec, w );
    // release the memory of w
    gsl_eigen_symmv_free( w );
    // sort eigenvalues in a descending order
    gsl_eigen_symmv_sort( eVal, eVec, GSL_EIGEN_SORT_VAL_DESC );

    // #ifdef MYDEBUG
    for ( unsigned int i = 0; i < dim; ++i ) {
	cerr << "Eigenvalue No. " << i << " = " << gsl_vector_get( eVal, i ) << endl;
	cerr << "Eigenvector No. " << i << endl;
	double length = 0.0;
	for ( unsigned int j = 0; j < dim; ++j ) {
	    cerr << gsl_matrix_get( eVec, i, j ) << " ";
	    length += gsl_matrix_get( eVec, i, j )*gsl_matrix_get( eVec, i, j );
	}
	cerr << " length = " << length << endl;
    }
    // #endif	// MYDEBUG

    // 0 1 2, 0 2 1, 
    

    Transformation3 map;
    map = 
	Transformation3( gsl_matrix_get(eVec,1,0), gsl_matrix_get(eVec,0,0), gsl_matrix_get(eVec,2,0),
			 gsl_matrix_get(eVec,1,1), gsl_matrix_get(eVec,0,1), gsl_matrix_get(eVec,2,1),
			 gsl_matrix_get(eVec,1,2), gsl_matrix_get(eVec,0,2), gsl_matrix_get(eVec,2,2) );
    if ( map.is_odd() ) {
	cerr << " Transformation matrix reflected" << endl;
	map = 
	    Transformation3( gsl_matrix_get(eVec,1,0), gsl_matrix_get(eVec,0,0), -gsl_matrix_get(eVec,2,0),
			     gsl_matrix_get(eVec,1,1), gsl_matrix_get(eVec,0,1), -gsl_matrix_get(eVec,2,1),
			     gsl_matrix_get(eVec,1,2), gsl_matrix_get(eVec,0,2), -gsl_matrix_get(eVec,2,2) );
    }

    for ( unsigned int i = 0; i < dim; ++i ) {
	cerr << "| ";
	for ( unsigned int j = 0; j < dim; ++j ) {
	    cerr << map.cartesian( i, j ) << " ";
	}
	cerr << "|" << endl;
    }

    transformMesh( poly, map );

    return;
}

//------------------------------------------------------------------------------
//	Normalize the 3D triangulated mesh with the display window
//------------------------------------------------------------------------------
void normalizeMesh( Polyhedron & poly, vector< Segment3 > & bone )
{
    Vector3 sum, ave;

    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	sum = sum + ( vi->point() - CGAL::ORIGIN );
    }
    ave = sum / ( double )poly.size_of_vertices();

    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	vi->point() = vi->point() - ave;
    }
    
    cerr << " ave = " << ave << endl;
    Transformation3 translate( CGAL::TRANSLATION, -ave );

    // fabs:absolute values-no negative values
    double sideMax = 0.0;
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	if ( fabs( vi->point().x() ) > sideMax ) sideMax = fabs( vi->point().x() );
	if ( fabs( vi->point().y() ) > sideMax ) sideMax = fabs( vi->point().y() );
	if ( fabs( vi->point().z() ) > sideMax ) sideMax = fabs( vi->point().z() );
    }
    
	// sideMax: the largest number
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	vi->point() = CGAL::ORIGIN + ( vi->point() - CGAL::ORIGIN ) / sideMax;
    }

    Transformation3 scale( CGAL::SCALING, 1.0/sideMax );

    Transformation3 composite = scale * translate;
	    
    for ( unsigned int k = 0; k < bone.size(); ++k ) {
	bone[ k ] = bone[ k ].transform( composite );
    }
}

void gnuplot_print_vertices(std::ostream& out, const Vertices& vertices)
{
    for (Vertex_iterator it = vertices.begin(); it != vertices.end(); ++it)
        out << (*it) << " " << it->unique_id() << endl;
}

//------------------------------------------------------------------------------
//	Align the mesh
//------------------------------------------------------------------------------
Vector3 principalAxis( Polyhedron & poly )
{
    int num = poly.size_of_vertices();

    // initialization here is very important!!
    Point3 ave( 0.0, 0.0, 0.0 );
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	ave = ave + ( vi->point() - CGAL::ORIGIN );
    }
    ave = CGAL::ORIGIN + ( ave - CGAL::ORIGIN )/( double )num;

    unsigned int dim	= 3;
    double * data	= new double [ num*dim ];

    int nPoints = 0;
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	data[ nPoints * dim + 0 ] = vi->point().x() - ave.x();
	data[ nPoints * dim + 1 ] = vi->point().y() - ave.y();
	data[ nPoints * dim + 2 ] = vi->point().z() - ave.z();
	nPoints++;
    }

    assert( nPoints == ( int )num );

    /*****************************************
      analyze the engenstructure of X^T X
     *****************************************/
    /* define the matrix X */
    gsl_matrix_view X = gsl_matrix_view_array( data, num, dim );

    /* memory reallocation */
    // proj = ( double * )realloc( proj, sizeof( double ) * PDIM * num );

    /* calculate the covariance matrix B */
    gsl_matrix * B              = gsl_matrix_alloc( dim, dim );
    gsl_blas_dgemm( CblasTrans, CblasNoTrans, 1.0,
                    &(X.matrix),
                    &(X.matrix), 0.0,
                    B );
       
    /* divided by the number of samples */
    gsl_matrix_scale( B, 1.0/(double)num );

    gsl_vector * eVal       = gsl_vector_alloc( dim );
    gsl_matrix * eVec       = gsl_matrix_alloc( dim, dim );
    gsl_eigen_symmv_workspace * w
                                = gsl_eigen_symmv_alloc( dim );

    // eigenanalysis of the matrix B
    gsl_eigen_symmv( B, eVal, eVec, w );
    // release the memory of w
    gsl_eigen_symmv_free( w );
    // sort eigenvalues in a descending order
    gsl_eigen_symmv_sort( eVal, eVec, GSL_EIGEN_SORT_VAL_DESC );

#ifdef MYDEBUG
    for ( unsigned int i = 0; i < dim; ++i ) {
	cerr << "Eigenvalue No. " << i << " = " << gsl_vector_get( eVal, i ) << endl;
	cerr << "Eigenvector No. " << i << endl;
	for ( unsigned int j = 0; j < dim; ++j ) {
	    length += gsl_matrix_get( eVec, i, j )*gsl_matrix_get( eVec, i, j );
	}
	cerr << " length = " << length << endl;
    }
#endif	// MYDEBUG

    Vector3 ref( gsl_matrix_get( eVec, 0, 0 ),
		 gsl_matrix_get( eVec, 0, 1 ),
		 gsl_matrix_get( eVec, 0, 2 ) );
    return ref;

#ifdef DEBUG
    gsl_vector_view eachVec = gsl_matrix_column( eigenVec, 0 );
    double cosRot = gsl_matrix_get( eigenVec, 0, 1 );
    double sinRot = gsl_matrix_get( eigenVec, 1, 1 );
#ifdef DEBUG
    cerr << " 2nd axis : " << cosRot << " , " << sinRot << endl;
#endif	// DEBUG

    Transformation2 rotate( CGAL::ROTATION, -sinRot, cosRot );
    for ( unsigned int i = 0; i < subpatch.size(); ++i ) {
	subpatch[ i ]->triangle() = subpatch[ i ]->triangle().transform( rotate );
    }
#endif	// DEBUG
}

void MSDM2_Component::Matching_Multires_Update(PolyhedronPtr m_PolyDegrad, Facet * _TabMatchedFacet)
	{
		
		
		int ind=0;
		for(Vertex_iterator	pVertex	= m_PolyDegrad->vertices_begin();
					pVertex	!= m_PolyDegrad->vertices_end();
					pVertex++)
		{
                        //Point3d Nearest=pVertex->match; // MT
			Facet* f_Nearest=&_TabMatchedFacet[ind];

			pVertex->tag(ind);

			ind++;

			//for debug
			//pVertex->point()=Nearest;

			///calculation of the nearest point curvature value using vertices of the Nearest triangle	
			//we use linear interpolation using barycentric coordinates
			Point3d x1=f_Nearest->halfedge()->vertex()->point();
			Point3d x2=f_Nearest->halfedge()->next()->vertex()->point();
			Point3d x3=f_Nearest->halfedge()->next()->next()->vertex()->point();

			double l1=sqrt((x3-x2)*(x3-x2));
			double l2=sqrt((x1-x3)*(x1-x3));
			double l3=sqrt((x1-x2)*(x1-x2));

			Vector v1=f_Nearest->halfedge()->vertex()->point()-pVertex->point();
			Vector v2=f_Nearest->halfedge()->next()->vertex()->point()-pVertex->point();
			Vector v3=f_Nearest->halfedge()->next()->next()->vertex()->point()-pVertex->point();

			double t1=sqrt(v1*v1);
			double t2=sqrt(v2*v2);
			double t3=sqrt(v3*v3);
			
			double p1=(l1+t2+t3)/2;
			double p2=(t1+l2+t3)/2;
			double p3=(t1+t2+l3)/2;

			double A1=(p1*(p1-l1)*(p1-t3)*(p1-t2));
			double A2=(p2*(p2-l2)*(p2-t3)*(p2-t1));
			double A3=(p3*(p3-l3)*(p3-t1)*(p3-t2));

			if(A1>0) A1=sqrt(A1); else	A1=0;
			if(A2>0) A2=sqrt(A2); else	A2=0;
			if(A3>0) A3=sqrt(A3); else	A3=0;
			
			double c1=f_Nearest->halfedge()->vertex()->KmaxCurv;
			double c2=f_Nearest->halfedge()->next()->vertex()->KmaxCurv;
			double c3=f_Nearest->halfedge()->next()->next()->vertex()->KmaxCurv;

			if((A1+A2+A3)>0)
				pVertex->curvmatch=(A1*c1+A2*c2+A3*c3)/(A1+A2+A3);
			else
				pVertex->curvmatch=(c1+c2+c3)/3;		


			

			
		}
		
}

//! Initialize an arrangement
void init_arr(Arrangement& arr, int verbose_level)
{
  // Initialize the data of the halfedges of the arrangement.
  Halfedge_iterator heit;
  for (heit = arr.halfedges_begin(); heit != arr.halfedges_end(); ++heit)
    heit->set_data(heit->curve().data() *
                   ((heit->direction() == CGAL::ARR_LEFT_TO_RIGHT) ? 1 : 2));

  // Initialize the data of the faces of the arrangement.
  Face_iterator fit;
  for (fit = arr.faces_begin(); fit != arr.faces_end(); ++fit) {
    unsigned int count = 0;

    // Outer ccb
    Outer_ccb_iterator ocit;
    for (ocit = fit->outer_ccbs_begin(); ocit != fit->outer_ccbs_end(); ++ocit) {
      Ccb_halfedge_circulator curr = *ocit;
      do count += curr->data() * 2;
      while (++curr != *ocit);
    }

    // Inner ccbs
    Inner_ccb_iterator icit;
    for (icit = fit->inner_ccbs_begin(); icit != fit->inner_ccbs_end(); ++icit) {
      Ccb_halfedge_circulator curr = *icit;
      do count += curr->data();
      while (++curr != *icit);
    }

    fit->set_data(count);
  }

  // Initialize the data of the vertices of the arrangement.
  Vertex_iterator vit;
  for (vit = arr.vertices_begin(); vit != arr.vertices_end(); ++vit) {
    unsigned int count = 0;
    if (vit->is_isolated()) count = vit->face()->data();
    else {
      Halfedge_around_vertex_const_circulator curr = vit->incident_halfedges();
      do count += curr->data();
      while (++curr != vit->incident_halfedges());
    }
    vit->set_data(count);
  }

  if (verbose_level > 0) std::cout << "Arrangement Input: " << std::endl;

  if (verbose_level > 1) {
    std::cout << "Halfedge Data: " << std::endl;
    Halfedge_iterator heit;
    for (heit = arr.halfedges_begin(); heit != arr.halfedges_end(); ++heit)
      std::cout << heit->source()->point() << " "
                << heit->target()->point() << " " << heit->data()
                << std::endl;
  }

  if (verbose_level > 0) {
    std::cout << "Face Data: " << std::endl;
    Face_iterator fit;
    for (fit = arr.faces_begin(); fit != arr.faces_end(); ++fit)
      std::cout << fit->data() << std::endl;
  }
}

//.........这里部分代码省略.........
									//	pointMap.insert( pair<Point,Point>( c, a1 ) );

									//	c = getOffsetPoint( a2, d2, this );
									//	newPoints.push_back( c );
									//	pointMap.insert( pair<Point,Point>( c, a2 ) );
									//}
									//else {

										c = getOffsetPoint( a, d, this );
										newPoints.push_back( c );
										pointMap.insert( pair<Point,Point>( c, a ) );
									//}
								}
								else if ( vhi->info().kill != TOO_CLOSE && dist2 < minSquareDistance( averageEdgeLength, vhi, vhj ) ) {

									// The higher-address endpoint
									vhj->info().kill = TOO_CLOSE;
								}
							}
						}
					}
				}
			}
		}
	}

	double bound[ 3 ][ 2 ] = {
		{ 1.0e30, -1.0e30 },
		{ 1.0e30, -1.0e30 },
		{ 1.0e30, -1.0e30 }
	};

	// Put in all the border vertices
	for ( Vertex_iterator it = newDT->finite_vertices_begin(); it != newDT->finite_vertices_end(); ++it ) {

		if ( it->point().x() < bound[ 0 ][ 0 ] ) {

			bound[ 0 ][ 0 ] = it->point().x();
		}
		else if ( it->point().x() > bound[ 0 ][ 1 ] ) {

			bound[ 0 ][ 1 ] = it->point().x();
		}

		if ( it->point().y() < bound[ 1 ][ 0 ] ) {

			bound[ 1 ][ 0 ] = it->point().y();
		}
		else if ( it->point().y() > bound[ 1 ][ 1 ] ) {

			bound[ 1 ][ 1 ] = it->point().y();
		}

		if ( it->point().z() < bound[ 2 ][ 0 ] ) {

			bound[ 2 ][ 0 ] = it->point().z();
		}
		else if ( it->point().z() > bound[ 2 ][ 1 ] ) {

			bound[ 2 ][ 1 ] = it->point().z();
		}

		if ( !live( it->point().x(), it->point().y(), it->point().z() ) ) {

			newPoints.push_back( it->point() );
			pointMap.insert( pair<Point,Point>( it->point(), it->point() ) );

	void MSDM2_Component::ProcessMSDM2_per_vertex( Vertex_iterator pVertex,double radius,std::vector<double> & TabDistance1,std::vector<double>& TabDistance2,std::vector<Point3d> &TabPoint1,std::vector<Point3d> &TabPoint2)
	{

		std::set<int> vertices ;

        std::stack<Vertex_iterator> S ;

		Point3d O = pVertex->point() ;

        S.push(pVertex) ;	
        vertices.insert(pVertex->tag()) ;
		
	
		TabDistance1.push_back(pVertex->KmaxCurv);
		TabPoint1.push_back(pVertex->point());

		TabDistance2.push_back(pVertex->curvmatch);
		TabPoint2.push_back(pVertex->match);

		int NbSommetInSphere=0;
                //double SommeDistance=0; // MT
	

        while(!S.empty())
		{
			Vertex_iterator v = S.top() ;
            S.pop() ;
            Point3d P = v->point() ;
            Halfedge_around_vertex_circulator h = v->vertex_begin();
			Halfedge_around_vertex_circulator pHalfedgeStart = h;
			CGAL_For_all(h,pHalfedgeStart)
			{
                Point3d p1 = h->vertex()->point();
				Point3d p2 = h->opposite()->vertex()->point();

				Point3d p1m = h->vertex()->match;
				Point3d p2m = h->opposite()->vertex()->match;

				Vector V = (p2-p1);
				Vector Vm = (p2m-p1m);

                if(v==pVertex || V * (P - O) > 0.0) 
				{
					double len_old = std::sqrt(V*V);
					bool isect = sphere_clip_vector_MSDM2(O, radius, P, V) ;
					double len_edge = std::sqrt(V*V);

					NbSommetInSphere++;
					
					
					double WeightedCurv1,WeightedCurv2;
					Point3d WeightedP1,WeightedP2;

					bool IsAlreadyIntegrated=false;
					if(!isect) 
					{
						
						Vertex_iterator w=h->opposite()->vertex();
                       if(vertices.find(w->tag()) == vertices.end())
						{
                            vertices.insert(w->tag()) ;
                            S.push(w) ;
                        }
					   else
						   IsAlreadyIntegrated=true;
                    }

					if (IsAlreadyIntegrated==false)
					{
						if(len_old!=0)
						{
							if(isect)
							{
								WeightedCurv1=(1-len_edge/len_old)*h->vertex()->KmaxCurv+len_edge/len_old*h->opposite()->vertex()->KmaxCurv;
								WeightedP1=p1+V;

								WeightedCurv2=(1-len_edge/len_old)*h->vertex()->curvmatch+len_edge/len_old*h->opposite()->vertex()->curvmatch;
								WeightedP2=p1m+(len_edge/len_old)*Vm;
							}
							else
							{
								WeightedCurv1=h->opposite()->vertex()->KmaxCurv;
								WeightedCurv2=h->opposite()->vertex()->curvmatch;
								WeightedP1=p2;
								WeightedP2=p2m;
							}
						}
						else
						{
							WeightedCurv1=h->opposite()->vertex()->KmaxCurv;
							WeightedCurv2=h->opposite()->vertex()->curvmatch;
							WeightedP1=p2;
							WeightedP2=p2m;
						}

						TabDistance1.push_back(WeightedCurv1);
						TabPoint1.push_back(WeightedP1);

						TabDistance2.push_back(WeightedCurv2);
						TabPoint2.push_back(WeightedP2);
//.........这里部分代码省略.........

/**
 * Sets the new vertex info for the vertices in the new mesh
 */
void StoneWeatherer::setVertexInfo() {

	// Clear everything
	for ( Vertex_iterator it = newDT->finite_vertices_begin(); it != newDT->finite_vertices_end(); ++it ) {

		it->info().clear();
	}

	int i;

	// Label the vertices to keep (those defining boundaries between materials, where "infinite" means "air")
	for ( All_Cell_iterator it = newDT->all_cells_begin(); it != newDT->all_cells_end(); ++it ) {

		if ( newDT->is_infinite( it ) ) {

			for ( i = 0; i < 4; ++i ) {

				it->vertex( i )->info().flag |= AIR;
			}
		}
		else {

			for ( i = 0; i < 4; ++i ) {

				it->vertex( i )->info().flag |= it->info();
			}
		}
	}

	// Simplify the labels for mesh simplification later
	for ( Vertex_iterator it = newDT->finite_vertices_begin(); it != newDT->finite_vertices_end(); ++it ) {

		switch ( it->info().flag ) {

			case AIR | DIRT | ROCK | MORE_ROCK:

				// border
			case AIR | DIRT | ROCK:

				// border

			case AIR | DIRT | MORE_ROCK:

				// border

			case AIR | ROCK | MORE_ROCK:

				// border

			case DIRT | ROCK | MORE_ROCK:

				// border

			case AIR | ROCK:

				// border

			case AIR | DIRT:

				// border

			case AIR | MORE_ROCK:

				// border

			case ROCK | DIRT:

				// border

			case DIRT | MORE_ROCK:

				// border

			case ROCK | MORE_ROCK:

				// border

				it->info().kill = BORDER;
				break;

			default:

				it->info().kill = FLOATER;
		}
	}

	int n;

	// Accumulate face normals for border faces
	for ( Cell_iterator it = newDT->finite_cells_begin(); it != newDT->finite_cells_end(); ++it ) {

		if ( it->info() != AIR ) {

			for ( n = 0; n < 4; ++n ) {

				if ( it->neighbor( n )->info() < it->info() || newDT->is_infinite( it->neighbor( n ) ) ) {

//.........这里部分代码省略.........

float Remesh::RunStep() {
	cur_iter++;
	// COMPUTE THE MESH DISPLACEMENT

	for (Vertex_iterator vi = data->mesh.p.vertices_begin(); vi!=data->mesh.p.vertices_end(); vi++) 
		vi->delta = Kernel::Vector_3(0,0,0);

	
	// from current to the closest destination point
	for (Vertex_iterator vi = data->mesh.p.vertices_begin(); vi!=data->mesh.p.vertices_end(); vi++) {
		Kernel::Point_3 closest_point = dst_mesh.closestPoint(vi->point());		
		//cerr << "closest_point (" << closest_point << ") has ID " << dst_mesh.vertex_mapping[closest_point]->id << endl;	
		Kernel::Vector_3 closest_normal = dst_mesh.vertex_mapping[closest_point]->normal();
		//Kernel::Vector_3 closest_normal = dst_mesh.vertexMapping(closest_point)->normal();

/*		
		vi->delta = closest_point - vi->point();
					
		bool is_outside = v_norm((closest_point + closest_normal*v_norm(vi->delta)) - vi->point()) < v_norm(vi->delta);
//		bool is_outside = v_norm(v_normalized(closest_normal) + v_normalized(vi->delta)) < 1.4142;
//		bool is_outside = v_angle(closest_normal, vi->delta) > PI/2;		
		double dist_sign = (is_outside?1:-1);
		vi->delta = vi->normal()*v_norm(vi->delta)*(-1)*dist_sign;
*/
		
		vi->delta = vi->normal()* (closest_normal*(closest_point-vi->point())); 

//		vi->delta == v_normalized(data->mesh.computeVectorComponent(vi->normal(),vi->delta,1))*v_norm(vi->delta);
//		vi->delta = v_normalized(vi->delta)*data->mesh.computeVertexStatistics(*vi,1)*0.05;
		//		vi->delta = vi->normal(); //*alg_dt
	}

	// NORMALIZE the movements
	float total_movement = 0;
	int total_elements = 0;
	double max_delta = data->mesh.edge_avg*alg_dt;
	for (Vertex_iterator vi = data->mesh.p.vertices_begin(); vi!=data->mesh.p.vertices_end(); vi++) {	
		//vi->delta = v_normalized(vi->delta)*data->mesh.computeVertexStatistics(*vi,1)*0.1;
//		double max_delta = data->mesh.computeVertexStatistics(*vi,1)*alg_dt;
//		double min_delta = data->mesh.edge_avg/5;

		
//		vi->delta = vi->delta;
//		vi->delta = vi->delta + v_normalized(vi->delta)*(RAND_MAX/2-std::rand())*1.0/RAND_MAX*data->mesh.edge_avg/2*alg_smoothing;
//		vi->delta = v_normalized(vi->delta)*max_delta;
		


		double the_norm = v_norm(vi->delta);		
		if (the_norm > max_delta) {
			 vi->delta = v_normalized(vi->delta)*max_delta;
		}
//		if (the_norm < min_delta) vi->delta = v_normalized(vi->delta)*min_delta;	
		the_norm = v_norm(vi->delta);
		if (the_norm > max_delta*0.5) {
			total_elements++;
			total_movement	+= v_norm(vi->delta);
		}
		
		//vi->delta = vi->delta + Kernel::Vector_3((RAND_MAX/2-std::rand())*1.0/RAND_MAX, (RAND_MAX/2-std::rand())*1.0/RAND_MAX, (RAND_MAX/2-std::rand())*1.0/RAND_MAX)*data->mesh.computeVertexStatistics(*vi,1)*alg_smoothing;
//		vi->delta = vi->delta + data->mesh.computeVectorComponent(vi->normal(),vi->laplacian()*alg_dt,0);		
		if (vi->border()==false) vi->delta = vi->delta + vi->laplacian()*alg_smoothing;	
	}
	// MOVE THE MESH
	OpenGLContext::mutex.lock();
	data->mesh.lock();
	
	for (Vertex_iterator vi = data->mesh.p.vertices_begin(); vi!=data->mesh.p.vertices_end(); vi++) {
		if (alg_keepVerticesConstant) vi->delta = vi->normal()*(vi->delta*vi->normal());	
		vi->prev_delta = vi->delta;
		vi->border()=false;
		vi->move ( vi->delta );
	}
	data->mesh.unlock();
	data->mesh.updateMeshData();
	OpenGLContext::mutex.unlock();
	
	//saveOutput
	if (alg_saveOutput) {
		char filename[300];
		sprintf(filename,"%s/output_%04d.off",alg_saveOutputPrefix,cur_iter);
		data->mesh.saveFormat(filename,"off");
		sprintf(filename,"%s/output_%04d.diff",alg_saveOutputPrefix,cur_iter);
		data->mesh.saveVectorField(filename);
		sprintf(filename,"%s/output_%04d.idx",alg_saveOutputPrefix,cur_iter);
		data->mesh.saveVertexIndices(filename);
		
	}
	
	emit stepFinished();
	return total_elements;
//	return (++cur_iter < iter);	

}

// TODO: fixing after quantization the potentially introduced degeneracies, such as removing the null surface facet
void Various_Processing_Component::CoordinateQuantization (PolyhedronPtr pMesh, int bitDepth)
{
	Vertex_iterator	pVertex;
	double quantizationLevel = std::pow(2.0,bitDepth);
	pVertex = pMesh->vertices_begin();
	Point3d point = pVertex->point();
	double xmax = double(point.x());
	double xmin = xmax;
	double ymax = double(point.y());
	double ymin = ymax;
	double zmax = double(point.z());
	double zmin = zmax;
	pVertex++;

	for (; pVertex != pMesh->vertices_end(); pVertex++)
	{
		point = pVertex->point();
		double x = double(point.x());
		double y = double(point.y());
		double z = double(point.z());
		if (x>xmax)
			xmax = x;
		if (x<xmin)
			xmin = x;
		if (y>ymax)
			ymax = y;
		if (y<ymin)
			ymin = y;
		if (z>zmax)
			zmax = z;
		if (z<zmin)
			zmin = z;
	}

	double xstep = (xmax-xmin)/quantizationLevel;
	double ystep = (ymax-ymin)/quantizationLevel;
	double zstep = (zmax-zmin)/quantizationLevel;

	for (pVertex = pMesh->vertices_begin(); pVertex != pMesh->vertices_end();pVertex++)
	{
		point = pVertex->point();
		double x = double(point.x());
		double y = double(point.y());
		double z = double(point.z());

		double xquantified, yquantified, zquantified;

		double xint = 1.0*std::floor((x-xmin)/xstep)*xstep + xmin;
		double xfrac = x - xint;
		if (xfrac<=(0.5*xstep))
			xquantified = xint;
		else
			xquantified = xint + xstep;

		double yint = 1.0*std::floor((y-ymin)/ystep)*ystep + ymin;
		double yfrac = y - yint;
		if (yfrac<=(0.5*ystep))
			yquantified = yint;
		else
			yquantified = yint +ystep;

		double zint = 1.0*std::floor((z-zmin)/zstep)*zstep + zmin;
		double zfrac = z - zint;
		if (zfrac<=(0.5*zstep))
			zquantified = zint;
		else
			zquantified = zint + zstep;

		pVertex->point() = Point3d(xquantified,yquantified,zquantified);
	}

	pMesh->compute_normals();
}

// this time, we add Gaussian-distributed additive noise to the mesh vertex coordinates
// the standard deviation of the Gaussian distribution is "noiseLevel = distancetoCentroid * noiseIntensity"
void Various_Processing_Component::NoiseAdditionGaussian (PolyhedronPtr pMesh, double noiseIntensity, bool preserveBoundaries)
{
	int numVertex = pMesh->size_of_vertices();;
	Vector centroid = Point3d(0,0,0) - CGAL::ORIGIN;
	double distancetoCentroid = 0.0;
	Vertex_iterator	pVertex;
	for (pVertex = pMesh->vertices_begin(); pVertex != pMesh->vertices_end(); pVertex++)
	{
		Vector vectemp = pVertex->point() - CGAL::ORIGIN;
		centroid = centroid + vectemp;
	}
	centroid = centroid/numVertex;

	for (pVertex = pMesh->vertices_begin(); pVertex!= pMesh->vertices_end(); pVertex++)
	{
		Vector vectemp = pVertex->point() - CGAL::ORIGIN;
		distancetoCentroid = distancetoCentroid + (double)std::sqrt((vectemp - centroid) * (vectemp - centroid));
	}
	distancetoCentroid = distancetoCentroid/numVertex;

	srand((unsigned)time(NULL));
	double noisex, noisey, noisez;
	double * gaussNumbers = new double[3];
	double noiseLevel = distancetoCentroid * noiseIntensity;
	for (pVertex = pMesh->vertices_begin(); pVertex!= pMesh->vertices_end(); pVertex++)
	{
		bool is_border_vertex = false;
		bool stopFlag = false;
		Halfedge_around_vertex_circulator hav = (*pVertex).vertex_begin();
		do
		{
			if (hav->is_border()==true)
			{
				is_border_vertex = true;
				stopFlag = true;
			}
			hav++;
		} while ((hav!=(*pVertex).vertex_begin())&&(stopFlag==false));

		if ((preserveBoundaries==true)&&(is_border_vertex==true))
			continue;

		// pseudo-random Gaussian-distributed numbers generation from uniformly-distributed pseudo-random numbers
		double x, y, r2;
		for (int i=0; i<3; i++)
		{
			do
			{
				x = -1.0 + 2.0 * 1.0*rand()/RAND_MAX;
				y = -1.0 + 2.0 * 1.0*rand()/RAND_MAX;
				r2 = x * x + y * y;
			} while ((r2>1.0)||(r2==0.0));
			gaussNumbers[i] = y * sqrt(-2.0 * log(r2) / r2);
		}

		noisex = noiseLevel * gaussNumbers[0];
		noisey = noiseLevel * gaussNumbers[1];
		noisez = noiseLevel * gaussNumbers[2];
		Vector temp = Point3d(noisex, noisey, noisez) - CGAL::ORIGIN;
		pVertex->point() = pVertex->point() + temp;
	}

	pMesh->compute_normals();

	delete [] gaussNumbers;
	gaussNumbers = 0;
}