Starting to implement the interface to the delaunay mesher

This commit is contained in:
Nico Kruithof 2006-07-27 13:35:44 +00:00
parent 970a0790ad
commit a52f20ae99
15 changed files with 720 additions and 135 deletions

1
.gitattributes vendored
View File

@ -1476,6 +1476,7 @@ Skin_surface_3/examples/Skin_surface_3/data/1IYE.pdb -text
Skin_surface_3/examples/Skin_surface_3/data/molecule_tunnel.cin -text
Skin_surface_3/examples/Skin_surface_3/dsrpdb/doc/doc_footer -text
Skin_surface_3/examples/Skin_surface_3/dsrpdb/doc/main_doc -text
Skin_surface_3/test/Skin_surface_3/data/inexact_tmc_fails.cin -text
Skin_surface_3/test/Skin_surface_3/data/test1.cin -text
Skin_surface_3/test/Skin_surface_3/data/test10.cin -text
Skin_surface_3/test/Skin_surface_3/data/test11.cin -text

View File

@ -1,2 +1,3 @@
TODO.txt
Doxyfile
test/Skin_surface_3/test_mixed_complex_point_location.cpp

View File

@ -44,7 +44,8 @@ all: \
skin_surface_simple$(EXE_EXT) \
skin_surface_subdiv$(EXE_EXT) \
skin_surface_subdiv_with_normals$(EXE_EXT) \
union_of_balls_simple$(EXE_EXT)
union_of_balls_simple$(EXE_EXT) \
skin_surface_with_surface_mesher$(EXE_EXT)
NGHK_skin_surface_simple$(EXE_EXT): NGHK_skin_surface_simple$(OBJ_EXT)
$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)NGHK_skin_surface_simple NGHK_skin_surface_simple$(OBJ_EXT) $(LDFLAGS)
@ -73,12 +74,16 @@ skin_surface_subdiv_with_normals$(EXE_EXT): skin_surface_subdiv_with_normals$(OB
union_of_balls_simple$(EXE_EXT): union_of_balls_simple$(OBJ_EXT)
$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)union_of_balls_simple union_of_balls_simple$(OBJ_EXT) $(LDFLAGS)
skin_surface_with_surface_mesher$(EXE_EXT): skin_surface_with_surface_mesher$(OBJ_EXT)
$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)skin_surface_with_surface_mesher skin_surface_with_surface_mesher$(OBJ_EXT) $(LDFLAGS)
clean: \
skin_surface_pdb_reader.clean \
skin_surface_simple.clean \
skin_surface_subdiv.clean \
skin_surface_subdiv_with_normals.clean \
union_of_balls_simple.clean \
skin_surface_with_surface_mesher.clean \
NGHK_skin_surface_simple.clean \
NGHK_skin_surface_subdiv.clean
@ -88,4 +93,3 @@ clean: \
.cpp$(OBJ_EXT):
$(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $<

View File

@ -0,0 +1,58 @@
#define CGAL_SURFACE_MESHER_VERBOSE 1
// file examples/Skin_surface_3/skin_surface_with_surface_mesher.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Skin_surface_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
#include <CGAL/IO/Complex_2_in_triangulation_3_file_writer.h>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Regular_triangulation_filtered_traits_3<Kernel> Traits;
typedef CGAL::Skin_surface_3<Traits> Skin_surface_3;
typedef Skin_surface_3::RT RT;
typedef Skin_surface_3::Weighted_point Weighted_point;
typedef Weighted_point::Point Bare_point;
typedef CGAL::Surface_mesh_vertex_base_3<Kernel> Vb;
typedef CGAL::Surface_mesh_cell_base_3<Kernel> Cb;
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> Tr;
typedef CGAL::Surface_mesh_complex_2_in_triangulation_3<Tr> C2t3;
typedef Kernel::Point_3 Point_3;
typedef Kernel::FT FT;
int main(int, char **) {
// Construct the Skin_surface_3 object:
std::list<Weighted_point> l;
RT shrinkfactor = 0.5;
l.push_front(Weighted_point(Bare_point(0,0,0), 1));
l.push_front(Weighted_point(Bare_point(0,1,0), 2));
l.push_front(Weighted_point(Bare_point(0,0,2), 1));
Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor);
// Construct the types needed for the surface mesher:
Tr tr; // 3D-Delaunay triangulation
C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation
// defining meshing criteria
CGAL::Surface_mesh_default_criteria_3<Tr> criteria(30., // angular bound
0.1, // radius bound
0.1); // distance bound
// meshing surface
make_surface_mesh(c2t3, skin_surface, criteria, CGAL::Non_manifold_tag());
std::ofstream out("delaunay_mesh.off");
CGAL::output_surface_facets_to_off(out, c2t3);
std::cout << "Final number of points: " << tr.number_of_vertices() << "\n";
}

