mirror of https://github.com/CGAL/cgal
Merge branch 'Surface_reconstruction_points_3-make_it_faster-afabri'
After approval of the release manager Conflicts: Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3_ref/Poisson_reconstruction_function.tex Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h
This commit is contained in:
Sébastien Loriot
commit
485d1bc5e0
|
|
@ -4291,6 +4291,7 @@ Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/contouri
|
|||
Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/eros.jpg -text svneol=unset#image/jpeg
|
||||
Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/fig/surface_reconstruction_points.jpg -text svneol=unset#image/jpeg
|
||||
Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/fig/surface_reconstruction_points_detail.png -text svneol=unset#image/png
|
||||
Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/fig/two-passes.png -text svneol=unset#image/png
|
||||
Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/flipped_normals.jpg -text svneol=unset#image/jpeg
|
||||
Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/holes_bad.jpg -text svneol=unset#image/jpeg
|
||||
Surface_reconstruction_points_3/doc_tex/Surface_reconstruction_points_3/holes_good.jpg -text svneol=unset#image/jpeg
|
||||
|
|
|
|||
|
|
@ -53,14 +53,13 @@ public:
|
|||
/// Create a square matrix initialized with zeros.
|
||||
Eigen_sparse_matrix(int dim, ///< Matrix dimension.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_matrix(dim,dim)
|
||||
: m_is_uptodate(false), m_matrix(dim,dim)
|
||||
{
|
||||
CGAL_precondition(dim > 0);
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
|
||||
// reserve memory for a regular 3D grid
|
||||
m_matrix.reserve(Eigen::VectorXi::Constant(m_matrix.outerSize(), 27));
|
||||
m_triplets.reserve(dim);
|
||||
}
|
||||
|
||||
/// Create a rectangular matrix initialized with zeros.
|
||||
|
|
@ -69,7 +68,7 @@ public:
|
|||
Eigen_sparse_matrix(int rows, ///< Number of rows.
|
||||
int columns, ///< Number of columns.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_matrix(rows,columns)
|
||||
: m_is_uptodate(false), m_matrix(rows,columns)
|
||||
{
|
||||
CGAL_precondition(rows > 0);
|
||||
CGAL_precondition(columns > 0);
|
||||
|
|
@ -78,9 +77,8 @@ public:
|
|||
}
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
|
||||
// reserve memory for a regular 3D grid
|
||||
m_matrix.reserve(Eigen::VectorXi::Constant(m_matrix.outerSize(), 27));
|
||||
m_triplets.reserve(rows);
|
||||
}
|
||||
|
||||
/// Delete this object and the wrapped TAUCS matrix.
|
||||
|
|
@ -112,9 +110,9 @@ public:
|
|||
|
||||
if (m_is_symmetric && (j > i))
|
||||
return;
|
||||
|
||||
if(new_coef) m_matrix.insert(i,j) = val;
|
||||
else m_matrix.coeffRef(i,j) = val;
|
||||
|
||||
m_triplets.push_back(Triplet(i,j,val));
|
||||
m_is_uptodate = false;
|
||||
}
|
||||
|
||||
/// Write access to a matrix coefficient: a_ij <- a_ij+val.
|
||||
|
|
@ -134,13 +132,19 @@ public:
|
|||
if (m_is_symmetric && (j > i))
|
||||
return;
|
||||
|
||||
m_matrix.coeffRef(i,j) += val;
|
||||
m_triplets.push_back(Triplet(i,j,val));
|
||||
m_is_uptodate = false;
|
||||
}
|
||||
|
||||
|
||||
|
||||
const EigenType& eigen_object() const
|
||||
{
|
||||
if(!m_is_uptodate)
|
||||
{
|
||||
m_matrix.setFromTriplets(m_triplets.begin(), m_triplets.end());
|
||||
m_is_uptodate = true;
|
||||
}
|
||||
// turns the matrix into compressed mode:
|
||||
// -> release some memory
|
||||
// -> required for some external solvers
|
||||
|
|
@ -157,6 +161,10 @@ private:
|
|||
|
||||
// Fields
|
||||
private:
|
||||
|
||||
mutable bool m_is_uptodate;
|
||||
typedef Eigen::Triplet<T,int> Triplet;
|
||||
mutable std::vector<Triplet> m_triplets;
|
||||
|
||||
mutable EigenType m_matrix;
|
||||
|
||||
|
|
|
|||
|
|
@ -36,10 +36,10 @@
|
|||
namespace CGAL {
|
||||
namespace Surface_mesher {
|
||||
|
||||
template <typename K2 = Exact_predicates_inexact_constructions_kernel>
|
||||
class Surface_mesh_default_triangulation_3_generator {
|
||||
|
||||
// traits class
|
||||
typedef Exact_predicates_inexact_constructions_kernel K2;
|
||||
typedef Robust_circumcenter_traits_3<K2> K;
|
||||
|
||||
// vertex and cell types
|
||||
|
|
@ -60,7 +60,7 @@ namespace CGAL {
|
|||
|
||||
} // end Surface_mesher
|
||||
|
||||
typedef Surface_mesher::Surface_mesh_default_triangulation_3_generator::Type
|
||||
typedef Surface_mesher::Surface_mesh_default_triangulation_3_generator<>::Type
|
||||
Surface_mesh_default_triangulation_3;
|
||||
|
||||
} // end namespace CGAL
|
||||
|
|
|
|||
|
|
@ -12,7 +12,7 @@
|
|||
#include <CGAL/Timer.h>
|
||||
#include <CGAL/Surface_mesh_default_triangulation_3.h>
|
||||
#include <CGAL/make_surface_mesh.h>
|
||||
#include <CGAL/Implicit_surface_3.h>
|
||||
#include <CGAL/Poisson_implicit_surface_3.h>
|
||||
#include <CGAL/IO/output_surface_facets_to_polyhedron.h>
|
||||
#include <CGAL/Poisson_reconstruction_function.h>
|
||||
#include <CGAL/compute_average_spacing.h>
|
||||
|
|
@ -37,13 +37,44 @@ typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_fun
|
|||
// Surface mesher
|
||||
typedef CGAL::Surface_mesh_default_triangulation_3 STr;
|
||||
typedef CGAL::Surface_mesh_complex_2_in_triangulation_3<STr> C2t3;
|
||||
typedef CGAL::Implicit_surface_3<Kernel, Poisson_reconstruction_function> Surface_3;
|
||||
typedef CGAL::Poisson_implicit_surface_3<Kernel, Poisson_reconstruction_function> Surface_3;
|
||||
|
||||
// AABB tree
|
||||
typedef CGAL::AABB_polyhedron_triangle_primitive<Kernel,Polyhedron> Primitive;
|
||||
typedef CGAL::AABB_traits<Kernel, Primitive> AABB_traits;
|
||||
typedef CGAL::AABB_tree<AABB_traits> AABB_tree;
|
||||
|
||||
struct Counter {
|
||||
int i, N;
|
||||
|
||||
Counter(int N)
|
||||
: i(0), N(N)
|
||||
{}
|
||||
|
||||
void operator()()
|
||||
{
|
||||
i++;
|
||||
if(i == N){
|
||||
std::cerr << "Counter reached " << N << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
struct InsertVisitor {
|
||||
|
||||
Counter& c;
|
||||
InsertVisitor(Counter& c)
|
||||
: c(c)
|
||||
{}
|
||||
|
||||
void before_insertion()
|
||||
{
|
||||
c();
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
// Poisson reconstruction method:
|
||||
// Reconstructs a surface mesh from a point set and returns it as a polyhedron.
|
||||
|
|
@ -73,6 +104,10 @@ Polyhedron* poisson_reconstruct(const Point_set& points,
|
|||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
Counter counter(std::distance(points.begin(), points.end()));
|
||||
InsertVisitor visitor(counter) ;
|
||||
|
||||
CGAL::Timer reconstruction_timer; reconstruction_timer.start();
|
||||
|
||||
//***************************************
|
||||
|
|
@ -89,7 +124,8 @@ Polyhedron* poisson_reconstruct(const Point_set& points,
|
|||
// The position property map can be omitted here as we use iterators over Point_3 elements.
|
||||
Poisson_reconstruction_function function(
|
||||
points.begin(), points.end(),
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()));
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()),
|
||||
visitor);
|
||||
|
||||
bool ok = false;
|
||||
#ifdef CGAL_TAUCS_ENABLED
|
||||
|
|
@ -129,6 +165,7 @@ Polyhedron* poisson_reconstruct(const Point_set& points,
|
|||
|
||||
// Computes the Poisson indicator function f()
|
||||
// at each vertex of the triangulation.
|
||||
|
||||
if ( ! ok )
|
||||
{
|
||||
std::cerr << "Error: cannot compute implicit function" << std::endl;
|
||||
|
|
@ -149,6 +186,7 @@ Polyhedron* poisson_reconstruct(const Point_set& points,
|
|||
FT average_spacing = CGAL::compute_average_spacing(points.begin(), points.end(),
|
||||
6 /* knn = 1 ring */);
|
||||
|
||||
|
||||
// Gets one point inside the implicit surface
|
||||
Point inner_point = function.get_inner_point();
|
||||
FT inner_point_value = function(inner_point);
|
||||
|
|
@ -183,6 +221,9 @@ Polyhedron* poisson_reconstruct(const Point_set& points,
|
|||
<< " dichotomy error=distance/"<<sm_distance*average_spacing/sm_dichotomy_error<<",\n"
|
||||
<< " Manifold_with_boundary_tag)\n";
|
||||
|
||||
|
||||
|
||||
|
||||
// Generates surface mesh with manifold option
|
||||
STr tr; // 3D Delaunay triangulation for surface mesh generation
|
||||
C2t3 c2t3(tr); // 2D complex in 3D Delaunay triangulation
|
||||
|
|
|
|||
Binary file not shown.
|
|
@ -143,7 +143,7 @@ The reconstruction algorithm expects a sufficiently dense point set. Although th
|
|||
|
||||
\subsubsection{Large Holes}
|
||||
|
||||
The reconstruction is devised to solve for an implicit function which is an approximate indicator function of an inferred solid. For this reason the contouring algorithm always extracts a closed surface mesh and hence is able to fill the small holes where data are missing due, e.g., to occlusions during acquisition. In case of large holes the algorithm still closes them all but sometimes in an unexpected manner. In addition the resulting piecewise linear implicit function may exhibit large triangle patches and sharp creases as the 3D Delaunay triangulation used for solving is very coarse where the holes are filled (see Figure~\ref{Surface_reconstruction_points_3-fig-holes_bad}).
|
||||
The reconstruction is devised to solve for an implicit function which is an approximate indicator function of an inferred solid. For this reason the contouring algorithm always extracts a closed surface mesh and hence is able to fill the small holes where data are missing due, e.g., to occlusions during acquisition(see Figure~\ref{Surface_reconstruction_points_3-fig-holes_bad}).
|
||||
|
||||
% Insert image holes_bad.jpg/.eps
|
||||
\begin{center}
|
||||
|
|
@ -163,6 +163,21 @@ The reconstruction is devised to solve for an implicit function which is an appr
|
|||
\end{figure}
|
||||
\end{center}
|
||||
|
||||
In case of large holes the algorithm still closes them all but the resulting piecewise linear implicit function may exhibit large triangle patches and sharp creases as the 3D Delaunay triangulation used for solving is very coarse where the holes are filled. This can be avoided by a two pass approach. The first pass for a subset of the points serves to get an approximation of the surface at the holes. This surface then serves to compute a smoother 3D Delaunay triangulation for the second pass with the full set of points.
|
||||
|
||||
\begin{center}
|
||||
\begin{ccTexOnly}
|
||||
\includegraphics[width=0.8\textwidth]{Surface_reconstruction_points_3/fig/two-passes.png}
|
||||
\end{ccTexOnly}
|
||||
\begin{ccHtmlOnly}
|
||||
<img style="max-width: 80%;" border=0 src="./fig/two-passes.png"><P>
|
||||
\end{ccHtmlOnly}
|
||||
\begin{figure}[h]
|
||||
\caption{Left: The wrist. Middle: one pass. Right: two passes.}
|
||||
\label{Surface_reconstruction_points_3-fig-two_passes}
|
||||
\end{figure}
|
||||
\end{center}
|
||||
|
||||
|
||||
\subsubsection{Wrongly Oriented Normals}
|
||||
|
||||
|
|
|
|||
Binary file not shown.
|
Before Width: | Height: | Size: 66 KiB After Width: | Height: | Size: 63 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 437 KiB |
|
|
@ -0,0 +1,7 @@
|
|||
|
||||
|
||||
\section{Design and Implementation History}
|
||||
|
||||
The initial implementation was essentially done by Laurent Saboret, guided by Pierre Alliez and Ga\"el Guennebaud.
|
||||
For later releases of the package Andreas Fabri worked on performance improvements, and Laurent Rineau added the
|
||||
two passes for dealing with holes.
|
||||
|
|
@ -13,3 +13,6 @@
|
|||
\input{Surface_reconstruction_points_3/contouring}
|
||||
\input{Surface_reconstruction_points_3/output}
|
||||
\input{Surface_reconstruction_points_3/case_studies}
|
||||
\input{Surface_reconstruction_points_3/performances}
|
||||
\input{Surface_reconstruction_points_3/implementation}
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,121 @@
|
|||
\section{Performances}
|
||||
\label{surface_reconstruction_section_performances}
|
||||
|
||||
We provide some performance numbers for scanning data. We measure the Poisson implicit function computation time, the contouring time for a range of approximation distances, the memory occupancy as well as the influence of the point set simplification. The machine used is a PC running Windows 7 64 bits with an Intel CPU Core 2 Duo processor clocked at 2.81 GHz and with 8 GB of RAM. The software is compiled with Visual C++ 2010 (VC9) compiler with the 03 option which maximizes speed. All measurements were done using the \ccThirdPartyEigen\ library.
|
||||
|
||||
|
||||
\subsection{Poisson implicit function}
|
||||
|
||||
The point set chosen for benchmarking the Poisson implicit function is the Bimba con Nastrino point set (1.6 million points) depicted by Figure~\ref{Surface_reconstruction_points_3-fig-poisson_bench}. We measure the Poisson implicit function computation (i.e., the call to \ccc{Poisson_reconstruction_function::compute_implicit_function()} denoted by Poisson solve hereafter) for this point set as well as for simplified versions obtained through random simplification. The following table provides Poisson solve computation times in seconds for an increasing number of points.
|
||||
|
||||
\begin{tabular}{|c|c|}
|
||||
\hline
|
||||
Number of points (x1000) & Poisson solve duration (in s) \\
|
||||
\hline
|
||||
60 & 15 \\
|
||||
100 & 25 \\
|
||||
250 & 96 \\
|
||||
500 & 150 \\
|
||||
1,000 & 249 \\
|
||||
1,800 & 478 \\
|
||||
\hline
|
||||
\end{tabular}
|
||||
|
||||
|
||||
|
||||
\subsection{Contouring}
|
||||
|
||||
The point set chosen for benchmarking the contouring stage is the Bimba con Nastrino point set simplified to 100k points. We measure the contouring (i.e. the call to \ccc{make_surface_mesh()}) duration and the reconstruction error for a range of approximation distances.
|
||||
The reconstruction error is expressed as the average distance from input points to the reconstructed surface in mm (the Bimba con Nastrino statue is 324 mm tall).
|
||||
|
||||
\begin{tabular}{|c|c|c|}
|
||||
\hline
|
||||
Approx. distance (*average spacing) & Contouring duration (in s) & Reconstruction error (mm) \\
|
||||
\hline
|
||||
0.1 & 19.2 & 0.055 \\
|
||||
0.25 & 6.9 & 0.106 \\
|
||||
0.5 & 3.2 & 0.18 \\
|
||||
1 & 1.65 & 0.36 \\
|
||||
2 & 0.8 & 0.76 \\
|
||||
\hline
|
||||
\end{tabular}
|
||||
|
||||
|
||||
% Insert image contouring_bench.jpg/.eps
|
||||
\begin{center}
|
||||
\begin{ccTexOnly}
|
||||
\includegraphics[width=1.0\textwidth]{Surface_reconstruction_points_3/contouring_bench}
|
||||
\end{ccTexOnly}
|
||||
\begin{ccHtmlOnly}
|
||||
<img style="max-width: 100%;" border=0 src="./contouring_bench.jpg"><P>
|
||||
\end{ccHtmlOnly}
|
||||
\begin{figure}[h]
|
||||
\caption{Contouring duration (in s) and reconstruction error (mm)
|
||||
against several approximation distance parameters
|
||||
for the Bimba con Nastrino point set simplified to 100k points.}
|
||||
\label{Surface_reconstruction_points_3-fig-contouring_bench}
|
||||
\end{figure}
|
||||
\end{center}
|
||||
|
||||
|
||||
|
||||
\subsection{Memory}
|
||||
|
||||
We measure the memory occupancy for the reconstruction of the full Bimba con Nastrino point set (1.8 millions points) as well as for simplified versions.\\
|
||||
The Poisson implicit function computation has a memory peak when solving the Poisson linear system using the sparse linear solver. \\
|
||||
|
||||
\begin{tabular}{|c|c|}
|
||||
\hline
|
||||
Number of points (x1000) & Memory occupancy (MBytes) \\
|
||||
\hline
|
||||
60 & 180 \\
|
||||
100 & 270 \\
|
||||
250 & 790 \\
|
||||
500 & 1300 \\
|
||||
1,000 & 2200 \\
|
||||
1,800 & 3800 \\
|
||||
\hline
|
||||
\end{tabular}
|
||||
|
||||
|
||||
|
||||
\subsection{Point Set Simplification}
|
||||
|
||||
Due to the memory limitations described above, we recommend to simplify the point sets captured by laser scanners.\\
|
||||
We measure the reconstruction error for the Bimba con Nastrino point set (1.6M points) as well as for simplified versions. All reconstructions use the recommended contouring parameter approximation distance = 0.25 * the input point set's average spacing.
|
||||
The reconstruction error is expressed as the average distance from input points to the reconstructed surface in mm (the Bimba con Nastrino statue is 324 mm tall).
|
||||
|
||||
\begin{tabular}{|c|c|}
|
||||
\hline
|
||||
Number of points (x1000) & Reconstruction error (mm) \\
|
||||
\hline
|
||||
60 & 0.27 \\
|
||||
120 & 0.15 \\
|
||||
250 & 0.11 \\
|
||||
500 & 0.079 \\
|
||||
1,000 & 0.066 \\
|
||||
1,500 & 0.061 \\
|
||||
1,600 & 0.06 \\
|
||||
\hline
|
||||
\end{tabular}
|
||||
|
||||
% Insert image simplification_bench.jpg/.eps
|
||||
\begin{center}
|
||||
\begin{ccTexOnly}
|
||||
\includegraphics[width=1.0\textwidth]{Surface_reconstruction_points_3/simplification_bench}
|
||||
\end{ccTexOnly}
|
||||
\begin{ccHtmlOnly}
|
||||
<img style="max-width: 100%;" border=0 src="./simplification_bench.jpg"><P>
|
||||
\end{ccHtmlOnly}
|
||||
\begin{figure}[h]
|
||||
\caption{Reconstruction error (mm) against number of points
|
||||
for the Bimba con Nastrino point set with 1.6M points
|
||||
as well as for simplified versions.}
|
||||
\label{Surface_reconstruction_points_3-fig-simplification_bench}
|
||||
\end{figure}
|
||||
\end{center}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -6,14 +6,8 @@ Given a set of 3D points with oriented normals (denoted oriented points in the s
|
|||
|
||||
\subsection{Interface}
|
||||
|
||||
The class template declaration is:
|
||||
|
||||
template$<$ \\
|
||||
class Gt$>$ \\
|
||||
class \ccc{Poisson_reconstruction_function};
|
||||
|
||||
with \\
|
||||
\ccc{Gt}: Geometric traits class.
|
||||
The class template declaration is \ccc{template<class Gt> class Poisson_reconstruction_function} where
|
||||
\ccc{Gt} is a geometric traits class.
|
||||
|
||||
Creation:
|
||||
|
||||
|
|
@ -36,9 +30,10 @@ The main operations are:
|
|||
Returns a sphere bounding the inferred surface.
|
||||
}
|
||||
|
||||
\ccFunction{bool compute_implicit_function();}
|
||||
\ccFunction{bool compute_implicit_function(bool smoother_hole_filling = false);}
|
||||
{
|
||||
The function \ccc{compute_implicit_function}() must be called after each insertion of oriented points. It computes the piecewise linear scalar function operator() by: applying Delaunay refinement, solving for operator() at each vertex of the triangulation with a sparse linear solver, and shifting and orienting operator() such that it is 0 at all input points and negative inside the inferred surface. \\
|
||||
The optional parameter \ccc{smoother_hole_filling} controls if the Delaunay refinement is done for the input points, or for an approximation of the surface obtained from a first pass of the algorithm on a sample of the points.
|
||||
Returns false if the linear solver fails.\\
|
||||
The sparse solver used for this step is a parameter of the function. We recommend to use the class \ccc{CGAL::Eigen_solver_traits<T>}
|
||||
instantiated with the iterative conjugate gradient solver \ccc{Eigen::ConjugateGradient} for \ccc{double} (which is the default when \ccThirdPartyEigen\ is available and \ccc{CGAL_EIGEN3_ENABLED} is defined).
|
||||
|
|
|
|||
|
|
@ -131,7 +131,7 @@ Creates a Poisson implicit function from the [first, beyond) range of points.
|
|||
Returns a sphere bounding the inferred surface.
|
||||
}
|
||||
\ccGlue
|
||||
\ccMethod{template<class SparseLinearAlgebraTraits_d> bool compute_implicit_function(SparseLinearAlgebraTraits_d solver);}
|
||||
\ccMethod{template<class SparseLinearAlgebraTraits_d> bool compute_implicit_function(SparseLinearAlgebraTraits_d solver, bool smoother_hole_filling = false);}
|
||||
{
|
||||
The function \ccc{compute_implicit_function}() must be called after the insertion of oriented points. It computes the piecewise linear scalar function operator() by: applying Delaunay refinement, solving for operator() at each vertex of the triangulation with a sparse linear solver, and shifting and orienting operator() such that it is 0 at all input points and negative inside the inferred surface.
|
||||
\ccCommentHeading{Template parameters} \\
|
||||
|
|
@ -140,7 +140,9 @@ If \eigen\ 3.1 (or greater) is available and \ccc{CGAL_EIGEN3_ENABLED} is define
|
|||
as default solver is provided.
|
||||
\ccCommentHeading{Returns} false if the linear solver fails.
|
||||
\ccCommentHeading{Parameters} \\
|
||||
\ccc{solver}: sparse linear solver.
|
||||
\ccc{solver}: sparse linear solver.\\
|
||||
\ccc{smoother_hole_filling} controls if the Delaunay refinement is done for the input points, or for an approximation of the surface obtained from a first pass of the algorithm on a sample of the points.\\
|
||||
Returns false if the linear solver fails.\\
|
||||
}
|
||||
\ccGlue
|
||||
\ccMethod{FT operator()(const Point& p) const;}
|
||||
|
|
|
|||
|
|
@ -19,7 +19,7 @@
|
|||
#include <CGAL/IO/Polyhedron_iostream.h>
|
||||
#include <CGAL/Surface_mesh_default_triangulation_3.h>
|
||||
#include <CGAL/make_surface_mesh.h>
|
||||
#include <CGAL/Implicit_surface_3.h>
|
||||
#include <CGAL/Poisson_implicit_surface_3.h>
|
||||
#include <CGAL/IO/output_surface_facets_to_polyhedron.h>
|
||||
#include <CGAL/Poisson_reconstruction_function.h>
|
||||
#include <CGAL/Point_with_normal_3.h>
|
||||
|
|
@ -61,7 +61,7 @@ typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_fun
|
|||
// Surface mesher
|
||||
typedef CGAL::Surface_mesh_default_triangulation_3 STr;
|
||||
typedef CGAL::Surface_mesh_complex_2_in_triangulation_3<STr> C2t3;
|
||||
typedef CGAL::Implicit_surface_3<Kernel, Poisson_reconstruction_function> Surface_3;
|
||||
typedef CGAL::Poisson_implicit_surface_3<Kernel, Poisson_reconstruction_function> Surface_3;
|
||||
|
||||
// AABB tree
|
||||
typedef CGAL::AABB_polyhedron_triangle_primitive<Kernel,Polyhedron> Primitive;
|
||||
|
|
@ -69,6 +69,37 @@ typedef CGAL::AABB_traits<Kernel, Primitive> AABB_traits;
|
|||
typedef CGAL::AABB_tree<AABB_traits> AABB_tree;
|
||||
|
||||
|
||||
struct Counter {
|
||||
int i, N;
|
||||
Counter(int N)
|
||||
: i(0), N(N)
|
||||
{}
|
||||
|
||||
void operator()()
|
||||
{
|
||||
i++;
|
||||
if(i == N){
|
||||
std::cerr << "Counter reached " << N << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
struct InsertVisitor {
|
||||
|
||||
Counter& c;
|
||||
InsertVisitor(Counter& c)
|
||||
: c(c)
|
||||
{}
|
||||
|
||||
void before_insertion()
|
||||
{
|
||||
c();
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// main()
|
||||
// ----------------------------------------------------------------------------
|
||||
|
|
@ -113,6 +144,8 @@ int main(int argc, char * argv[])
|
|||
std::string solver_name = "no_solver_available";
|
||||
#endif
|
||||
#endif
|
||||
double approximation_ratio = 0.02;
|
||||
double average_spacing_ratio = 5;
|
||||
|
||||
// decode parameters
|
||||
std::string input_filename = argv[1];
|
||||
|
|
@ -125,6 +158,10 @@ int main(int argc, char * argv[])
|
|||
sm_distance = atof(argv[++i]);
|
||||
else if (std::string(argv[i])=="-solver")
|
||||
solver_name = argv[++i];
|
||||
else if (std::string(argv[i])=="-approx")
|
||||
approximation_ratio = atof(argv[++i]);
|
||||
else if (std::string(argv[i])=="-ratio")
|
||||
average_spacing_ratio = atof(argv[++i]);
|
||||
else {
|
||||
std::cerr << "Error: invalid option " << argv[i] << "\n";
|
||||
return EXIT_FAILURE;
|
||||
|
|
@ -215,6 +252,11 @@ int main(int argc, char * argv[])
|
|||
|
||||
CGAL::Timer reconstruction_timer; reconstruction_timer.start();
|
||||
|
||||
|
||||
Counter counter(std::distance(points.begin(), points.end()));
|
||||
InsertVisitor visitor(counter) ;
|
||||
|
||||
|
||||
//***************************************
|
||||
// Computes implicit function
|
||||
//***************************************
|
||||
|
|
@ -227,7 +269,8 @@ int main(int argc, char * argv[])
|
|||
// The position property map can be omitted here as we use iterators over Point_3 elements.
|
||||
Poisson_reconstruction_function function(
|
||||
points.begin(), points.end(),
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()));
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()),
|
||||
visitor);
|
||||
|
||||
#ifdef CGAL_TAUCS_ENABLED
|
||||
if (solver_name == "taucs")
|
||||
|
|
@ -249,7 +292,10 @@ int main(int argc, char * argv[])
|
|||
|
||||
// Computes the Poisson indicator function f()
|
||||
// at each vertex of the triangulation.
|
||||
if ( ! function.compute_implicit_function(solver) )
|
||||
|
||||
if ( ! function.compute_implicit_function(solver, visitor,
|
||||
approximation_ratio,
|
||||
average_spacing_ratio) )
|
||||
{
|
||||
std::cerr << "Error: cannot compute implicit function" << std::endl;
|
||||
return EXIT_FAILURE;
|
||||
|
|
@ -262,7 +308,10 @@ int main(int argc, char * argv[])
|
|||
if (solver_name == "eigen")
|
||||
{
|
||||
std::cerr << "Use Eigen 3\n";
|
||||
if ( ! function.compute_implicit_function() )
|
||||
CGAL::Eigen_solver_traits<Eigen::ConjugateGradient<CGAL::Eigen_sparse_symmetric_matrix<double>::EigenType> > solver;
|
||||
if ( ! function.compute_implicit_function(solver, visitor,
|
||||
approximation_ratio,
|
||||
average_spacing_ratio) )
|
||||
{
|
||||
std::cerr << "Error: cannot compute implicit function" << std::endl;
|
||||
return EXIT_FAILURE;
|
||||
|
|
|
|||
|
|
@ -28,6 +28,38 @@ typedef CGAL::Surface_mesh_default_triangulation_3 STr;
|
|||
typedef CGAL::Surface_mesh_complex_2_in_triangulation_3<STr> C2t3;
|
||||
typedef CGAL::Implicit_surface_3<Kernel, Poisson_reconstruction_function> Surface_3;
|
||||
|
||||
|
||||
struct Counter {
|
||||
int i, N;
|
||||
Counter(int N)
|
||||
: i(0), N(N)
|
||||
{}
|
||||
|
||||
void operator()()
|
||||
{
|
||||
i++;
|
||||
if(i == N){
|
||||
std::cerr << "Counter reached " << N << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
struct InsertVisitor {
|
||||
|
||||
Counter& c;
|
||||
InsertVisitor(Counter& c)
|
||||
: c(c)
|
||||
{}
|
||||
|
||||
void before_insertion()
|
||||
{
|
||||
c();
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
int main(void)
|
||||
{
|
||||
// Poisson options
|
||||
|
|
@ -51,17 +83,20 @@ int main(void)
|
|||
return EXIT_FAILURE;
|
||||
}
|
||||
|
||||
Counter counter(std::distance(points.begin(), points.end()));
|
||||
InsertVisitor visitor(counter) ;
|
||||
// Creates implicit function from the read points using the default solver.
|
||||
|
||||
// Note: this method requires an iterator over points
|
||||
// + property maps to access each point's position and normal.
|
||||
// The position property map can be omitted here as we use iterators over Point_3 elements.
|
||||
Poisson_reconstruction_function function(
|
||||
points.begin(), points.end(),
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()));
|
||||
Poisson_reconstruction_function function(points.begin(), points.end(),
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()),
|
||||
visitor);
|
||||
|
||||
// Computes the Poisson indicator function f()
|
||||
// at each vertex of the triangulation.
|
||||
if ( ! function.compute_implicit_function() )
|
||||
if ( ! function.compute_implicit_function() )
|
||||
return EXIT_FAILURE;
|
||||
|
||||
// Computes average spacing
|
||||
|
|
|
|||
|
|
@ -0,0 +1,143 @@
|
|||
// Copyright (c) 2006-2007 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) : Laurent RINEAU
|
||||
|
||||
#ifndef CGAL_POISSON_IMPLICIT_SURFACE_3_H
|
||||
#define CGAL_POISSON_IMPLICIT_SURFACE_3_H
|
||||
|
||||
#include <CGAL/make_surface_mesh.h>
|
||||
#include <CGAL/Surface_mesher/Poisson_implicit_surface_oracle_3.h>
|
||||
|
||||
#include <functional>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
template<
|
||||
typename GT,
|
||||
typename Function_
|
||||
>
|
||||
class Poisson_implicit_surface_3
|
||||
{
|
||||
public:
|
||||
typedef GT Geom_traits;
|
||||
typedef typename Geom_traits::Sphere_3 Sphere_3;
|
||||
typedef typename Geom_traits::FT FT;
|
||||
typedef typename Geom_traits::Point_3 Point;
|
||||
typedef Function_ Function;
|
||||
typedef Poisson_implicit_surface_3<Geom_traits, Function> Self;
|
||||
|
||||
Function& function() { return func; }
|
||||
|
||||
typedef Surface_mesher::Poisson_implicit_surface_oracle_3<
|
||||
Geom_traits,
|
||||
Self> Surface_mesher_traits_3;
|
||||
|
||||
Poisson_implicit_surface_3(Function f,
|
||||
const Sphere_3 bounding_sphere,
|
||||
const FT error_bound = FT(1e-3),
|
||||
Geom_traits gt = Geom_traits())
|
||||
: func(f),
|
||||
sphere(bounding_sphere),
|
||||
gt(gt)
|
||||
{
|
||||
squared_error = error_bound * error_bound;
|
||||
squared_error = squared_error *
|
||||
gt.compute_squared_radius_3_object()(bounding_sphere);
|
||||
}
|
||||
|
||||
FT operator()(Point p) const
|
||||
{
|
||||
return func(p);
|
||||
}
|
||||
|
||||
const FT& squared_error_bound() const
|
||||
{
|
||||
return squared_error;
|
||||
}
|
||||
|
||||
const Sphere_3& bounding_sphere() const
|
||||
{
|
||||
return sphere;
|
||||
}
|
||||
|
||||
const Sphere_3& bounding_sphere_squared_radius() const
|
||||
{
|
||||
return gt.compute_squared_radius_3_object()(sphere);
|
||||
}
|
||||
|
||||
template <typename Vertex_handle>
|
||||
bool vertices_not_on_same_surface_patch(const Vertex_handle& v1,
|
||||
const Vertex_handle& v2,
|
||||
const Vertex_handle& v3) const
|
||||
{
|
||||
return
|
||||
v1->point().element_index() != v2->point().element_index() ||
|
||||
v1->point().element_index() != v3->point().element_index();
|
||||
}
|
||||
|
||||
const Function& function() const
|
||||
{
|
||||
return func;
|
||||
}
|
||||
|
||||
private:
|
||||
Function func;
|
||||
Sphere_3 sphere;
|
||||
FT squared_error;
|
||||
Geom_traits gt;
|
||||
}; // end Poisson_implicit_surface_3
|
||||
|
||||
|
||||
template <typename GT, typename Function>
|
||||
Poisson_implicit_surface_3<GT, Function>
|
||||
make_implicit_surface_3(GT, Function f,
|
||||
typename GT::Sphere_3 sphere,
|
||||
typename GT::FT error_bound)
|
||||
{
|
||||
typedef Poisson_implicit_surface_3<GT, Function> surface;
|
||||
return surface(f, sphere, error_bound);
|
||||
}
|
||||
|
||||
// template <typename GT, typename Function>
|
||||
// struct Surface_mesh_traits_generator_3<Poisson_implicit_surface_3<GT, Function> >
|
||||
// {
|
||||
// typedef Poisson_implicit_surface_3<GT, Function> Surface_type;
|
||||
// typedef typename Surface_mesher::Poisson_implicit_surface_oracle_3<GT,
|
||||
// Surface_type> Type;
|
||||
// typedef Type type; // Boost meta-programming compatibility
|
||||
// };
|
||||
|
||||
// non documented class
|
||||
template <typename FT, typename Point>
|
||||
class Poisson_implicit_function_wrapper : public std::unary_function<Point, FT>
|
||||
{
|
||||
typedef FT (*Poisson_implicit_function)(FT, FT, FT);
|
||||
|
||||
Poisson_implicit_function function;
|
||||
|
||||
public:
|
||||
Poisson_implicit_function_wrapper(Poisson_implicit_function f) : function(f) {}
|
||||
|
||||
FT operator()(Point p) const
|
||||
{
|
||||
return function(p.x(), p.y(), p.z());
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_POISSON_IMPLICIT_SURFACE_3_H
|
||||
|
|
@ -25,11 +25,39 @@
|
|||
|
||||
namespace CGAL {
|
||||
|
||||
namespace internal {
|
||||
namespace Poisson {
|
||||
|
||||
template <class Tr>
|
||||
class Constant_sizing_field {
|
||||
double sq_radius_bound;
|
||||
public:
|
||||
double cell_radius_bound() const { return CGAL::sqrt(sq_radius_bound); }
|
||||
|
||||
Constant_sizing_field(double sq_radius_bound = 0.)
|
||||
: sq_radius_bound(sq_radius_bound) {}
|
||||
|
||||
template <typename Point>
|
||||
double operator()(const Point&) const { return sq_radius_bound; }
|
||||
}; // end class Constant_sizing_field
|
||||
|
||||
} // end namespace Poisson
|
||||
} // end namespace internal
|
||||
|
||||
template <
|
||||
class Tr,
|
||||
typename Sizing_field = internal::Poisson::Constant_sizing_field<Tr>,
|
||||
typename Second_sizing_field = internal::Poisson::Constant_sizing_field<Tr>
|
||||
>
|
||||
class Poisson_mesh_cell_criteria_3
|
||||
{
|
||||
double squared_radius_bound_;
|
||||
Sizing_field sizing_field;
|
||||
Second_sizing_field second_sizing_field;
|
||||
double radius_edge_bound_;
|
||||
|
||||
typedef Poisson_mesh_cell_criteria_3<Tr,
|
||||
Sizing_field,
|
||||
Second_sizing_field> Self;
|
||||
public:
|
||||
struct Cell_quality : public std::pair<double, double>
|
||||
{
|
||||
|
|
@ -57,25 +85,34 @@ public:
|
|||
}
|
||||
};
|
||||
|
||||
inline
|
||||
double squared_radius_bound() const
|
||||
{
|
||||
return squared_radius_bound_;
|
||||
}
|
||||
// inline
|
||||
// double squared_radius_bound() const
|
||||
// {
|
||||
// return squared_radius_bound_;
|
||||
// }
|
||||
|
||||
typedef typename Tr::Cell_handle Cell_handle;
|
||||
|
||||
Poisson_mesh_cell_criteria_3(const double radius_edge_bound = 2, //< radius edge ratio bound (ignored if zero)
|
||||
const double radius_bound = 0) //< cell radius bound (ignored if zero)
|
||||
: squared_radius_bound_(radius_bound*radius_bound),
|
||||
: sizing_field(radius_bound*radius_bound),
|
||||
second_sizing_field(),
|
||||
radius_edge_bound_(radius_edge_bound)
|
||||
{}
|
||||
|
||||
inline
|
||||
void set_squared_radius_bound(const double squared_radius_bound)
|
||||
{
|
||||
squared_radius_bound_ = squared_radius_bound;
|
||||
}
|
||||
Poisson_mesh_cell_criteria_3(const double radius_edge_bound, //< radius edge ratio bound (ignored if zero)
|
||||
const Sizing_field& sizing_field,
|
||||
const Second_sizing_field second_sizing_field = Second_sizing_field())
|
||||
: sizing_field(sizing_field),
|
||||
second_sizing_field(second_sizing_field),
|
||||
radius_edge_bound_(radius_edge_bound)
|
||||
{}
|
||||
|
||||
// inline
|
||||
// void set_squared_radius_bound(const double squared_radius_bound)
|
||||
// {
|
||||
// squared_radius_bound_ = squared_radius_bound;
|
||||
// }
|
||||
|
||||
inline
|
||||
double radius_edge_bound() const
|
||||
|
|
@ -92,15 +129,12 @@ public:
|
|||
class Is_bad
|
||||
{
|
||||
protected:
|
||||
const double radius_edge_bound_;
|
||||
const double squared_radius_bound_;
|
||||
const Self* cell_criteria;
|
||||
public:
|
||||
typedef typename Tr::Point Point_3;
|
||||
|
||||
Is_bad(const double radius_edge_bound,
|
||||
const double squared_radius_bound)
|
||||
: radius_edge_bound_(radius_edge_bound),
|
||||
squared_radius_bound_(squared_radius_bound) {}
|
||||
Is_bad(const Self* cell_criteria)
|
||||
: cell_criteria(cell_criteria) {}
|
||||
|
||||
bool operator()(const Cell_handle& c,
|
||||
Cell_quality& qual) const
|
||||
|
|
@ -114,23 +148,30 @@ public:
|
|||
typedef typename Geom_traits::Compute_squared_radius_3 Radius;
|
||||
typedef typename Geom_traits::Compute_squared_distance_3 Distance;
|
||||
|
||||
Radius radius = Geom_traits().compute_squared_radius_3_object();
|
||||
Radius sq_radius = Geom_traits().compute_squared_radius_3_object();
|
||||
Distance distance = Geom_traits().compute_squared_distance_3_object();
|
||||
|
||||
double size = to_double(radius(p, q, r, s));
|
||||
double size = to_double(sq_radius(p, q, r, s));
|
||||
|
||||
if( squared_radius_bound_ != 0 )
|
||||
{
|
||||
qual.second = size / squared_radius_bound_;
|
||||
// normalized by size bound to deal
|
||||
// with size field
|
||||
// If the tetrahedron is small enough according to the second sizing
|
||||
// field, then let's say the tetrahedron is OK according to the
|
||||
// sizing fields.
|
||||
if( size > cell_criteria->second_sizing_field(p) )
|
||||
{
|
||||
// first sizing field
|
||||
const double size_bound = cell_criteria->sizing_field(p);
|
||||
if(size_bound > 0.) {
|
||||
qual.second = size / size_bound;
|
||||
// normalized by size bound to deal
|
||||
// with size field
|
||||
if( qual.sq_size() > 1 )
|
||||
{
|
||||
qual.first = 1; // (do not compute aspect)
|
||||
return true;
|
||||
}
|
||||
{
|
||||
qual.first = 1; // (do not compute aspect)
|
||||
return true;
|
||||
}
|
||||
}
|
||||
if( radius_edge_bound_ == 0 )
|
||||
}
|
||||
if( cell_criteria->radius_edge_bound_ == 0 )
|
||||
{
|
||||
qual = Cell_quality(0,1);
|
||||
return false;
|
||||
|
|
@ -145,14 +186,14 @@ public:
|
|||
|
||||
qual.first = size / min_sq_length;
|
||||
|
||||
return (qual.first > radius_edge_bound_);
|
||||
return (qual.first > cell_criteria->radius_edge_bound_);
|
||||
}
|
||||
|
||||
}; // end Is_bad
|
||||
|
||||
|
||||
Is_bad is_bad_object() const
|
||||
{ return Is_bad(radius_edge_bound_, squared_radius_bound_); }
|
||||
{ return Is_bad(this); }
|
||||
|
||||
}; // end Poisson_mesh_cell_criteria_3
|
||||
|
||||
|
|
|
|||
|
|
@ -20,6 +20,12 @@
|
|||
#ifndef CGAL_POISSON_RECONSTRUCTION_FUNCTION_H
|
||||
#define CGAL_POISSON_RECONSTRUCTION_FUNCTION_H
|
||||
|
||||
#ifndef CGAL_DIV_NORMALIZED
|
||||
# ifndef CGAL_DIV_NON_NORMALIZED
|
||||
# define CGAL_DIV_NON_NORMALIZED 1
|
||||
# endif
|
||||
#endif
|
||||
|
||||
#include <vector>
|
||||
#include <deque>
|
||||
#include <algorithm>
|
||||
|
|
@ -38,14 +44,57 @@
|
|||
#include <CGAL/centroid.h>
|
||||
#include <CGAL/property_map.h>
|
||||
#include <CGAL/surface_reconstruction_points_assertions.h>
|
||||
#include <CGAL/Memory_sizer.h>
|
||||
#include <CGAL/poisson_refine_triangulation.h>
|
||||
#include <CGAL/Robust_circumcenter_filtered_traits_3.h>
|
||||
#include <CGAL/compute_average_spacing.h>
|
||||
|
||||
#include <boost/shared_ptr.hpp>
|
||||
#include <boost/array.hpp>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace internal {
|
||||
template <class RT>
|
||||
bool
|
||||
invert(
|
||||
const RT& a0, const RT& a1, const RT& a2,
|
||||
const RT& a3, const RT& a4, const RT& a5,
|
||||
const RT& a6, const RT& a7, const RT& a8,
|
||||
RT& i0, RT& i1, RT& i2,
|
||||
RT& i3, RT& i4, RT& i5,
|
||||
RT& i6, RT& i7, RT& i8)
|
||||
{
|
||||
// Compute the adjoint.
|
||||
i0 = a4*a8 - a5*a7;
|
||||
i1 = a2*a7 - a1*a8;
|
||||
i2 = a1*a5 - a2*a4;
|
||||
i3 = a5*a6 - a3*a8;
|
||||
i4 = a0*a8 - a2*a6;
|
||||
i5 = a2*a3 - a0*a5;
|
||||
i6 = a3*a7 - a4*a6;
|
||||
i7 = a1*a6 - a0*a7;
|
||||
i8 = a0*a4 - a1*a3;
|
||||
|
||||
RT det = a0*i0 + a1*i3 + a2*i6;
|
||||
|
||||
if(det != 0) {
|
||||
RT idet = (RT(1.0))/det;
|
||||
i0 *= idet;
|
||||
i1 *= idet;
|
||||
i2 *= idet;
|
||||
i3 *= idet;
|
||||
i4 *= idet;
|
||||
i5 *= idet;
|
||||
i6 *= idet;
|
||||
i7 *= idet;
|
||||
i8 *= idet;
|
||||
return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/// Given a set of 3D points with oriented normals sampled on the boundary of a 3D solid,
|
||||
/// the Poisson Surface Reconstruction method [Kazhdan06] solves for an approximate indicator function
|
||||
|
|
@ -62,6 +111,45 @@ namespace CGAL {
|
|||
/// @heading Parameters:
|
||||
/// @param Gt Geometric traits class.
|
||||
|
||||
|
||||
struct Poisson_visitor {
|
||||
void before_insertion() const
|
||||
{}
|
||||
};
|
||||
|
||||
struct Poisson_skip_vertices {
|
||||
double ratio;
|
||||
Random& m_random;
|
||||
Poisson_skip_vertices(const double ratio, Random& random)
|
||||
: ratio(ratio), m_random(random) {}
|
||||
|
||||
template <typename Iterator>
|
||||
bool operator()(Iterator) const {
|
||||
return m_random.get_double() < ratio;
|
||||
}
|
||||
};
|
||||
|
||||
// Given f1 and f2, two sizing fields, that functor wrapper returns
|
||||
// max(f1, f2*f2)
|
||||
// The wrapper stores only pointers to the two functors.
|
||||
template <typename F1, typename F2>
|
||||
struct Special_wrapper_of_two_functions_keep_pointers {
|
||||
F1 *f1;
|
||||
F2 *f2;
|
||||
Special_wrapper_of_two_functions_keep_pointers(F1* f1, F2* f2)
|
||||
: f1(f1), f2(f2) {}
|
||||
|
||||
template <typename X>
|
||||
double operator()(const X& x) const {
|
||||
return (std::max)((*f1)(x), CGAL::square((*f2)(x)));
|
||||
}
|
||||
|
||||
template <typename X>
|
||||
double operator()(const X& x) {
|
||||
return (std::max)((*f1)(x), CGAL::square((*f2)(x)));
|
||||
}
|
||||
}; // end struct Special_wrapper_of_two_functions_keep_pointers<F1, F2>
|
||||
|
||||
template <class Gt>
|
||||
class Poisson_reconstruction_function
|
||||
{
|
||||
|
|
@ -69,6 +157,9 @@ class Poisson_reconstruction_function
|
|||
public:
|
||||
|
||||
typedef Gt Geom_traits; ///< Geometric traits class
|
||||
typedef Reconstruction_triangulation_3<Robust_circumcenter_filtered_traits_3<Gt> >
|
||||
Triangulation;
|
||||
typedef typename Triangulation::Cell_handle Cell_handle;
|
||||
|
||||
// Geometric types
|
||||
typedef typename Geom_traits::FT FT; ///< typedef to Geom_traits::FT
|
||||
|
|
@ -81,8 +172,6 @@ private:
|
|||
|
||||
/// Internal 3D triangulation, of type Reconstruction_triangulation_3.
|
||||
// Note: poisson_refine_triangulation() requires a robust circumcenter computation.
|
||||
typedef Reconstruction_triangulation_3<Robust_circumcenter_filtered_traits_3<Gt> >
|
||||
Triangulation;
|
||||
|
||||
// Repeat Triangulation types
|
||||
typedef typename Triangulation::Triangulation_data_structure Triangulation_data_structure;
|
||||
|
|
@ -91,7 +180,6 @@ private:
|
|||
typedef typename Geom_traits::Segment_3 Segment;
|
||||
typedef typename Geom_traits::Triangle_3 Triangle;
|
||||
typedef typename Geom_traits::Tetrahedron_3 Tetrahedron;
|
||||
typedef typename Triangulation::Cell_handle Cell_handle;
|
||||
typedef typename Triangulation::Vertex_handle Vertex_handle;
|
||||
typedef typename Triangulation::Cell Cell;
|
||||
typedef typename Triangulation::Vertex Vertex;
|
||||
|
|
@ -120,10 +208,17 @@ private:
|
|||
// the Poisson equation Laplacian(f) = divergent(normals field).
|
||||
boost::shared_ptr<Triangulation> m_tr;
|
||||
|
||||
mutable boost::shared_ptr<std::vector<boost::array<double,9> > > m_Bary;
|
||||
mutable std::vector<Point> Dual;
|
||||
mutable std::vector<Vector> Normal;
|
||||
|
||||
// contouring and meshing
|
||||
Point m_sink; // Point with the minimum value of operator()
|
||||
mutable Cell_handle m_hint; // last cell found = hint for next search
|
||||
|
||||
FT average_spacing;
|
||||
|
||||
|
||||
// Public methods
|
||||
public:
|
||||
|
||||
|
|
@ -138,14 +233,17 @@ public:
|
|||
// This variant requires all parameters.
|
||||
template <typename InputIterator,
|
||||
typename PointPMap,
|
||||
typename NormalPMap
|
||||
typename NormalPMap,
|
||||
typename Visitor
|
||||
>
|
||||
Poisson_reconstruction_function(
|
||||
InputIterator first, ///< iterator over the first input point.
|
||||
InputIterator beyond, ///< past-the-end iterator over the input points.
|
||||
PointPMap point_pmap, ///< property map to access the position of an input point.
|
||||
NormalPMap normal_pmap) ///< property map to access the *oriented* normal of an input point.
|
||||
: m_tr(new Triangulation)
|
||||
NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point.
|
||||
Visitor visitor)
|
||||
: m_tr(new Triangulation), m_Bary(new std::vector<boost::array<double,9> > )
|
||||
, average_spacing(CGAL::compute_average_spacing(first, beyond, 6))
|
||||
{
|
||||
CGAL::Timer task_timer; task_timer.start();
|
||||
CGAL_TRACE_STREAM << "Creates Poisson triangulation...\n";
|
||||
|
|
@ -154,24 +252,27 @@ public:
|
|||
m_tr->insert(
|
||||
first,beyond,
|
||||
point_pmap,
|
||||
normal_pmap);
|
||||
normal_pmap,
|
||||
visitor);
|
||||
|
||||
// Prints status
|
||||
CGAL_TRACE_STREAM << "Creates Poisson triangulation: " << task_timer.time() << " seconds, "
|
||||
<< (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated"
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
/// @cond SKIP_IN_MANUAL
|
||||
// This variant creates a default point property map = Dereference_property_map.
|
||||
template <typename InputIterator,
|
||||
typename NormalPMap
|
||||
typename NormalPMap,
|
||||
typename Visitor
|
||||
>
|
||||
Poisson_reconstruction_function(
|
||||
InputIterator first, ///< iterator over the first input point.
|
||||
InputIterator beyond, ///< past-the-end iterator over the input points.
|
||||
NormalPMap normal_pmap) ///< property map to access the *oriented* normal of an input point.
|
||||
: m_tr(new Triangulation)
|
||||
NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point.
|
||||
Visitor visitor)
|
||||
: m_tr(new Triangulation), m_Bary(new std::vector<boost::array<double,9> > )
|
||||
, average_spacing(CGAL::compute_average_spacing(first, beyond, 6))
|
||||
{
|
||||
CGAL::Timer task_timer; task_timer.start();
|
||||
CGAL_TRACE_STREAM << "Creates Poisson triangulation...\n";
|
||||
|
|
@ -179,11 +280,11 @@ public:
|
|||
// Inserts points in triangulation
|
||||
m_tr->insert(
|
||||
first,beyond,
|
||||
normal_pmap);
|
||||
normal_pmap,
|
||||
visitor);
|
||||
|
||||
// Prints status
|
||||
CGAL_TRACE_STREAM << "Creates Poisson triangulation: " << task_timer.time() << " seconds, "
|
||||
<< (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated"
|
||||
<< std::endl;
|
||||
}
|
||||
/// @endcond
|
||||
|
|
@ -191,7 +292,11 @@ public:
|
|||
/// Returns a sphere bounding the inferred surface.
|
||||
Sphere bounding_sphere() const
|
||||
{
|
||||
return m_tr->input_points_bounding_sphere();
|
||||
return m_tr->bounding_sphere();
|
||||
}
|
||||
|
||||
const Triangulation& tr() const {
|
||||
return *m_tr;
|
||||
}
|
||||
|
||||
/// The function compute_implicit_function() must be called
|
||||
|
|
@ -208,9 +313,13 @@ public:
|
|||
/// @return false if the linear solver fails.
|
||||
|
||||
// This variant requires all parameters.
|
||||
template <class SparseLinearAlgebraTraits_d>
|
||||
template <class SparseLinearAlgebraTraits_d,
|
||||
class Visitor>
|
||||
bool compute_implicit_function(
|
||||
SparseLinearAlgebraTraits_d solver = SparseLinearAlgebraTraits_d()) ///< sparse linear solver
|
||||
SparseLinearAlgebraTraits_d solver,// = SparseLinearAlgebraTraits_d(),
|
||||
Visitor visitor,
|
||||
double approximation_ratio = 0,
|
||||
double average_spacing_ratio = 5)
|
||||
{
|
||||
CGAL::Timer task_timer; task_timer.start();
|
||||
CGAL_TRACE_STREAM << "Delaunay refinement...\n";
|
||||
|
|
@ -221,16 +330,102 @@ public:
|
|||
const FT enlarge_ratio = 1.5;
|
||||
const FT radius = sqrt(bounding_sphere().squared_radius()); // get triangulation's radius
|
||||
const FT cell_radius_bound = radius/5.; // large
|
||||
unsigned int nb_vertices_added = delaunay_refinement(radius_edge_ratio_bound,cell_radius_bound,max_vertices,enlarge_ratio);
|
||||
|
||||
internal::Poisson::Constant_sizing_field<Triangulation> sizing_field(CGAL::square(cell_radius_bound));
|
||||
|
||||
std::vector<int> NB;
|
||||
|
||||
NB.push_back( delaunay_refinement(radius_edge_ratio_bound,sizing_field,max_vertices,enlarge_ratio));
|
||||
|
||||
while(m_tr->insert_fraction(visitor)){
|
||||
|
||||
NB.push_back( delaunay_refinement(radius_edge_ratio_bound,sizing_field,max_vertices,enlarge_ratio));
|
||||
}
|
||||
|
||||
if(approximation_ratio > 0. &&
|
||||
approximation_ratio * std::distance(m_tr->input_points_begin(),
|
||||
m_tr->input_points_end()) > 20) {
|
||||
|
||||
// Add a pass of Delaunay refinement.
|
||||
//
|
||||
// In that pass, the sizing field, of the refinement process of the
|
||||
// triangulation, is based on the result of a poisson function with a
|
||||
// sample of the input points. The ratio is 'approximation_ratio'.
|
||||
//
|
||||
// For optimization reasons, the cell criteria of the refinement
|
||||
// process uses two sizing fields:
|
||||
//
|
||||
// - the minimum of the square of 'coarse_poisson_function' and the
|
||||
// square of the constant field equal to 'average_spacing',
|
||||
//
|
||||
// - a second sizing field that is constant, and equal to:
|
||||
//
|
||||
// average_spacing*average_spacing_ratio
|
||||
//
|
||||
// If a given cell is smaller than the constant second sizing field,
|
||||
// then the cell is considered as small enough, and the first sizing
|
||||
// field, more costly, is not evaluated.
|
||||
|
||||
typedef Filter_iterator<typename Triangulation::Input_point_iterator,
|
||||
Poisson_skip_vertices> Some_points_iterator;
|
||||
//make it deterministic
|
||||
Random random(0);
|
||||
Poisson_skip_vertices skip(1.-approximation_ratio,random);
|
||||
|
||||
CGAL_TRACE_STREAM << "SPECIAL PASS that uses an approximation of the result (approximation ratio: "
|
||||
<< approximation_ratio << ")" << std::endl;
|
||||
CGAL::Timer approximation_timer; approximation_timer.start();
|
||||
|
||||
CGAL::Timer sizing_field_timer; sizing_field_timer.start();
|
||||
Poisson_reconstruction_function<Geom_traits>
|
||||
coarse_poisson_function(Some_points_iterator(m_tr->input_points_end(),
|
||||
skip,
|
||||
m_tr->input_points_begin()),
|
||||
Some_points_iterator(m_tr->input_points_end(),
|
||||
skip),
|
||||
Normal_of_point_with_normal_pmap<Geom_traits>(),
|
||||
Poisson_visitor());
|
||||
coarse_poisson_function.compute_implicit_function(solver, Poisson_visitor(),
|
||||
0.);
|
||||
internal::Poisson::Constant_sizing_field<Triangulation>
|
||||
min_sizing_field(CGAL::square(average_spacing));
|
||||
internal::Poisson::Constant_sizing_field<Triangulation>
|
||||
sizing_field_ok(CGAL::square(average_spacing*average_spacing_ratio));
|
||||
|
||||
Special_wrapper_of_two_functions_keep_pointers<
|
||||
internal::Poisson::Constant_sizing_field<Triangulation>,
|
||||
Poisson_reconstruction_function<Geom_traits> > sizing_field2(&min_sizing_field,
|
||||
&coarse_poisson_function);
|
||||
|
||||
sizing_field_timer.stop();
|
||||
std::cerr << "Construction time of the sizing field: " << sizing_field_timer.time()
|
||||
<< " seconds" << std::endl;
|
||||
|
||||
NB.push_back( delaunay_refinement(radius_edge_ratio_bound,
|
||||
sizing_field2,
|
||||
max_vertices,
|
||||
enlarge_ratio,
|
||||
sizing_field_ok) );
|
||||
approximation_timer.stop();
|
||||
CGAL_TRACE_STREAM << "SPECIAL PASS END (" << approximation_timer.time() << " seconds)" << std::endl;
|
||||
}
|
||||
|
||||
|
||||
// Prints status
|
||||
CGAL_TRACE_STREAM << "Delaunay refinement: " << "added " << nb_vertices_added << " Steiner points, "
|
||||
<< task_timer.time() << " seconds, "
|
||||
<< (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated"
|
||||
<< std::endl;
|
||||
CGAL_TRACE_STREAM << "Delaunay refinement: " << "added ";
|
||||
for(std::size_t i = 0; i < NB.size()-1; i++){
|
||||
CGAL_TRACE_STREAM << NB[i] << " + ";
|
||||
}
|
||||
CGAL_TRACE_STREAM << NB.back() << " Steiner points, "
|
||||
<< task_timer.time() << " seconds, "
|
||||
<< std::endl;
|
||||
task_timer.reset();
|
||||
|
||||
#ifdef CGAL_DIV_NON_NORMALIZED
|
||||
CGAL_TRACE_STREAM << "Solve Poisson equation with non-normalized divergence...\n";
|
||||
#else
|
||||
CGAL_TRACE_STREAM << "Solve Poisson equation with normalized divergence...\n";
|
||||
#endif
|
||||
|
||||
// Computes the Poisson indicator function operator()
|
||||
// at each vertex of the triangulation.
|
||||
|
|
@ -248,40 +443,68 @@ public:
|
|||
|
||||
// Prints status
|
||||
CGAL_TRACE_STREAM << "Solve Poisson equation: " << task_timer.time() << " seconds, "
|
||||
<< (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated"
|
||||
<< std::endl;
|
||||
task_timer.reset();
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class SparseLinearAlgebraTraits_d>
|
||||
bool compute_implicit_function(SparseLinearAlgebraTraits_d solver, bool smoother_hole_filling = false)
|
||||
{
|
||||
if (smoother_hole_filling)
|
||||
return compute_implicit_function<SparseLinearAlgebraTraits_d,Poisson_visitor>(solver,Poisson_visitor(),0.02,5);
|
||||
else
|
||||
return compute_implicit_function<SparseLinearAlgebraTraits_d,Poisson_visitor>(solver,Poisson_visitor());
|
||||
}
|
||||
|
||||
#ifdef CGAL_EIGEN3_ENABLED
|
||||
#ifdef CGAL_EIGEN3_ENABLED
|
||||
/// @cond SKIP_IN_MANUAL
|
||||
// This variant provides the default sparse linear traits class = Eigen_solver_traits.
|
||||
bool compute_implicit_function()
|
||||
bool compute_implicit_function(bool smoother_hole_filling = false)
|
||||
{
|
||||
return compute_implicit_function<
|
||||
Eigen_solver_traits<Eigen::ConjugateGradient<Eigen_sparse_symmetric_matrix<double>::EigenType> >
|
||||
>();
|
||||
typedef Eigen_solver_traits<Eigen::ConjugateGradient<Eigen_sparse_symmetric_matrix<double>::EigenType> > Solver;
|
||||
return compute_implicit_function<Solver>(Solver(), smoother_hole_filling);
|
||||
}
|
||||
/// @endcond
|
||||
#else
|
||||
#else
|
||||
#ifdef CGAL_TAUCS_ENABLED
|
||||
/// @cond SKIP_IN_MANUAL
|
||||
// This variant provides the default sparse linear traits class = Taucs_symmetric_solver_traits.
|
||||
bool compute_implicit_function()
|
||||
bool compute_implicit_function(bool smoother_hole_filling = false)
|
||||
{
|
||||
return compute_implicit_function< Taucs_symmetric_solver_traits<double> >();
|
||||
typedef Taucs_symmetric_solver_traits<double> Solver;
|
||||
return compute_implicit_function<Solver>(Solver(), smoother_hole_filling);
|
||||
}
|
||||
/// @endcond
|
||||
#endif
|
||||
#endif
|
||||
#endif
|
||||
|
||||
boost::tuple<FT, Cell_handle, bool> special_func(const Point& p) const
|
||||
{
|
||||
m_hint = m_tr->locate(p ,m_hint ); // no hint when we use hierarchy
|
||||
|
||||
if(m_tr->is_infinite(m_hint)) {
|
||||
int i = m_hint->index(m_tr->infinite_vertex());
|
||||
return boost::make_tuple(m_hint->vertex((i+1)&3)->f(),
|
||||
m_hint, true);
|
||||
}
|
||||
|
||||
FT a,b,c,d;
|
||||
barycentric_coordinates(p,m_hint,a,b,c,d);
|
||||
return boost::make_tuple(a * m_hint->vertex(0)->f() +
|
||||
b * m_hint->vertex(1)->f() +
|
||||
c * m_hint->vertex(2)->f() +
|
||||
d * m_hint->vertex(3)->f(),
|
||||
m_hint, false);
|
||||
}
|
||||
|
||||
>>>>>>> Surface_reconstruction_points_3-make_it_faster-afabri
|
||||
|
||||
/// 'ImplicitFunction' interface: evaluates the implicit function at a given 3D query point.
|
||||
FT operator()(const Point& p) const
|
||||
{
|
||||
m_hint = m_tr->locate(p,m_hint);
|
||||
m_hint = m_tr->locate(p ,m_hint);
|
||||
|
||||
if(m_tr->is_infinite(m_hint)) {
|
||||
int i = m_hint->index(m_tr->infinite_vertex());
|
||||
|
|
@ -296,6 +519,84 @@ public:
|
|||
d * m_hint->vertex(3)->f();
|
||||
}
|
||||
|
||||
void initialize_cell_indices()
|
||||
{
|
||||
int i=0;
|
||||
for(Finite_cells_iterator fcit = m_tr->finite_cells_begin();
|
||||
fcit != m_tr->finite_cells_end();
|
||||
++fcit){
|
||||
fcit->info()= i++;
|
||||
}
|
||||
}
|
||||
|
||||
void initialize_barycenters() const
|
||||
{
|
||||
m_Bary->resize(m_tr->number_of_cells());
|
||||
|
||||
for(std::size_t i=0; i< m_Bary->size();i++){
|
||||
(*m_Bary)[i][0]=-1;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
void initialize_cell_normals() const
|
||||
{
|
||||
Normal.resize(m_tr->number_of_cells());
|
||||
int i = 0;
|
||||
int N = 0;
|
||||
for(Finite_cells_iterator fcit = m_tr->finite_cells_begin();
|
||||
fcit != m_tr->finite_cells_end();
|
||||
++fcit){
|
||||
Normal[i] = cell_normal(fcit);
|
||||
if(Normal[i] == NULL_VECTOR){
|
||||
N++;
|
||||
}
|
||||
++i;
|
||||
}
|
||||
std::cerr << N << " out of " << i << " cells have NULL_VECTOR as normal" << std::endl;
|
||||
}
|
||||
|
||||
void initialize_duals() const
|
||||
{
|
||||
Dual.resize(m_tr->number_of_cells());
|
||||
int i = 0;
|
||||
for(Finite_cells_iterator fcit = m_tr->finite_cells_begin();
|
||||
fcit != m_tr->finite_cells_end();
|
||||
++fcit){
|
||||
Dual[i++] = m_tr->dual(fcit);
|
||||
}
|
||||
}
|
||||
|
||||
void clear_duals() const
|
||||
{
|
||||
Dual.clear();
|
||||
}
|
||||
|
||||
void clear_normals() const
|
||||
{
|
||||
Normal.clear();
|
||||
}
|
||||
|
||||
void initialize_matrix_entry(Cell_handle ch) const
|
||||
{
|
||||
boost::array<double,9> & entry = (*m_Bary)[ch->info()];
|
||||
const Point& pa = ch->vertex(0)->point();
|
||||
const Point& pb = ch->vertex(1)->point();
|
||||
const Point& pc = ch->vertex(2)->point();
|
||||
const Point& pd = ch->vertex(3)->point();
|
||||
|
||||
Vector va = pa - pd;
|
||||
Vector vb = pb - pd;
|
||||
Vector vc = pc - pd;
|
||||
|
||||
internal::invert(va.x(), va.y(), va.z(),
|
||||
vb.x(), vb.y(), vb.z(),
|
||||
vc.x(), vc.y(), vc.z(),
|
||||
entry[0],entry[1],entry[2],entry[3],entry[4],entry[5],entry[6],entry[7],entry[8]);
|
||||
}
|
||||
|
||||
/// Returns a point located inside the inferred surface.
|
||||
Point get_inner_point() const
|
||||
{
|
||||
|
|
@ -310,18 +611,30 @@ private:
|
|||
/// bad means badly shaped or too big). The normal of
|
||||
/// Steiner points is set to zero.
|
||||
/// Returns the number of vertices inserted.
|
||||
|
||||
template <typename Sizing_field>
|
||||
unsigned int delaunay_refinement(FT radius_edge_ratio_bound, ///< radius edge ratio bound (ignored if zero)
|
||||
FT cell_radius_bound, ///< cell radius bound (ignored if zero)
|
||||
Sizing_field sizing_field, ///< cell radius bound (ignored if zero)
|
||||
unsigned int max_vertices, ///< number of vertices bound
|
||||
FT enlarge_ratio) ///< bounding box enlarge ratio
|
||||
{
|
||||
CGAL_TRACE("Calls delaunay_refinement(radius_edge_ratio_bound=%lf, cell_radius_bound=%lf, max_vertices=%u, enlarge_ratio=%lf)\n",
|
||||
radius_edge_ratio_bound, cell_radius_bound, max_vertices, enlarge_ratio);
|
||||
return delaunay_refinement(radius_edge_ratio_bound,
|
||||
sizing_field,
|
||||
max_vertices,
|
||||
enlarge_ratio,
|
||||
internal::Poisson::Constant_sizing_field<Triangulation>());
|
||||
}
|
||||
|
||||
template <typename Sizing_field,
|
||||
typename Second_sizing_field>
|
||||
unsigned int delaunay_refinement(FT radius_edge_ratio_bound, ///< radius edge ratio bound (ignored if zero)
|
||||
Sizing_field sizing_field, ///< cell radius bound (ignored if zero)
|
||||
unsigned int max_vertices, ///< number of vertices bound
|
||||
FT enlarge_ratio, ///< bounding box enlarge ratio
|
||||
Second_sizing_field second_sizing_field)
|
||||
{
|
||||
Sphere elarged_bsphere = enlarged_bounding_sphere(enlarge_ratio);
|
||||
unsigned int nb_vertices_added = poisson_refine_triangulation(*m_tr,radius_edge_ratio_bound,cell_radius_bound,max_vertices,elarged_bsphere);
|
||||
|
||||
CGAL_TRACE("End of delaunay_refinement()\n");
|
||||
unsigned int nb_vertices_added = poisson_refine_triangulation(*m_tr,radius_edge_ratio_bound,sizing_field,second_sizing_field,max_vertices,elarged_bsphere);
|
||||
|
||||
return nb_vertices_added;
|
||||
}
|
||||
|
|
@ -343,39 +656,46 @@ private:
|
|||
double duration_assembly = 0.0;
|
||||
double duration_solve = 0.0;
|
||||
|
||||
CGAL_TRACE(" Creates matrix...\n");
|
||||
|
||||
initialize_cell_indices();
|
||||
initialize_barycenters();
|
||||
|
||||
// get #variables
|
||||
unsigned int nb_variables = m_tr->index_unconstrained_vertices();
|
||||
constrain_one_vertex_on_convex_hull();
|
||||
m_tr->index_unconstrained_vertices();
|
||||
unsigned int nb_variables = m_tr->number_of_vertices()-1;
|
||||
|
||||
// at least one vertex must be constrained
|
||||
if(nb_variables == m_tr->number_of_vertices())
|
||||
{
|
||||
constrain_one_vertex_on_convex_hull();
|
||||
nb_variables = m_tr->index_unconstrained_vertices();
|
||||
}
|
||||
CGAL_TRACE(" Number of variables: %ld\n", (long)(nb_variables));
|
||||
|
||||
// Assemble linear system A*X=B
|
||||
typename SparseLinearAlgebraTraits_d::Matrix A(nb_variables); // matrix is symmetric definite positive
|
||||
typename SparseLinearAlgebraTraits_d::Vector X(nb_variables), B(nb_variables);
|
||||
|
||||
initialize_duals();
|
||||
#ifndef CGAL_DIV_NON_NORMALIZED
|
||||
initialize_cell_normals();
|
||||
#endif
|
||||
Finite_vertices_iterator v, e;
|
||||
for(v = m_tr->finite_vertices_begin(),
|
||||
e = m_tr->finite_vertices_end();
|
||||
v != e;
|
||||
++v)
|
||||
{
|
||||
if(!v->constrained())
|
||||
{
|
||||
if(!m_tr->is_constrained(v)) {
|
||||
#ifdef CGAL_DIV_NON_NORMALIZED
|
||||
B[v->index()] = div(v); // rhs -> divergent
|
||||
#else // not defined(CGAL_DIV_NORMALIZED)
|
||||
B[v->index()] = div_normalized(v); // rhs -> divergent
|
||||
#endif // not defined(CGAL_DIV_NORMALIZED)
|
||||
assemble_poisson_row<SparseLinearAlgebraTraits_d>(A,v,B,lambda);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
clear_duals();
|
||||
clear_normals();
|
||||
duration_assembly = (clock() - time_init)/CLOCKS_PER_SEC;
|
||||
CGAL_TRACE(" Creates matrix: done (%.2lf s)\n", duration_assembly);
|
||||
|
||||
CGAL_TRACE(" %ld Mb allocated\n", long(CGAL::Memory_sizer().virtual_size()>>20));
|
||||
CGAL_TRACE(" Solve sparse linear system...\n");
|
||||
|
||||
// Solve "A*X = B". On success, solution is (1/D) * X.
|
||||
|
|
@ -391,10 +711,9 @@ private:
|
|||
// copy function's values to vertices
|
||||
unsigned int index = 0;
|
||||
for (v = m_tr->finite_vertices_begin(), e = m_tr->finite_vertices_end(); v!= e; ++v)
|
||||
if(!v->constrained())
|
||||
if(!m_tr->is_constrained(v))
|
||||
v->f() = X[index++];
|
||||
|
||||
CGAL_TRACE(" %ld Mb allocated\n", long(CGAL::Memory_sizer().virtual_size()>>20));
|
||||
CGAL_TRACE("End of solve_poisson()\n");
|
||||
|
||||
return true;
|
||||
|
|
@ -453,16 +772,43 @@ private:
|
|||
FT& c,
|
||||
FT& d) const
|
||||
{
|
||||
const Point& pa = cell->vertex(0)->point();
|
||||
const Point& pb = cell->vertex(1)->point();
|
||||
const Point& pc = cell->vertex(2)->point();
|
||||
const Point& pd = cell->vertex(3)->point();
|
||||
|
||||
// const Point& pa = cell->vertex(0)->point();
|
||||
// const Point& pb = cell->vertex(1)->point();
|
||||
// const Point& pc = cell->vertex(2)->point();
|
||||
const Point& pd = cell->vertex(3)->point();
|
||||
#if 1
|
||||
//Vector va = pa - pd;
|
||||
//Vector vb = pb - pd;
|
||||
//Vector vc = pc - pd;
|
||||
Vector vp = p - pd;
|
||||
|
||||
//FT i00, i01, i02, i10, i11, i12, i20, i21, i22;
|
||||
//internal::invert(va.x(), va.y(), va.z(),
|
||||
// vb.x(), vb.y(), vb.z(),
|
||||
// vc.x(), vc.y(), vc.z(),
|
||||
// i00, i01, i02, i10, i11, i12, i20, i21, i22);
|
||||
const boost::array<double,9> & i = (*m_Bary)[cell->info()];
|
||||
if(i[0]==-1){
|
||||
initialize_matrix_entry(cell);
|
||||
}
|
||||
// UsedBary[cell->info()] = true;
|
||||
a = i[0] * vp.x() + i[3] * vp.y() + i[6] * vp.z();
|
||||
b = i[1] * vp.x() + i[4] * vp.y() + i[7] * vp.z();
|
||||
c = i[2] * vp.x() + i[5] * vp.y() + i[8] * vp.z();
|
||||
d = 1 - ( a + b + c);
|
||||
#else
|
||||
FT v = volume(pa,pb,pc,pd);
|
||||
a = CGAL::abs(volume(pb,pc,pd,p) / v);
|
||||
b = CGAL::abs(volume(pa,pc,pd,p) / v);
|
||||
c = CGAL::abs(volume(pb,pa,pd,p) / v);
|
||||
d = CGAL::abs(volume(pb,pc,pa,p) / v);
|
||||
a = std::fabs(volume(pb,pc,pd,p) / v);
|
||||
b = std::fabs(volume(pa,pc,pd,p) / v);
|
||||
c = std::fabs(volume(pb,pa,pd,p) / v);
|
||||
d = std::fabs(volume(pb,pc,pa,p) / v);
|
||||
|
||||
std::cerr << "_________________________________\n";
|
||||
std::cerr << aa << " " << bb << " " << cc << " " << dd << std::endl;
|
||||
std::cerr << a << " " << b << " " << c << " " << d << std::endl;
|
||||
|
||||
#endif
|
||||
}
|
||||
|
||||
FT find_sink()
|
||||
|
|
@ -506,31 +852,29 @@ private:
|
|||
|
||||
Vertex_handle any_vertex_on_convex_hull()
|
||||
{
|
||||
// TODO: return NULL if none and assert
|
||||
std::vector<Vertex_handle> vertices;
|
||||
vertices.reserve(32);
|
||||
m_tr->adjacent_vertices(m_tr->infinite_vertex(),std::back_inserter(vertices));
|
||||
typename std::vector<Vertex_handle>::iterator it = vertices.begin();
|
||||
return *it;
|
||||
Cell_handle ch = m_tr->infinite_vertex()->cell();
|
||||
return ch->vertex( (ch->index( m_tr->infinite_vertex())+1)%4);
|
||||
}
|
||||
|
||||
|
||||
void constrain_one_vertex_on_convex_hull(const FT value = 0.0)
|
||||
{
|
||||
Vertex_handle v = any_vertex_on_convex_hull();
|
||||
v->constrained() = true;
|
||||
m_tr->constrain(v);
|
||||
v->f() = value;
|
||||
}
|
||||
|
||||
// divergent
|
||||
FT div(Vertex_handle v)
|
||||
// TODO: Some entities are computed too often
|
||||
// - nn and area should not be computed for the face and its opposite face
|
||||
//
|
||||
// divergent
|
||||
FT div_normalized(Vertex_handle v)
|
||||
{
|
||||
std::vector<Cell_handle> cells;
|
||||
cells.reserve(32);
|
||||
m_tr->incident_cells(v,std::back_inserter(cells));
|
||||
if(cells.size() == 0)
|
||||
return 0.0;
|
||||
|
||||
FT div = 0.0;
|
||||
|
||||
FT div = 0;
|
||||
typename std::vector<Cell_handle>::iterator it;
|
||||
for(it = cells.begin(); it != cells.end(); it++)
|
||||
{
|
||||
|
|
@ -539,35 +883,84 @@ private:
|
|||
continue;
|
||||
|
||||
// compute average normal per cell
|
||||
Vector n = cell_normal(cell);
|
||||
Vector n = get_cell_normal(cell);
|
||||
|
||||
// zero normal - no need to compute anything else
|
||||
if(n == CGAL::NULL_VECTOR)
|
||||
continue;
|
||||
|
||||
|
||||
// compute n'
|
||||
int index = cell->index(v);
|
||||
const Point& x = cell->vertex(index)->point();
|
||||
const Point& a = cell->vertex((index+1)%4)->point();
|
||||
const Point& b = cell->vertex((index+2)%4)->point();
|
||||
const Point& c = cell->vertex((index+3)%4)->point();
|
||||
Vector nn = (index%2==0) ? CGAL::cross_product(b-a,c-a) : CGAL::cross_product(c-a,b-a);
|
||||
div += n * nn;
|
||||
nn = nn / std::sqrt(nn*nn); // normalize
|
||||
Vector p = a - x;
|
||||
Vector q = b - x;
|
||||
Vector r = c - x;
|
||||
FT p_n = std::sqrt(p*p);
|
||||
FT q_n = std::sqrt(q*q);
|
||||
FT r_n = std::sqrt(r*r);
|
||||
FT solid_angle = p*(CGAL::cross_product(q,r));
|
||||
solid_angle = std::abs(solid_angle / (p_n*q_n*r_n + (p*q)*r_n + (q*r)*p_n + (r*p)*q_n));
|
||||
|
||||
FT area = std::sqrt(squared_area(a,b,c));
|
||||
FT length = p_n + q_n + r_n;
|
||||
div += n * nn * area / length ;
|
||||
}
|
||||
return div * FT(3.0);
|
||||
}
|
||||
|
||||
FT div(Vertex_handle v)
|
||||
{
|
||||
std::vector<Cell_handle> cells;
|
||||
cells.reserve(32);
|
||||
m_tr->incident_cells(v,std::back_inserter(cells));
|
||||
|
||||
FT div = 0.0;
|
||||
typename std::vector<Cell_handle>::iterator it;
|
||||
for(it = cells.begin(); it != cells.end(); it++)
|
||||
{
|
||||
Cell_handle cell = *it;
|
||||
if(m_tr->is_infinite(cell))
|
||||
continue;
|
||||
|
||||
const int index = cell->index(v);
|
||||
const Point& a = cell->vertex(m_tr->vertex_triple_index(index, 0))->point();
|
||||
const Point& b = cell->vertex(m_tr->vertex_triple_index(index, 1))->point();
|
||||
const Point& c = cell->vertex(m_tr->vertex_triple_index(index, 2))->point();
|
||||
const Vector nn = CGAL::cross_product(b-a,c-a);
|
||||
|
||||
div+= nn * (//v->normal() +
|
||||
cell->vertex((index+1)%4)->normal() +
|
||||
cell->vertex((index+2)%4)->normal() +
|
||||
cell->vertex((index+3)%4)->normal());
|
||||
}
|
||||
return div;
|
||||
}
|
||||
|
||||
Vector cell_normal(Cell_handle cell)
|
||||
Vector get_cell_normal(Cell_handle cell)
|
||||
{
|
||||
return Normal[cell->info()];
|
||||
}
|
||||
|
||||
Vector cell_normal(Cell_handle cell) const
|
||||
{
|
||||
const Vector& n0 = cell->vertex(0)->normal();
|
||||
const Vector& n1 = cell->vertex(1)->normal();
|
||||
const Vector& n2 = cell->vertex(2)->normal();
|
||||
const Vector& n3 = cell->vertex(3)->normal();
|
||||
Vector n = n0 + n1 + n2 + n3;
|
||||
FT sq_norm = n*n;
|
||||
if(sq_norm != 0.0)
|
||||
return n / std::sqrt(sq_norm); // normalize
|
||||
else
|
||||
return CGAL::NULL_VECTOR;
|
||||
if(n != NULL_VECTOR){
|
||||
FT sq_norm = n*n;
|
||||
if(sq_norm != 0.0){
|
||||
return n / std::sqrt(sq_norm); // normalize
|
||||
}
|
||||
}
|
||||
return NULL_VECTOR;
|
||||
}
|
||||
|
||||
// cotan formula as area(voronoi face) / len(primal edge)
|
||||
|
|
@ -593,11 +986,12 @@ private:
|
|||
Cell_circulator circ = m_tr->incident_cells(edge);
|
||||
Cell_circulator done = circ;
|
||||
std::vector<Point> voronoi_points;
|
||||
voronoi_points.reserve(9);
|
||||
do
|
||||
{
|
||||
Cell_handle cell = circ;
|
||||
if(!m_tr->is_infinite(cell))
|
||||
voronoi_points.push_back(m_tr->dual(cell));
|
||||
voronoi_points.push_back(Dual[cell->info()]);
|
||||
else // one infinite tet, switch to another calculation
|
||||
return area_voronoi_face_boundary(edge);
|
||||
circ++;
|
||||
|
|
@ -618,8 +1012,7 @@ private:
|
|||
{
|
||||
const Point& b = voronoi_points[i];
|
||||
const Point& c = voronoi_points[i+1];
|
||||
Triangle triangle(a,b,c);
|
||||
area += std::sqrt(triangle.squared_area());
|
||||
area += std::sqrt(squared_area(a,b,c));
|
||||
}
|
||||
return area;
|
||||
}
|
||||
|
|
@ -644,13 +1037,14 @@ private:
|
|||
if(!m_tr->is_infinite(cell))
|
||||
{
|
||||
// circumcenter of cell
|
||||
Point c = m_tr->dual(cell);
|
||||
Point c = Dual[cell->info()];
|
||||
Tetrahedron tet = m_tr->tetrahedron(cell);
|
||||
|
||||
int i = cell->index(vi);
|
||||
int j = cell->index(vj);
|
||||
int k = -1, l = -1;
|
||||
other_two_indices(i,j, &k,&l);
|
||||
int k = Triangulation_utils_3::next_around_edge(i,j);
|
||||
int l = Triangulation_utils_3::next_around_edge(j,i);
|
||||
|
||||
Vertex_handle vk = cell->vertex(k);
|
||||
Vertex_handle vl = cell->vertex(l);
|
||||
|
||||
|
|
@ -668,11 +1062,8 @@ private:
|
|||
Point ck = CGAL::circumcenter(pi,pj,pk);
|
||||
Point cl = CGAL::circumcenter(pi,pj,pl);
|
||||
|
||||
Triangle mcck(m,c,ck);
|
||||
Triangle mccl(m,c,cl);
|
||||
|
||||
area += std::sqrt(mcck.squared_area());
|
||||
area += std::sqrt(mccl.squared_area());
|
||||
area += std::sqrt(squared_area(m,c,ck));
|
||||
area += std::sqrt(squared_area(m,c,cl));
|
||||
}
|
||||
circ++;
|
||||
}
|
||||
|
|
@ -680,32 +1071,6 @@ private:
|
|||
return area;
|
||||
}
|
||||
|
||||
// Gets indices different from i and j
|
||||
void other_two_indices(int i, int j, int* k, int* l)
|
||||
{
|
||||
CGAL_surface_reconstruction_points_assertion(i != j);
|
||||
bool k_done = false;
|
||||
bool l_done = false;
|
||||
for(int index=0;index<4;index++)
|
||||
{
|
||||
if(index != i && index != j)
|
||||
{
|
||||
if(!k_done)
|
||||
{
|
||||
*k = index;
|
||||
k_done = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
*l = index;
|
||||
l_done = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
CGAL_surface_reconstruction_points_assertion(k_done);
|
||||
CGAL_surface_reconstruction_points_assertion(l_done);
|
||||
}
|
||||
|
||||
/// Assemble vi's row of the linear system A*X=B
|
||||
///
|
||||
/// @commentheading Template parameters:
|
||||
|
|
@ -722,39 +1087,51 @@ private:
|
|||
|
||||
double diagonal = 0.0;
|
||||
|
||||
for(typename std::vector<Edge>::iterator it = edges.begin();
|
||||
for(typename std::vector<Edge>::iterator it = edges.begin();
|
||||
it != edges.end();
|
||||
it++)
|
||||
{
|
||||
Vertex_handle vj = it->first->vertex(it->third);
|
||||
if(vj == vi){
|
||||
vj = it->first->vertex(it->second);
|
||||
{
|
||||
Vertex_handle vj = it->first->vertex(it->third);
|
||||
if(vj == vi){
|
||||
vj = it->first->vertex(it->second);
|
||||
}
|
||||
if(m_tr->is_infinite(vj))
|
||||
continue;
|
||||
|
||||
// get corresponding edge
|
||||
Edge edge( it->first, it->first->index(vi), it->first->index(vj));
|
||||
if(vi->index() < vj->index()){
|
||||
std::swap(edge.second, edge.third);
|
||||
}
|
||||
|
||||
double cij = cotan_geometric(edge);
|
||||
|
||||
if(m_tr->is_constrained(vj)){
|
||||
if(! is_valid(vj->f())){
|
||||
std::cerr << "vj->f() = " << vj->f() << " is not valid" << std::endl;
|
||||
}
|
||||
B[vi->index()] -= cij * vj->f(); // change rhs
|
||||
if(! is_valid( B[vi->index()])){
|
||||
std::cerr << " B[vi->index()] = " << B[vi->index()] << " is not valid" << std::endl;
|
||||
}
|
||||
|
||||
} else {
|
||||
if(! is_valid(cij)){
|
||||
std::cerr << "cij = " << cij << " is not valid" << std::endl;
|
||||
}
|
||||
A.set_coef(vi->index(),vj->index(), -cij, true /*new*/); // off-diagonal coefficient
|
||||
}
|
||||
|
||||
diagonal += cij;
|
||||
}
|
||||
if(m_tr->is_infinite(vj))
|
||||
continue;
|
||||
|
||||
// get corresponding edge
|
||||
Edge edge( it->first, it->first->index(vi), it->first->index(vj));
|
||||
if(vi->index() < vj->index()){
|
||||
std::swap(edge.second, edge.third);
|
||||
}
|
||||
|
||||
double cij = cotan_geometric(edge);
|
||||
if(vj->constrained())
|
||||
B[vi->index()] -= cij * vj->f(); // change rhs
|
||||
else
|
||||
A.set_coef(vi->index(),vj->index(), -cij, true /*new*/); // off-diagonal coefficient
|
||||
|
||||
diagonal += cij;
|
||||
}
|
||||
// diagonal coefficient
|
||||
if (vi->type() == Triangulation::INPUT)
|
||||
if (vi->type() == Triangulation::INPUT){
|
||||
A.set_coef(vi->index(),vi->index(), diagonal + lambda, true /*new*/) ;
|
||||
else
|
||||
} else{
|
||||
A.set_coef(vi->index(),vi->index(), diagonal, true /*new*/);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/// Computes enlarged geometric bounding sphere of the embedded triangulation.
|
||||
Sphere enlarged_bounding_sphere(FT ratio) const
|
||||
|
|
|
|||
|
|
@ -28,10 +28,14 @@
|
|||
#include <CGAL/surface_reconstruction_points_assertions.h>
|
||||
|
||||
#include <CGAL/Delaunay_triangulation_3.h>
|
||||
#include <CGAL/Triangulation_cell_base_with_info_3.h>
|
||||
#include <CGAL/Min_sphere_of_spheres_d.h>
|
||||
#include <CGAL/Min_sphere_of_spheres_d_traits_3.h>
|
||||
#include <CGAL/Min_sphere_of_points_d_traits_3.h>
|
||||
#include <CGAL/centroid.h>
|
||||
|
||||
#include <boost/random.hpp>
|
||||
#include <boost/random/linear_congruential.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
|
|
@ -90,26 +94,21 @@ private:
|
|||
public:
|
||||
|
||||
Reconstruction_vertex_base_3()
|
||||
: Vb(), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0)
|
||||
: Vb(), m_f(FT(0.0)), m_type(0), m_index(0)
|
||||
{}
|
||||
|
||||
Reconstruction_vertex_base_3(const Point_with_normal& p)
|
||||
: Vb(p), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0)
|
||||
: Vb(p), m_f(FT(0.0)), m_type(0), m_index(0)
|
||||
{}
|
||||
|
||||
Reconstruction_vertex_base_3(const Point_with_normal& p, Cell_handle c)
|
||||
: Vb(p,c), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0)
|
||||
: Vb(p,c), m_f(FT(0.0)), m_type(0), m_index(0)
|
||||
{}
|
||||
|
||||
Reconstruction_vertex_base_3(Cell_handle c)
|
||||
: Vb(c), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0)
|
||||
: Vb(c), m_f(FT(0.0)), m_type(0), m_index(0)
|
||||
{}
|
||||
|
||||
/// Is vertex constrained, i.e.
|
||||
/// does it contribute to the right or left member of the linear system?
|
||||
/// Default value is false.
|
||||
bool constrained() const { return m_constrained; }
|
||||
bool& constrained() { return m_constrained; }
|
||||
|
||||
/// Gets/sets the value of the implicit function.
|
||||
/// Default value is 0.0.
|
||||
|
|
@ -170,7 +169,7 @@ struct Reconstruction_triangulation_default_geom_traits_3 : public BaseGt
|
|||
|
||||
template <class BaseGt,
|
||||
class Gt = Reconstruction_triangulation_default_geom_traits_3<BaseGt>,
|
||||
class Tds_ = Triangulation_data_structure_3<Reconstruction_vertex_base_3<Gt> > >
|
||||
class Tds_ = Triangulation_data_structure_3<Reconstruction_vertex_base_3<Gt>, Triangulation_cell_base_with_info_3<int,Gt> > >
|
||||
class Reconstruction_triangulation_3 : public Delaunay_triangulation_3<Gt,Tds_>
|
||||
{
|
||||
// Private types
|
||||
|
|
@ -237,8 +236,8 @@ public:
|
|||
|
||||
/// Point type
|
||||
enum Point_type {
|
||||
INPUT, ///< Input point.
|
||||
STEINER ///< Steiner point created by Delaunay refinement.
|
||||
INPUT=0, ///< Input point.
|
||||
STEINER=1 ///< Steiner point created by Delaunay refinement.
|
||||
};
|
||||
|
||||
/// Iterator over input vertices.
|
||||
|
|
@ -249,13 +248,21 @@ public:
|
|||
typedef Iterator_project<Input_vertices_iterator,
|
||||
Project_point<Vertex> > Input_point_iterator;
|
||||
|
||||
// Public methods
|
||||
mutable Sphere sphere;
|
||||
std::vector<Point_with_normal> points;
|
||||
std::size_t fraction;
|
||||
std::list<double> fractions;
|
||||
Vertex_handle constrained_vertex;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
/// Default constructor.
|
||||
Reconstruction_triangulation_3()
|
||||
{
|
||||
}
|
||||
{}
|
||||
|
||||
~Reconstruction_triangulation_3()
|
||||
{}
|
||||
|
||||
// Default copy constructor and operator =() are fine.
|
||||
|
||||
|
|
@ -295,62 +302,55 @@ public:
|
|||
return Input_point_iterator(input_vertices_end());
|
||||
}
|
||||
|
||||
/// Gets the bounding sphere of all points.
|
||||
/// Gets the bounding sphere of input points.
|
||||
|
||||
|
||||
Sphere bounding_sphere() const
|
||||
{
|
||||
typedef Min_sphere_of_spheres_d_traits_3<Gt,FT> Traits;
|
||||
typedef Min_sphere_of_spheres_d<Traits> Min_sphere;
|
||||
typedef typename Traits::Sphere Traits_sphere;
|
||||
|
||||
// Represents *all* points by a set of spheres with 0 radius
|
||||
std::vector<Traits_sphere> spheres;
|
||||
spheres.reserve(number_of_vertices());
|
||||
|
||||
for (Point_iterator it=points_begin(), eit=points_end();
|
||||
it != eit; ++it)
|
||||
spheres.push_back(Traits_sphere(*it,0));
|
||||
|
||||
// Computes min sphere
|
||||
Min_sphere ms(spheres.begin(),spheres.end());
|
||||
typename Min_sphere::Cartesian_const_iterator coord = ms.center_cartesian_begin();
|
||||
FT cx = *coord++;
|
||||
FT cy = *coord++;
|
||||
FT cz = *coord++;
|
||||
return Sphere(Point(cx,cy,cz), ms.radius()*ms.radius());
|
||||
return sphere;
|
||||
}
|
||||
|
||||
/// Gets the bounding sphere of input points.
|
||||
Sphere input_points_bounding_sphere() const
|
||||
void initialize_bounding_sphere() const
|
||||
{
|
||||
typedef Min_sphere_of_spheres_d_traits_3<Gt,FT> Traits;
|
||||
typedef Min_sphere_of_points_d_traits_3<Gt,FT> Traits;
|
||||
typedef Min_sphere_of_spheres_d<Traits> Min_sphere;
|
||||
typedef typename Traits::Sphere Traits_sphere;
|
||||
|
||||
// Represents *input* points by a set of spheres with 0 radius
|
||||
std::vector<Traits_sphere> spheres;
|
||||
for (Input_point_iterator it=input_points_begin(), eit=input_points_end();
|
||||
it != eit;
|
||||
++it)
|
||||
spheres.push_back(Traits_sphere(*it,0));
|
||||
|
||||
|
||||
// Computes min sphere
|
||||
Min_sphere ms(spheres.begin(),spheres.end());
|
||||
Min_sphere ms(points.begin(),points.end());
|
||||
|
||||
typename Min_sphere::Cartesian_const_iterator coord = ms.center_cartesian_begin();
|
||||
FT cx = *coord++;
|
||||
FT cy = *coord++;
|
||||
FT cz = *coord++;
|
||||
return Sphere(Point(cx,cy,cz), ms.radius()*ms.radius());
|
||||
FT cz = *coord;
|
||||
|
||||
sphere = Sphere(Point(cx,cy,cz), ms.radius()*ms.radius());
|
||||
}
|
||||
|
||||
/// Insert point in the triangulation.
|
||||
/// Default type is INPUT.
|
||||
template <typename Visitor>
|
||||
Vertex_handle insert(const Point_with_normal& p,
|
||||
Point_type type = INPUT,
|
||||
Cell_handle start = Cell_handle())
|
||||
Point_type type,// = INPUT,
|
||||
Cell_handle start,// = Cell_handle(),
|
||||
Visitor visitor)
|
||||
{
|
||||
Vertex_handle v = Base::insert(p, start);
|
||||
|
||||
if(type == INPUT){
|
||||
visitor.before_insertion();
|
||||
}
|
||||
if(this->dimension() < 3){
|
||||
Vertex_handle v = Base::insert(p, start);
|
||||
v->type() = type;
|
||||
return v;
|
||||
}
|
||||
typename Base::Locate_type lt;
|
||||
int li, lj;
|
||||
Cell_handle ch = Base::locate(p, lt, li, lj, start);
|
||||
|
||||
Vertex_handle v = Base::insert(p, lt, ch, li, lj);
|
||||
v->type() = type;
|
||||
return v;
|
||||
|
||||
}
|
||||
|
||||
/// Insert the [first, beyond) range of points in the triangulation using a spatial sort.
|
||||
|
|
@ -367,57 +367,91 @@ public:
|
|||
// This variant requires all parameters.
|
||||
template <typename InputIterator,
|
||||
typename PointPMap,
|
||||
typename NormalPMap
|
||||
typename NormalPMap,
|
||||
typename Visitor
|
||||
>
|
||||
int insert(
|
||||
InputIterator first, ///< iterator over the first input point.
|
||||
InputIterator beyond, ///< past-the-end iterator over the input points.
|
||||
PointPMap point_pmap, ///< property map to access the position of an input point.
|
||||
NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point.
|
||||
Point_type type = INPUT)
|
||||
Visitor visitor)
|
||||
{
|
||||
int n = number_of_vertices();
|
||||
|
||||
if(! points.empty()){
|
||||
std::cerr << "WARNING: not all points inserted yet" << std::endl;
|
||||
}
|
||||
// Convert input points to Point_with_normal_3
|
||||
std::vector<Point_with_normal> points;
|
||||
//std::vector<Point_with_normal> points;
|
||||
for (InputIterator it = first; it != beyond; ++it)
|
||||
{
|
||||
Point_with_normal pwn(get(point_pmap,it), get(normal_pmap,it));
|
||||
points.push_back(pwn);
|
||||
}
|
||||
int n = points.size();
|
||||
|
||||
// Spatial sorting
|
||||
std::random_shuffle (points.begin(), points.end());
|
||||
spatial_sort (points.begin(), points.end(), geom_traits());
|
||||
initialize_bounding_sphere();
|
||||
|
||||
// Insert in triangulation
|
||||
boost::rand48 random;
|
||||
boost::random_number_generator<boost::rand48> rng(random);
|
||||
std::random_shuffle (points.begin(), points.end(), rng);
|
||||
fraction = 0;
|
||||
|
||||
fractions.clear();
|
||||
fractions.push_back(1.0);
|
||||
|
||||
double m = n;
|
||||
|
||||
while(m > 500){
|
||||
m /= 2;
|
||||
fractions.push_front(m/n);
|
||||
}
|
||||
|
||||
insert_fraction(visitor);
|
||||
return 0;
|
||||
}
|
||||
|
||||
template <typename Visitor>
|
||||
bool insert_fraction(Visitor visitor)
|
||||
{
|
||||
if(fractions.empty()){
|
||||
points.clear();
|
||||
return false;
|
||||
}
|
||||
double frac = fractions.front();
|
||||
fractions.pop_front();
|
||||
std::size_t more = (std::size_t)(points.size() * frac) - fraction;
|
||||
if((fraction+more) > points.size()){
|
||||
more = points.size() - fraction;
|
||||
}
|
||||
Cell_handle hint;
|
||||
for (typename std::vector<Point_with_normal>::const_iterator p = points.begin();
|
||||
p != points.end(); ++p)
|
||||
spatial_sort (points.begin()+fraction, points.begin()+fraction+more, geom_traits());
|
||||
for (typename std::vector<Point_with_normal>::const_iterator p = points.begin()+fraction;
|
||||
p != points.begin()+fraction+more; ++p)
|
||||
{
|
||||
Vertex_handle v = insert(*p, type, hint);
|
||||
Vertex_handle v = insert(*p, INPUT, hint, visitor);
|
||||
hint = v->cell();
|
||||
}
|
||||
|
||||
return number_of_vertices() - n;
|
||||
fraction += more;
|
||||
return true;
|
||||
}
|
||||
|
||||
/// @cond SKIP_IN_MANUAL
|
||||
// This variant creates a default point property map = Dereference_property_map.
|
||||
template <typename InputIterator,
|
||||
typename NormalPMap
|
||||
typename NormalPMap,
|
||||
typename Visitor
|
||||
>
|
||||
int insert(
|
||||
InputIterator first, ///< iterator over the first input point.
|
||||
InputIterator beyond, ///< past-the-end iterator over the input points.
|
||||
NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point.
|
||||
Point_type type = INPUT)
|
||||
Visitor visitor)
|
||||
{
|
||||
return insert(
|
||||
first,beyond,
|
||||
make_dereference_property_map(first),
|
||||
normal_pmap,
|
||||
type);
|
||||
visitor);
|
||||
}
|
||||
|
||||
/// Delaunay refinement callback:
|
||||
|
|
@ -443,12 +477,26 @@ public:
|
|||
v!= e;
|
||||
++v)
|
||||
{
|
||||
if(!v->constrained())
|
||||
if(! is_constrained(v))
|
||||
v->index() = index++;
|
||||
}
|
||||
return index;
|
||||
}
|
||||
|
||||
/// Is vertex constrained, i.e.
|
||||
/// does it contribute to the right or left member of the linear system?
|
||||
|
||||
bool is_constrained(Vertex_handle v) const
|
||||
{
|
||||
return v == constrained_vertex;
|
||||
}
|
||||
|
||||
void constrain(Vertex_handle v)
|
||||
{
|
||||
constrained_vertex = v;
|
||||
}
|
||||
|
||||
|
||||
}; // end of Reconstruction_triangulation_3
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,523 @@
|
|||
// Copyright (c) 2003-2007 INRIA Sophia-Antipolis (France).
|
||||
// Copyright (c) 2008 GeometryFactory, 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) : Steve OUDOT, Laurent RINEAU
|
||||
|
||||
#ifndef CGAL_SURFACE_MESHER_POISSON_IMPLICIT_SURFACE_ORACLE_3_H
|
||||
#define CGAL_SURFACE_MESHER_POISSON_IMPLICIT_SURFACE_ORACLE_3_H
|
||||
|
||||
#include <CGAL/Surface_mesher/Null_oracle_visitor.h>
|
||||
#include <CGAL/point_generators_3.h>
|
||||
#include <CGAL/Surface_mesher/Sphere_oracle_3.h>
|
||||
#include <CGAL/Real_embeddable_traits.h>
|
||||
|
||||
#include <boost/mpl/has_xxx.hpp>
|
||||
|
||||
#include <queue>
|
||||
|
||||
#include <CGAL/Surface_mesher/Profile_timer.h>
|
||||
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_IMPLICIT_ORACLE
|
||||
# define CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
|
||||
#endif
|
||||
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
|
||||
# include <string>
|
||||
# include <sstream>
|
||||
# include <boost/format.hpp>
|
||||
#endif
|
||||
|
||||
// NB: this oracle requires that the user provide a function that can
|
||||
// compute the value of the potential in any point of space
|
||||
namespace CGAL {
|
||||
|
||||
#ifdef BOOST_MPL_CFG_NO_HAS_XXX
|
||||
# error "BOOST_MPL_HAS_XXX_TRAIT_DEF is needed!"
|
||||
#else
|
||||
BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(has_set_on_surface, Set_on_surface_tag, true)
|
||||
#endif
|
||||
|
||||
namespace Surface_mesher {
|
||||
|
||||
/**
|
||||
Type and functions required in GT:
|
||||
- Compute_scalar_product_3 (from rev. 29646)
|
||||
- Compute_squared_distance_3
|
||||
- Compute_squared_radius_3
|
||||
- Construct_center_3
|
||||
- Construct_midpoint_3
|
||||
- Construct_point_on_3 (from rev. 29646)
|
||||
- Construct_scaled_vector_3
|
||||
- Construct_segment_3 (from rev. 29646)
|
||||
- Construct_translated_point_3
|
||||
- Construct_vector_3
|
||||
- FT
|
||||
- Has_on_bounded_side_3 (from rev. 29646)
|
||||
- Line_3
|
||||
- Point_3
|
||||
- Ray_3
|
||||
- Segment_3
|
||||
- Sphere_3
|
||||
- Vector_3 (from rev. 29646)
|
||||
- compute_scalar_product_3_object (from rev. 29646)
|
||||
- compute_squared_distance_3_object
|
||||
- compute_squared_radius_3_object
|
||||
- construct_center_3_object
|
||||
- construct_line_3_object (no longer, since rev. 29646)
|
||||
- construct_midpoint_3_object
|
||||
- construct_point_on_3_object (from rev. 29646)
|
||||
- construct_scaled_vector_3_object
|
||||
- construct_segment_3_object (from rev. 29646)
|
||||
- construct_translated_point_3_object
|
||||
- construct_vector_3_object
|
||||
- has_on_bounded_side_3_object (from rev. 29646)
|
||||
(Computed by use of:
|
||||
perl -ne '/GT\(\)\.([a-zA-Z_0-9]+)/
|
||||
&& print "$1\n";' {implicit_surface_oracle_3.h,Sphere_oracle_3.h} | sort -u
|
||||
perl -ne '/GT::([a-zA-Z_0-9]+)/
|
||||
&& print "$1\n";' {implicit_surface_oracle_3.h,Sphere_oracle_3.h} | sort -u
|
||||
)
|
||||
*/
|
||||
|
||||
namespace {
|
||||
|
||||
template <typename T>
|
||||
struct Return_min : std::binary_function<T, T, T>
|
||||
{
|
||||
T operator()(const T& a, const T& b) const
|
||||
{
|
||||
return (std::min)(a, b);
|
||||
}
|
||||
};
|
||||
|
||||
} // end anonymous namespace
|
||||
|
||||
template <
|
||||
class GT,
|
||||
class Surface,
|
||||
class Transform_functor_ =
|
||||
typename Real_embeddable_traits<typename Surface::FT>::Sgn,
|
||||
class Surface_identifiers_generator_ =
|
||||
Return_min<typename Transform_functor_::result_type>,
|
||||
class Point_creator = Creator_uniform_3<typename GT::FT,
|
||||
typename GT::Point_3>,
|
||||
class Visitor = Null_oracle_visitor
|
||||
>
|
||||
class Poisson_implicit_surface_oracle_3
|
||||
{
|
||||
// private types
|
||||
typedef Poisson_implicit_surface_oracle_3<GT,
|
||||
Surface,
|
||||
Transform_functor_,
|
||||
Surface_identifiers_generator_,
|
||||
Point_creator,
|
||||
Visitor> Self;
|
||||
|
||||
typedef 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 Transform_functor_ Transform_functor;
|
||||
typedef Surface_identifiers_generator_ Surface_identifiers_generator;
|
||||
|
||||
typedef Point_3 Intersection_point;
|
||||
|
||||
typedef Surface Surface_3;
|
||||
|
||||
typedef typename Surface::FT Surface_value_type;
|
||||
typedef typename Transform_functor::result_type result_type;
|
||||
|
||||
private:
|
||||
// Private members
|
||||
Visitor visitor; // a visitor that can modify a point, before returning
|
||||
// it.
|
||||
Transform_functor transform_functor;
|
||||
Surface_identifiers_generator surface_identifiers_generator;
|
||||
|
||||
public:
|
||||
|
||||
// Constructors
|
||||
// Poisson_implicit_surface_oracle_3 (Visitor visitor_ = Visitor()) :
|
||||
// visitor(visitor_),
|
||||
// transform_functor(),
|
||||
// surface_identifiers_generator()
|
||||
// {
|
||||
// #ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS
|
||||
// std::cerr << "CONS: Poisson_implicit_surface_oracle_3\n";
|
||||
// #endif
|
||||
// }
|
||||
|
||||
Poisson_implicit_surface_oracle_3 (Transform_functor transform_functor
|
||||
= Transform_functor(),
|
||||
Surface_identifiers_generator
|
||||
surface_identifiers_generator
|
||||
= Surface_identifiers_generator(),
|
||||
Visitor visitor_ = Visitor()) :
|
||||
visitor(visitor_),
|
||||
transform_functor(transform_functor),
|
||||
surface_identifiers_generator(surface_identifiers_generator)
|
||||
{
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS
|
||||
std::cerr << "CONS: Poisson_implicit_surface_oracle_3\n";
|
||||
#endif
|
||||
}
|
||||
|
||||
class Intersect_3
|
||||
{
|
||||
Visitor visitor;
|
||||
Transform_functor transform_functor;
|
||||
Surface_identifiers_generator surface_identifiers_generator;
|
||||
|
||||
// The two following template and specialization are used in
|
||||
// intersect_clipped_segment(...).
|
||||
template <typename Point, typename Identifier, bool b>
|
||||
struct Set_on_surface {
|
||||
void operator()(Point &, const Identifier&) const
|
||||
{
|
||||
}
|
||||
};
|
||||
|
||||
template <typename Point, typename Identifier>
|
||||
struct Set_on_surface<Point, Identifier, true> {
|
||||
void operator()(Point & p, const Identifier& id) const
|
||||
{
|
||||
p.set_on_surface(id);
|
||||
}
|
||||
};
|
||||
|
||||
Set_on_surface<Point,
|
||||
typename Surface_identifiers_generator::result_type,
|
||||
CGAL::has_set_on_surface<Point>::value> set_on_surface;
|
||||
|
||||
public:
|
||||
Intersect_3(Visitor visitor, Transform_functor transform_functor,
|
||||
Surface_identifiers_generator surface_identifiers_generator)
|
||||
: visitor(visitor),
|
||||
transform_functor(transform_functor),
|
||||
surface_identifiers_generator(surface_identifiers_generator)
|
||||
{
|
||||
}
|
||||
|
||||
Object operator()(const Surface_3& surface, Segment_3 s)
|
||||
// s is passed by value, because it is clipped below
|
||||
{
|
||||
CGAL_SURFACE_MESHER_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Segment_3)");
|
||||
CGAL_SURFACE_MESHER_TIME_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Segment_3)");
|
||||
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);
|
||||
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_IMPLICIT_ORACLE
|
||||
std::cerr <<
|
||||
boost::format("** Poisson_implicit_surface_oracle_3::operator()( (%1%), (%2%) )\n")
|
||||
% a % b;
|
||||
#endif
|
||||
// if both extremities are on the same side of the surface, return
|
||||
// no intersection
|
||||
if(surf_equation(surface, a) == surf_equation(surface, b))
|
||||
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
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_IMPLICIT_ORACLE
|
||||
std::cerr << "No intersection\n";
|
||||
#endif
|
||||
return Object();
|
||||
} // end operator()(Surface_3, Segment_3)
|
||||
|
||||
Object operator()(const Surface_3& surface, const Ray_3& r) {
|
||||
CGAL_SURFACE_MESHER_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Ray_3)");
|
||||
CGAL_SURFACE_MESHER_TIME_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Ray_3)");
|
||||
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) {
|
||||
CGAL_SURFACE_MESHER_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Line_3)");
|
||||
CGAL_SURFACE_MESHER_TIME_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Line_3)");
|
||||
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
|
||||
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()))
|
||||
<< ", value=" << surf_equation(surface, p)
|
||||
<< ")";
|
||||
return s.str();
|
||||
}
|
||||
|
||||
result_type surf_equation (const Surface_3& surface,
|
||||
const Point& p)
|
||||
{
|
||||
return transform_functor(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)
|
||||
{
|
||||
CGAL_SURFACE_MESHER_PROFILER("intersect_clipped_segment");
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
|
||||
std::cerr << "clipped_segment( (" << p1 << "), (" << p2 << ") )\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();
|
||||
|
||||
typedef typename Surface::Function::Cell_handle Cell_handle;
|
||||
|
||||
// Cannot be const: those values are modified below.
|
||||
FT value_at_p1, value_at_p2;
|
||||
Cell_handle c1, c2;
|
||||
bool c1_is_inf, c2_is_inf;
|
||||
|
||||
boost::tie(value_at_p1, c1, c1_is_inf) = surface.function().special_func(p1);
|
||||
boost::tie(value_at_p2, c2, c2_is_inf) = surface.function().special_func(p2);
|
||||
|
||||
// If both extremities are in the same volume component, returns
|
||||
// no intersection.
|
||||
if(transform_functor(value_at_p1) == transform_functor(value_at_p2))
|
||||
return Object();
|
||||
|
||||
#ifdef CGAL_SURFACE_MESHER_PROFILE
|
||||
int steps = 0;
|
||||
#endif
|
||||
while(true)
|
||||
{
|
||||
if(c1 == c2) {
|
||||
if(c1_is_inf) {
|
||||
return Object();
|
||||
} else {
|
||||
return make_object(Point(ORIGIN+ ((value_at_p2 * (p1 - ORIGIN)) - (value_at_p1 * (p2 - ORIGIN))) / (value_at_p2 - value_at_p1)));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#ifdef CGAL_SURFACE_MESHER_PROFILE
|
||||
++steps;
|
||||
#endif
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_CLIPPED_SEGMENT
|
||||
std::cerr << debug_point(surface, p1) << ", "
|
||||
<< debug_point(surface, p2) << "\n";
|
||||
#endif
|
||||
// mid cannot be const, because it is modified below.
|
||||
Point mid = midpoint(p1, p2);
|
||||
Cell_handle c_at_mid;
|
||||
FT value_at_mid;
|
||||
bool c_is_inf;
|
||||
boost::tie(value_at_mid, c_at_mid, c_is_inf) = surface.function().special_func(mid);
|
||||
|
||||
if ( 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
|
||||
// the following function conditionnally call
|
||||
// mid.set_on_surface(...) if mid has such a function.
|
||||
set_on_surface(mid,
|
||||
surface_identifiers_generator(transform_functor(value_at_p1),
|
||||
transform_functor(value_at_p2)));
|
||||
|
||||
visitor.new_point(mid);
|
||||
CGAL_SURFACE_MESHER_HISTOGRAM_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Segment_3) bissection steps", steps)
|
||||
return make_object(mid);
|
||||
}
|
||||
|
||||
// Else we must go on
|
||||
if ( transform_functor(value_at_p1) != transform_functor(value_at_mid) )
|
||||
{
|
||||
p2 = mid;
|
||||
value_at_p2 = value_at_mid;
|
||||
c2 = c_at_mid;
|
||||
c2_is_inf = c_is_inf;
|
||||
}
|
||||
else
|
||||
{
|
||||
p1 = mid;
|
||||
value_at_p1 = value_at_mid;
|
||||
c1 = c_at_mid;
|
||||
c1_is_inf = c_is_inf;
|
||||
}
|
||||
}
|
||||
} // end intersect_clipped_segment
|
||||
|
||||
}; // end nested class Intersect_3
|
||||
|
||||
class Construct_initial_points
|
||||
{
|
||||
const Self& oracle;
|
||||
public:
|
||||
Construct_initial_points(const Self& oracle) : oracle(oracle)
|
||||
{
|
||||
}
|
||||
|
||||
// Random points
|
||||
template <typename OutputIteratorPoints>
|
||||
OutputIteratorPoints operator() (const Surface_3& surface,
|
||||
OutputIteratorPoints out,
|
||||
int n = 20) // WARNING: why 20?
|
||||
{
|
||||
const Sphere_3& sphere = surface.bounding_sphere();
|
||||
const Point initial_center = GT().construct_center_3_object()(sphere);
|
||||
const FT squared_radius =
|
||||
GT().compute_squared_radius_3_object()(sphere);
|
||||
const FT radius = CGAL::sqrt(squared_radius);
|
||||
typename Self::Intersect_3 intersect = oracle.intersect_3_object();
|
||||
|
||||
Random rng(0);
|
||||
typename CGAL::Random_points_on_sphere_3<Point,
|
||||
Point_creator> random_point_on_sphere(CGAL::to_double(radius), rng);
|
||||
typename CGAL::Random_points_in_sphere_3<Point,
|
||||
Point_creator> random_point_in_sphere(CGAL::to_double(radius), rng);
|
||||
typename GT::Construct_segment_3 segment_3 =
|
||||
GT().construct_segment_3_object();
|
||||
typename GT::Construct_vector_3 vector_3 =
|
||||
GT().construct_vector_3_object();
|
||||
typename GT::Construct_translated_point_3 translate =
|
||||
GT().construct_translated_point_3_object();
|
||||
|
||||
Point center = initial_center;
|
||||
while (n>0)
|
||||
{
|
||||
const Point p = translate(*random_point_on_sphere++,
|
||||
vector_3(CGAL::ORIGIN, initial_center));
|
||||
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_INITIAL_POINTS
|
||||
std::cerr << "test "
|
||||
<< intersect.debug_point(surface, center)
|
||||
<< ", " << intersect.debug_point(surface, p);
|
||||
#endif
|
||||
|
||||
Object o = intersect(surface, segment_3(center, p));
|
||||
if (const Point* intersection = object_cast<Point>(&o)) {
|
||||
*out++= *intersection;
|
||||
--n;
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_INITIAL_POINTS
|
||||
std::cerr << " = "
|
||||
<< intersect.debug_point(surface, *intersection)
|
||||
<< "\n";
|
||||
#endif
|
||||
}
|
||||
else
|
||||
{
|
||||
center = translate(*random_point_in_sphere++,
|
||||
vector_3(CGAL::ORIGIN, initial_center));
|
||||
#ifdef CGAL_SURFACE_MESHER_VERBOSE
|
||||
std::cerr << "Warning: new center " << center << "\n";
|
||||
#endif
|
||||
#ifdef CGAL_SURFACE_MESHER_DEBUG_INITIAL_POINTS
|
||||
std::cerr << " = null\n";
|
||||
#endif
|
||||
}
|
||||
}
|
||||
return out;
|
||||
}
|
||||
}; // end nested class Construct_initial_points
|
||||
|
||||
Construct_initial_points construct_initial_points_object() const
|
||||
{
|
||||
return Construct_initial_points(*this);
|
||||
}
|
||||
|
||||
Intersect_3 intersect_3_object() const
|
||||
{
|
||||
return Intersect_3(visitor,
|
||||
transform_functor,
|
||||
surface_identifiers_generator);
|
||||
}
|
||||
|
||||
bool is_in_volume(const Surface_3& surface, const Point& p) const
|
||||
{
|
||||
return transform_functor(surface(p)) != 0;
|
||||
}
|
||||
|
||||
}; // end Poisson_implicit_surface_oracle_3
|
||||
|
||||
template <typename FT>
|
||||
FT approximate_sqrt(const FT x) {
|
||||
return FT (CGAL_NTS sqrt(CGAL_NTS to_double(x)));
|
||||
}
|
||||
|
||||
} // namespace Surface_mesher
|
||||
|
||||
} // namespace CGAL
|
||||
|
||||
|
||||
#endif // CGAL_SURFACE_MESHER_POISSON_IMPLICIT_SURFACE_ORACLE_3_H
|
||||
|
|
@ -1,98 +0,0 @@
|
|||
// Copyright (c) 2007-08 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// 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) : Pierre Alliez and Laurent Saboret
|
||||
|
||||
#ifndef CGAL_K_NEIGHBOR_NEIGHBOR_H
|
||||
#define CGAL_K_NEIGHBOR_NEIGHBOR_H
|
||||
|
||||
#include <CGAL/Orthogonal_k_neighbor_search.h>
|
||||
#include <CGAL/Search_traits_vertex_handle_3.h>
|
||||
#include <list>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
|
||||
/// Wrapper around Orthogonal_k_neighbor_search for Triangulation_3 Vertex_handles.
|
||||
template <class Gt, class Vertex_handle>
|
||||
class K_nearest_neighbor
|
||||
{
|
||||
public:
|
||||
typedef Gt Geom_traits;
|
||||
|
||||
typedef typename Geom_traits::FT FT;
|
||||
typedef typename Geom_traits::Point_3 Point;
|
||||
typedef CGAL::Point_vertex_handle_3<Vertex_handle> Point_vertex_handle_3;
|
||||
typedef Search_traits_vertex_handle_3<Vertex_handle> Traits;
|
||||
typedef Euclidean_distance_vertex_handle_3<Vertex_handle> KDistance;
|
||||
typedef Orthogonal_k_neighbor_search<Traits,KDistance> Neighbor_search;
|
||||
typedef typename Neighbor_search::Tree Tree;
|
||||
typedef typename Neighbor_search::iterator Search_iterator;
|
||||
|
||||
private:
|
||||
Tree m_tree;
|
||||
|
||||
public:
|
||||
K_nearest_neighbor() {}
|
||||
|
||||
/// @commentheading Precondition:
|
||||
/// InputIterator value_type must be convertible to Point_vertex_handle_3.
|
||||
template <class InputIterator>
|
||||
K_nearest_neighbor(InputIterator first, InputIterator beyond)
|
||||
{
|
||||
m_tree = Tree(first, beyond);
|
||||
}
|
||||
|
||||
/// Default copy constructor, operator =() and destructor are fine.
|
||||
|
||||
/// @commentheading Precondition:
|
||||
/// InputIterator value_type must be convertible to Point_vertex_handle_3.
|
||||
template <class InputIterator>
|
||||
void insert(InputIterator first, InputIterator beyond)
|
||||
{
|
||||
m_tree = Tree(first, beyond);
|
||||
}
|
||||
|
||||
/// Empty KD-tree.
|
||||
void clear()
|
||||
{
|
||||
m_tree = Tree();
|
||||
}
|
||||
|
||||
/// Search 'nb' nearest_neighbors of 'query' point.
|
||||
bool get_k_nearest_neighbors(const Point& query,
|
||||
const unsigned int nb,
|
||||
std::list<Vertex_handle>& nearest_neighbors)
|
||||
{
|
||||
Point_vertex_handle_3 point_wrapper(query.x(), query.y(), query.z(), NULL);
|
||||
Neighbor_search search(m_tree, point_wrapper, nb); // only nb nearest neighbors
|
||||
Search_iterator it = search.begin();
|
||||
for(unsigned int i=0;i<nb;i++,it++)
|
||||
{
|
||||
if(it == search.end())
|
||||
return false;
|
||||
nearest_neighbors.push_back((Vertex_handle)it->first);
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
}; // K_nearest_neighbor
|
||||
|
||||
|
||||
} //namespace CGAL
|
||||
|
||||
#endif // CGAL_K_NEIGHBOR_NEIGHBOR_H
|
||||
|
|
@ -68,7 +68,6 @@ public:
|
|||
Poisson_mesher_level_impl_base(Tr& t, Criteria crit, unsigned int max_vert, Surface& surface, Oracle& oracle)
|
||||
: Base(t, crit, surface, oracle),
|
||||
max_vertices(max_vert) ///< number of vertices bound (ignored if zero)
|
||||
|
||||
{}
|
||||
|
||||
protected:
|
||||
|
|
@ -192,31 +191,49 @@ public:
|
|||
/// @commentheading Template Parameters:
|
||||
/// @param Tr 3D Delaunay triangulation.
|
||||
/// @param Surface Sphere_3 or Iso_cuboid_3.
|
||||
/// @param Sizing_field A sizing field functor type
|
||||
/// @param Second_sizing_field A sizing field functor type
|
||||
///
|
||||
/// @commentheading Sizing fields
|
||||
/// - The first sizing field is the real sizing field that is targeted by
|
||||
/// the refinement process. It may be costly to use.
|
||||
/// - The second sizing field is supposed to be a sizing field that is less
|
||||
/// costly to use (such as a constant sizing field). If a cell has a radius
|
||||
/// that is smaller than the value of the second sizing field at its
|
||||
/// center, then the cell is considered as small enough and the first
|
||||
/// sizing field is not evaluated for that cell.
|
||||
template <typename Tr,
|
||||
typename Surface>
|
||||
typename Surface,
|
||||
typename Sizing_field,
|
||||
typename Second_sizing_field>
|
||||
unsigned int poisson_refine_triangulation(
|
||||
Tr& tr,
|
||||
double radius_edge_ratio_bound, ///< radius edge ratio bound (>= 1.0)
|
||||
double cell_radius_bound, ///< cell radius bound (ignored if zero)
|
||||
const Sizing_field& sizing_field, ///< sizing field for cell radius bound
|
||||
const Second_sizing_field& second_sizing_field, ///< second sizing field for cell radius bound
|
||||
unsigned int max_vertices, ///< number of vertices bound (ignored if zero)
|
||||
Surface& enlarged_bbox) ///< new bounding sphere or box
|
||||
{
|
||||
CGAL_TRACE("Calls poisson_refine_triangulation()\n");
|
||||
|
||||
// Convergence is guaranteed if radius_edge_ratio_bound >= 1.0
|
||||
CGAL_surface_reconstruction_points_precondition(radius_edge_ratio_bound >= 1.0);
|
||||
|
||||
// Mesher_level types
|
||||
typedef Poisson_mesh_cell_criteria_3<Tr> Tets_criteria;
|
||||
typedef Poisson_mesh_cell_criteria_3<
|
||||
Tr
|
||||
, Sizing_field
|
||||
, Second_sizing_field
|
||||
> Tets_criteria;
|
||||
typedef typename CGAL::Surface_mesh_traits_generator_3<Surface>::type Oracle;
|
||||
typedef Poisson_mesher_level<Tr, Tets_criteria, Surface, Oracle, Null_mesher_level> Refiner;
|
||||
|
||||
CGAL_TRACE(" Creates queue\n");
|
||||
|
||||
int nb_vertices = tr.number_of_vertices(); // get former #vertices
|
||||
|
||||
// Delaunay refinement
|
||||
Tets_criteria tets_criteria(radius_edge_ratio_bound*radius_edge_ratio_bound, cell_radius_bound);
|
||||
Tets_criteria tets_criteria(radius_edge_ratio_bound*radius_edge_ratio_bound,
|
||||
sizing_field,
|
||||
second_sizing_field);
|
||||
Oracle oracle;
|
||||
Null_mesher_level null_mesher_level;
|
||||
Refiner refiner(tr, tets_criteria, max_vertices, enlarged_bbox, oracle, null_mesher_level);
|
||||
|
|
@ -225,7 +242,6 @@ unsigned int poisson_refine_triangulation(
|
|||
|
||||
int nb_vertices_added = tr.number_of_vertices() - nb_vertices;
|
||||
|
||||
CGAL_TRACE("End of poisson_refine_triangulation()\n");
|
||||
|
||||
return (unsigned int) nb_vertices_added;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -201,7 +201,8 @@ int main(int argc, char * argv[])
|
|||
// The position property map can be omitted here as we use iterators over Point_3 elements.
|
||||
Poisson_reconstruction_function function(
|
||||
points.begin(), points.end(),
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()));
|
||||
CGAL::make_normal_of_point_with_normal_pmap(points.begin()),
|
||||
CGAL::Poisson_visitor());
|
||||
|
||||
// Computes the Poisson indicator function f()
|
||||
// at each vertex of the triangulation.
|
||||
|
|
|
|||
Loading…
Reference in New Issue