CGAL users discussion list

[weihy1984 <>]

Hello, everyone,

   I have a problem with the output of the implicit surface mesh data to a txt file. And my CGAL version is 3.3.1 of windows XP. At first, I want to get the number and coordinates of the vertices. Furthermore, I want to know the number of the triangles on the implicit surface and the indexes of three vertices of these triangles. So I add some output codes to 3d implicit surface mesh generation example in  examples/Surface_mesher/mesh_an_implicit_function.cpp.

    But the result is not what I want. Many triangles I got were not on the surface. And I have spent four days to solve it, but there is no result. And I still don't know the reason why . I just be a cgal novice  and I need your help very much, and I will appreciate it for your help very much.

    The following is my code, you can check the output part.

  
#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
#include <fstream>

// default triangulation for Surface_mesher
typedef CGAL::Surface_mesh_default_triangulation_3 Tr;

// c2t3
typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;

typedef Tr::Geom_traits GT;
typedef GT::Sphere_3 Sphere_3;
typedef GT::Point_3 Point_3;
typedef GT::FT FT;
typedef Tr::Finite_facets_iterator Fi;
typedef Tr::Finite_cells_iterator Fc;

//typedef C2t3::Facet_iterator Fi;
typedef Tr::Finite_vertices_iterator Vi;
typedef C2t3::Vertex_handle Vertex_handle;

typedef FT (*Function)(Point_3);

typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;

FT sphere_function (Point_3 p) {
  const FT x2=p.x()*p.x(), y2=p.y()*p.y(), z2=p.z()*p.z();
  return x2+y2+z2-1;
}

int main() {

Tr tr; // 3D-Delaunay triangulation

C2t3 c2t3(tr); // 2D-complex in 3D-Delaunay triangulation

// defining the surface
Surface_3 surface(sphere_function, // pointer to function
                  Sphere_3(CGAL::ORIGIN, 2.)); // bounding sphere
// defining meshing criteria
CGAL::Surface_mesh_default_criteria_3<Tr> criteria(30., // angular bound
                                                   0.1, // radius bound
                                                   0.1); // distance bound
// meshing surface
CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Non_manifold_tag());

std::cout << "Final number of points: " << tr.number_of_vertices() << "\n";

 

// the following is the output part
std::ofstream os("out.mesh");
int n = tr.number_of_vertices();  
int m =c2t3.number_of_facets();

 

// outputs  number of vertices and the number of the triangles on the implicit surface
os << "OFF"<<std::endl<<n <<" "<< m<<" "<<0<<std::endl;

 

std::vector<Vertex_handle> TV(n);
int i = 0;

// write the vertices

for (Vi it=tr.finite_vertices_begin(); it!=tr.finite_vertices_end(); ++it)
    TV[i++] = it;

std::map<Vertex_handle, int > V;

for (i=0; i < n; i++) {
    os << *TV[i];
    V[TV[i]] = i;
   os << std::endl;

}

// These part is aimed to write the indexes of three vertices of the triangles on the implicit surface

for(Fi it =tr.finite_facets_begin(); it !=tr.finite_facets_end(); ++it) {

   // here I want to judge whether or not the facet is on the surface, is it right?
 if(c2t3.is_in_complex(*it)){
    int k=-1;
    for(int j = 0; j < 4; j++){

//we know that the type of a Facet is std::pair<Cell handle, int>, (c,j) is the facet of c opposite   

//to the vertex of index j, where c is a Cell_handle.   

     if ((*it).second!=j){
        os << V.find((*it).first->vertex(j))->second;
        k++;
        if ( k==2 ){
           os << std::endl;
          break;
        }

       else
         os <<  ' ';
     }
   }

  }

}
  os.close();
  return 0;
}

Best Wishes

Huayi

---

穿越地震带 纪念汶川地震一周年