View File

@ -0,0 +1,94 @@
// Copyright (c) 2005 Rijksuniversiteit Groningen (Netherlands)
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Nico Kruithof <Nico@cs.rug.nl>
#ifndef CGAL_MIXED_COMPLEX_TRAITS_3_H
#define CGAL_MIXED_COMPLEX_TRAITS_3_H
#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
// Contains the weighted converter:
#include <CGAL/Regular_triangulation_filtered_traits_3.h>
CGAL_BEGIN_NAMESPACE
/** Input: a list of n weighted points p_1...p_n and a query point x.
There is a plane separating the mixed cell defined by p_1...p_n-1
and the mixed cell defined by p_1...p_n. The predicate tests
whether x lies on the same side of this plane as the mixed cell
defined by p_1...p_n-1 (NEGATIVE), on the plane (ZERO) or on the
opposite side (POSITIVE).
**/
template <class K>
class Side_of_mixed_cell {
public:
typedef typename K::FT FT;
typedef typename K::Bare_point Bare_point;
typedef typename K::Weighted_point Weighted_point;
Side_of_mixed_cell(const FT &shrink) : s(shrink) {}
typedef CGAL::Arity_tag< 5 > Arity;
typedef CGAL::Sign result_type;
Test_add() : _t(1) {}
Test_add(FT t) : _t(t) {}
result_type operator()(const Weighted_point &p1,
const Weighted_point &p2,
const Bare_point &x) const {
return side_of_mixed_cellC3(p1.x(),p1.y(),p1.z(),p1.weight(),
p2.x(),p2.y(),p2.z(),p2.weight(),
x.x(),x.y(),x.z(),
s);
}
result_type operator()(const Weighted_point &p1,
const Weighted_point &p2,
const Weighted_point &p3,
const Bare_point &x) const {
return side_of_mixed_cellC3(p1.x(),p1.y(),p1.z(),p1.weight(),
p2.x(),p2.y(),p2.z(),p2.weight(),
p3.x(),p3.y(),p3.z(),p3.weight(),
x.x(),x.y(),x.z(),
s);
}
result_type operator()(const Weighted_point &p1,
const Weighted_point &p2,
const Weighted_point &p3,
const Weighted_point &p4,
const Bare_point &x) const {
return side_of_mixed_cellC3(p1.x(),p1.y(),p1.z(),p1.weight(),
p2.x(),p2.y(),p2.z(),p2.weight(),
p3.x(),p3.y(),p3.z(),p3.weight(),
p4.x(),p4.y(),p4.z(),p4.weight(),
x.x(),x.y(),x.z(),
s);
}
private:
FT s;
}
template <class K>
class Mixed_complex_traits_3_base {
};
CGAL_END_NAMESPACE
#endif // CGAL_MIXED_COMPLEX_TRAITS_3_H

View File

