songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Optimizing the new implicit mixed complex code 

next: do the same for the subdivision
This commit is contained in:
Nico Kruithof 2006-08-17 09:31:06 +00:00
parent beedd15017
commit 9fb81a174d
8 changed files with 106 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_simple.cpp
View File

@ -1,21 +1,25 @@
// examples/Skin_surface_3/NGHK_skin_surface_simple.C
//#define CGAL_PROFILE
// examples/Skin_surface_3/skin_surface_simple.C
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Skin_surface_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/mesh_skin_surface_3.h>
#include <CGAL/subdivide_skin_surface_mesh_3.h>
#include <CGAL/Skin_surface_polyhedral_items_3.h>
#include <list>
#include <fstream>
#include "skin_surface_writer.h"
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_euclidean_traits_3<K> Traits;
typedef CGAL::Mixed_complex_traits_3<K> 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::Skin_surface_polyhedral_items_3<Skin_surface_3> Poly_items;
typedef CGAL::Polyhedron_3<K,Poly_items> Polyhedron;
typedef CGAL::Polyhedron_3<
CGAL::Simple_cartesian<double>,
CGAL::Skin_surface_polyhedral_items_3<Skin_surface_3> > Polyhedron;
int main(int argc, char *argv[]) {
if (argc < 2) {
@ -28,14 +32,14 @@ int main(int argc, char *argv[]) {
Weighted_point wp;
std::ifstream in(argv[1]);
while (in >> wp) l.push_front(wp);
while (in >> wp) l.push_back(wp);
Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor, false);
Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor);
Polyhedron p;
CGAL::mesh_skin_surface_3(skin_surface, p);
CGAL::subdivide_skin_surface_mesh_3(p, skin_surface,1);
// CGAL::subdivide_skin_surface_mesh_3(p, skin_surface,1);
std::cout << "Is closed: " << (p.is_closed() ? "Yes" : "No") << std::endl;

1
Skin_surface_3/examples/Skin_surface_3/makefile
View File

@ -86,6 +86,7 @@ clean: \
skin_surface_with_surface_mesher.clean \
NGHK_skin_surface_simple.clean \
NGHK_skin_surface_subdiv.clean
rm -f $(OBJ)
#---------------------------------------------------------------------#
# suffix rules

2
Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.cpp
View File

@ -1,5 +1,3 @@
#define CGAL_PROFILE
// examples/Skin_surface_3/skin_surface_simple.C
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Skin_surface_3.h>

2
Skin_surface_3/include/CGAL/Marching_tetrahedra_traits_skin_surface_3.h
View File

@ -49,7 +49,7 @@ public:
HDS_point intersection(Cell_iterator const ch, int i, int j) const {
// Precondition: ch is not an infinite cell: their surface is not set
Skin_point p;
ss_3.intersect(ch->vertex(i), ch->vertex(j), p);
ss_3.intersect(ch->vertex(i), ch->vertex(j), ch->mixed_cell(), p);
return
Cartesian_converter<typename Skin_point::R, typename HDS_point::R>()(p);

76
Skin_surface_3/include/CGAL/Skin_surface_3.h
View File

@ -353,19 +353,29 @@ public:
return construct_surface(vh->first, K()).sign(p);
}
// Trivial caching: check wether the surface is the same as the previous:
mutable Skin_surface_quadratic_surface_3<
Simple_cartesian<Interval_nt_advanced> > previous_sign_surface;
mutable Simplex previous_sign_simplex;
Sign sign(const Simplex &sim, const Bare_point &p) const {
if (previous_sign_simplex != sim) {
previous_sign_simplex = sim;
previous_sign_surface =
construct_surface(sim,
Simple_cartesian<Interval_nt_advanced>());
}
try
{
CGAL_PROFILER(std::string("NGHK: calls to : ") + std::string(CGAL_PRETTY_FUNCTION));
CGAL_PROFILER(std::string("NGHK: calls to : ") +
std::string(CGAL_PRETTY_FUNCTION));
Protect_FPU_rounding<true> P;
Sign result = construct_surface
(sim,
Exact_predicates_inexact_constructions_kernel()).sign(p);
Sign result = previous_sign_surface.sign(p);
if (! is_indeterminate(result))
return result;
}
catch (Interval_nt_advanced::unsafe_comparison) {}
CGAL_PROFILER(std::string("NGHK: failures of : ") + std::string(CGAL_PRETTY_FUNCTION));
CGAL_PROFILER(std::string("NGHK: failures of : ") +
std::string(CGAL_PRETTY_FUNCTION));
Protect_FPU_rounding<false> P(CGAL_FE_TONEAREST);
return construct_surface
(sim,
@ -400,6 +410,15 @@ public:
Simplex s2 = vh2->first;
intersect(p1,p2, s1,s2, p);
}
void intersect(const CMCT_Vertex_handle vh1,
const CMCT_Vertex_handle vh2,
const Simplex &s,
Bare_point &p) const {
Bare_point p1 = mc_triangulator->location(vh1, gt);
Bare_point p2 = mc_triangulator->location(vh2, gt);
Simplex sp = s;
intersect(p1,p2, sp,sp, p);
}
void intersect(Bare_point &p1, Bare_point &p2,
Simplex &s1, Simplex &s2,
@ -432,20 +451,15 @@ public:
p = midpoint(p1, p2);
}
template <class Point>
void intersect_with_transversal_segment(Point &p) const {
typedef typename Point::R Traits;
typedef typename Traits::Plane_3 Plane;
typedef typename Traits::Line_3 Line;
Mixed_complex_traits_3<Traits> point_traits(gt.get_shrink());
Cartesian_converter<Traits,
typename Geometric_traits::Bare_point::R> converter;
void intersect_with_transversal_segment(Bare_point &p) const {
typedef typename Geometric_traits::Kernel::Plane_3 Plane;
typedef typename Geometric_traits::Kernel::Line_3 Line;
Simplex sim = locate_mixed(converter(p));
Simplex sim = locate_mixed(p);
CMCT_Cell tet = locate_tet(p, sim);
// get transversal segment:
Point p1, p2;
Bare_point p1, p2;
// Compute signs on vertices and sort them:
int nIn = 0;
@ -460,12 +474,12 @@ public:
Object obj;
if (nIn==1) {
p1 = mc_triangulator->location(tet.vertex(sortedV[0]), point_traits);
p1 = mc_triangulator->location(tet.vertex(sortedV[0]), gt);
obj = CGAL::intersection(
Plane
(mc_triangulator->location(tet.vertex(sortedV[1]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[2]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[3]), point_traits)),
(mc_triangulator->location(tet.vertex(sortedV[1]), gt),
mc_triangulator->location(tet.vertex(sortedV[2]), gt),
mc_triangulator->location(tet.vertex(sortedV[3]), gt)),
Line(p1, p));
if ( !assign(p2, obj) ) {
CGAL_assertion_msg(false,"intersection: no intersection.");
@ -473,33 +487,33 @@ public:
} else if (nIn==2) {
obj = CGAL::intersection(
Plane
(mc_triangulator->location(tet.vertex(sortedV[2]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[3]), point_traits),
(mc_triangulator->location(tet.vertex(sortedV[2]), gt),
mc_triangulator->location(tet.vertex(sortedV[3]), gt),
p),
Line(
mc_triangulator->location(tet.vertex(sortedV[0]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[1]), point_traits)));
mc_triangulator->location(tet.vertex(sortedV[0]), gt),
mc_triangulator->location(tet.vertex(sortedV[1]), gt)));
if ( !assign(p1, obj) ) {
CGAL_assertion_msg(false,"intersection: no intersection.");
}
obj = CGAL::intersection(
Plane
(mc_triangulator->location(tet.vertex(sortedV[0]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[1]), point_traits),
(mc_triangulator->location(tet.vertex(sortedV[0]), gt),
mc_triangulator->location(tet.vertex(sortedV[1]), gt),
p),
Line(
mc_triangulator->location(tet.vertex(sortedV[2]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[3]), point_traits)));
mc_triangulator->location(tet.vertex(sortedV[2]), gt),
mc_triangulator->location(tet.vertex(sortedV[3]), gt)));
if ( !assign(p2, obj) ) {
CGAL_assertion_msg(false,"intersection: no intersection.");
}
} else if (nIn==3) {
p2 = mc_triangulator->location(tet.vertex(sortedV[3]), point_traits);
p2 = mc_triangulator->location(tet.vertex(sortedV[3]), gt);
obj = CGAL::intersection(
Plane(
mc_triangulator->location(tet.vertex(sortedV[0]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[1]), point_traits),
mc_triangulator->location(tet.vertex(sortedV[2]), point_traits)),
mc_triangulator->location(tet.vertex(sortedV[0]), gt),
mc_triangulator->location(tet.vertex(sortedV[1]), gt),
mc_triangulator->location(tet.vertex(sortedV[2]), gt)),
Line(p2, p));
if ( !assign(p1, obj) ) {
CGAL_assertion_msg(false,"intersection: no intersection.");

5
Skin_surface_3/include/CGAL/Skin_surface_polyhedral_items_3.h
View File

@ -28,10 +28,9 @@ CGAL_BEGIN_NAMESPACE
template <class Refs, class SkinSurface3>
struct Skin_Surface_polyhedral_face : public CGAL::HalfedgeDS_face_base<Refs> {
typedef SkinSurface3 Skin_surface;
// typedef typename Skin_surface::Triangulated_mixed_complex TMC;
// typedef typename TMC::Cell_handle Triang_Cell_handle;
typedef typename SkinSurface3::Simplex Simplex;
// Triang_Cell_handle triang_ch;
Simplex sim;
};
template < class SkinSurface3 >

17
Skin_surface_3/include/CGAL/Skin_surface_refinement_traits_3.h
View File

@ -48,13 +48,24 @@ public:
}
P_point to_surface(P_vertex_handle vh) {
P_point result = vh->point();
typename Skin_surface::Bare_point result =
Cartesian_converter<P_traits,
typename Skin_surface::Geometric_traits::Kernel>()(vh->point());
ss_3.intersect_with_transversal_segment(result);
return result;
return
Cartesian_converter
<typename Skin_surface::Geometric_traits::Kernel, P_traits>()( result );
}
P_vector normal(P_vertex_handle vh) {
return ss_3.normal(vh->point());
// Convert to and from the skin surface kernel
return
Cartesian_converter
<typename Skin_surface::Geometric_traits::Kernel,
P_traits>()( ss_3.normal
(Cartesian_converter<P_traits,
typename Skin_surface::Geometric_traits::Kernel>()
(vh->point())));
}
protected:

36
Skin_surface_3/include/CGAL/triangulate_mixed_complex_3.h
View File

@ -37,16 +37,18 @@ class CMCT_Cell {
public:
typedef CMCT_Cell<CMCT> Self;
typedef typename CMCT::Vertex_handle Vertex_handle;
typedef typename CMCT::Rt_Simplex Rt_Simplex;
CMCT_Cell() {
}
CMCT_Cell(Vertex_handle vs[]) {
CMCT_Cell(Vertex_handle vs[], Rt_Simplex &sim) : sim_(sim) {
for (int i=0; i<4; i++) _vs[i] = vs[i];
}
CMCT_Cell(const Vertex_handle &vh0,
const Vertex_handle &vh1,
const Vertex_handle &vh2,
const Vertex_handle &vh3)
const Vertex_handle &vh3,
Rt_Simplex &sim) : sim_(sim)
{
_vs[0] = vh0;
_vs[1] = vh1;
@ -61,6 +63,9 @@ public:
CGAL_assertion((0 <= i) && (i<4));
return _vs[i];
}
const Rt_Simplex &mixed_cell() const {
return sim_;
}
bool operator==(const Self &other) const {
return ((_vs[0] == other._vs[0]) &&
(_vs[1] == other._vs[1]) &&
@ -70,6 +75,7 @@ public:
bool operator!=(const Self &other) const { return !operator==(other); }
private:
Vertex_handle _vs[4];
Rt_Simplex sim_;
};
template < class CMCT >
@ -501,6 +507,7 @@ private:
template <class OutputIteratorCells>
void add_cell(Vertex_handle vh[],
int orient,
Rt_Simplex &simplex,
OutputIteratorCells cells) const;
Vertex_handle get_vertex(Rt_Simplex &sDel, Rt_Simplex &sVor) const;
@ -941,7 +948,8 @@ construct_0_cell(Rt_Vertex_handle rt_vh, OutputIteratorCells cells) const
CGAL_assertion(sVor_e != sVor_f);
CGAL_assertion(sVor_f != sVor_c);
add_cell(vh,(index + (j==(i%3+1)? 1:0))&1, cells);
add_cell(vh,(index + (j==(i%3+1)? 1:0))&1,
simplex, cells);
}
}
}
@ -998,15 +1006,15 @@ construct_1_cell(const Rt_Edge &e,
orient = (index1 + (fi==1))&1;
}
// vh: dimension are (01,11,12,13)
add_cell(vh,orient,cells);
add_cell(vh,orient,mixed_cell_simplex,cells);
vh[1] = get_vertex(sDel_v, sVor_f);
// vh: dimension are (01,02,12,13)
add_cell(vh,1-orient,cells);
add_cell(vh,1-orient,mixed_cell_simplex,cells);
vh[2] = get_vertex(sDel_v, sVor_c);
// vh: dimension are (01,02,03,13)
add_cell(vh,orient,cells);
add_cell(vh,orient,mixed_cell_simplex,cells);
}
}
}
@ -1071,13 +1079,13 @@ construct_2_cell(const Rt_Facet &f,
orient = (index1 + (((4+index0-fi)&3)==1))&1;
}
add_cell(vh,orient,cells);
add_cell(vh,orient,simplex,cells);
vh[2] = get_vertex(sDel_e, sVor_c);
add_cell(vh,1-orient,cells);
add_cell(vh,1-orient,simplex,cells);
vh[1] = get_vertex(sDel_v, sVor_c);
add_cell(vh,orient,cells);
add_cell(vh,orient,simplex,cells);
}
}
}
@ -1135,7 +1143,7 @@ construct_3_cell(Rt_Cell_handle rt_ch,
} else {
orient = (index1 + (vi==3))&1;
}
add_cell(vh, orient, cells);
add_cell(vh, orient, simplex, cells);
}
}
}
@ -1181,7 +1189,9 @@ template <class MixedComplexTraits_3>
template <class OutputIteratorCells>
void
Combinatorial_mixed_complex_triangulator_3<MixedComplexTraits_3>::
add_cell(Vertex_handle vh[], int orient, OutputIteratorCells cells) const {
add_cell(Vertex_handle vh[], int orient,
Rt_Simplex &simplex,
OutputIteratorCells cells) const {
assert((orient==0) || (orient==1));
assert(*vh[0] != Vertex()); assert(*vh[1] != Vertex());
assert(*vh[2] != Vertex()); assert(*vh[3] != Vertex());
@ -1189,9 +1199,9 @@ add_cell(Vertex_handle vh[], int orient, OutputIteratorCells cells) const {
assert(*vh[1] != *vh[2]); assert(*vh[1] != *vh[3]); assert(*vh[2] != *vh[3]);
if (orient) {
*cells++ = Cell(vh[0], vh[1], vh[2], vh[3]);
*cells++ = Cell(vh[0], vh[1], vh[2], vh[3], simplex);
} else {
*cells++ = Cell(vh[0], vh[1], vh[3], vh[2]);
*cells++ = Cell(vh[0], vh[1], vh[3], vh[2], simplex);
}
}