diff --git a/.gitattributes b/.gitattributes index 3d84586e007..c4ef0e57ff5 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1002,6 +1002,8 @@ Mesh_2/doc_tex/Mesh_2_ref/part_of_a_cluster.gif -text svneol=unset#unset Mesh_2/test/Mesh_2/fish-and-rectangle.poly -text Mesh_2/test/Mesh_2/fish.edg -text Mesh_2/test/Mesh_2/fish.poly -text +Mesh_3/applications/GNUmakefile -text +Mesh_3/applications/medit_to_cgal.C -text Min_sphere_of_spheres_d/test_extensive/stability/maple/balls-on-boundary-2.mws -text Min_sphere_of_spheres_d/test_extensive/stability/maple/balls-on-boundary-3.mws -text Min_sphere_of_spheres_d/web/figs/heuristic/excess.eps -text diff --git a/.gitignore b/.gitignore index befeea9f1ee..be97ec6fe1d 100644 --- a/.gitignore +++ b/.gitignore @@ -73,6 +73,24 @@ Mesh_2/test/Mesh_2/test_double_map Mesh_2/test/Mesh_2/test_filtred_container Mesh_2/test/Mesh_2/test_meshing Mesh_3/*.tags.xml +Mesh_3/applications/*.cgal +Mesh_3/applications/*.exe +Mesh_3/applications/*.maillage +Mesh_3/applications/*.mesh +Mesh_3/applications/*.png +Mesh_3/applications/*.surface* +Mesh_3/applications/cgal_to_medit +Mesh_3/applications/depends +Mesh_3/applications/display_distribution +Mesh_3/applications/lanteri +Mesh_3/applications/lanteri_output_tet_mesh +Mesh_3/applications/medit_to_cgal +Mesh_3/applications/my_makefile +Mesh_3/applications/off_to_ghs +Mesh_3/applications/off_to_medit +Mesh_3/applications/read_mesh +Mesh_3/applications/slivers_exuder +Mesh_3/applications/stat_mesh Mesh_3/doxygen Mesh_3/examples/Mesh_3/chair-after.mesh Mesh_3/examples/Mesh_3/chair-after.png diff --git a/Mesh_3/applications/Distribution_displayer.h b/Mesh_3/applications/Distribution_displayer.h new file mode 100644 index 00000000000..7b790b5aa7c --- /dev/null +++ b/Mesh_3/applications/Distribution_displayer.h @@ -0,0 +1,18 @@ +#ifndef DISTRIBUTION_DISPLAYER_H +#define DISTRIBUTION_DISPLAYER_H +#include + +struct Distribution_displayer +{ + virtual void fill_rectangle(double x1, double y1, + double x2, double y2, + CGAL::Color c) = 0; + + virtual void segment(double x1, double y1, + double x2, double y2, + CGAL::Color c) = 0; + + virtual ~Distribution_displayer() {}; +}; + +#endif diff --git a/Mesh_3/applications/GNUmakefile b/Mesh_3/applications/GNUmakefile new file mode 100644 index 00000000000..d44133c8f8a --- /dev/null +++ b/Mesh_3/applications/GNUmakefile @@ -0,0 +1,28 @@ +include makefile + +#---------------------------------------------------------------------# +# +# dependencies +# +# if you want deps, create a file my_makefile that contains +# -include depends +# +#---------------------------------------------------------------------# + +-include my_makefile + +dep: + rm depends; $(MAKE) depends + +.PHONY: dep + +depends: *.C + cat /dev/null > depends + for f in *.C; do \ + echo >> depends; \ + echo >> depends; \ + echo "$${f%.C}$(OBJ_EXT): \\" >> depends; \ + $(CGAL_CXX) $(CXXFLAGS) -M -MG $$f \ + | grep '\.\..*/include/CGAL' >> depends; \ + done; \ + test -f depends diff --git a/Mesh_3/applications/Gd_displayer.C b/Mesh_3/applications/Gd_displayer.C new file mode 100644 index 00000000000..05bad871ec3 --- /dev/null +++ b/Mesh_3/applications/Gd_displayer.C @@ -0,0 +1,79 @@ +#include "Gd_displayer.h" + +Gd_displayer::Gd_displayer(int sx, int sy) +{ + im = gdImageCreate(sx, sy); + gdImageColorAllocate(im, 255, 255, 255); + set_window(0., 1., 0., 1.); +} + +Gd_displayer::Gd_displayer(gdImagePtr gd) +{ + im = gd; + set_window(0., 1., 0., 1.); +} + +Gd_displayer::~Gd_displayer() +{ + gdImageDestroy(im); +} + +void Gd_displayer::fill_rectangle(double x1, double y1, + double x2, double y2, + CGAL::Color c) +{ + gdImageFilledRectangle(im, + x_pixel(x1), y_pixel(y2), + x_pixel(x2), y_pixel(y1), + gd_color(c)); +} + +void Gd_displayer::segment(double x1, double y1, + double x2, double y2, + CGAL::Color c) +{ + gdImageLine(im, + x_pixel(x1), y_pixel(y1), + x_pixel(x2), y_pixel(y2), + gd_color(c)); +} + +void Gd_displayer::set_window(const double x_min, const double x_max, + const double y_min, const double y_max, + int x, int y, int sx, int sy) +{ + xmin = x_min; + xmax = x_max; + ymin = y_min; + ymax = y_max; + xclip = x; + yclip = y; + sxclip = sx; + syclip = sy; + gdImageSetClip(im, xclip, yclip, xclip+sxclip, yclip+syclip); + xscal = sxclip / (xmax-xmin); + yscal = syclip / (ymax-ymin); +} + +int Gd_displayer::x_pixel(double x) const +{ + return( xclip + static_cast((x-xmin)*xscal) ); +} + +int Gd_displayer::y_pixel(double y) const +{ + return( yclip - static_cast((y-ymax)*yscal) ); +} + +bool Gd_displayer::save_png(const char* filename) +{ + FILE *pngout = fopen(filename, "wb"); + if (pngout != NULL) + { + gdImagePng(im, pngout); + fclose(pngout); + return true; + } + else + return false; +} diff --git a/Mesh_3/applications/Gd_displayer.h b/Mesh_3/applications/Gd_displayer.h new file mode 100644 index 00000000000..553fcb996dd --- /dev/null +++ b/Mesh_3/applications/Gd_displayer.h @@ -0,0 +1,104 @@ +#include +#include "Distribution_displayer.h" + +struct Gd_displayer : public Distribution_displayer +{ + // --- CONSTRUCTORS --- + + /** Create a Gd_displayer using with a new GD image. + /param sx width of the image + /param sy height of the image + */ + Gd_displayer(int sx, int sy); + + /** Create a Gd_displayer using an existing GD Image. */ + Gd_displayer(gdImagePtr gd); + + ~Gd_displayer(); + + + /** \name VIRTUAL FUNCTIONS FROM Distribution_displayer */ + //@{ + void fill_rectangle(double x1, double y1, + double x2, double y2, + CGAL::Color c); + + void segment(double x1, double y1, + double x2, double y2, + CGAL::Color c); + //@} + + /** \name FUNCTIONS SPECIFIC TO Gd_displayer */ + //@{ + /** Set world coordinates of the image.*/ + inline + void set_window(const double x_min, + const double x_max, + const double y_min, + const double y_max) + { + set_window(x_min, x_max, y_min, y_max, + 0, 0, gdImageSX(im), gdImageSY(im)); + } + + /** Set world coordinates and a clip zone. + \param x x coordinate (in pixel) of the upper left point of + the clip zone + \param y y coordinate (in pixel) of the upper left point of + the clip zone + \param sx width of the clip zone + \param sy height of the clip zone + */ + void set_window(const double x_min, + const double x_max, + const double y_min, + const double y_max, + int x, int y, + int sx, int sy); + + /** Returns the x coordinate, in pixel, of the world-coordinate x. */ + int x_pixel(double x) const; + + /** Returns the y coordinate, in pixel, of the world-coordinate y. */ + int y_pixel(double y) const; + + /** Save the GD image to a PNG file. */ + bool save_png(const char* filename); + + /** Returns the GD image. */ + inline gdImagePtr image() + { + return im; + } + + /** Returns the index of the color c in the image palette. */ + inline + int gd_color(CGAL::Color c) + { + int i = gdImageColorExact(im, + c.red(), + c.green(), + c.blue()); + if( i < 0 ) + return gdImageColorAllocate(im, + c.red(), + c.green(), + c.blue()); + else + return i; + } + //@} + // end specific functions +private: + /** \name real dimensions + @{ */ + double xmin, xmax, ymin, ymax; //@} + /** \name clip zone + @{ */ + int xclip, yclip, sxclip, syclip; //@} + /** \name scaling factors + @{ */ + double xscal, yscal; //@} + + gdImagePtr im; /**< the GD image */ +}; diff --git a/Mesh_3/applications/Qt_widget_displayer.C b/Mesh_3/applications/Qt_widget_displayer.C new file mode 100644 index 00000000000..1526edfe798 --- /dev/null +++ b/Mesh_3/applications/Qt_widget_displayer.C @@ -0,0 +1,30 @@ +#include +#include + +#include + +#include "Qt_widget_displayer.h" + +typedef CGAL::Simple_cartesian Simple_kernel; +typedef Simple_kernel::Iso_rectangle_2 Rectangle_2; +typedef Simple_kernel::Segment_2 Segment_2; +typedef Simple_kernel::Point_2 Point_2; + +Qt_widget_displayer::Qt_widget_displayer(CGAL::Qt_widget* w) : widget(w) {} + +void Qt_widget_displayer::fill_rectangle(double x1, double y1, + double x2, double y2, + CGAL::Color color) +{ + *widget << CGAL::FillColor(color) + << color + << Rectangle_2(Point_2(x1, y1), Point_2(x2, y2)); +} + +void Qt_widget_displayer::segment(double x1, double y1, + double x2, double y2, + CGAL::Color color) +{ + *widget << color << Segment_2(Point_2(x1, y1), + Point_2(x2, y2)); +} diff --git a/Mesh_3/applications/Qt_widget_displayer.h b/Mesh_3/applications/Qt_widget_displayer.h new file mode 100644 index 00000000000..3f5224c8887 --- /dev/null +++ b/Mesh_3/applications/Qt_widget_displayer.h @@ -0,0 +1,20 @@ +#include "Distribution_displayer.h" + +namespace CGAL { + class Qt_widget; +} + +struct Qt_widget_displayer : public Distribution_displayer +{ + Qt_widget_displayer(CGAL::Qt_widget* widget); + + void fill_rectangle(double x1, double y1, + double x2, double y2, + CGAL::Color c); + + void segment(double x1, double y1, + double x2, double y2, + CGAL::Color c); +private: + CGAL::Qt_widget* widget; +}; diff --git a/Mesh_3/applications/README b/Mesh_3/applications/README new file mode 100644 index 00000000000..56c662164e1 --- /dev/null +++ b/Mesh_3/applications/README @@ -0,0 +1,193 @@ +This directory is a repository of tools, related to Mesh_3, but which +are not demos nor examples. Some of them (like slivers_exuder) will +perhaps be moved to demo/Mesh_3/ in the future. + +Here is the list of tools: + cgal_to_medit + medit_to_cgal + display_distribution + filter_remove_tets_from_medit + lanteri + lanteri_output_tet_mesh + off_to_ghs + read_mesh + slivers_exuder + stat_mesh + +Details follows... + +Files formats: +------------- +[too be written] + +filter_remove_tets_from_medit: +----------------------------- +This is a filter, written in Perl. + +Use it like that: + ./filter_remove_tets_from_medit < file.mesh > file-notets.mesh + +It takes a medit file as standard input, and returns to standard +output the same medit file, without its tetrahedra. + + +cgal_to_medit: +------------- +Usage: + ./cgal_to_medit INPUT OUPUT + + INPUT must be a file of format produced by the output operator of + CGAL::Triangulation_3, with points + CGAL::Weighted_point_with_surface_index and cells + CGAL::Mesh_3::Complex_2_in_triangulation_cell_base_3. + + OUTPUT will be a medit file. + +medit_to_cgal: +------------- +Usage: + ./medit_to_cgal INPUT OUPUT + + INPUT must be a medit file. + + OUTPUT will be a file of format produced by the output operator of + CGAL::Triangulation_3, with points + CGAL::Weighted_point_with_surface_index and cells + CGAL::Mesh_3::Complex_2_in_triangulation_cell_base_3. + +lanteri: +-------- + This tool is designed specificaly for answering to a request of + Stephane Lanteri, INRIA. It uses several files: + - utils.h contains templated functions + 'display_faces_counts', + 'display_facets_by_surface_indices_statistics', + 'display_vertices_by_surface_indices_statistics', + 'display_cells_by_volume_indices_statistics'. + - distribution.h and distribution.C declare and define functions + 'compute_distribution', + 'display_distribution'. + They use Distribution_displayer.h. + - Distribution_displayer.h defines a base class + 'Distribution_displayer', with two pure virtual functions: + 'fill_rectangle', + 'segment' + - Gd_displayer.h Qt_widget_displayer.h + Gd_displayer.C Qt_widget_displayer.C + implementes two different classes derived of + 'Distribution_displayer' (two implementation for displaying). + - lanteri_process_results.C defines three functions (declared in + lanteri.C): + 'process_cells', + 'process_volume_edges', + 'process_surface_edges'. + - lanteri_utils.h defines two templated functions: + 'scan_edges_and_process', + 'scan_cells_and_process'. + These two functions use functions of lanteri_process_results.C. + +lanteri_output_tet_mesh: +----------------------- + This tool is designed specificaly for answering to a request of + Stephane Lanteri, INRIA. + +display_distribution: +-------------------- +Usage: + ./display_distribution [--scale x] (--tets|--noboite|--mesh) FILE +Options: + --scale SCALE SCALE is a real number, which will be the + the vertical scaling of the histogram. + --tets FILE + --noboite FILE + --mesh FILE + Read input file FILE. + FILE must a .tets file, or a .noboite file, + or a .mesh file. + +This tool displays the radius-radius ratio distribution of a mesh file, +and save it in a png file named FILE-scale-SCALE.png. + +off_to_ghs: +---------- +Usage: +./off_to_ghs [FILE_BASENAME] + will convert FILE_BASENAME.off to FILE_BASENAME.faces + and FILE_BASENAME.points. + +off_to_medit: +---------- +This tool creates a medit files with 0 tetrahedra, from an .off file. + +Usage: +./off_to_medit [FILE_BASENAME] + will convert FILE_BASENAME.off to FILE_BASENAME.mesh. + +read_mesh: +--------- +Usage: + ./read_mesh FILE + + FILE must be a file of format produced by the output operator of + CGAL::Triangulation_3, with points + CGAL::Weighted_point_with_surface_index and cells + CGAL::Mesh_3::Complex_2_in_triangulation_cell_base_3. + +This tools read a .mesh file and displays faces counts, such as in: + Vertices: 3247 + Facets on surface: 6474 + Cells: 19670 + Cells in volume: 18335 + +This file uses utils.h, and types.h. + +stat_mesh: +--------- +Usage: + ./stat_mesh FILE + + FILE must be a file of format produced by the output operator of + CGAL::Triangulation_3, with points + CGAL::Weighted_point_with_surface_index and cells + CGAL::Mesh_3::Complex_2_in_triangulation_cell_base_3. + +Same as read_mesh, but displays more informations: + Vertices: 3247 + Facets on surface: 6474 + Cells: 19670 + Cells in volume: 18335 + + Statistics: +(vertices) + Vertices with index #0: 0 + Vertices with index #1: 450 + Vertices with index #2: 563 + Vertices with index #3: 705 + Vertices with index #4: 863 + Vertices with index #5: 666 +(facets) + Hybrid facets: 0 + Facets on surface #0: 0 + Facets on surface #1: 896 + Facets on surface #2: 1122 + Facets on surface #3: 1406 + Facets on surface #4: 1722 + Facets on surface #5: 1328 + +slivers_exuder: +-------------- +This tools exudes slivers from a cgal triangulation. + +Usage: + ./slivers_exuder INPUT OUTPUT bound + + INPUT must be a file of format produced by the output operator of + CGAL::Triangulation_3, with points + CGAL::Weighted_point_with_surface_index and cells + CGAL::Mesh_3::Complex_2_in_triangulation_cell_base_3. + + OUPUT is the name of the ouput file. It will be the same format + + 0<=bound<=1.0 is a bound on pumping process. + +Uses utils.h and weighted_types.h. diff --git a/Mesh_3/applications/cgal_to_medit.C b/Mesh_3/applications/cgal_to_medit.C new file mode 100644 index 00000000000..788f73affb1 --- /dev/null +++ b/Mesh_3/applications/cgal_to_medit.C @@ -0,0 +1,46 @@ +#include "weighted_types.h" + +// ios +#include +#include +#include +#include + +#include "utils.h" + +int main(int , char** argv) +{ + Tr tr; + C2T3 c2t3(tr); + + std::ifstream ifs(argv[1]); + std::ofstream ofs(argv[2]); + if( !ifs || !ofs ) + { + std::cerr << "Usage:\n" + << " " << argv[0] << " INPUT OUPUT\n" << +"\n" +" INPUT must be " << format_cgal_description << +"\n" +" OUTPUT will be a medit file.\n" +"\n"; + return EXIT_FAILURE; + } + + std::cout << " Reading " << argv[1] << std::endl; + + if( CGAL::Mesh_3::input_mesh(ifs, c2t3, + true, + &std::cerr) ) +// if( CGAL::input_pslg_from_medit(ifs, +// c2t3, +// true, // debug +// &std::cout) ) // debug to cout + { + std::cout << " Writing " << argv[2] << std::endl; + CGAL::output_to_medit(ofs, c2t3); + return EXIT_SUCCESS; + } + else + return EXIT_FAILURE; +} diff --git a/Mesh_3/applications/display_distribution.C b/Mesh_3/applications/display_distribution.C new file mode 100644 index 00000000000..64621ead75e --- /dev/null +++ b/Mesh_3/applications/display_distribution.C @@ -0,0 +1,412 @@ +#include +#include +#include +#include + +#include + +#ifndef CGAL_USE_QT +int main() +{ + std::cerr << "This tool requires Qt.\n"; + exit(1); +} +#else +#include +#include +#include + +struct K : public CGAL::Simple_cartesian {}; +typedef K::Point_3 Point_3; +typedef K::Tetrahedron_3 Tetrahedron_3; +typedef K::Point_2 Point_2; +typedef K::Segment_2 Segment_2; +typedef CGAL::Iso_rectangle_2 Rectangle_2; + +/// Global variables +typedef std::map String_options; +typedef std::map Double_options; + +String_options string_options; +Double_options double_options; + +void init_options() +{ + string_options["tets"] = ""; + string_options["noboite"] = ""; + string_options["mesh"] = ""; + double_options["scale"] = 1; +} + +void usage(char *argv0, std::string error = "") +{ + if( error != "" ) + std:: cerr << "Error: " << error << std::endl; + std::cerr << "Usage:\n " + << argv0 + << " [--scale x] (--tets|--noboite|--mesh) FILE\n" + << "Options:\n" + << " --scale SCALE " + << "SCALE is a real number, which will be the\n" + << " " + << "the vertical scaling of the histogram.\n" + << " --tets FILE\n" + << " --noboite FILE\n" + << " --mesh FILE\n" + << " " + << "Read input file FILE.\n" + << " " + << "FILE must a .tets file, or a .noboite file,\n" + << " or a .mesh file." + << std::endl; + exit(1); +} + +void parse_argv(int argc, char** argv, int extra_args = 0) +{ + if (argc >=(2 + extra_args)) + { + std::string arg = argv[1+extra_args]; + if( arg.substr(0, 2) == "--" ) + { + Double_options::iterator opt_it = + double_options.find(arg.substr(2, arg.length()-2)); + if( opt_it != double_options.end() ) + { + if( argc < (3 + extra_args) ) + usage(argv[0], + (arg + " must be followed by a double!").c_str()); + std::stringstream s; + double val; + s << argv[extra_args + 2]; + s >> val; + if( !s ) + usage(argv[0], ("Bad double after " + arg + "!").c_str()); + opt_it->second = val; + parse_argv(argc, argv, extra_args + 2); + } + else + { + String_options::iterator opt_it = + string_options.find(arg.substr(2, arg.length()-2)); + if( opt_it != string_options.end() ) + { + if( argc < (3 + extra_args) ) + usage(argv[0], + (arg + " must be followed by a string!").c_str()); + std::string s = argv[extra_args + 2]; + opt_it->second = s; + parse_argv(argc, argv, extra_args + 2); + } + else + usage(argv[0], ("Invalid option: " + arg + "!").c_str()); + } + } + } +} // end parse_argv + +void output_distribution_to_png(std::vector& elements, + double max, + const int number_of_classes, + std::string filename) +{ + const int number_of_cells = elements.size(); + + std::vector distribution(number_of_classes); + + for(int j=0;j(dratio*(double)number_of_classes); + } + distribution[index]++; + } + +// const int max_occurrence = *std::max_element(distribution.begin(), +// distribution.end()); + + CGAL::Qt_widget *widget = new CGAL::Qt_widget(); + qApp->setMainWidget(widget); + widget->resize(400, 400); + widget->set_window(0, 1, 0, 1, true); // x_min, x_max, + // y_min, y_max. + widget->show(); + + widget->lock(); + *widget << CGAL::FillColor(CGAL::Color(200, 200, 200)) + << CGAL::Color(200, 200, 200) + << Rectangle_2(Point_2(0, 0), Point_2(1,1)); + + if( number_of_classes == 0 ) return; + const double width = 1.0 / number_of_classes; + + const double scale = double_options["scale"]; + + *widget << CGAL::FillColor(CGAL::BLACK); + // *widget << Segment_2(Point_2(0., 0.), Point_2(1., 0.)); + for(int k=0;k0) + { + double height = ( (distribution[k]+0.)/number_of_cells ) * scale; + *widget << CGAL::BLACK + << Rectangle_2(Point_2(k*width, 0), + Point_2((k+1)*width, height)); + } + else + *widget << CGAL::RED << Segment_2(Point_2(k*width, 0), + Point_2((k+1)*width, 0)); + + widget->unlock(); + if( widget->get_pixmap().save( QString(filename.c_str()), + "PNG") ) + std::cerr << "Distribution saved to file " << filename + << std::endl; + else + { + std::cerr << "Error: cannot save distribution to file " + << filename << std::endl; + exit(1); + } + qApp->exec(); +} + +bool failed(const char* error) +{ + std::cerr << error << std::endl; + return false; +} + +bool read_tets(std::vector& elements, std::istream& in) +{ + // read header + int nb_vertices = 0; + int nb_tets = 0; + std::string head; + + in >> nb_vertices >> head; + if( !in || head != "vertices" ) + return false; + + in >> nb_tets >> head; + if( !in || head != "tets" ) + return false; + + std::vector points; + points.reserve(nb_vertices); + + // read points + for(int i=0;i> x >> y >> z; + if( !in ) + return false; + points.push_back(Point_3(x,y,z)); + } + + // read tets + for(int i=0;i> dummy >> i0 >> i1 >> i2 >> i3; + if( dummy != 4 || !in ) + return false; + Tetrahedron_3 t = Tetrahedron_3(points[i0], + points[i1], + points[i2], + points[i3]); + elements.push_back(radius_ratio(t)); + } + return true; +} + +bool read_mesh(std::vector& elements, std::istream& in) +{ + // Header. + std::string head; + int version; + in >> head >> version; + if( head != "MeshVersionFormatted" || + version != 1 || + ! in) + return failed("read_mesh: bad version"); + + int dimension; + in >> head >> dimension; + if( head != "Dimension" || + dimension!= 3 || + ! in) + return failed("read_mesh: bad dimension"); + + // Vertices + int number_of_vertices; + in >> head >> number_of_vertices; + if( head != "Vertices" || ! in ) + return failed("read_mesh: bad file (missing Vertices)"); + std::vector points; + points.reserve(number_of_vertices); + for(int i = 0; i < number_of_vertices; ++i) + { + int dummy_i; + in >> points[i] >> dummy_i; + if( !in ) + return failed("read_mesh: bad file (reading of vertices)"); + } + + // Facets + int number_of_facets_on_surface; + in >> head >> number_of_facets_on_surface; + if( !in || head != "Triangles" ) + return failed("read_mesh: bad file (missing Triangles)"); + for(int i = 0; i < 4 * number_of_facets_on_surface; ++i) + { + double dummy; + in >> dummy; + } + + // Tetrahedra + int number_of_cells; + in >> head >> number_of_cells; + if( !in || head != "Tetrahedra") + return failed("read_mesh: bad file (missing Tetrahedra)"); + for(int i = 0; i < number_of_cells; ++i) + { + int i0, i1, i2, i3, dummy; + in >> i0 >> i1 >> i2 >> i3 >> dummy; + if( !in ) + return failed("read_mesh: bad file (reading of cells)"); + const Tetrahedron_3 t = Tetrahedron_3(points[i0-1], + points[i1-1], + points[i2-1], + points[i3-1]); + elements.push_back(radius_ratio(t)); + } + in >> head; + if ( !in || head != "End") + return failed("read_mesh: bad file (missing End)"); + else + return true; +} + +bool read_noboite(std::vector& elements, std::istream& in) +{ + int nb_vertices = 0; + int nb_tets = 0; + int nb_input_points; + int dummy; + + in >> nb_tets >> nb_vertices >> nb_input_points + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy + >> dummy; + + if( !in ) + return false; + + std::vector points; + points.reserve(nb_vertices); + + elements.clear(); + elements.reserve(nb_tets); + + // read tets + std::vector tets; + tets.reserve(4 * nb_tets); + for(int i=0;i> i0 >> i1 >> i2 >> i3; + if( !in ) + return false; + tets.push_back(i0-1); + tets.push_back(i1-1); + tets.push_back(i2-1); + tets.push_back(i3-1); + } + + + // read points + for(int i=0;i> x >> y >> z; + if( !in ) + return false; + points.push_back(Point_3(x,y,z)); + } + + // compute elements + for(int i = 0; i < 4 * nb_tets; i += 4) + { + const Tetrahedron_3 t = Tetrahedron_3(points[tets[i]], + points[tets[i+1]], + points[tets[i+2]], + points[tets[i+3]]); + elements.push_back(radius_ratio(t)); + } + return true; +} + +int main(int argc, char** argv) +{ + QApplication app(argc, argv); + init_options(); + parse_argv(argc, argv, 0); + + bool ghs = false; + bool tets = false; + bool mesh = false; + std::string filename = ""; + std::string ghs_filename = string_options["noboite"]; + std::string tets_filename = string_options["tets"]; + std::string mesh_filename = string_options["mesh"]; + + if(ghs_filename != "") + { + ghs = true; + filename = ghs_filename; + } + if(mesh_filename != "") + { + mesh = true; + filename = mesh_filename; + } + if(tets_filename != "") + { + tets = true; + filename = tets_filename; + } + + std::vector elements; + std::ifstream in_file(filename.c_str()); + if(tets) + tets = read_tets(elements, in_file); + if(mesh) + mesh = read_mesh(elements, in_file); + if(ghs) + ghs = read_noboite(elements, in_file); + std::stringstream png_name; + png_name << filename << "-scale-" << double_options["scale"] + << ".png"; + if(tets || mesh || ghs) + output_distribution_to_png(elements, 1., 100, png_name.str()); + else + usage(argv[0], "Error: cannot read file " + filename + "!"); +} +#endif // CGAL_USE_QT diff --git a/Mesh_3/applications/distribution.C b/Mesh_3/applications/distribution.C new file mode 100644 index 00000000000..42fb5678dde --- /dev/null +++ b/Mesh_3/applications/distribution.C @@ -0,0 +1,43 @@ +#include "distribution.h" + +void compute_distribution(const Qualities& qualities, + const double max_quality, + Distribution& distribution) +{ + const int number_of_classes = distribution.size(); + + const int qualities_size = qualities.size(); + + for(int j = 0; j < qualities_size; ++j) + { + if(qualities[j] < max_quality) + { + ++distribution[static_cast((qualities[j]/max_quality) + *number_of_classes)]; + } + } +} + +void display_distribution(Distribution_displayer* display, + const Distribution& distribution, + const double echelle) +{ + const int number_of_classes = distribution.size(); + + if( number_of_classes == 0 ) return; + const double width = 1.0 / number_of_classes; + + display->fill_rectangle(0., 0., 1., 1., CGAL::Color(200, 200, 200)); + for(int k = 0; k < number_of_classes; ++k) + if(distribution[k]>0) + { + const double height = ( distribution[k]+0. ) * echelle; + display->fill_rectangle(k * width, 0, + (k+1)* width, height, + CGAL::BLACK); + } + else + display->segment(k * width, 0., + (k+1) * width, 0., + CGAL::RED); +} diff --git a/Mesh_3/applications/distribution.h b/Mesh_3/applications/distribution.h new file mode 100644 index 00000000000..de300cba15d --- /dev/null +++ b/Mesh_3/applications/distribution.h @@ -0,0 +1,13 @@ +#include +#include "Distribution_displayer.h" + +typedef std::vector Qualities; +typedef std::vector Distribution; + +void compute_distribution(const Qualities& qualities, + const double max_quality, + Distribution& distribution); + +void display_distribution(Distribution_displayer* display, + const Distribution& distribution, + const double echelle); diff --git a/Mesh_3/applications/dump_to_ghs_points b/Mesh_3/applications/dump_to_ghs_points new file mode 100755 index 00000000000..073c83af814 --- /dev/null +++ b/Mesh_3/applications/dump_to_ghs_points @@ -0,0 +1,11 @@ +#!/usr/bin/perl -w + +my $n = <>; +chop $n; +print "$n\n"; + +while(<>) +{ + chop; + print; print " 0\n"; +} diff --git a/Mesh_3/applications/filter_remove_tets_from_medit b/Mesh_3/applications/filter_remove_tets_from_medit new file mode 100755 index 00000000000..76e35ae41f2 --- /dev/null +++ b/Mesh_3/applications/filter_remove_tets_from_medit @@ -0,0 +1,14 @@ +#!/usr/bin/perl -W -T + +$stop = 0; + +while( $stop == 0 ) +{ + $line = <>; + print $line; + if($line =~ /Tetrahedra/) + { + $stop = 1; + print "0\n" + } +} diff --git a/Mesh_3/applications/lanteri.C b/Mesh_3/applications/lanteri.C new file mode 100644 index 00000000000..bb9b7d57476 --- /dev/null +++ b/Mesh_3/applications/lanteri.C @@ -0,0 +1,229 @@ +#include + +// regular +#include +#include +#include + +// IO.h must be included before vertex and cell bases. +#include + +// vertex and cell bases +#include +#include +#include + +// c2t3 +#include + +// traits class for reading meshes +#include + +// radius_ratio +#include + +#include + +// traits class +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Regular_triangulation_euclidean_traits_3 Regular_traits; +typedef CGAL::Weighted_point_with_surface_index_geom_traits My_traits; + +// vertex and cell types +typedef CGAL::Surface_mesh_vertex_base_3 Vb; +typedef CGAL::Triangulation_cell_base_3 Cb1; +typedef CGAL::Surface_mesh_cell_base_3 Cb2; +typedef CGAL::Volume_mesher_cell_base_3 Cb3; + +template > +class Cell_with_volume_index : public Cb +{ +private: + int volume; + +public: + typedef typename GT::Point_3 Point; + typedef typename Cb::Triangulation_data_structure Tds; + typedef typename Tds::Vertex_handle Vertex_handle; + typedef typename Tds::Cell_handle Cell_handle; + + template < class TDS3 > + struct Rebind_TDS { + typedef typename Cb::template Rebind_TDS::Other Cb3; + typedef Cell_with_volume_index Other; + }; + + // Constructors + Cell_with_volume_index() : Cb(), volume(-1) + { + } + Cell_with_volume_index (Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2, + Vertex_handle v3) + : Cb (v0, v1, v2, v3), volume(-1) + { + } + Cell_with_volume_index (Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2, + Vertex_handle v3, + Cell_handle n0, + Cell_handle n1, + Cell_handle n2, + Cell_handle n3) + : Cb (v0, v1, v2, v3, n0, n1, n2, n3), volume(-1) + { + } + + // volume index + int volume_index() const + { + return volume; + } + + void set_volume_index(const int i) + { + volume = i; + } + +#ifdef CGAL_MESH_3_IO_H + static + std::string io_signature() + { + return CGAL::Get_io_signature()(); + } +#endif + +}; // end template class Cell_with_volume_index + +template +inline +std::istream& +operator>>(std::istream &is, Cell_with_volume_index& c) +{ + return is >> static_cast(c); +} + +template +inline +std::ostream& +operator<<(std::ostream &os, const Cell_with_volume_index& c) +{ + return os << static_cast(c); +} + +typedef Cell_with_volume_index Cb; + +// triangulation +typedef CGAL::Triangulation_data_structure_3 Tds; + +typedef CGAL::Regular_triangulation_3 Tr; + +// c2t3 +typedef CGAL::Complex_2_in_triangulation_3 C2T3; + +// ios +#include +#include +#include +#include +#include + +#include "utils.h" +#include "distribution.h" + +bool process_cells(const std::vector& volume_cells_quality, + const std::string filename_prefix); +bool process_volume_edges(const std::vector& volume_edges_length, + const std::vector& length_bounds, + const std::string filename_prefix); +bool process_surface_edges(const std::vector& surface_edges_length, + const std::vector& length_bounds, + const std::string filename_prefix); + +#include "lanteri_utils.h" + +typedef Tr::Vertex_handle Vertex_handle; + +int main(int , char**) +{ + Tr tr; + C2T3 c2t3(tr); + + double r1, r2, r3, r4, r5; + std::vector size_bounds(5); + std::vector radii(5); + + std::cout << "Input r1, r2, r3, r4, r5:" << std::endl; + std::cin >> r1 >> r2 >> r3 >> r4 >> r5; + std::cout << "Input the corresponding 5 size bounds:" << std::endl; + std::cin >> size_bounds[0] + >> size_bounds[1] + >> size_bounds[2] + >> size_bounds[3] + >> size_bounds[4]; + if(!std::cin) + return EXIT_FAILURE; + + std::string filename; + std::cout << "Input filename (with extension):" << std::endl; + std::cin >> filename; + std::ifstream ifs((filename+".cgal").c_str()); + if( !ifs || !std::cin) + { + return EXIT_FAILURE; + } + + std::cout << " Reading " << (filename+".cgal") << std::endl; + if( ! CGAL::Mesh_3::input_mesh(ifs, + c2t3, + true, // debug + &std::cerr) ) // debug to cerr + return EXIT_FAILURE; + + display_faces_counts(tr, " ", &std::cout); + + std::cout << "\n Statistics:\n"; + + std::cout << "(vertices)\n"; + display_vertices_by_surface_indices_statistics(tr, " ", &std::cout); + + std:: cout << "(facets)\n"; + display_facets_by_surface_indices_statistics(c2t3, " ", &std::cout); + + // sets volume indices + for(Tr::Finite_cells_iterator cit = tr.finite_cells_begin(); + cit != tr.finite_cells_end(); + ++cit) + if(cit->is_in_domain()) + { + const double sq_r = + squared_distance(K::Point_3(0, 0, 0), + static_cast(tr.dual(cit))); + + if( sq_r < r1*r1 ) + cit->set_volume_index(1); + else if( sq_r < r2*r2 ) + cit->set_volume_index(2); + else if( sq_r < r3*r3 ) + cit->set_volume_index(3); + else if( sq_r < r4*r4 ) + cit->set_volume_index(4); + else if( sq_r < r5*r5 ) + cit->set_volume_index(5); + } + + std::cout << "(cells)\n"; + display_cells_by_volume_indices_statistics(tr, " ", &std::cout); + + std::cout << "\n" + << " Edges lengths:\n"; + if(!scan_edges_and_process(tr, size_bounds, filename, " ", &std::cout)) + return EXIT_FAILURE; + + if(!scan_cells_and_process(tr, filename)) + return EXIT_FAILURE; + + return EXIT_SUCCESS; +} diff --git a/Mesh_3/applications/lanteri_output_tet_mesh.C b/Mesh_3/applications/lanteri_output_tet_mesh.C new file mode 100644 index 00000000000..034ea30a8a9 --- /dev/null +++ b/Mesh_3/applications/lanteri_output_tet_mesh.C @@ -0,0 +1,315 @@ +#include + +// regular +#include +#include +#include + +// IO.h must be included before vertex and cell bases. +#include + +// vertex and cell bases +#include +#include +#include + +// c2t3 +#include + +// traits class for reading meshes +#include + +// radius_ratio +#include + +// Point_traits +#include + +#include + +// traits class +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Regular_triangulation_euclidean_traits_3 Regular_traits; +typedef CGAL::Weighted_point_with_surface_index_geom_traits My_traits; + +// vertex and cell types +typedef CGAL::Surface_mesh_vertex_base_3 Vb; +typedef CGAL::Regular_triangulation_cell_base_3 Cb1; +typedef CGAL::Surface_mesh_cell_base_3 Cb2; +typedef CGAL::Volume_mesher_cell_base_3 Cb3; + +template > +class Cell_with_volume_index : public Cb +{ +private: + int volume; + +public: + typedef typename GT::Point_3 Point; + typedef typename Cb::Triangulation_data_structure Tds; + typedef typename Tds::Vertex_handle Vertex_handle; + typedef typename Tds::Cell_handle Cell_handle; + + template < class TDS3 > + struct Rebind_TDS { + typedef typename Cb::template Rebind_TDS::Other Cb3; + typedef Cell_with_volume_index Other; + }; + + // Constructors + Cell_with_volume_index() : Cb(), volume(0) + { + } + Cell_with_volume_index (Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2, + Vertex_handle v3) + : Cb (v0, v1, v2, v3), volume(0) + { + } + Cell_with_volume_index (Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2, + Vertex_handle v3, + Cell_handle n0, + Cell_handle n1, + Cell_handle n2, + Cell_handle n3) + : Cb (v0, v1, v2, v3, n0, n1, n2, n3), volume(0) + { + } + + // volume index + int volume_index() const + { + return volume; + } + + void set_volume_index(const int i) + { + volume = i; + } + +#ifdef CGAL_MESH_3_IO_H + static + std::string io_signature() + { + return CGAL::Get_io_signature()(); + } +#endif +}; // end template class Cell_with_volume_index + +template +inline +std::istream& +operator>>(std::istream &is, Cell_with_volume_index& c) +{ + return is >> static_cast(c); +} + +template +inline +std::ostream& +operator<<(std::ostream &os, const Cell_with_volume_index& c) +{ + return os << static_cast(c); +} + +typedef Cell_with_volume_index Cb; + +// triangulation +typedef CGAL::Triangulation_data_structure_3 Tds; +typedef CGAL::Regular_triangulation_3 Tr; + +// c2t3 +typedef CGAL::Complex_2_in_triangulation_3 C2T3; + +// ios +#include +#include +#include +#include +#include +#include + + +int main(int , char**) +{ + Tr tr; + C2T3 c2t3(tr); + + double r1, r2, r3, r4, r5; + std::vector size_bounds(5); + std::vector radii(5); + + std::cout << "Input r1, r2, r3, r4, r5:" << std::endl; + std::cin >> r1 >> r2 >> r3 >> r4 >> r5; + std::cout << "Input the corresponding 5 size bounds:" << std::endl; + std::cin >> size_bounds[0] + >> size_bounds[1] + >> size_bounds[2] + >> size_bounds[3] + >> size_bounds[4]; + if(!std::cin) + return EXIT_FAILURE; + + std::string file_prefix; + std::cout << "Input filename prefix:" << std::endl; + std::cin >> file_prefix; + + // input file + std::ifstream ifs((file_prefix + ".cgal").c_str()); + + // ouput tets mesh file + std::ofstream ofs_maillage((file_prefix + ".maillage").c_str()); + + // output five surface meshes files + std::vector ofs_surfaces(6); + // ofs_surfaces[0] will not be used. + + for(int i = 1; i <= 5; ++i) + { + std::stringstream str_stream; + str_stream << file_prefix << ".surface" << i; + + ofs_surfaces[i] = new std::ofstream(str_stream.str().c_str()); + } + + std::cout << " Reading file " << (file_prefix + ".cgal") << std::endl; + + if( CGAL::Mesh_3::input_mesh(ifs, c2t3, + true, + &std::cerr) ) + { + std::cout << " Writing file " << (file_prefix + ".maillage") << std::endl + << " and files " << (file_prefix + ".surface*") << std::endl; + + typedef C2T3::Triangulation Tr; + typedef Tr::Finite_cells_iterator Finite_cells_iterator; + typedef Tr::Finite_facets_iterator Finite_facets_iterator; + typedef Tr::Finite_vertices_iterator Finite_vertices_iterator; + typedef Tr::Vertex_handle Vertex_handle; + typedef Tr::Point Point; + typedef CGAL::Point_traits P_traits; + typedef P_traits::Bare_point Bare_point; + + // sets volume indices + for(Tr::Finite_cells_iterator cit = tr.finite_cells_begin(); + cit != tr.finite_cells_end(); + ++cit) + if(cit->is_in_domain()) + { + const double sq_r = + squared_distance(K::Point_3(0, 0, 0), + static_cast(tr.dual(cit))); + + if( sq_r < r1*r1 ) + cit->set_volume_index(1); // brain + else if( sq_r < r2*r2 ) + cit->set_volume_index(2); // LCR + else if( sq_r < r3*r3 ) + cit->set_volume_index(3); // head skull + else if( sq_r < r4*r4 ) + cit->set_volume_index(4); // skin + else if( sq_r < r5*r5 ) + cit->set_volume_index(-1);// air + } + + const Tr& tr = c2t3.triangulation(); + + // Headers + + ofs_maillage << tr.number_of_vertices() << " " + << CGAL::number_of_cells_in_domain(tr) << " " + << CGAL::number_of_facets_on_surface_with_index(c2t3, 5) + << std::endl; + + for(int i = 1; i <= 5; ++i) + { + *ofs_surfaces[i] << tr.number_of_vertices() << " " + << CGAL::number_of_facets_on_surface_with_index(c2t3, i) + << std::endl; + } + + // precision + ofs_maillage << std::setprecision(20); + for(int i = 1; i <= 5; ++i) + { + *ofs_surfaces[i] << std::setprecision(20); + } + + // Vertices + + std::map V; // vertices are counter from 1 + int inum = 1; + for( Finite_vertices_iterator vit = tr.finite_vertices_begin(); + vit != tr.finite_vertices_end(); + ++vit) + { + V[vit] = inum++; + Point p = static_cast(vit->point()); + ofs_maillage << p.x() << " " << p.y() << " " << p.z() + << std::endl; + for(int i = 1; i <= 5; ++i) + { + *ofs_surfaces[i] << p.x() << " " << p.y() << " " << p.z() + << std::endl; + } + } + + // Tetrahedra + + for( Finite_cells_iterator cit = tr.finite_cells_begin(); + cit != tr.finite_cells_end(); ++cit) + if( cit->is_in_domain() ) + { + for (int i=0; i<4; i++) + ofs_maillage << V[cit->vertex(i)] << " "; + + ofs_maillage << cit->volume_index() + << std::endl; + } + + // Facets + for( Finite_facets_iterator fit = tr.finite_facets_begin(); + fit != tr.finite_facets_end(); ++fit) + { + int surface_index = 0; + if (c2t3.face_status(fit->first,fit->second) + != C2T3::NOT_IN_COMPLEX) + { + for (int i=0; i<4; i++) + if (i != (*fit).second) + { + const Vertex_handle& vh = (*fit).first->vertex(i); + surface_index = vh->point().surface_index(); + if(surface_index == 5) + { + ofs_maillage << V[vh] << " "; + } + *ofs_surfaces[surface_index] << V[vh] << " "; + } + if(surface_index == 5) + ofs_maillage << "4" + << std::endl; + *ofs_surfaces[surface_index] << surface_index + << std::endl; + } + } + + for(int i = 1; i <= 5; ++i) + if( ofs_surfaces[i]->bad() ) + return EXIT_FAILURE; + + if( ofs_maillage.good() ) + { + for(int i = 1; i <= 5; ++i) + delete ofs_surfaces[i]; + return EXIT_SUCCESS; + } + else + return EXIT_FAILURE; + } + else + return EXIT_FAILURE; +} diff --git a/Mesh_3/applications/lanteri_process_results.C b/Mesh_3/applications/lanteri_process_results.C new file mode 100644 index 00000000000..6b676faa817 --- /dev/null +++ b/Mesh_3/applications/lanteri_process_results.C @@ -0,0 +1,117 @@ +#include "distribution.h" +#include "Gd_displayer.h" + +#include +#include + +const int global_number_of_classes = 20; + +bool process_aux_1(const std::vector& qualities, + const std::string filename, + const int number_of_classes) +{ + Gd_displayer main_display(1000, 200); + + main_display.set_window(-1.1, 1.1, -1.1, 1.1, + 0, 0, 200, 200); + + std::vector displays(5); + std::vector distributions(5); + + for(int i = 0; i < 5; ++i) + { + distributions[i].resize(number_of_classes); + + compute_distribution(qualities[i+1], + 1., + distributions[i]); + + displays[i] = new Gd_displayer(main_display.image()); + displays[i]->set_window(-.1, 1.1, -.1, 1.1, + i * 200, 0, 200, 200); + + const int max = *(std::max_element(distributions[i].begin(), + distributions[i].end())); + + display_distribution(displays[i], + distributions[i], + 1. / max); + } + + return main_display.save_png(filename.c_str()); +} + +bool process_aux_2(const std::vector& qualities, + const std::vector& length_bounds, + const std::string filename, + const int number_of_classes) +{ + Gd_displayer main_display(1000, 200); + + main_display.set_window(-1.1, 1.1, -1.1, 1.1, + 0, 0, 200, 200); + + std::vector displays(5); + std::vector distributions(5); + + for(int i = 0; i < 5; ++i) + { + distributions[i].resize(number_of_classes); + + double max_quality = *(std::max_element(qualities[i+1].begin(), + qualities[i+1].end())); + + max_quality = std::max(max_quality, length_bounds[i]); + + compute_distribution(qualities[i+1], + max_quality, + distributions[i]); + + displays[i] = new Gd_displayer(main_display.image()); + displays[i]->set_window(-.1, 1.1, -.1, 1.1, + i * 200, 0, 200, 200); + + const int max = *(std::max_element(distributions[i].begin(), + distributions[i].end())); + + display_distribution(displays[i], + distributions[i], + 1. / max); + + double x_position_of_length_bound = length_bounds[i] / max_quality; + + displays[i]->segment(x_position_of_length_bound, 0.0, + x_position_of_length_bound, -0.05, + CGAL::BLUE); + } + + return main_display.save_png(filename.c_str()); +} + +bool process_cells(const std::vector& volume_cells_quality, + const std::string filename_prefix) +{ + return process_aux_1(volume_cells_quality, + filename_prefix + "_cells_radius_radius_ratio.png", + global_number_of_classes); +} + +bool process_volume_edges(const std::vector& volume_edges_length, + const std::vector& length_bounds, + const std::string filename_prefix) +{ + return process_aux_2(volume_edges_length, + length_bounds, + filename_prefix + "_volume_edges_lengths.png", + global_number_of_classes); +} + +bool process_surface_edges(const std::vector& surface_edges_length, + const std::vector& length_bounds, + const std::string filename_prefix) +{ + return process_aux_2(surface_edges_length, + length_bounds, + filename_prefix + "_surface_edges_lengths.png", + global_number_of_classes); +} diff --git a/Mesh_3/applications/lanteri_utils.h b/Mesh_3/applications/lanteri_utils.h new file mode 100644 index 00000000000..b40c95abf9b --- /dev/null +++ b/Mesh_3/applications/lanteri_utils.h @@ -0,0 +1,119 @@ +#include + +template +bool +scan_edges_and_process(const Tr& tr, + std::vector length_bounds, + std::string filename_prefix, + std::string prefix = "", + // prefix to each line output + std::ostream* out_stream = + &std::cout + // output stream + ) +{ + // reminder: Qualities is std::vector (voir "distribution.h") + std::vector surface_edges_length; + std::vector volume_edges_length; + + std::vector surface_edges_length_max; + std::vector volume_edges_length_max; + + for(typename Tr::Finite_edges_iterator fit = tr.finite_edges_begin(); + fit!=tr.finite_edges_end(); + ++fit) + { + const typename Tr::Vertex_handle& va = fit->first->vertex(fit->second); + const typename Tr::Vertex_handle& vb = fit->first->vertex(fit->third); + + const double length = + CGAL::sqrt(CGAL::to_double(squared_distance(va->point(), + vb->point()))); + + const unsigned int& index_a = va->point().surface_index(); + const unsigned int& index_b = vb->point().surface_index(); + + if( index_a != 0 && index_a == index_b ) // surface edge + { + if( surface_edges_length.size() <= index_a ) + { + surface_edges_length.resize(index_a+1); + surface_edges_length_max.resize(index_a+1); + } + if( length > surface_edges_length_max[index_a] ) + surface_edges_length_max[index_a] = length; + + surface_edges_length[index_a].push_back(length); + } + else // volume edge + { + const int index = fit->first->volume_index(); + + if(index >=0) + { + const unsigned int positive_index = index; + if( volume_edges_length.size() <= positive_index ) + { + volume_edges_length.resize(positive_index+1); + volume_edges_length_max.resize(positive_index+1); + } + if( length > volume_edges_length_max[positive_index] ) + volume_edges_length_max[positive_index] = length; + volume_edges_length[positive_index].push_back(length); + } + } + } + + const typename Qualities::size_type surface_vector_size = + surface_edges_length_max.size(); + for(unsigned int i = 0; i < surface_vector_size; ++i) + { + *out_stream << prefix + << "Maximum length for edges on surface #" << i << ": " + << surface_edges_length_max[i] << std::endl; + } + + const typename Qualities::size_type volume_vector_size = + volume_edges_length_max.size(); + for(unsigned int i = 0; i < volume_vector_size; ++i) + { + *out_stream << prefix + << "Maximum length for edges in volume #" << i << ": " + << volume_edges_length_max[i] << std::endl; + } + + return process_surface_edges(surface_edges_length, + length_bounds, + filename_prefix) && + process_volume_edges(volume_edges_length, + length_bounds, + filename_prefix); +} + +template +bool +scan_cells_and_process(const Tr& tr, std::string filename_prefix) +{ + std::vector volume_cells_quality; + + for(typename Tr::Finite_cells_iterator cit = tr.finite_cells_begin(); + cit != tr.finite_cells_end(); + ++cit) + if(cit->is_in_domain()) + { + const double quality = + CGAL::to_double(radius_ratio(tr.tetrahedron(cit))); + // radius ratio is in common namespace, in Slivers_exuder.h + int index = cit->volume_index(); + if(index < 0) + index = 0; + + const unsigned int positive_index = index; + + if( volume_cells_quality.size() <= positive_index ) + volume_cells_quality.resize(positive_index+1); + volume_cells_quality[positive_index].push_back(quality); + } + + return process_cells(volume_cells_quality, filename_prefix); +} diff --git a/Mesh_3/applications/makefile b/Mesh_3/applications/makefile new file mode 100644 index 00000000000..cde99223c76 --- /dev/null +++ b/Mesh_3/applications/makefile @@ -0,0 +1,116 @@ +# Created by the script cgal_create_makefile +# This is the makefile for compiling a CGAL application. + +#---------------------------------------------------------------------# +# include platform specific settings +#---------------------------------------------------------------------# +# Choose the right include file from the /make directory. + +# CGAL_MAKEFILE = ENTER_YOUR_INCLUDE_MAKEFILE_HERE +include $(CGAL_MAKEFILE) + +#---------------------------------------------------------------------# +# compiler flags +#---------------------------------------------------------------------# + +CXXFLAGS = -I../../Triangulation_3/include \ + -I../../Mesh_2/include \ + -I../include \ + -I../../Data_structure_for_queries_3/include \ + -I../../Surface_mesher/include \ + $(CGAL_CXXFLAGS) \ + $(LONG_NAME_PROBLEM_CXXFLAGS) \ + $(OTHER_CXXFLAGS) + +#---------------------------------------------------------------------# +# linker flags +#---------------------------------------------------------------------# + +LIBPATH = \ + $(CGAL_LIBPATH) + +LDFLAGS = \ + $(LONG_NAME_PROBLEM_LDFLAGS) \ + $(CGAL_QT_LDFLAGS) + +#---------------------------------------------------------------------# +# target entries +#---------------------------------------------------------------------# + +all: \ + display_distribution$(EXE_EXT) \ + off_to_ghs$(EXE_EXT) \ + off_to_medit$(EXE_EXT) \ + lanteri$(EXE_EXT) \ + lanteri_output_tet_mesh$(EXE_EXT) \ + read_mesh$(EXE_EXT) \ + slivers_exuder$(EXE_EXT) \ + cgal_to_medit$(EXE_EXT) \ + medit_to_cgal$(EXE_EXT) \ + stat_mesh$(EXE_EXT) + +cgal_to_medit$(EXE_EXT): cgal_to_medit$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)cgal_to_medit cgal_to_medit$(OBJ_EXT) $(LDFLAGS) + +medit_to_cgal$$(EXE_EXT): medit_to_cgal$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)medit_to_cgal medit_to_cgal$(OBJ_EXT) $(LDFLAGS) + +display_distribution$(EXE_EXT): display_distribution$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)display_distribution display_distribution$(OBJ_EXT) $(LDFLAGS) + +off_to_ghs$(EXE_EXT): off_to_ghs$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)off_to_ghs off_to_ghs$(OBJ_EXT) $(LDFLAGS) + +off_to_medit$(EXE_EXT): off_to_medit$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)off_to_medit off_to_medit$(OBJ_EXT) $(LDFLAGS) + +lanteri_output_tet_mesh$(EXE_EXT): lanteri_output_tet_mesh$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)lanteri_output_tet_mesh lanteri_output_tet_mesh$(OBJ_EXT) $(LDFLAGS) + +lanteri$(EXE_EXT): lanteri$(OBJ_EXT) \ + distribution$(OBJ_EXT) \ + Qt_widget_displayer$(OBJ_EXT) \ + Gd_displayer$(OBJ_EXT) \ + lanteri_process_results$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)lanteri \ + lanteri$(OBJ_EXT) \ + distribution$(OBJ_EXT) \ + Qt_widget_displayer$(OBJ_EXT) \ + Gd_displayer$(OBJ_EXT) \ + lanteri_process_results$(OBJ_EXT) \ + $(LDFLAGS) -lgd + +read_mesh$(EXE_EXT): read_mesh$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)read_mesh read_mesh$(OBJ_EXT) $(LDFLAGS) + +slivers_exuder$(OBJ_EXT): ../include/CGAL/Mesh_3/Slivers_exuder.h + +slivers_exuder$(EXE_EXT): slivers_exuder$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)slivers_exuder slivers_exuder$(OBJ_EXT) $(LDFLAGS) + +stat_mesh$(EXE_EXT): stat_mesh$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)stat_mesh stat_mesh$(OBJ_EXT) $(LDFLAGS) + +clean: \ + distribution.clean \ + Gd_displayer.clean \ + cgal_to_medit.clean \ + medit_to_cgal.clean \ + lanteri.clean \ + lanteri_output_tet_mesh.clean \ + lanteri_process_results.clean \ + Qt_widget_displayer.clean \ + read_mesh.clean \ + slivers_exuder.clean \ + stat_mesh.clean \ + display_distribution.clean \ + off_to_ghs.clean \ + off_to_medit.clean + + +#---------------------------------------------------------------------# +# suffix rules +#---------------------------------------------------------------------# + +.C$(OBJ_EXT): + $(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $< diff --git a/Mesh_3/applications/medit_to_cgal.C b/Mesh_3/applications/medit_to_cgal.C new file mode 100644 index 00000000000..a48fcb5af6f --- /dev/null +++ b/Mesh_3/applications/medit_to_cgal.C @@ -0,0 +1,43 @@ +#include "weighted_types.h" + +// ios +#include +#include +#include +#include + +#include "utils.h" + +int main(int , char** argv) +{ + Tr tr; + C2T3 c2t3(tr); + + std::ifstream ifs(argv[1]); + std::ofstream ofs(argv[2]); + if( !ifs || !ofs ) + { + std::cerr << "Usage:\n" + << " " << argv[0] << " INPUT OUPUT\n" << +"\n" +" INPUT must be a medit file.\n" +"\n" +" OUTPUT will be " << format_cgal_description << +"\n"; + return EXIT_FAILURE; + } + + std::cout << " Reading " << argv[1] << std::endl; + + if( CGAL::input_from_medit(ifs, + c2t3, + true, // debug + &std::cout) ) // debug to cout + { + std::cout << " Writing " << argv[2] << std::endl; + CGAL::Mesh_3::output_mesh(ofs, c2t3); + return EXIT_SUCCESS; + } + else + return EXIT_FAILURE; +} diff --git a/Mesh_3/applications/off_to_ghs.C b/Mesh_3/applications/off_to_ghs.C new file mode 100644 index 00000000000..b8a4d9f4ddc --- /dev/null +++ b/Mesh_3/applications/off_to_ghs.C @@ -0,0 +1,71 @@ +#include +#include +#include + +void usage(char* name) +{ + std::cerr << "Usage:\n" + << name << " [FILE_BASENAME]\n" + << " will convert FILE_BASENAME.off to FILE_BASENAME.faces \n" + << " and FILE_BASENAME.points.\n"; + exit(1); +} + +int main(int argc, char** argv) +{ + if(argc < 2) + usage(argv[0]); + + std::string file_base = argv[1]; + + std::ifstream off((file_base + ".off").c_str()); + std::ofstream points((file_base + ".points").c_str()); + std::ofstream faces((file_base + ".faces").c_str()); + + std::string dummy_s; + + int number_of_points; + int number_of_faces; + int dummy_i; + + off >> dummy_s; // OFF + off >> number_of_points >> number_of_faces >> dummy_i; + points << number_of_points << "\n"; + faces << number_of_faces << "\n"; + + while( number_of_points > 0 && ! off.eof() ) + { + std::string x, y, z; + off >> x >> y >> z; + if(off) + { + points << x << " " << y << " " << z << " 0\n"; + --number_of_points; + } + } + if( number_of_points != 0) + { + std::cerr << "Error: bad number of points in OFF file " + << file_base << ".off\n"; + exit(2); + } + points.close(); + while( ! off.eof() ) + { + int i0, i1, i2, i3; + off >> i0 >> i1 >> i2 >> i3; + if(off) + { + faces << i0 << " " << (i1 + 1) << " " + << (i2 + 1) << " " << (i3 + 1) << " 0 0 0 0\n"; + --number_of_faces; + } + } + faces.close(); + if( number_of_faces != 0) + { + std::cerr << "Error: bad number of faces in OFF file " + << file_base << ".off\n"; + exit(2); + } +} diff --git a/Mesh_3/applications/off_to_medit.C b/Mesh_3/applications/off_to_medit.C new file mode 100644 index 00000000000..712509eb229 --- /dev/null +++ b/Mesh_3/applications/off_to_medit.C @@ -0,0 +1,80 @@ +#include +#include +#include + +void usage(char* name) +{ + std::cerr << "Usage:\n" + << name << " [FILE_BASENAME]\n" + << " will convert FILE_BASENAME.off to FILE_BASENAME.mesh.\n"; + exit(1); +} + +int main(int argc, char** argv) +{ + if(argc < 2) + usage(argv[0]); + + std::string file_base = argv[1]; + + std::ifstream off((file_base + ".off").c_str()); + std::ofstream mesh((file_base + ".mesh").c_str()); + + std::string dummy_s; + + int number_of_points; + int number_of_faces; + int dummy_i; + + mesh << "MeshVersionFormatted 1\n" + << "Dimension 3\n" + << "Vertices\n"; + + off >> dummy_s; // OFF + off >> number_of_points >> number_of_faces >> dummy_i; + mesh << number_of_points << "\n"; + + while( number_of_points > 0 && ! off.eof() ) + { + std::string x, y, z; + off >> x >> y >> z; + if(off) + { + mesh << x << " " << y << " " << z << " 1\n"; + --number_of_points; + } + } + + if( number_of_points != 0) + { + std::cerr << "Error: bad number of points in OFF file " + << file_base << ".off\n"; + exit(2); + } + + mesh << "Triangles" << "\n" + << number_of_faces << "\n"; + + while( ! off.eof() ) + { + int i0, i1, i2, i3; + off >> i0 >> i1 >> i2 >> i3; + if(off) + { + mesh << (i1 + 1) << " " + << (i2 + 1) << " " + << (i3 + 1) << " 0\n"; + --number_of_faces; + } + } + + mesh << "Tetrahedra\n0\nEnd\n"; + + mesh.close(); + if( number_of_faces != 0) + { + std::cerr << "Error: bad number of faces in OFF file " + << file_base << ".off\n"; + exit(2); + } +} diff --git a/Mesh_3/applications/read_mesh.C b/Mesh_3/applications/read_mesh.C new file mode 100644 index 00000000000..e6d19b0406b --- /dev/null +++ b/Mesh_3/applications/read_mesh.C @@ -0,0 +1,37 @@ +#include "weighted_types.h" + +// ios +#include +#include +#include + +#include "utils.h" + +int main(int , char** argv) +{ + Tr tr; + C2T3 c2t3(tr); + + std::ifstream ifs(argv[1]); + if( !ifs ) + { + std::cerr << "Usage:\n" + << " " << argv[0] << " FILE\n" + << "\n" + << " FILE must be " << format_cgal_description + << "\n"; + return EXIT_FAILURE; + } + + std::cout << " Reading " << argv[1] << std::endl; + + if( CGAL::Mesh_3::input_mesh(ifs, c2t3, + true, + &std::cerr) ) + { + display_faces_counts(tr, " ", &std::cout); + return EXIT_SUCCESS; + } + else + return EXIT_FAILURE; +} diff --git a/Mesh_3/applications/slivers_exuder.C b/Mesh_3/applications/slivers_exuder.C new file mode 100644 index 00000000000..8c36f853265 --- /dev/null +++ b/Mesh_3/applications/slivers_exuder.C @@ -0,0 +1,58 @@ +#include "weighted_types.h" + +// ios +#include +#include +#include +#include + +#include + +#include "utils.h" + +int main(int argc, char** argv) +{ + Tr tr; + C2T3 c2t3(tr); + + std::ifstream ifs(argv[1]); + std::ofstream ofs(argv[2]); + if( argc != 4 || !ifs || !ofs ) + { + std::cerr << +"Usage:\n" << +" " << argv[0] << " INPUT OUTPUT bound\n" << +"\n" << +" INPUT must be " << format_cgal_description << +"\n" << +" OUPUT is the name of the ouput file. It will be in the same format.\n" << +"\n" +" 0<=bound<=1.0 is a bound on pumping process.\n"; + return EXIT_FAILURE; + } + + std::stringstream str_stream(argv[3]); + + double pumping_bound; + + str_stream >> pumping_bound; + + std::cout << " Reading " << argv[1] << std::endl; + if(! CGAL::Mesh_3::input_mesh(ifs, + c2t3, + true, // debug + &std::cerr) ) // debug to cerr + return EXIT_FAILURE; + ifs.close(); + + CGAL::Mesh_3::Slivers_exuder exuder(tr); + + std::cout << " Pumping" << std::endl; + exuder.init(pumping_bound); + exuder.pump_vertices(pumping_bound); + + std::cout << " Writing " << argv[2] << std::endl; + CGAL::Mesh_3::output_mesh(ofs, c2t3); + + return EXIT_SUCCESS; +} diff --git a/Mesh_3/applications/stat_mesh.C b/Mesh_3/applications/stat_mesh.C new file mode 100644 index 00000000000..1d7ef8f59ff --- /dev/null +++ b/Mesh_3/applications/stat_mesh.C @@ -0,0 +1,51 @@ +#include "weighted_types.h" + +// ios +#include +#include +#include +#include + +#include "utils.h" + +int main(int , char** argv) +{ + Tr tr; + C2T3 c2t3(tr); + + std::ifstream ifs(argv[1]); + if( !ifs ) + { + std::cerr << "Usage:\n" + << " " << argv[0] << " FILE\n" + << "\n" + << " FILE must be " << format_cgal_description + << "\n"; + return EXIT_FAILURE; + } + + std::cout << " Reading " << argv[1] << std::endl; + + if( CGAL::Mesh_3::input_mesh(ifs, c2t3, + true, + &std::cerr) ) +// if( CGAL::input_pslg_from_medit(ifs, +// c2t3, +// true, // debug +// &std::cout) ) // debug to cout + { + display_faces_counts(tr, " ", &std::cout); + + std::cout << "\n Statistics:\n"; + + std::cout << "(vertices)\n"; + display_vertices_by_surface_indices_statistics(tr, " ", &std::cout); + + std:: cout << "(facets)\n"; + display_facets_by_surface_indices_statistics(c2t3, " ", &std::cout); + + return EXIT_SUCCESS; + } + else + return EXIT_FAILURE; +} diff --git a/Mesh_3/applications/types.h b/Mesh_3/applications/types.h new file mode 100644 index 00000000000..197f3a54919 --- /dev/null +++ b/Mesh_3/applications/types.h @@ -0,0 +1,32 @@ +#include + +// vertex and cell bases +#include +#include +#include + +// c2t3 +#include + +// delaunay +#include + +// traits class for reading meshes +#include + +// traits class +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Point_with_surface_index_geom_traits My_traits; + +// vertex and cell types +typedef CGAL::Surface_mesh_vertex_base_3 Vb; +typedef CGAL::Regular_triangulation_cell_base_3 Cb1; +typedef CGAL::Surface_mesh_cell_base_3 Cb2; +typedef CGAL::Volume_mesher_cell_base_3 Cb; + +// triangulation +typedef CGAL::Triangulation_data_structure_3 Tds; +typedef CGAL::Delaunay_triangulation_3 Tr; + +// c2t3 +typedef CGAL::Complex_2_in_triangulation_3 C2T3; diff --git a/Mesh_3/applications/utils.h b/Mesh_3/applications/utils.h new file mode 100644 index 00000000000..bb6a9542325 --- /dev/null +++ b/Mesh_3/applications/utils.h @@ -0,0 +1,153 @@ +#include +#include +#include + +template +void display_faces_counts(const Tr& tr, + std::string prefix = "", + // prefix to each line output + std::ostream* out_stream = &std::cout + // output stream + ) +{ + *out_stream << prefix << "Vertices: " << tr.number_of_vertices() << std::endl + << prefix << "Facets on surface: " + << CGAL::number_of_facets_on_surface(tr) << std::endl + << prefix << "Cells: " + << tr.number_of_cells() << std::endl + << prefix << "Cells in volume: " + << CGAL::number_of_cells_in_domain(tr) << std::endl; +} + +template +void +display_facets_by_surface_indices_statistics(C2T3& c2t3, + std::string prefix = "", + // prefix to each line output + std::ostream* out_stream = + &std::cout + // output stream + ) +{ + typedef typename C2T3::Triangulation Tr; + + const Tr& tr = c2t3.triangulation(); + + std::vector number_of_facets_by_surface_index; + int number_of_hybrid_facets = 0; + + for(typename Tr::Finite_facets_iterator fit = tr.finite_facets_begin(); + fit != tr.finite_facets_end(); + ++fit) + { + const typename Tr::Cell_handle& cell = fit->first; + const int index = fit->second; + + if(c2t3.face_status(cell, index) != C2T3::NOT_IN_COMPLEX) + { + const typename Tr::Vertex_handle& va = cell->vertex((index+1)&3); + const typename Tr::Vertex_handle& vb = cell->vertex((index+1)&3); + const typename Tr::Vertex_handle& vc = cell->vertex((index+1)&3); + + const unsigned int va_index = va->point().surface_index(); + const unsigned int vb_index = vb->point().surface_index(); + const unsigned int vc_index = vc->point().surface_index(); + + if(va_index != vb_index || va_index != vc_index ) + ++number_of_hybrid_facets; + else + { + if(number_of_facets_by_surface_index.size() <= va_index) + number_of_facets_by_surface_index.resize(va_index+1); + ++number_of_facets_by_surface_index[va_index]; + } + } + } + + *out_stream << prefix << "Hybrid facets: " + << number_of_hybrid_facets << std::endl; + + std::vector::size_type vector_size = + number_of_facets_by_surface_index.size(); + for(unsigned int i = 0; i < vector_size; ++i) + { + *out_stream << prefix << "Facets on surface #" << i << ": " + << number_of_facets_by_surface_index[i] << std::endl; + } +} + +template +void +display_vertices_by_surface_indices_statistics(const Tr& tr, + std::string prefix = "", + // prefix to each line output + std::ostream* out_stream = + &std::cout + // output stream + ) +{ + std::vector number_of_vertices_by_surface_index; + + for(typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin(); + vit != tr.finite_vertices_end(); + ++vit) + { + const unsigned int index = vit->point().surface_index(); + + if(number_of_vertices_by_surface_index.size() <= index) + number_of_vertices_by_surface_index.resize(index+1); + ++number_of_vertices_by_surface_index[index]; + } + + std::vector::size_type vector_size = + number_of_vertices_by_surface_index.size(); + for(unsigned int i = 0; i < vector_size; ++i) + { + *out_stream << prefix << "Vertices with index #" << i << ": " + << number_of_vertices_by_surface_index[i] << std::endl; + } +} + +template +void +display_cells_by_volume_indices_statistics(const Tr& tr, + std::string prefix = "", + // prefix to each line output + std::ostream* out_stream = + &std::cout + // output stream + ) +{ + std::vector number_of_cells_by_volume_index; + int number_of_cells_with_default_volume_index = 0; + + for(typename Tr::Finite_cells_iterator cit = tr.finite_cells_begin(); + cit != tr.finite_cells_end(); + ++cit) + if(cit->is_in_domain()) + { + const int index = cit->volume_index(); + + if(index < 0) + ++number_of_cells_with_default_volume_index; + else + { + const unsigned int positive_index = index; + + if(number_of_cells_by_volume_index.size() <= positive_index) + number_of_cells_by_volume_index.resize(index+1); + ++number_of_cells_by_volume_index[index]; + } + } + + *out_stream << prefix << "Cells with default index: " + << number_of_cells_with_default_volume_index << std::endl; + + std::vector::size_type vector_size = + number_of_cells_by_volume_index.size(); + for(unsigned int i = 0; i < vector_size; ++i) + { + *out_stream << prefix << "Cells with index #" << i << ": " + << number_of_cells_by_volume_index[i] << std::endl; + } +} diff --git a/Mesh_3/applications/weighted_types.h b/Mesh_3/applications/weighted_types.h new file mode 100644 index 00000000000..2c8e9b27bd8 --- /dev/null +++ b/Mesh_3/applications/weighted_types.h @@ -0,0 +1,44 @@ +#include + +// regular +#include +#include +#include + +// IO.h must be included before vertex and cell bases. +#include + +// vertex and cell bases +#include +#include +#include + +// c2t3 +#include + +// traits class for reading meshes +#include + +// traits class +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Regular_triangulation_euclidean_traits_3 Regular_traits; +typedef CGAL::Weighted_point_with_surface_index_geom_traits My_traits; + +// vertex and cell types +typedef CGAL::Surface_mesh_vertex_base_3 Vb; +typedef CGAL::Regular_triangulation_cell_base_3 Cb1; +typedef CGAL::Surface_mesh_cell_base_3 Cb2; +typedef CGAL::Volume_mesher_cell_base_3 Cb; + +// triangulation +typedef CGAL::Triangulation_data_structure_3 Tds; +typedef CGAL::Regular_triangulation_3 Tr; + +// c2t3 +typedef CGAL::Complex_2_in_triangulation_3 C2T3; + +const std::string format_cgal_description = + "a file of format produced by the output operator of\n" +" CGAL::Triangulation_3, with points\n" +" CGAL::Weighted_point_with_surface_index and cells\n" +" CGAL::Volume_mesher_cell_base_3.\n";