@ -31,10 +31,14 @@
#include <CGAL/triangulate_mixed_complex_3.h>
// Needed for the (Delaunay) surface mesher
#include <CGAL/Skin_surface_mesher_oracle_3.h>
#include <CGAL/Triangulation_simplex_3.h>
CGAL_BEGIN_NAMESPACE
template < class GT,
class QuadraticSurface,
class SkinSurface_3,
class Cb = Triangulation_cell_base_3<GT> >
class Triangulated_mixed_complex_cell_3 : public Cb
{
@ -43,12 +47,13 @@ public:
typedef typename Triangulation_data_structure::Vertex_handle Vertex_handle;
typedef typename Triangulation_data_structure::Cell_handle Cell_handle;
typedef QuadraticSurface Quadratic_surface;
typedef typename SkinSurface_3::Quadratic_surface Quadratic_surface;
typedef typename SkinSurface_3::Simplex Simplex;
template < class TDS2 >
struct Rebind_TDS {
typedef typename Cb::template Rebind_TDS<TDS2>::Other Cb2;
typedef Triangulated_mixed_complex_cell_3<GT, Quadratic_surface, Cb2>
typedef Triangulated_mixed_complex_cell_3<GT, SkinSurface_3, Cb2>
Other;
};
@ -64,6 +69,7 @@ public:
// return surf->sign(p);
// }
Quadratic_surface *surf;
Simplex simp;
};
template < class GT,
@ -92,35 +98,41 @@ public:
template <class SkinSurfaceTraits_3>
class Skin_surface_3 {
typedef SkinSurfaceTraits_3 Gt;
typedef Skin_surface_3<Gt> Self;
public:
typedef SkinSurfaceTraits_3 Geometric_traits;
typedef typename Gt::Weighted_point Weighted_point;
typedef typename Weighted_point::Weight RT;
// NGHK:: added for the Delaunay mesher
typedef typename Gt::Sphere_3 Sphere_3;
private:
typedef typename Weighted_point::Point Bare_point;
typedef Regular_triangulation_3<Gt> Regular;
public:
// NGHK: remove later?
typedef Triangulation_simplex_3<Regular> Simplex;
// defining the triangulated mixed complex:
typedef Exact_predicates_exact_constructions_kernel TMC_traits;
#ifdef CGAL_SKIN_SURFACE_USE_EXACT_IMPLICIT_SURFACE
typedef Skin_surface_quadratic_surface_3<TMC_traits> Quadratic_surface;
#else
typedef Skin_surface_quadratic_surface_3<Simple_cartesian<double> >
Quadratic_surface;
#endif // CGAL_SKIN_SURFACE_USE_EXACT_IMPLICIT_SURFACE
typedef Triangulation_3<
TMC_traits,
Triangulation_data_structure_3
< Triangulated_mixed_complex_vertex_3<TMC_traits>,
Triangulated_mixed_complex_cell_3<TMC_traits,Quadratic_surface> >
Triangulated_mixed_complex_cell_3<TMC_traits,Self> >
> Triangulated_mixed_complex;
typedef typename Triangulated_mixed_complex::Vertex_handle TMC_Vertex_handle;
typedef typename Triangulated_mixed_complex::Cell_handle TMC_Cell_handle;
// NGHK: added for the (Delaunay) surface mesher, document
typedef Exact_predicates_inexact_constructions_kernel Mesher_Gt;
typedef Skin_surface_mesher_oracle_3<Mesher_Gt,Self> Surface_mesher_traits_3;
private:
typedef typename TMC_traits::Point_3 TMC_Point;
@ -153,7 +165,7 @@ public:
}
// Construct the triangulated mixed complex:
triangulate_mixed_complex_3(regular, shrink, _tmc,verbose);
triangulate_mixed_complex_3(regular, shrink, _tmc, verbose);
CGAL_assertion(_tmc.is_valid());
if (verbose) {
@ -174,15 +186,45 @@ public:
}
TMC_Cell_handle locate(const TMC_Point &p) const{
return _tmc.locate(p);
last_ch = _tmc.locate(p, last_ch);
return last_ch;
}
// NGHK: added for the (Delaunay) surface mesher, document
Sphere_3 bounding_sphere() const {
return _bounding_sphere;
}
RT squared_error_bound() const {
return .01;
}
Sign operator()(const Bare_point &p) const {
Cartesian_converter<typename Bare_point::R, TMC_traits > converter;
TMC_Point p_tmc = converter(p);
TMC_Cell_handle ch = locate(p_tmc);
if (_tmc.is_infinite(ch)) {
// Infinite cells do not have a pointer to a surface
return NEGATIVE;
}
return ch->surf->sign(p_tmc);
}
typename Mesher_Gt::FT
get_density(const typename Mesher_Gt::Point_3 &p) const {
// NGHK: Make adaptive
return 1;
}
private:
// Used to optimize the point location in TMC:
mutable TMC_Cell_handle last_ch;
void construct_bounding_box(Regular &regular);
Gt &gt;
RT shrink;
Triangulated_mixed_complex _tmc;
bool verbose;
Sphere_3 _bounding_sphere;
};
template <class SkinSurfaceTraits_3>
@ -226,6 +268,9 @@ construct_bounding_box(Regular &regular)
Bare_point(mid.x(),mid.y(),bbox.zmax()+(dx+dy+dr)/shrink),-1));
regular.insert(Weighted_point(
Bare_point(mid.x(),mid.y(),bbox.zmin()-(dx+dy+dr)/shrink),-1));
// Set the bounding sphere for the Delaunay mesher
_bounding_sphere = Sphere_3(mid, dr*dr+1);
}
}

