int main(int argc,char** argv)
{

	setCameraDistance(30.f);

#define TRISIZE 10.f
#ifdef DEBUG_MESH
	SimdVector3 vert0(-TRISIZE ,0,TRISIZE );
	SimdVector3 vert1(TRISIZE ,10,TRISIZE );
	SimdVector3 vert2(TRISIZE ,0,-TRISIZE );
	meshData.AddTriangle(vert0,vert1,vert2);
	SimdVector3 vert3(-TRISIZE ,0,TRISIZE );
	SimdVector3 vert4(TRISIZE ,0,-TRISIZE );
	SimdVector3 vert5(-TRISIZE ,0,-TRISIZE );
	meshData.AddTriangle(vert3,vert4,vert5);
#else
#ifdef ODE_MESH
	SimdVector3 Size = SimdVector3(15.f,15.f,12.5f);
	
  gVertices[0][0] = -Size[0];
  gVertices[0][1] = Size[2];
  gVertices[0][2] = -Size[1];
  
  gVertices[1][0] = Size[0];
  gVertices[1][1] = Size[2];
  gVertices[1][2] = -Size[1];
  
  gVertices[2][0] = Size[0];
  gVertices[2][1] = Size[2];
  gVertices[2][2] = Size[1];  

  gVertices[3][0] = -Size[0];
  gVertices[3][1] = Size[2];
  gVertices[3][2] = Size[1];
  
  gVertices[4][0] = 0;
  gVertices[4][1] = 0;
  gVertices[4][2] = 0;
  
  gIndices[0] = 0;
  gIndices[1] = 1;
  gIndices[2] = 4;
  
  gIndices[3] = 1;
  gIndices[4] = 2;
  gIndices[5] = 4;
  
  gIndices[6] = 2;
  gIndices[7] = 3;
  gIndices[8] = 4;
  
  gIndices[9] = 3;
  gIndices[10] = 0;
  gIndices[11] = 4;

  int vertStride = sizeof(SimdVector3);
  int indexStride = 3*sizeof(int);

	TriangleIndexVertexArray* indexVertexArrays = new TriangleIndexVertexArray(NUM_TRIANGLES,
		gIndices,
		indexStride,
		NUM_VERTICES,(float*) &gVertices[0].x(),vertStride);

	//shapePtr[4] = new TriangleMeshShape(indexVertexArrays);
	shapePtr[4] = new BvhTriangleMeshShape(indexVertexArrays);
#else

	int vertStride = sizeof(SimdVector3);
	int indexStride = 3*sizeof(int);

	const int NUM_VERTS_X = 50;
	const int NUM_VERTS_Y = 50;
	const int totalVerts = NUM_VERTS_X*NUM_VERTS_Y;
	
	const int totalTriangles = 2*(NUM_VERTS_X-1)*(NUM_VERTS_Y-1);

	SimdVector3*	gVertices = new SimdVector3[totalVerts];
	int*	gIndices = new int[totalTriangles*3];

	int i;

	for ( i=0;i<NUM_VERTS_X;i++)
	{
		for (int j=0;j<NUM_VERTS_Y;j++)
		{
			gVertices[i+j*NUM_VERTS_X].setValue((i-NUM_VERTS_X*0.5f)*10.f,2.f*sinf((float)i)*cosf((float)j),(j-NUM_VERTS_Y*0.5f)*10.f);
		}
	}

	int index=0;
	for ( i=0;i<NUM_VERTS_X-1;i++)
	{
		for (int j=0;j<NUM_VERTS_Y-1;j++)
		{
			gIndices[index++] = j*NUM_VERTS_X+i;
			gIndices[index++] = j*NUM_VERTS_X+i+1;
			gIndices[index++] = (j+1)*NUM_VERTS_X+i+1;

			gIndices[index++] = j*NUM_VERTS_X+i;
			gIndices[index++] = (j+1)*NUM_VERTS_X+i+1;
//.........这里部分代码省略.........

void PoissonSurfaceReconstruction::reconstruct(std::vector<Eigen::Vector3d> &points, std::vector<Eigen::Vector3d> &normals, TriangleMesh &mesh){

	assert(points.size() == normals.size());

  std::cout << "creating points with normal..." << std::endl;
  std::vector<Point_with_normal> points_with_normal;
  points_with_normal.resize((int)points.size());
  for(int i=0; i<(int)points.size(); i++){
    Vector vec(normals[i][0], normals[i][1], normals[i][2]);
    //Point_with_normal pwn(points[i][0], points[i][1], points[i][2], vec);
    //points_with_normal[i]  = pwn;
    points_with_normal[i] = Point_with_normal(points[i][0], points[i][1], points[i][2], vec);
  }

  std::cout << "constructing poisson reconstruction function..." << std::endl;
  Poisson_reconstruction_function function(points_with_normal.begin(), points_with_normal.end(),
                                           CGAL::make_normal_of_point_with_normal_pmap(PointList::value_type()));

  std::cout << "computing implicit function..." << std::endl;
  if( ! function.compute_implicit_function() ) {
  	std::cout << "compute implicit function is failure" << std::endl;
  	return;
  }
  	//return EXIT_FAILURE;

  // Computes average spacing
  std::cout << "compute average spacing..." << std::endl;
  FT average_spacing = CGAL::compute_average_spacing(points_with_normal.begin(), points_with_normal.end(),
                                                       6 /* knn = 1 ring */);
  // Gets one point inside the implicit surface
  // and computes implicit function bounding sphere radius.
  Point inner_point = function.get_inner_point();
  Sphere bsphere = function.bounding_sphere();
  FT radius = std::sqrt(bsphere.squared_radius());

  // Defines the implicit surface: requires defining a
  // conservative bounding sphere centered at inner point.
  FT sm_sphere_radius = 5.0 * radius;
  FT sm_dichotomy_error = distance_criteria*average_spacing/1000.0; // Dichotomy error must be << sm_distance
  //FT sm_dichotomy_error = distance_criteria*average_spacing/10.0; // Dichotomy error must be << sm_distance
  std::cout << "reconstructed surface" << std::endl;
  Surface_3 reconstructed_surface(function,
                    Sphere(inner_point,sm_sphere_radius*sm_sphere_radius),
                    sm_dichotomy_error/sm_sphere_radius);

  // Defines surface mesh generation criteria    
  CGAL::Surface_mesh_default_criteria_3<STr> criteria(angle_criteria,  // Min triangle angle (degrees)
                                                        radius_criteria*average_spacing,  // Max triangle size
                                                        distance_criteria*average_spacing); // Approximation error
  
  std::cout << "generating surface mesh..." << std::endl;
  // Generates surface mesh with manifold option
  STr tr; // 3D Delaunay triangulation for surface mesh generation
  C2t3 c2t3(tr); // 2D complex in 3D Delaunay triangulation
  CGAL::make_surface_mesh(c2t3,                                 // reconstructed mesh
                            reconstructed_surface,                              // implicit surface
                            criteria,                             // meshing criteria
                            CGAL::Manifold_tag());  // require manifold mesh

  if(tr.number_of_vertices() == 0){
  	std::cout << "surface mesh generation is failed" << std::endl;
  	return;
  }

  Polyhedron surface;
  CGAL::output_surface_facets_to_polyhedron(c2t3, surface);

  // convert CGAL::surface to TriangleMesh //
  std::cout << "converting CGA::surface to TriangleMesh..." << std::endl;
	std::vector<Eigen::Vector3d> pts;
	std::vector<std::vector<int> > faces;
	pts.resize(surface.size_of_vertices());
	faces.resize(surface.size_of_facets());

	Polyhedron::Point_iterator pit;
	int index = 0;
	for(pit=surface.points_begin(); pit!=surface.points_end(); ++pit){
		pts[index][0] = pit->x();
		pts[index][1] = pit->y();
		pts[index][2] = pit->z();				
		index ++;
	}
	index = 0;
	Polyhedron::Face_iterator fit;
	for(fit=surface.facets_begin(); fit!=surface.facets_end(); ++fit){
		std::vector<int > face(3);
		Halfedge_facet_circulator j = fit->facet_begin();
		int f_index = 0;
		do {
			face[f_index] = std::distance(surface.vertices_begin(), j->vertex());
			f_index++;
    } while ( ++j != fit->facet_begin());    

    faces[index] = face;
		index++;
	}

	mesh.createFromFaceVertex(pts, faces);  
}

void simplifyMesh(TriangleMesh &src, TriangleMesh &dst){
	namespace SMS = CGAL::Surface_mesh_simplification; 

	Surface surface;
  PolyhedronBuilder<HalfedgeDS> builder(src.points, src.facets);
	surface.delegate(builder);


	// --- before mesh simplification, to maintain mesh quality high poisson surface reconstruction is applied ---


	// This is a stop predicate (defines when the algorithm terminates).
	// In this example, the simplification stops when the number of undirected edges
	// left in the surface drops below the specified number (1000) 
	// SMS::Count_stop_predicate<Surface> stop(3000);
   //SMS::Count_ratio_stop_predicate<Surface> stop(0.1);
	SMS::Count_stop_predicate<Surface> stop(100);	
  // This the actual call to the simplification algorithm.
  // The surface and stop conditions are mandatory arguments. 
  // The index maps are needed because the vertices and edges 
  // of this surface lack an "id()" field.
  int r = SMS::edge_collapse
             (surface
             ,stop
             ,CGAL::vertex_index_map(boost::get(CGAL::vertex_external_index, surface))
             .edge_index_map(boost::get(CGAL::edge_external_index, surface))
             .get_cost(CGAL::Surface_mesh_simplification::Edge_length_cost<Polyhedron>())
             .get_placement(CGAL::Surface_mesh_simplification::Midpoint_placement<Polyhedron>())
             //.get_cost(CGAL::Surface_mesh_simplification::LindstromTurk_cost<Polyhedron>())            
             //.get_placement(CGAL::Surface_mesh_simplification::LindstromTurk_placement<Polyhedron>())
             );
  std::cout << "\nFinished...\n" << r << " edges removed.\n"
           << (surface.size_of_halfedges()/2) << " final edges.\n" ;

	std::vector<Eigen::Vector3d> points;
	std::vector<std::vector<int> > faces;
	points.resize(surface.size_of_vertices());
	faces.resize(surface.size_of_facets());

	Polyhedron::Point_iterator pit;
	int index = 0;
	for(pit=surface.points_begin(); pit!=surface.points_end(); ++pit){
		points[index][0] = pit->x();
		points[index][1] = pit->y();
		points[index][2] = pit->z();				
		index ++;
	}
	index = 0;
	Polyhedron::Face_iterator fit;
	for(fit=surface.facets_begin(); fit!=surface.facets_end(); ++fit){
		std::vector<int > face(3);
		Halfedge_facet_circulator j = fit->facet_begin();
		int f_index = 0;
		do {
			face[f_index] = std::distance(surface.vertices_begin(), j->vertex());
			f_index++;
    } while ( ++j != fit->facet_begin());    

    faces[index] = face;
		index++;
	}

	dst.createFromFaceVertex(points, faces);
}

//.........这里部分代码省略.........
    {
        if (m_radius < 0) m_radius = 0.5f*extend.length_2d();
        m_shape = new btCylinderShape(0.5f*extend);
        break;
    }
    case MP_CYLINDER_X:
    {
        if (m_radius < 0)
            m_radius = 0.5f*sqrt(extend.getY()*extend.getY() +
                                 extend.getZ()*extend.getZ());
        m_shape = new btCylinderShapeX(0.5f*extend);
        break;
    }
    case MP_CYLINDER_Z:
    {
        if (m_radius < 0)
            m_radius = 0.5f*sqrt(extend.getX()*extend.getX() +
                                 extend.getY()*extend.getY());
        m_shape = new btCylinderShapeZ(0.5f*extend);
        break;
    }
    case MP_SPHERE:
    {
        if(m_radius<0)
        {
            m_radius =      std::max(extend.getX(), extend.getY());
            m_radius = 0.5f*std::max(m_radius,      extend.getZ());
        }
        m_shape = new btSphereShape(m_radius);
        break;
    }
    case MP_EXACT:
    {
        TriangleMesh* triangle_mesh = new TriangleMesh();

        // In case of readonly materials we have to get the material from
        // the mesh, otherwise from the node. This is esp. important for
        // water nodes, which only have the material defined in the node,
        // but not in the mesh at all!
        bool is_readonly_material = false;

        scene::IMesh* mesh = NULL;
        switch (presentation->getNode()->getType())
        {
            case scene::ESNT_MESH          :
            case scene::ESNT_WATER_SURFACE :
            case scene::ESNT_OCTREE        :
                {
                    scene::IMeshSceneNode *node =
                        (scene::IMeshSceneNode*)presentation->getNode();
                    mesh = node->getMesh();
                    is_readonly_material = node->isReadOnlyMaterials();
                    break;
                }
            case scene::ESNT_ANIMATED_MESH :
                {
                    // for now just use frame 0
                    scene::IAnimatedMeshSceneNode *node =
                      (scene::IAnimatedMeshSceneNode*)presentation->getNode();
                    mesh = node->getMesh()->getMesh(0);
                    is_readonly_material = node->isReadOnlyMaterials();
                    break;
                }
            default:
                Log::warn("PhysicalObject", "Unknown object type, "
                                        "cannot create exact collision body!");

//////////////////////////////////////////////////////////
//
// MeshSweeper implementation
// ===========
//
TriangleMesh*
MeshSweeper::makeCylinder(
  const Polyline& circle,
  const vec3& path,
  const mat4& m)
//[]----------------------------------------------------[]
//|  Make cylinder                                       |
//[]----------------------------------------------------[]
{
  int np = circle.getNumberOfVertices();
  int nv = np * 2; // number of vertices
  int nb = np - 2; // number of triangles of the base
  int nt = nv + 2 * nb; // number of triangles
  TriangleMesh::Arrays data;

  data.vertices = new vec3[data.numberOfVertices = nv];
  data.triangles = new TriangleMesh::Triangle[data.numberOfTriangles = nt];

  vec3 c(0, 0, 0);

  if (true)
  {
    Polyline::VertexIterator vit(circle.getVertexIterator());

    for (int i = 0; i < np; i++)
    {
      const vec3& p = vit++.position;

      c += p;
      data.vertices[i + np] = m.transform3x4(p);
      data.vertices[i] = m.transform3x4(p + path);
    }
    c *= Math::inverse(REAL(np));
  }

  TriangleMesh::Triangle* triangle = data.triangles;

  for (int i = 0; i < np; i++)
  {
    int j = (i + np);
    int k = (i + 1) % np;

    triangle->setVertices(i, j, k);
    triangle[1].setVertices(j, k + np, k);
    triangle += 2;
  }

  int v0 = 0;
  int v1 = 1;
  int v2 = 2;

  for (int i = 0; i < nb; i++)
  {
    triangle->setVertices(v0, v1, v2);
    triangle[nb].setVertices(v0 + np, v2 + np, v1 + np);
    triangle++;
    v2 = ((v1 = (v0 = v2) + 1) + 1) % np;
  }

  TriangleMesh* mesh = new TriangleMesh(data);

  mesh->computeNormals();
  return mesh;
}

void
makeBunny20Scene()
{
	g_camera = new Camera;
	g_scene = new Scene;
	g_image = new Image;

	g_image->resize(512, 512);

	// set up the camera
	g_camera->setBGColor(Vector3(0.0f, 0.0f, 0.2f));
	g_camera->setEye(Vector3(0, 5, 15));
	g_camera->setLookAt(Vector3(0, 0, 0));
	g_camera->setUp(Vector3(0, 1, 0));
	g_camera->setFOV(45);

	// create and place a point light source
	PointLight * light = new PointLight;
	light->setPosition(Vector3(10, 20, 10));
	light->setColor(Vector3(1, 1, 1));
	light->setWattage(4.0 * PI * 1000);
	g_scene->addLight(light);

	TriangleMesh * mesh;
	Material* material = new Lambert(Vector3(1.0f));
	Matrix4x4 xform;
	Matrix4x4 xform2;
	xform2 *= rotate(110, 0, 1, 0);
	xform2 *= scale(.6, 1, 1.1);


	// bunny 1
	xform.setIdentity();
	xform *= scale(0.3, 2.0, 0.7);
	xform *= translate(-1, .4, .3);
	xform *= rotate(25, .3, .1, .6);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 2
	xform.setIdentity();
	xform *= scale(.6, 1.2, .9);
	xform *= translate(7.6, .8, .6);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 3
	xform.setIdentity();
	xform *= translate(.7, 0, -2);
	xform *= rotate(120, 0, .6, 1);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 4
	xform.setIdentity();
	xform *= translate(3.6, 3, -1);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 5
	xform.setIdentity();
	xform *= translate(-2.4, 2, 3);
	xform *= scale(1, .8, 2);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 6
	xform.setIdentity();
	xform *= translate(5.5, -.5, 1);
	xform *= scale(1, 2, 1);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 7
	xform.setIdentity();
	xform *= rotate(15, 0, 0, 1);
	xform *= translate(-4, -.5, -6);
	xform *= scale(1, 2, 1);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 8
	xform.setIdentity();
	xform *= rotate(60, 0, 1, 0);
	xform *= translate(5, .1, 3);
	mesh = new TriangleMesh;
	mesh->load("bunny.obj", xform);
	addMeshTrianglesToScene(mesh, material);

	// bunny 9
	xform.setIdentity();
	xform *= translate(-3, .4, 6);
	xform *= rotate(-30, 0, 1, 0);
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
//  BoundaryMap private method definitions
//-----------------------------------------------------------------------------
void BoundaryMap::computeSignedDistanceField(SparseVoxelMap<float*>& map, 
        const TriangleMesh& mesh)
{
    const float* vertexList = mesh.getVertexList();
    const unsigned int* faceList = mesh.getFaceList();
    const float *v1, *v2, *v3;
    unsigned int cMin[3];
    unsigned int cMax[3];
    float x0 = mDomain.getV1().getX();
    float y0 = mDomain.getV1().getY();
    float z0 = mDomain.getV1().getZ();
    float x[3];
    float normal[3];
    float dist;

    //
    // compute distance
    //

    // for each triangle of the mesh
    for (unsigned int m = 0; m < mesh.getNumFaces(); m++) 
    {
        std::cout << m << " of " << mesh.getNumFaces() << std::endl;

        //get triangle vertices
        v1 = &vertexList[3*faceList[3*m + 0]];
        v2 = &vertexList[3*faceList[3*m + 1]];
        v3 = &vertexList[3*faceList[3*m + 2]];

        // compute bounding box of the triangle
        this->computeGridCoordinates(cMin, cMax, v1, v2, v3);

        //std::cout << "cMin[0] = " << cMin[0] << std::endl;
        //std::cout << "cMin[1] = " << cMin[1] << std::endl;
        //std::cout << "cMin[2] = " << cMin[2] << std::endl;
        //std::cout << "cMax[0] = " << cMax[0] << std::endl;
        //std::cout << "cMax[1] = " << cMax[1] << std::endl;
        //std::cout << "cMax[2] = " << cMax[2] << std::endl;

        // for each coordinate within the bounding box
        for (unsigned int k = cMin[2]; k <= cMax[2]; k++) 
        {
            for (unsigned int j = cMin[1]; j <= cMax[1]; j++) 
            {
                for (unsigned int i = cMin[0]; i <= cMax[0]; i++) 
                {
                    // translate coordinate to world coordinates
                    x[0] = x0 + i*mDx;
                    x[1] = y0 + j*mDx;
                    x[2] = z0 + k*mDx;

                    // compute distance to the triangle
                    compute_distance_point_triangle(dist, normal, x, v1, v2, v3);

                    // update sparse voxel map
                    if (dist <= mMaxDist)
                    {
                        Coordinate coord(i, j, k);

                        // if the map already contains distance information
                        // at the current coordinate
                        if (map.contains(coord)) 
                        {
                            // check if the previously computed distance is
                            // is greater than [dist].
                            float prev;
                            float* nodeContent;
                            map.get(nodeContent, coord);
                            prev = nodeContent[NC_DISTANCE];
                            if (prev > dist)
                            {
                                // if yes, update the map
                                delete[] nodeContent;
                                nodeContent = new float[NC_NUM_ELEMENTS];
                                nodeContent[NC_DISTANCE] = dist;
                                nodeContent[NC_NORMAL_X] = normal[0];
                                nodeContent[NC_NORMAL_Y] = normal[1];
                                nodeContent[NC_NORMAL_Z] = normal[2];
                                map.add(coord, nodeContent);
                            }
                        }
                        // if not yet contained in the map
                        else
                        {
                            float* nodeContent = new float[NC_NUM_ELEMENTS];
                            nodeContent[NC_DISTANCE] = dist;
                            nodeContent[NC_NORMAL_X] = normal[0];
                            nodeContent[NC_NORMAL_Y] = normal[1];
                            nodeContent[NC_NORMAL_Z] = normal[2];
                            //std::cout << dist << std::endl;
                            // just add val to the map
                            map.add(coord, nodeContent);
                        }
                    }
                }
            }
        }
//.........这里部分代码省略.........

int
main(const int argc, const char **argv)
{
    // parse all command line options into a DynamicConfig
    ConfigDef config_def;
    config_def.merge(cli_config_def);
    config_def.merge(print_config_def);
    DynamicConfig config(&config_def);
    t_config_option_keys input_files;
    config.read_cli(argc, argv, &input_files);
    
    // apply command line options to a more handy CLIConfig
    CLIConfig cli_config;
    cli_config.apply(config, true);
    
    DynamicPrintConfig print_config;
    
    // load config files supplied via --load
    for (std::vector<std::string>::const_iterator file = cli_config.load.values.begin();
        file != cli_config.load.values.end(); ++file) {
        if (!boost::filesystem::exists(*file)) {
            std::cout << "No such file: " << *file << std::endl;
            exit(1);
        }
        
        DynamicPrintConfig c;
        try {
            c.load(*file);
        } catch (std::exception &e) {
            std::cout << "Error while reading config file: " << e.what() << std::endl;
            exit(1);
        }
        c.normalize();
        print_config.apply(c);
    }
    
    // apply command line options to a more specific DynamicPrintConfig which provides normalize()
    // (command line options override --load files)
    print_config.apply(config, true);
    print_config.normalize();
    
    // write config if requested
    if (!cli_config.save.value.empty()) print_config.save(cli_config.save.value);
    
    // read input file(s) if any
    std::vector<Model> models;
    for (t_config_option_keys::const_iterator it = input_files.begin(); it != input_files.end(); ++it) {
        if (!boost::filesystem::exists(*it)) {
            std::cout << "No such file: " << *it << std::endl;
            exit(1);
        }
        
        Model model;
        // TODO: read other file formats with Model::read_from_file()
        try {
            Slic3r::IO::STL::read(*it, &model);
        } catch (std::exception &e) {
            std::cout << *it << ": " << e.what() << std::endl;
            exit(1);
        }
        
        if (model.objects.empty()) {
            printf("Error: file is empty: %s\n", it->c_str());
            continue;
        }
        
        model.add_default_instances();
        
        // apply command line transform options
        for (ModelObjectPtrs::iterator o = model.objects.begin(); o != model.objects.end(); ++o) {
            if (cli_config.scale_to_fit.is_positive_volume())
                (*o)->scale_to_fit(cli_config.scale_to_fit.value);
            
            (*o)->scale(cli_config.scale.value);
            (*o)->rotate(cli_config.rotate.value, Z);
        }
        
        // TODO: handle --merge
        models.push_back(model);
    }
    
    for (std::vector<Model>::iterator model = models.begin(); model != models.end(); ++model) {
        if (cli_config.info) {
            // --info works on unrepaired model
            model->print_info();
        } else if (cli_config.export_obj) {
            std::string outfile = cli_config.output.value;
            if (outfile.empty()) outfile = model->objects.front()->input_file + ".obj";
    
            TriangleMesh mesh = model->mesh();
            mesh.repair();
            Slic3r::IO::OBJ::write(mesh, outfile);
            printf("File exported to %s\n", outfile.c_str());
        } else if (cli_config.export_pov) {
            std::string outfile = cli_config.output.value;
            if (outfile.empty()) outfile = model->objects.front()->input_file + ".pov";
    
            TriangleMesh mesh = model->mesh();
            mesh.repair();
            Slic3r::IO::POV::write(mesh, outfile);
//.........这里部分代码省略.........

int main(int argc, const char * argv[])
{

  Electrostatics e;

  TriangleMesh *teststl = new TriangleMesh();
  teststl->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere5120.stl", MeshFileFormat::MFF_STL);
  teststl->write("/Users/phaedon/github/bem-laplace-simple/meshes/sphere5120.obj", MeshFileFormat::MFF_OBJ);
  
  TriangleMesh *s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h0cm_s3.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);
  
  s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h10cm.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);

  s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h30cm.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);
  
  s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h50cm.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);
  
  s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h80cm.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);
  
  s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h90cm.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);
  
  s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h95cm.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);
  
  
  s = new TriangleMesh();
  s->read("/Users/phaedon/github/bem-laplace-simple/meshes/sphere_5mm_h995mm_s3.obj", MeshFileFormat::MFF_OBJ);
  e.addBubble(s);
  
  TriangleMesh *plane = new TriangleMesh();
  plane->read("/Users/phaedon/github/bem-laplace-simple/meshes/plane_h1_s50mm_t128.obj", MeshFileFormat::MFF_OBJ);

  e.setSurface(plane);
  std::vector<double> caps;
  e.capacitance(caps);
  for (size_t i = 0; i < caps.size(); i++) {
    std::cout << "capacitance of " << i << "th bubble: " << caps[i] << std::endl;
  }
  
  std::cout << "Finished!" << std::endl;
  
  return 0;
}

void GenerateBezierPatch(const BezierPatchDescriptor& desc, TriangleMesh& mesh)
{
    auto idxOffset = mesh.vertices.size();

    const auto segsHorz = std::max(1u, desc.segments.x);
    const auto segsVert = std::max(1u, desc.segments.y);

    const auto invHorz  = Gs::Real(1) / static_cast<Gs::Real>(segsHorz);
    const auto invVert  = Gs::Real(1) / static_cast<Gs::Real>(segsVert);

    auto AddQuad = [&](unsigned int u, unsigned int v, VertexIndex i0, VertexIndex i1, VertexIndex i2, VertexIndex i3)
    {
        if (desc.backFacing)
            AddTriangulatedQuad(mesh, desc.alternateGrid, u, v, i1, i0, i3, i2, idxOffset);
        else
            AddTriangulatedQuad(mesh, desc.alternateGrid, u, v, i0, i1, i2, i3, idxOffset);
    };

    /* Generate vertices */
    static const Gs::Real delta = Gs::Real(0.01);

    Gs::Vector3 coord, normal;
    Gs::Vector2 texCoord;

    for (unsigned int i = 0; i <= segsVert; ++i)
    {
        for (unsigned int j = 0; j <= segsHorz; ++j)
        {
            /* Compute coordinate and texture-coordinate */
            texCoord.x = static_cast<Gs::Real>(j) * invHorz;
            texCoord.y = static_cast<Gs::Real>(i) * invVert;

            coord = desc.bezierPatch(texCoord.x, texCoord.y);
            
            /* Sample bezier patch to approximate normal */
            auto uOffset = (desc.bezierPatch(texCoord.x + delta, texCoord.y) - coord);
            auto vOffset = (desc.bezierPatch(texCoord.x, texCoord.y + delta) - coord);
            normal = Gs::Cross(uOffset, vOffset).Normalized();

            /* Add vertex */
            if (!desc.backFacing)
            {
                texCoord.y = Gs::Real(1) - texCoord.y;
                normal = -normal;
            }

            mesh.AddVertex(coord, normal, texCoord);
        }
    }

    /* Generate indices */
    const auto strideHorz = segsHorz + 1;

    for (unsigned int v = 0; v < segsVert; ++v)
    {
        for (unsigned int u = 0; u < segsHorz; ++u)
        {
            AddQuad(
                u, v,
                (  v   *strideHorz + u   ),
                ( (v+1)*strideHorz + u   ),
                ( (v+1)*strideHorz + u+1 ),
                (  v   *strideHorz + u+1 )
            );
        }
    }
}

			virtual void ConvexDecompResult(ConvexDecomposition::ConvexResult &result)
			{

				TriangleMesh* trimesh = new TriangleMesh();

				SimdVector3 localScaling(6.f,6.f,6.f);

				//export data to .obj
				printf("ConvexResult\n");
				if (mOutputFile)
				{
					fprintf(mOutputFile,"## Hull Piece %d with %d vertices and %d triangles.\r\n", mHullCount, result.mHullVcount, result.mHullTcount );

					fprintf(mOutputFile,"usemtl Material%i\r\n",mBaseCount);
					fprintf(mOutputFile,"o Object%i\r\n",mBaseCount);

					for (unsigned int i=0; i<result.mHullVcount; i++)
					{
						const float *p = &result.mHullVertices[i*3];
						fprintf(mOutputFile,"v %0.9f %0.9f %0.9f\r\n", p[0], p[1], p[2] );
					}

					//calc centroid, to shift vertices around center of mass
					centroid.setValue(0,0,0);
					if ( 1 )
					{
						const unsigned int *src = result.mHullIndices;
						for (unsigned int i=0; i<result.mHullTcount; i++)
						{
							unsigned int index0 = *src++;
							unsigned int index1 = *src++;
							unsigned int index2 = *src++;
							SimdVector3 vertex0(result.mHullVertices[index0*3], result.mHullVertices[index0*3+1],result.mHullVertices[index0*3+2]);
							SimdVector3 vertex1(result.mHullVertices[index1*3], result.mHullVertices[index1*3+1],result.mHullVertices[index1*3+2]);
							SimdVector3 vertex2(result.mHullVertices[index2*3], result.mHullVertices[index2*3+1],result.mHullVertices[index2*3+2]);
							vertex0 *= localScaling;
							vertex1 *= localScaling;
							vertex2 *= localScaling;
							centroid += vertex0;
							centroid += vertex1;
							centroid += vertex2;
							
						}
					}

					centroid *= 1.f/(float(result.mHullTcount) * 3);

					if ( 1 )
					{
						const unsigned int *src = result.mHullIndices;
						for (unsigned int i=0; i<result.mHullTcount; i++)
						{
							unsigned int index0 = *src++;
							unsigned int index1 = *src++;
							unsigned int index2 = *src++;


							SimdVector3 vertex0(result.mHullVertices[index0*3], result.mHullVertices[index0*3+1],result.mHullVertices[index0*3+2]);
							SimdVector3 vertex1(result.mHullVertices[index1*3], result.mHullVertices[index1*3+1],result.mHullVertices[index1*3+2]);
							SimdVector3 vertex2(result.mHullVertices[index2*3], result.mHullVertices[index2*3+1],result.mHullVertices[index2*3+2]);
							vertex0 *= localScaling;
							vertex1 *= localScaling;
							vertex2 *= localScaling;
							
							vertex0 -= centroid;
							vertex1 -= centroid;
							vertex2 -= centroid;

							trimesh->AddTriangle(vertex0,vertex1,vertex2);

							index0+=mBaseCount;
							index1+=mBaseCount;
							index2+=mBaseCount;
							
							fprintf(mOutputFile,"f %d %d %d\r\n", index0+1, index1+1, index2+1 );
						}
					}

					bool isDynamic = true;
					float mass = 1.f;
					CollisionShape* convexShape = new ConvexTriangleMeshShape(trimesh);
					SimdTransform trans;
					trans.setIdentity();
					trans.setOrigin(centroid);
					m_convexDemo->LocalCreatePhysicsObject(isDynamic, mass, trans,convexShape);

					mBaseCount+=result.mHullVcount; // advance the 'base index' counter.


				}
			}

void
makeTeapotScene()
{
	g_camera = new Camera;
	g_scene = new Scene;
	g_image = new Image;

	g_image->resize(512, 512);

	// set up the camera
	g_camera->setBGColor(Vector3(0.0f, 0.0f, 0.2f));
	g_camera->setEye(Vector3(0, 6, 2));
	g_camera->setLookAt(Vector3(0, 0, 0));
	g_camera->setUp(Vector3(0, 1, 0));
	g_camera->setFOV(45);

	// create and place a point light source
	PointLight * light = new PointLight;
	light->setPosition(Vector3(-2.5, 10, -2.5));
	light->setColor(Vector3(1, 1, 1));
	light->setWattage(4.0 * PI * 700);
	
	PointLight * light1 = new PointLight;
	light1->setPosition(Vector3(10, 15, 10));
	light1->setColor(Vector3(1, 1, 1));
	light1->setWattage(4.0 * PI * 700);
	
	g_scene->addLight(light);
	//g_scene->addLight(light1); //for test shadow

	Material* lambert = new Lambert(Vector3(1.0f));
	Material* mirror = new Mirror(Vector3(0.8f));
	Material* glass = new Glass(Vector3(0.95f), 1.5f);
	Material* brdf = new BPhong(Vector3(0.4f,0.4f,0.0f), Vector3(0.6f,0.6f,0.0f), 50);

	TriangleMesh * teapot = new TriangleMesh;
	Matrix4x4 xform;

	teapot->load("teapot.obj");
	addMeshTrianglesToScene(teapot, brdf);

	xform.setIdentity();
	xform *= translate(-5, 0, 0);
	teapot = new TriangleMesh;
	teapot->load("teapot.obj",xform);
	addMeshTrianglesToScene(teapot, lambert);

	xform.setIdentity();
	xform *= translate(-5, 0, -5);
	teapot = new TriangleMesh;
	teapot->load("teapot.obj", xform);
	addMeshTrianglesToScene(teapot, mirror);



	// create the floor triangle
	TriangleMesh * floor = new TriangleMesh;
	floor->createSingleTriangle();
	floor->setV1(Vector3(-100, 0, -100));
	floor->setV2(Vector3(  0, 0,  100));
	floor->setV3(Vector3( 100, 0, -100));
	floor->setN1(Vector3(0, 1, 0));
	floor->setN2(Vector3(0, 1, 0));
	floor->setN3(Vector3(0, 1, 0));
    
	Triangle* t = new Triangle;
	t->setIndex(0);
	t->setMesh(floor);
	t->setMaterial(lambert); 
	g_scene->addObject(t);




	// let objects do pre-calculations if needed
	g_scene->preCalc();
}

void ConvexDecompositionDemo::initPhysics(const char* filename)
{
	ConvexDecomposition::WavefrontObj wo;

	tcount = wo.loadObj(filename);

	CollisionDispatcher* dispatcher = new	CollisionDispatcher();


	SimdVector3 worldAabbMin(-10000,-10000,-10000);
	SimdVector3 worldAabbMax(10000,10000,10000);

	OverlappingPairCache* broadphase = new AxisSweep3(worldAabbMin,worldAabbMax);
	//OverlappingPairCache* broadphase = new SimpleBroadphase();

	m_physicsEnvironmentPtr = new CcdPhysicsEnvironment(dispatcher,broadphase);
	m_physicsEnvironmentPtr->setDeactivationTime(2.f);

	m_physicsEnvironmentPtr->setGravity(0,-10,0);

	SimdTransform startTransform;
	startTransform.setIdentity();
	startTransform.setOrigin(SimdVector3(0,-4,0));

	LocalCreatePhysicsObject(false,0,startTransform,new BoxShape(SimdVector3(30,2,30)));

	class MyConvexDecomposition : public ConvexDecomposition::ConvexDecompInterface
	{

		ConvexDecompositionDemo*	m_convexDemo;
		public:

		MyConvexDecomposition (FILE* outputFile,ConvexDecompositionDemo* demo)
			:m_convexDemo(demo),
				mBaseCount(0),
			mHullCount(0),
			mOutputFile(outputFile)

		{
		}
		
			virtual void ConvexDecompResult(ConvexDecomposition::ConvexResult &result)
			{

				TriangleMesh* trimesh = new TriangleMesh();

				SimdVector3 localScaling(6.f,6.f,6.f);

				//export data to .obj
				printf("ConvexResult\n");
				if (mOutputFile)
				{
					fprintf(mOutputFile,"## Hull Piece %d with %d vertices and %d triangles.\r\n", mHullCount, result.mHullVcount, result.mHullTcount );

					fprintf(mOutputFile,"usemtl Material%i\r\n",mBaseCount);
					fprintf(mOutputFile,"o Object%i\r\n",mBaseCount);

					for (unsigned int i=0; i<result.mHullVcount; i++)
					{
						const float *p = &result.mHullVertices[i*3];
						fprintf(mOutputFile,"v %0.9f %0.9f %0.9f\r\n", p[0], p[1], p[2] );
					}

					//calc centroid, to shift vertices around center of mass
					centroid.setValue(0,0,0);
					if ( 1 )
					{
						const unsigned int *src = result.mHullIndices;
						for (unsigned int i=0; i<result.mHullTcount; i++)
						{
							unsigned int index0 = *src++;
							unsigned int index1 = *src++;
							unsigned int index2 = *src++;
							SimdVector3 vertex0(result.mHullVertices[index0*3], result.mHullVertices[index0*3+1],result.mHullVertices[index0*3+2]);
							SimdVector3 vertex1(result.mHullVertices[index1*3], result.mHullVertices[index1*3+1],result.mHullVertices[index1*3+2]);
							SimdVector3 vertex2(result.mHullVertices[index2*3], result.mHullVertices[index2*3+1],result.mHullVertices[index2*3+2]);
							vertex0 *= localScaling;
							vertex1 *= localScaling;
							vertex2 *= localScaling;
							centroid += vertex0;
							centroid += vertex1;
							centroid += vertex2;
							
						}
					}

					centroid *= 1.f/(float(result.mHullTcount) * 3);

					if ( 1 )
					{
						const unsigned int *src = result.mHullIndices;
						for (unsigned int i=0; i<result.mHullTcount; i++)
						{
							unsigned int index0 = *src++;
							unsigned int index1 = *src++;
							unsigned int index2 = *src++;


							SimdVector3 vertex0(result.mHullVertices[index0*3], result.mHullVertices[index0*3+1],result.mHullVertices[index0*3+2]);
							SimdVector3 vertex1(result.mHullVertices[index1*3], result.mHullVertices[index1*3+1],result.mHullVertices[index1*3+2]);
//.........这里部分代码省略.........

void
makeBallScene()
{
	g_camera = new Camera;
	g_scene = new Scene;
	g_image = new Image;


	g_image->resize(512, 512);

	Material* lambert = new Lambert(Vector3(1.0f,0.5f,0.0f));
	Material* lambert_ground = new Lambert(Vector3(0.9f,0.9f,0.9f));
	Material* mirror = new Mirror(Vector3(0.6f));
	Material* glass = new Glass(Vector3(0.95f),1.5f);
	Material* brdf = new BPhong(Vector3(0.4f),Vector3(0.6f),200);

	// set up the camera
	g_camera->setBGColor(Vector3(0.0f, 0.0f, 0.0f));
	g_camera->setEye(Vector3( 10, 3, 10));
	g_camera->setLookAt(Vector3(0, 0, 0));
	g_camera->setUp(Vector3(0, 1, 0));
	g_camera->setFOV(45);

	// create and place a point light source
	PointLight * light = new PointLight;
	light->setPosition(Vector3(5, 15, 5));
	light->setColor(Vector3(1, 1, 1));
	light->setWattage(4.0 * PI * 700);

	g_scene->addLight(light);

	// create the floor triangle
	TriangleMesh * floor = new TriangleMesh;
	floor->createSingleTriangle();
	floor->setV1(Vector3(-100, 0, -100));
	floor->setV2(Vector3(0, 0, 100));
	floor->setV3(Vector3(100, 0, 0));
	floor->setN1(Vector3(0, 1, 0));
	floor->setN2(Vector3(0, 1, 0));
	floor->setN3(Vector3(0, 1, 0));


	TriangleMesh * right = new TriangleMesh;
	right->createSingleTriangle();
	right->setV1(Vector3(10, 0, -10));
	right->setV2(Vector3(0, 100, 0));
	right->setV3(Vector3(-10, 0, 10));
	right->setN1(Vector3(1, 0, 1));
	right->setN2(Vector3(1, 0, 1));
	right->setN3(Vector3(1, 0, 1));

	TriangleMesh * left = new TriangleMesh;
	left->createSingleTriangle();
	left->setV1(Vector3(0, 0, 0));
	left->setV2(Vector3(100, 0, 0));
	left->setV3(Vector3(0, 100, 0));
	left->setN1(Vector3(0, 0, 1));
	left->setN2(Vector3(0, 0, 1));
	left->setN3(Vector3(0, 0, 1));

	Triangle* t = new Triangle;
	t->setIndex(0);
	t->setMesh(floor);
	t->setMaterial(lambert_ground);

	Triangle* t2 = new Triangle;
	t2->setIndex(0);
	t2->setMesh(right);
	t2->setMaterial(lambert_ground);

	Triangle* t3 = new Triangle;
	t3->setIndex(0);
	t3->setMesh(left);
	t3->setMaterial(lambert_ground);
	

	g_scene->addObject(t);
	g_scene->addObject(t2);
	//g_scene->addObject(t3);

	//create two balls
	Sphere * sphere1 = new Sphere;
	sphere1->setCenter(Vector3(0.5, 1, 1));
	sphere1->setRadius(0.2);
	sphere1->setMaterial(lambert);
	Sphere * sphere2 = new Sphere;
	sphere2->setCenter(Vector3(1, 1, 0.5));
	sphere2->setRadius(0.2);
	sphere2->setMaterial(lambert);

	Sphere * sphere3= new Sphere;
	sphere3->setCenter(Vector3(4, 1, 2));
	sphere3->setRadius(1);
	sphere3->setMaterial(glass);

	Sphere * sphere4 = new Sphere;
	sphere4->setCenter(Vector3(2, 1, 4));
	sphere4->setRadius(1);
	sphere4->setMaterial(brdf);

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

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

	for(i = 0; i < numNodesSurf; i++)
	{
		for(j = 0; j < numElemsTet; j++)
		{
			g_Element elem;
			set<int>::iterator triangleIt;
			for(triangleIt = trianglesOnSurface[j].begin(); triangleIt != trianglesOnSurface[j].end(); triangleIt++)
				elem.node(volumetricMeshNodes[*triangleIt]);

			// g_Vector closestPoint = computeClosestPoint(rigidlyAlignedSurfaceTemplate->nodes()[i], &elem);
			// dist = closestPoint.DistanceTo(rigidlyAlignedSurfaceTemplate->nodes()[i]->coordinate());
		        g_Vector closestPoint;
		        dist = computeClosestPoint(*rigidlyAlignedSurfaceTemplate->nodes()[i], elem, closestPoint );
			dist = sqrt(dist);

			if((j == 0) || (dist < minDist))
			{
				minDist = dist;

				k = 0;
				for(triangleIt = trianglesOnSurface[j].begin(); triangleIt != trianglesOnSurface[j].end(); triangleIt++, k++)
					mappingIndicesSurf[3*i+k] = *triangleIt;

				vector<double> barycentrics = computeBarycentrics(rigidlyAlignedSurfaceTemplate->nodes()[i]->coordinate(), closestPoint, &elem);
				mappingWeightsSurf[3*i] = barycentrics[0];
				mappingWeightsSurf[3*i+1] = barycentrics[1];
				mappingWeightsSurf[3*i+2] = barycentrics[2];
				offsetsSurf[i] = barycentrics[3];
			}


			//---------------------------
			if(i == 0 && contactPoint != NULL){
			        g_Vector closestContactPoint;
				distContactPoint = computeClosestPoint(*contactPoint, elem, closestContactPoint );
				distContactPoint = sqrt(distContactPoint);
				// g_Vector closestContactPoint = computeClosestPoint(contactPoint, &elem);
				// distContactPoint = closestContactPoint.DistanceTo(contactPoint->coordinate());
				if((j == 0) || (distContactPoint < minDistContactPoint))
				{
					minDistContactPoint = distContactPoint;

					k = 0;
					for(triangleIt = trianglesOnSurface[j].begin(); triangleIt != trianglesOnSurface[j].end(); triangleIt++, k++)
						contactTriangleId[3*i+k] = *triangleIt;

					vector<double> barycentricsContactPoint = computeBarycentrics(contactPoint->coordinate()