View File

@ -0,0 +1,329 @@
// Copyright (c) 2003-2006 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Nico Kruithof
#ifndef CGAL_SKIN_SURFACE_MESHER_ORACLE_3_H
#define CGAL_SKIN_SURFACE_MESHER_ORACLE_3_H
#include <CGAL/point_generators_3.h>
//#include <CGAL/Surface_mesher/Sphere_oracle_3.h>
#include <CGAL/Surface_mesher/Oracles/Sphere_oracle_3.h>
#include <queue>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/mesh_skin_surface_3.h>
#ifdef CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
#include <string>
#endif
#include <sstream>
#include <iostream>
// NGHK: remove later, for debugging
#include <CGAL/IO/Polyhedron_iostream.h>
CGAL_BEGIN_NAMESPACE
template <
class GT,
class Surface,
class Point_creator = Creator_uniform_3<typename GT::FT,
typename GT::Point_3>
>
class Skin_surface_mesher_oracle_3
{
// private types
typedef Skin_surface_mesher_oracle_3<GT, Surface, Point_creator> Self;
typedef typename Surface_mesher::Sphere_oracle_3<GT, Point_creator> Sphere_oracle;
typedef typename GT::Point_3 Point;
typedef typename GT::FT FT;
typedef typename GT::Sphere_3 Sphere_3;
public:
// Public types
typedef GT Geom_traits;
typedef typename GT::Point_3 Point_3;
typedef typename GT::Segment_3 Segment_3;
typedef typename GT::Ray_3 Ray_3;
typedef typename GT::Line_3 Line_3;
typedef Surface Surface_3;
public:
// Constructors
Skin_surface_mesher_oracle_3 ()
{
#ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS
std::cerr << "CONS: Skin_surface_mesher_oracle_3\n";
#endif
}
class Intersect_3
{
public:
Intersect_3()
{
}
Object operator()(const Surface_3& surface, Segment_3 s) const
// s is passed by value, because it is clipped below
{
typename GT::Construct_point_on_3 point_on =
GT().construct_point_on_3_object();
typename Sphere_oracle::Intersect_3 clip =
Sphere_oracle().intersect_3_object();
const Sphere_3& sphere = surface.bounding_sphere();
Point_3 a = point_on(s, 0);
Point_3 b = point_on(s, 1);
// if both extremities are on the same side of the surface, return
// no intersection
if(surf_equation(surface, a) * surf_equation(surface, b) > 0)
return Object();
// Code for surfaces with boundaries
// First rescale segment s = [a, b]
if( clip.clip_segment(sphere, a, b) )
return intersect_clipped_segment(surface,
a,
b,
surface.squared_error_bound());
// else
return Object();
} // end operator()(Surface_3, Segment_3)
Object operator()(const Surface_3& surface, const Ray_3& r) const {
typename Sphere_oracle::Intersect_3 clip =
Sphere_oracle().intersect_3_object();
const Sphere_3& sphere = surface.bounding_sphere();
Point_3 a, b;
if(clip.clip_ray(sphere, r, a, b))
{
return intersect_clipped_segment(surface,
a,
b,
surface.squared_error_bound());
}
// else
return Object();
} // end operator()(Surface_3, Ray_3)
Object operator()(const Surface_3& surface, const Line_3& l) const {
typename Sphere_oracle::Intersect_3 clip =
Sphere_oracle().intersect_3_object();
const Sphere_3& sphere = surface.bounding_sphere();
Point_3 a, b;
if(clip.clip_line(sphere, l, a, b))
{
return intersect_clipped_segment(surface,
a,
b,
surface.squared_error_bound());
}
else
return Object();
}; // end operator()(Surface_3, Line_3)
// debug function
static std::string debug_point(const Surface_3& surface,
const Point& p)
{
std::stringstream s;
s << p << " (distance="
<< CGAL::sqrt(CGAL::squared_distance(p,
surface.bounding_sphere().center()))
<< ", sign=" << surf_equation(surface, p)
<< ")";
return s.str();
}
static CGAL::Sign surf_equation (Surface_3 surface,
const Point& p)
{
return CGAL::sign(surface(p));
} // @TODO, @WARNING: we use x(), y() and z()
private:
// Private functions
Object intersect_clipped_segment(const Surface_3& surface,
Point p1,
Point p2,
const FT& squared_distance_bound) const
{
#ifdef CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
std::cerr << "clipped_segment\n";
#endif
typename GT::Compute_squared_distance_3 squared_distance =
GT().compute_squared_distance_3_object();
typename GT::Construct_midpoint_3 midpoint =
GT().construct_midpoint_3_object();
Sign sign_at_p1 = surf_equation(surface, p1);
Sign sign_at_p2 = surf_equation(surface, p2);
if( sign_at_p1 == ZERO )
{
return make_object(p1);
}
if( sign_at_p2 == ZERO )
{
return make_object(p2);
}
// if both extremities are on the same side of the surface, return
// no intersection
if(sign_at_p1 * sign_at_p2 > 0)
return Object();
while(true)
{
#ifdef CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
std::cerr << debug_point(surface, p1) << ", "
<< debug_point(surface, p2) << "\n";
#endif
Point mid = midpoint(p1, p2);
const Sign sign_at_mid = surf_equation(surface, mid);
if ( sign_at_mid == ZERO ||
squared_distance(p1, p2) < squared_distance_bound )
// If the two points are close, then we must decide
{
#ifdef CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
std::cerr << "=" << debug_point(surface, mid) << "\n";
#endif
return make_object(mid);
}
// Else we must go on
if ( sign_at_p1 * sign_at_mid < 0 )
{
p2 = mid;
sign_at_p2 = sign_at_mid;
}
else
{
p1 = mid;
sign_at_p1 = sign_at_mid;
}
}
} // end intersect_clipped_segment
}; // end nested class Intersect_3
class Construct_initial_points
{
Self& oracle;
public:
Construct_initial_points(Self& oracle) : oracle(oracle)
{
}
// Random points
template <typename OutputIteratorPoints>
OutputIteratorPoints operator() (const Surface_3& surface,
OutputIteratorPoints out,
int n = 20) // WARNING: why 20?
{
typedef CGAL::Polyhedron_3<GT> Polyhedron;
typedef typename Polyhedron::Vertex_iterator V_iterator;
typedef typename Polyhedron::Halfedge_around_vertex_circulator
HAV_circulator;
Polyhedron coarse_mesh;
mesh_skin_surface_3(surface, coarse_mesh);
for (V_iterator vit = coarse_mesh.vertices_begin();
vit != coarse_mesh.vertices_end(); ) {
HAV_circulator hav_start, hav_it;
hav_start = hav_it = vit->vertex_begin();
bool ok = true;
FT dist = surface.get_density(vit->point());
do {
ok = (squared_distance(vit->point(),
hav_it->vertex()->point()) < dist);
} while (ok && (hav_start != (++hav_it)));
vit++;
if (!ok) {
std::cout << "remove Vertex" << std::endl;
V_iterator tmp = vit; tmp--;
CGAL_assertion(squared_distance(tmp->point(),
hav_it->vertex()->point()) < dist);
if (tmp->degree() == 3) {
CGAL_assertion(hav_it->opposite()->vertex() == tmp);
coarse_mesh.erase_center_vertex(hav_it->opposite());
continue;
}
if (hav_it->vertex()->degree() == 3) {
coarse_mesh.erase_center_vertex(hav_it);
continue;
}
coarse_mesh.join_vertex(hav_it);
}
}
std::copy(coarse_mesh.points_begin(), coarse_mesh.points_end(), out);
{
std::ofstream os("delaunay_coarsed_mesh.off");
os << coarse_mesh;
}
return out;
}
}; // end nested class Construct_initial_points
Construct_initial_points construct_initial_points_object()
{
return Construct_initial_points(*this);
}
Intersect_3 intersect_3_object()
{
return Intersect_3();
}
bool is_in_volume(const Surface_3& surface, const Point& p)
{
return Intersect_3::surf_equation(surface, p)<0.;
}
}; // end Skin_surface_mesher_oracle_3
template <typename FT>
FT approximate_sqrt(const FT x) {
return FT (CGAL_NTS sqrt(CGAL_NTS to_double(x)));
}
CGAL_END_NAMESPACE
#endif // CGAL_SKIN_SURFACE_MESHER_ORACLE_3_H

View File

@ -41,12 +41,14 @@ public:
RT value(Input_point const &x) const {
typedef Cartesian_converter<typename Input_point::R, K> Converter;
Vector v = Converter()(x) - p;
RT vx = Converter()(x.x()) - p.x();
RT vy = Converter()(x.y()) - p.y();
RT vz = Converter()(x.z()) - p.z();
return
v.x()*(Q[0]*v.x()) +
v.y()*(Q[1]*v.x()+Q[2]*v.y()) +
v.z()*(Q[3]*v.x()+Q[4]*v.y()+Q[5]*v.z()) +
vx*(Q[0]*vx) +
vy*(Q[1]*vx+Q[2]*vy) +
vz*(Q[3]*vx+Q[4]*vy+Q[5]*vz) +
c;
}
template <class Input_point>

View File

@ -84,7 +84,7 @@ public:
Rt_Edge e;
Rt_Facet f;
Rt_Cell_handle ch;
switch (s.dimension()) {
case 0:
{
@ -156,6 +156,7 @@ public:
}
}
ch->surf = surf;
ch->simplex = s;
}
Surface_RT shrink;

View File

@ -0,0 +1,97 @@
// Copyright (c) 1997 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Nico Kruithof
#ifndef CGAL_PREDICATES_FOR_MIXED_COMPLEX_3_H
#define CGAL_PREDICATES_FOR_MIXED_COMPLEX_3_H
#include <CGAL/determinant.h>
#include <CGAL/enum.h>
CGAL_BEGIN_NAMESPACE
/* q: circumcenter(p1,p2)
* <p1-x,p1-q> - s*<p1-q,p1-q>
* = <p1-x-s(p1-q),p1-q>
* = <(1-s)p1+s*q-x,p1-q>
*/
template < class FT>
Sign
side_of_mixed_cellC3(const FT &p1x, const FT &p1y, const FT &p1z, const FT &p1w,
const FT &p2x, const FT &p2y, const FT &p2z, const FT &p2w,
const FT &xx, const FT &xy, const FT &xz,
const FT &s) {
FT x, y, z, t;
t = 1-s;
weighted_circumcenterC3(p1x, p1y, p1z, p1w,
p2x, p2y, p2z, p2w,
x,y,z);
return CGAL_NTS sign((t*p1x+s*x-xx)*(p1x-x) +
(t*p1y+s*y-xy)*(p1y-y) +
(t*p1z+s*z-xz)*(p1z-z));
}
template < class FT>
Sign
side_of_mixed_cellC3(const FT &p1x, const FT &p1y, const FT &p1z, const FT &p1w,
const FT &p2x, const FT &p2y, const FT &p2z, const FT &p2w,
const FT &p3x, const FT &p3y, const FT &p3z, const FT &p3w,
const FT &xx, const FT &xy, const FT &xz,
const FT &s) {
FT x1, y1, z1, x2, y2, z2, t;
t = 1-s;
weighted_circumcenterC3(p1x, p1y, p1z, p1w,
p2x, p2y, p2z, p2w,
x1,y1,z1);
weighted_circumcenterC3(p1x, p1y, p1z, p1w,
p2x, p2y, p2z, p2w,
p3x, p3y, p3z, p3w,
x2,y2,z2);
return CGAL_NTS sign((t*x1+s*x2-xx)*(x1-x2) +
(t*y1+s*y2-xy)*(y1-y2) +
(t*z1+s*z2-xz)*(z1-z2));
}
template < class FT>
Sign
side_of_mixed_cellC3(const FT &p1x, const FT &p1y, const FT &p1z, const FT &p1w,
const FT &p2x, const FT &p2y, const FT &p2z, const FT &p2w,
const FT &p3x, const FT &p3y, const FT &p3z, const FT &p3w,
const FT &p4x, const FT &p4y, const FT &p4z, const FT &p4w,
const FT &xx, const FT &xy, const FT &xz,
const FT &s) {
FT x1, y1, z1, x2, y2, z2, t;
t = 1-s;
weighted_circumcenterC3(p1x, p1y, p1z, p1w,
p2x, p2y, p2z, p2w,
p3x, p3y, p3z, p3w,
x1,y1,z1);
weighted_circumcenterC3(p1x, p1y, p1z, p1w,
p2x, p2y, p2z, p2w,
p3x, p3y, p3z, p3w,
p4x, p4y, p4z, p4w,
x2,y2,z2);
return CGAL_NTS sign((t*x1+s*x2-xx)*(x1-x2) +
(t*y1+s*y2-xy)*(y1-y2) +
(t*z1+s*z2-xz)*(z1-z2));
}
CGAL_END_NAMESPACE
#endif // CGAL_PREDICATES_FOR_MIXED_COMPLEX_3_H

View File

@ -48,7 +48,7 @@ public:
typedef RegularTriangulation_3 Regular;
typedef TriangulatedMixedComplex_3 Triangulated_mixed_complex;
typedef TriangulatedMixedComplexObserver_3
Triangulated_mixed_complex_observer;
Triangulated_mixed_complex_observer;
private:
typedef typename Regular::Vertex_handle Rt_Vertex_handle;
typedef typename Regular::Edge Rt_Edge;

View File

@ -0,0 +1,4 @@
260.94 149.49 119.62 0
258.23 146.88 117.85 0
264.54 146.22 118.93 0
262.02 143.15 118.82 0

View File

@ -15,110 +15,51 @@ typedef Skin_surface_3::Weighted_point Weighted_point;
typedef Weighted_point::Point Bare_point;
typedef CGAL::Polyhedron_3<K> Polyhedron;
bool test(char * filename, double shrink) {
std::list<Weighted_point> l;
std::ifstream in(filename);
CGAL_assertion(in.is_open());
Weighted_point wp;
while (in >> wp) l.push_front(wp);
class Test_file {
public:
Test_file(double shrink) : s(shrink) {
}
void operator()(char * filename) {
std::cout << filename << std::endl;
Skin_surface_3 skin_surface(l.begin(), l.end(), shrink);
std::list<Weighted_point> l;
std::ifstream in(filename);
CGAL_assertion(in.is_open());
Weighted_point wp;
while (in >> wp) l.push_front(wp);
Skin_surface_3 skin_surface(l.begin(), l.end(), s);
Polyhedron p;
CGAL::mesh_skin_surface_3(skin_surface, p);
CGAL_assertion(p.is_valid() && p.is_closed());
}
private:
double s;
};
Polyhedron p;
CGAL::mesh_skin_surface_3(skin_surface, p);
return (p.is_valid() && p.is_closed());
}
int main(int argc, char *argv[]) {
bool result;
char *filename;
filename = "data/caffeine.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/ball.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/degenerate.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test1.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test2.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test3.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test4.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test5.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test6.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test7.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test8.cin";
CGAL_assertion(result);
result = test(filename, .5);
filename = "data/test9.cin";
CGAL_assertion(result);
result = test(filename, .5);
filename = "data/test10.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/test11.cin";
result = test(filename, .5);
CGAL_assertion(result);
filename = "data/caffeine.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/ball.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/degenerate.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test1.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test2.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test3.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test4.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test5.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test6.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test7.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test8.cin";
CGAL_assertion(result);
result = test(filename, .85);
filename = "data/test9.cin";
CGAL_assertion(result);
result = test(filename, .85);
filename = "data/test10.cin";
result = test(filename, .85);
CGAL_assertion(result);
filename = "data/test11.cin";
result = test(filename, .85);
CGAL_assertion(result);
std::vector<char *> filenames;
filenames.push_back("data/caffeine.cin");
filenames.push_back("data/ball.cin");
filenames.push_back("data/degenerate.cin");
filenames.push_back("data/test1.cin");
filenames.push_back("data/test2.cin");
filenames.push_back("data/test3.cin");
filenames.push_back("data/test4.cin");
filenames.push_back("data/test5.cin");
filenames.push_back("data/test6.cin");
filenames.push_back("data/test7.cin");
filenames.push_back("data/test8.cin");
filenames.push_back("data/test9.cin");
filenames.push_back("data/test10.cin");
filenames.push_back("data/test11.cin");
std::for_each(filenames.begin(), filenames.end(), Test_file(.5));
std::for_each(filenames.begin(), filenames.end(), Test_file(.85));
return 0;
}

View File

@ -32,18 +32,38 @@ Skin_surface_3 create_skin_surface(char * filename, double shrink) {
}
template < class Skin_surface_3, class Polyhedron>
bool construct_and_subdivide_mesh(Skin_surface_3 &skin_surface,
void construct_and_subdivide_mesh(Skin_surface_3 &skin_surface,
Polyhedron &polyhedron)
{
CGAL::mesh_skin_surface_3(skin_surface, polyhedron);
CGAL_assertion(polyhedron.is_valid() && polyhedron.is_closed());
CGAL::subdivide_skin_surface_mesh_3(polyhedron, skin_surface);
return (polyhedron.is_valid() && polyhedron.is_closed());
CGAL_assertion(polyhedron.is_valid() && polyhedron.is_closed());
}
class Test_file {
public:
Test_file(double shrink) : s(shrink) {
}
void operator()(char * filename) {
std::cout << filename << std::endl;
Skin_surface_3 skin_surface = create_skin_surface(filename, .5);
Polyhedron p;
construct_and_subdivide_mesh(skin_surface, p);
p.clear();
Polyhedron_skin p_skin;
construct_and_subdivide_mesh(skin_surface, p_skin);
p_skin.clear();
}
private:
double s;
};
int main(int argc, char *argv[]) {
bool result;
std::vector<char *> filenames;
filenames.push_back("data/caffeine.cin");
@ -61,20 +81,8 @@ int main(int argc, char *argv[]) {
filenames.push_back("data/test10.cin");
filenames.push_back("data/test11.cin");
Polyhedron p;
Polyhedron_skin p_skin;
for (std::vector<char *>::size_type i=0; i<filenames.size(); i++) {
Skin_surface_3 skin_surface = create_skin_surface(filenames[i], .5);
Polyhedron p;
result = construct_and_subdivide_mesh(skin_surface, p);
CGAL_assertion(result);
p.clear();
result = construct_and_subdivide_mesh(skin_surface, p_skin);
CGAL_assertion(result);
p_skin.clear();
}
for_each(filenames.begin(), filenames.end(), Test_file(.5));
for_each(filenames.begin(), filenames.end(), Test_file(.25));
return 0;
}