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

CGAL_pre and post condition added

This commit is contained in:
Marc Pouget 2006-08-24 12:54:27 +00:00
parent c8502e1c5a
commit 3b4f578ae6
15 changed files with 474 additions and 387 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitattributes vendored
View File

@ -1546,7 +1546,7 @@ Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges_mesh.pdf -text svneol=unset#applicati
Ridges_3/examples/Ridges_3/data/bezierpb3-0.007812.off -text svneol=unset#application/octet-stream
Ridges_3/examples/Ridges_3/data/ellipe0.003.off -text svneol=unset#application/octet-stream
Ridges_3/examples/Ridges_3/data/ellipsoid_u_0.02.off -text
Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500-off -text
Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500.off -text
Ridges_3/examples/Ridges_3/data/venus.off -text svneol=unset#application/octet-stream
Robustness/demo/Robustness/help/index.html svneol=native#text/html
Scripts/developer_scripts/create_macosx_installer -text

2
Ridges_3/demo/Ridges_3/README
View File

@ -12,3 +12,5 @@ it displays an OpenGL view of both the mesh and the ridge lines
./introspect-qt.exe ../../examples/Ridges_3/data/poly2x^2+y^2-0.062500.off ../../examples/Ridges_3/data_poly2x^2+y^2-0.062500.off.4ogl.txt

51
Ridges_3/demo/Ridges_3/SketchSample.cpp
View File

@ -6,6 +6,9 @@
#include "SketchSample.h"
extern double strength_threshold;
extern double sharpness_threshold;
SketchSample::SketchSample(Mesh* mesh, DS* ridge_data) {
highlight = false;
p_mesh = mesh;
@ -19,8 +22,8 @@ void SketchSample::buildDisplayList(GLuint surf) {
glNewList(surf, GL_COMPILE);
glEnable(GL_LIGHTING);
// glShadeModel(GL_SMOOTH);
glShadeModel(GL_FLAT);
glShadeModel(GL_SMOOTH);
//glShadeModel(GL_FLAT);
//mesh glMaterialfv is def in the mesh with offset to see the lines z_buffered
p_mesh->gl_draw_facets(true);
@ -29,8 +32,11 @@ void SketchSample::buildDisplayList(GLuint surf) {
glColor3d(1., 1., 0.);
for (DS_iterator it=p_ridge_data->begin();it!=p_ridge_data->end();it++)
draw_one_ridge(*it);
{
if ( (*it)->strength >= strength_threshold &&
(*it)->sharpness >= sharpness_threshold )
draw_one_ridge(*it);
}
glEnable(GL_LIGHTING);
//additional objects that may be displayed
@ -43,6 +49,26 @@ void SketchSample::buildDisplayList(GLuint surf) {
glEndList();
}
void SketchSample::draw_one_ridge(data_line* line)
{
if (line->ridge_type == CGAL::BLUE_ELLIPTIC_RIDGE) glColor3f(0.,0.,1.);
if (line->ridge_type == CGAL::BLUE_HYPERBOLIC_RIDGE) glColor3f(0.,1.,0.);
if (line->ridge_type == CGAL::BLUE_CREST) glColor3f(0.,0.,1.);
if (line->ridge_type == CGAL::RED_ELLIPTIC_RIDGE) glColor3f(1.,0.,0.);
if (line->ridge_type == CGAL::RED_HYPERBOLIC_RIDGE) glColor3f(1.,1.,0.);
if (line->ridge_type == CGAL::RED_CREST) glColor3f(1.,0.,0.);
std::vector<Point>::iterator iter = line->ridge_points.begin(),
ite = line->ridge_points.end();
glLineWidth(3 );
glBegin(GL_LINE_STRIP);
for (;iter!=ite;iter++) glVertex3d(iter->x(), iter->y(), iter->z());
glEnd();
}
void SketchSample::buildPicking(GLuint pickSurf) {
int red, green, blue;
int i = 1;
@ -99,20 +125,3 @@ const double* SketchSample::rcoord() {
}
void SketchSample::draw_one_ridge(data_line* line)
{
if (line->ridge_type == BE) glColor3f(0.,0.,1.);
if (line->ridge_type == BH) glColor3f(0.,1.,0.);
if (line->ridge_type == BC) glColor3f(0.,0.,1.);
if (line->ridge_type == RE) glColor3f(1.,0.,0.);
if (line->ridge_type == RH) glColor3f(1.,1.,0.);
if (line->ridge_type == RC) glColor3f(1.,0.,0.);
std::vector<Point>::iterator iter = line->ridge_points.begin(),
ite = line->ridge_points.end();
glLineWidth(3 );
glBegin(GL_LINE_STRIP);
for (;iter!=ite;iter++) glVertex3d(iter->x(), iter->y(), iter->z());
glEnd();
}

7
Ridges_3/demo/Ridges_3/SketchSample.h
View File

@ -7,11 +7,16 @@
#include "Sketcher.h"
#include <iostream>
#include <CGAL/basic.h>
#include <CGAL/Cartesian.h>
#include <CGAL/Ridges.h>
//marc
#include "visu_poly.h"
#include "enriched_polyhedron.h"
enum Ridge_type {NONE=0, BLUE_RIDGE, RED_RIDGE, CREST, BE, BH, BC, RE, RH, RC};
/* extern CGAL::Ridge_type NONE, BLUE_RIDGE, RED_RIDGE, CREST, */
/* BLUE_ELLIPTIC_RIDGE, BLUE_HYPERBOLIC_RIDGE, BLUE_CREST, */
/* RED_ELLIPTIC_RIDGE, RED_HYPERBOLIC_RIDGE, RED_CREST; */
//Data structure for one line of the input file .4ogl.txt
struct data_line{

2
Ridges_3/demo/Ridges_3/makefile
View File

@ -75,7 +75,5 @@ depend:
# DO NOT DELETE
SketchSample.o: SketchSample.h visu_poly.h enriched_polyhedron.h
SketchSample.o: visu_poly.h enriched_polyhedron.h
visu.o: SketchSample.h visu_poly.h enriched_polyhedron.h
visu_poly.o: enriched_polyhedron.h

77
Ridges_3/demo/Ridges_3/visu.cpp
View File

@ -20,8 +20,9 @@
Mesh m_mesh;
DS ridge_data;
double strength_threshold = 0;
double sharpness_threshold = 0;
void clean_DS(DS& L)
{
@ -57,7 +58,7 @@ read_line(std::ifstream& stream_res, int& ridge_type,
>> strength
>> sharpness;
//debug
std::cout <<ridge_type << " ";
// std::cout << strength << " ";
while ( issline.good() ) {
Point p;//there are at least 2 points in a line
@ -91,7 +92,15 @@ void load_data_from_file(DS& l, const char* file_res)
void load_geom(int argc, char* argv[])
{
assert(argc==3);
if (argc != 5)
{ std::cout << "Usage : prog + " << std::endl
<< " file.off " << std::endl
<< " file....4ogl.txt" << std::endl
<< " strength_threshold "<< std::endl
<< " sharpness_threshold " << std::endl;
exit(0);
}
const char* file_off = argv[1];
const char* file_res = argv[2];
@ -107,6 +116,9 @@ void load_geom(int argc, char* argv[])
ridge_data.clear();
load_data_from_file(ridge_data, file_res);
//init thresholds
strength_threshold = atof(argv[3]);
sharpness_threshold = atof(argv[4]);
}
@ -137,42 +149,43 @@ int main(int argc, char** argv) {
main.show();
//debug visu with a jvx file
std::cout << ridge_data.size() << std::endl;
DS_iterator iter_lines = ridge_data.begin(), iter_end = ridge_data.end();
std::ofstream out_jvx("debug.jvx");
CGAL::Javaview_writer<std::ofstream> jvw(out_jvx);
jvw.set_title("ridges");
jvw.write_header();
// //debug visu with a jvx file
// std::cout << ridge_data.size() << std::endl;
// DS_iterator iter_lines = ridge_data.begin(), iter_end = ridge_data.end();
// //first the polysurf
// jvw.set_geometry_name("polysurf");
// jvw.begin_geometry();
// polyhedron_javaview_writer(jvw, P);
// jvw.end_geometry();
// std::ofstream out_jvx("debug.jvx");
// CGAL::Javaview_writer<std::ofstream> jvw(out_jvx);
// jvw.set_title("ridges");
// jvw.write_header();
int compt = 0;
// // //first the polysurf
// // jvw.set_geometry_name("polysurf");
// // jvw.begin_geometry();
// // polyhedron_javaview_writer(jvw, P);
// // jvw.end_geometry();
for (;iter_lines!=iter_end;iter_lines++) {
compt++;
// create the name of the ridge
std::ostringstream str;
str << "ridge " << compt;
jvw.set_geometry_name(str.str());
// int compt = 0;
//color
if ((*iter_lines)->ridge_type ==BC) jvw.set_color(CGAL::BLUE);
else jvw.set_color(CGAL::RED);
// for (;iter_lines!=iter_end;iter_lines++) {
// compt++;
// // create the name of the ridge
// std::ostringstream str;
// str << "ridge " << compt;
// jvw.set_geometry_name(str.str());
//lines
jvw.begin_geometry();
polyline_javaview_writer(jvw, (*iter_lines)->ridge_points.begin(),
(*iter_lines)->ridge_points.end());
jvw.end_geometry();
}
// //color
// if ((*iter_lines)->ridge_type == CGAL::BLUE_CREST) jvw.set_color(CGAL::BLUE);
// else jvw.set_color(CGAL::RED);
jvw.write_footer();
// //lines
// jvw.begin_geometry();
// polyline_javaview_writer(jvw, (*iter_lines)->ridge_points.begin(),
// (*iter_lines)->ridge_points.end());
// jvw.end_geometry();
// }
// jvw.write_footer();

160
Ridges_3/doc_tex/Ridges_3/Ridges_3_user.tex
View File

@ -1,49 +1,57 @@
\newtheorem{definition}{Definition.}
This chapter describes the \cgal's package for the extraction of
ridges on meshes. Given a smooth surface, a ridge is a curve along
which one of the principal curvatures has an extremum along its
curvature line. Ridges are curves of {\em extremal} curvature and
ridges and umbilics on meshes. Given a smooth surface, a ridge is a
curve along which one of the principal curvatures has an extremum
along its curvature line. An umbilic is a point at which both
principal curvatures are equal. Umbilics are special points on the
ridge lines. Ridges are curves of {\em extremal} curvature and
therefore encode important informations used in segmentation,
registration, matching and surface analysis. Based on the results of
the article \cite{rr}, we propose several algorithms to identify and
extract different parts of this singular ridge curve and some of its
singularites on a surface given as a mesh.
the article
\cite{rr}, we propose several algorithms to identify and extract
different parts of this singular ridge curve as well as umbilics on a
surface given as a triangular mesh.
\subsection{Overview}
%%%%%%%%%%%%%%%%%%%%%%
2 main versions
2 main independent classes
\begin{itemize}
\item
non-topologically oriented: we output several lists of ridges which
are polylines, each list associated to a type: crest min/max,
ellip/hyper, red/blue.
each line has a weight enabling interactive filtering for the visu
\ccc{Umbilic_approximation}
\item
topologically oriented : we output a list of umbilics and some lists
of ridges outside umbilic patches, turning points?
\ccc{Ridge_approximation}
\end{itemize}
And two ``container'' classes.
%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Ridges of a smooth surface}
\subsection{Ridges and umbilics of a smooth surface}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A comprehensive description of ridges can be found in
\cite{ip-nss-71,hgyg-ttdpf-99,ip-gd-01}, and in the sequel, we just
\cite{hgyg-ttdpf-99,ip-gd-01}, and in the sequel, we just
introduce the basic notions so as to explain our algorithms.
Consider a smooth embedded surface, and denote $k_1$ and $k_2$ the
principal curvatures, with $k_1\geq k_2$. Umbilics are the points
where $k_1=k_2$. For any non umbilical point, the corresponding
principal directions of curvature are well defined, and we denote them
$d_1$ and $d_2$. In local coordinates, we denote $\langle , \rangle$
$d_1$ and $d_2$.
%%
Anything related to the maximal (minimal) curvature is qualified blue
(red), for example we shall speak of the blue curvature for $k_1$ or
the red direction for $d_2$.
%%
In local coordinates, we denote $\langle , \rangle$
the inner product induced by the ambient Euclidean space, and $dk_1$,
$dk_2$ the gradients of the principal curvatures. Ridges, illustrated
on Figs \ref{pict:ellipsoid_ridges} and \ref{fig:ridges_ellipsoid},
@ -69,71 +77,67 @@ dk_2,d_2 \rangle$ vanishes, i.e. $b_3=0$.
%\rangle$, and we call these functions the {\em extremality
%coefficients}.
%%
Notice that, as the principal curvatures are not differentiable at
umbilics, the extremality coefficients are not defined at such
points. In addition, the sign of the extremality coefficients is not
defined since the principal directions can be oriented by two opposite
unit vectors. Apart from umbilics, special points on ridges are {\em
purple} points, which correspond to intersections between red and a
blue ridges. The previous characterization of ridges involves
third-order differential properties. Using fourth-order differential
quantities, a ridge point can further be qualified as {\em elliptic}
if it corresponds to a maximum of $k_1$ or a minimum of $k_2$, or {\em
hyperbolic} otherwise. Ridges of a given color change from elliptic to
hyperbolic at special points called {\em turning} points.
The previous characterization of ridges involves third-order
differential properties. Using fourth-order differential quantities, a
ridge point can further be qualified as {\em elliptic} if it
corresponds to a maximum of $k_1$ or a minimum of $k_2$, or {\em
hyperbolic} otherwise. Hence we end with four types of ridges, namely
: \ccc{BLUE_ELLIPTIC_RIDGE}, \ccc{BLUE_HYPERBOLIC_RIDGE}, \ccc{RED_ELLIPTIC_RIDGE},
\ccc{RED_HYPERBOLIC_RIDGE}.
In addition, a subset of elliptic ridges, called the crest lines,
which can be seen as the visually most salient curves on a surface are
also of interest. A crest line is an elliptic ridge which is a maximum
of $\max(|k_1|,|k_2|)$. Hence we provide two additional ridge types :
\ccc{BLUE_CREST} and \ccc{RED_CREST}.
The calculation of ridges poses difficult problems, which are of three
kinds.
Weigths are associated to each line
\begin{itemize}
\item
the {\bf strength} which is the integral of the curvature along the
line.
\item
the {\bf sharpness} which is the integral of the absolute value of
the second derivative of the curvature along the ridge.
\end{itemize}
\paragraph{Topological difficulties.}
Ridges of a smooth surface form a singular curve on the surface, with
self-intersections at umbilics (more precisely at so-called 3-ridges
umbilics), and purple points. Moreover, ridges have complex
interactions with curvature lines at turning points. From the
application standpoint, reporting ridges of a surface faithfully
requires reporting umbilics, purple points and turning points.
%\begin{minipage}[c]{.45\linewidth}
\paragraph{Numerical difficulties.}
As pointed out above, ridges are characterized and qualified through
third and fourth order derivatives of the surface. Estimating them
depends on the particular type of surface processed ---implicitly
defined, parameterized, discretized by a mesh--- and is numerically a
difficult task.
\paragraph{Orientation difficulties.}
Since coefficients $b_0$ and $b_3$ depend on a given orientation of
the principal directions, their sign is not well defined.
Practically, tracking the sign change of functions whose sign depends
on the particular orientation of the frame in which they are expressed
poses a problem. In particular, tracking a zero-crossing of $b_0$ or
$b_3$ from sign changes along a curve segment on the surface imposes
to find a coherent orientation of the principal frame at the
endpoints. Given two principal directions at these endpoints, one way
to find a local orientation consists of choosing two vectors so that
they make an acute angle, whence the name {\em Acute Rule}.
%%
\begin{minipage}[c]{.45\linewidth}
\begin{figure}[H]
\begin{figure}[!ht]
\begin{ccTexOnly}
\centerline{
\includegraphics[height=5cm]{Ridges_3/ellipsoid_ridges_mesh}}
\end{ccTexOnly}
\begin{ccHtmlOnly}
<CENTER> <img border=0 src="./ellipsoid_ridges_mesh.jpg" width=400>
</CENTER>
\end{ccHtmlOnly}
\caption{Umbilics, ridges, and principal blue foliation on the
ellipsoid (10k points)}
\label{pict:ellipsoid_ridges}
\end{figure}
\end{minipage}
%\end{minipage}
%%
\hfill
%%
\begin{minipage}[c]{.45\linewidth}
%\begin{minipage}[c]{.45\linewidth}
\begin{figure}[H]
\begin{ccTexOnly}
\centerline{
\includegraphics[height=5cm]{Ridges_3/ellipsoid_ridges}}
\end{ccTexOnly}
\begin{ccHtmlOnly}
<CENTER> <img border=0 src="./ellipsoid_ridges.jpg" width=400>
</CENTER>
\end{ccHtmlOnly}
\caption{Schematic view of the umbilics and the ridges. Max of $k_1$:
blue; Min of $k_1$: green; Min of $k_2$: red; Max of $k_2$: yellow}
\label{fig:ridges_ellipsoid}
\end{figure}
\end{minipage}
%\end{minipage}
@ -141,17 +145,43 @@ blue; Min of $k_1$: green; Min of $k_2$: red; Max of $k_2$: yellow}
\section{Software Design}
%%%%%%%%%%%%%%%%%%%%%%%
usage of pm, triangular meshes.
\subsection{Options and interface specifications}
%%%%%%%%%%%%%%%%%%%%%
ridges: approx and container class
ridges compute with para type and tag
Umbilics : approx and container,
umbilic compute with para
size
\subsection{Template parameters}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
On a des concepts: Vertex2FTPropertyMap et Vertex2VectorPropertyMap
qui sont specialisent le concept de Propety Map de boost, avec
\ccc{key_type CGAL::Polyhedron_3::Vertex_handle} et de \ccc{value_type
CGAL::Polyhedron_3::Traits::FT ou Vector_3 }. On a donc les fct get,
put, []; preciser LvaluePropertyMap ? Utilisation stockage
d'information scalaire ou vectorielle pour les vertex d'un polyhedron,
peut-etre dans le vertex lui-meme ou externe avec une std::map par
exemple.
Poly est un concept qui a pour model \ccc{CGAL::Polyhedron_3} (faut-il
detaille les requirements?) TriangularPolyhedralSurface?
OutputIt est un concept de stl Output Iterator avec \ccc{value_type
CGAL::Ridge_line*} .\\
ou Outputit est un Output Iterator de la stl sur des Umbilic*
\subsection{Output}
%%%%%%%%%%%%%%%%%%%%
Classes \ccc{Umbilic} and \ccc{Ridge_line}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Examples}

2
Ridges_3/doc_tex/Ridges_3/main.tex
View File

@ -1,4 +1,4 @@
\ccUserChapter{Ridge extraction on meshes}
\ccUserChapter{Ridges and Umbilics Extraction on Meshes}
\label{chap:Ridges_3}
\ccChapterAuthor{Marc Pouget and Frédéric Cazals}

3
Ridges_3/examples/Ridges_3/README
View File

@ -28,6 +28,7 @@ Note : if the nb of collected points is less than the required min number of
blind -f data/models/aim@shape/brain.off -d4 -m4 -a3 -t4
@ -39,6 +40,8 @@ Note : if the nb of collected points is less than the required min number of
visu with:
ridge_introspect-qt data/models/mailleurNew/surfaces/off/blend_c.00005.off data/data_models_mailleurNew_surfaces_off_blend_c.00005.offRIDGES-d4-m4-t3-a3-p0.4ogl.txt
*
introspect-qt ../../examples/Ridges_3/data/poly2x^2+y^2-0.062500.off ../../examples/Ridges_3/data_poly2x^2+y^2-0.062500.off.4ogl.txt

37
Ridges_3/examples/Ridges_3/blind.cpp
View File

@ -1,3 +1,4 @@
#include <CGAL/basic.h>
#include <CGAL/Cartesian.h>
@ -11,6 +12,7 @@
#include <CGAL/Ridges.h>
#include <CGAL/Umbilic.h>
// #include <CGAL/Monge_via_jet_fitting.h> //does not work since not in the release yet
// #include <CGAL/Lapack/Linear_algebra_lapack.h>
#include "../../../Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h"
@ -159,7 +161,7 @@ void compute_differential_quantities(PolyhedralSurf& P, Poly_rings& poly_rings)
//exit if the nb of points is too small
if ( in_points.size() < min_nb_points )
{std::cerr << "Too few points to perform the fitting" << std::endl; exit(0);}
{std::cerr << "Too few points to perform the fitting" << std::endl; exit(1);}
//For Ridges we need at least 3rd order info
assert( d_monge >= 3);
@ -302,10 +304,10 @@ int main(int argc, char *argv[])
back_insert_iterator<std::vector<Ridge_line*> > ii(ridge_lines);
//Find BLUE_RIDGE, RED_RIDGE, CREST or all ridges
// ridge_approximation.compute_ridges(P, CGAL::BLUE_RIDGE, ii, tag_order);
// ridge_approximation.compute_ridges(P, CGAL::RED_RIDGE, ii, tag_order);
// ridge_approximation.compute_ridges(P, CGAL::CREST, ii, tag_order);
ridge_approximation.compute_all_ridges(P, ii, tag_order);
// ridge_approximation.compute_ridges(CGAL::BLUE_RIDGE, ii, tag_order);
// ridge_approximation.compute_ridges(CGAL::RED_RIDGE, ii, tag_order);
//ridge_approximation.compute_ridges(CGAL::CREST, ii, tag_order);
ridge_approximation.compute_all_ridges(ii, tag_order);
std::vector<Ridge_line*>::iterator iter_lines = ridge_lines.begin(),
@ -321,21 +323,18 @@ int main(int argc, char *argv[])
//---------------------------------------------------------------------------
// UMBILICS
//--------------------------------------------------------------------------
// Umbilic_approximation umbilic_approximation(P,
// vertex2k1_pm, vertex2k2_pm,
// vertex2d1_pm, vertex2d2_pm);
// std::vector<Umbilic*> umbilics;
// back_insert_iterator<std::vector<Umbilic*> > umb_it(umbilics);
// umbilic_approximation.compute(P, umb_it, umb_size);
Umbilic_approximation umbilic_approximation(P,
vertex2k1_pm, vertex2k2_pm,
vertex2d1_pm, vertex2d2_pm);
std::vector<Umbilic*> umbilics;
back_insert_iterator<std::vector<Umbilic*> > umb_it(umbilics);
umbilic_approximation.compute(umb_it, umb_size);
// std::vector<Umbilic*>::iterator iter_umb = umbilics.begin(),
// iter_umb_end = umbilics.end();
// // output
// std::cout << "nb of umbilics " << umbilics.size() << std::endl;
// for (;iter_umb!=iter_umb_end;iter_umb++)
// {
// std::cout << "umbilic type " << (*iter_umb)->umb_type << std::endl;
// }
std::vector<Umbilic*>::iterator iter_umb = umbilics.begin(),
iter_umb_end = umbilics.end();
// output
std::cout << "nb of umbilics " << umbilics.size() << std::endl;
for (;iter_umb!=iter_umb_end;iter_umb++) std::cout << **iter_umb;
return 1;
}

0
Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500-off → Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500.off
View File

5
Ridges_3/examples/Ridges_3/makefile
View File

@ -7,7 +7,7 @@
# Choose the right include file from the <cgalroot>/make directory.
include $(CGAL_MAKEFILE)
PROFOPT= -g
PROFOPT= -O3 #-g
CFLAGS=${PROFOPT}
#---------------------------------------------------------------------#
# compiler flags
@ -21,8 +21,7 @@ CXXFLAGS = \
-I../../../Jet_fitting_3/include \
$(CGAL_CXXFLAGS) \
$(LONG_NAME_PROBLEM_CXXFLAGS) \
$(LAPACK_INCDIRS) \
${CFLAGS}
$(LAPACK_INCDIRS)
#---------------------------------------------------------------------#
# linker flags

199
Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h
View File

@ -2,6 +2,7 @@
#define _POLYHEDRALSURF_NEIGHBORS_H_
#include <queue>
#include <algorithm>
#include <CGAL/basic.h>
CGAL_BEGIN_NAMESPACE
@ -14,38 +15,40 @@ CGAL_BEGIN_NAMESPACE
template < class TPoly > class T_Gate
{
public:
typedef typename TPoly::Traits::FT FT;
typedef typename TPoly::Traits::FT FT;
typedef typename TPoly::Traits::Vector_3 Vector_3;
typedef typename TPoly::Traits::Point_3 Point_3;
typedef typename TPoly::Vertex_handle Vertex_handle;
typedef typename TPoly::Halfedge_handle Halfedge_handle;
typedef typename TPoly::Traits::Point_3 Point_3;
typedef typename TPoly::Vertex_handle Vertex_handle;
typedef typename TPoly::Halfedge_handle Halfedge_handle;
T_Gate( Vertex_handle v, Halfedge_handle he);
FT& d() { return m_d;}
const FT d() const { return m_d;}
Halfedge_handle he() { return m_he;}
private:
double m_d;
FT m_d;
Halfedge_handle m_he;
public:
T_Gate( Vertex_handle v, Halfedge_handle he)
{
m_he = he;
Point_3 p0 = v->point(),
p1 = he->vertex()->point(),
p2 = he->next()->vertex()->point(),
p3 = he->prev()->vertex()->point();
Vector_3 p0p1 = p0 - p1,
p0p2 = p0 - p2,
p0p3 = p0 - p3;
FT d1 = p0p1*p0p1,
d2 = p0p2*p0p2,
d3 = p0p3*p0p3;
m_d = CGAL::sqrt( (std::max)( (std::max)(d1,d2), d3) );
}
FT& d() {return m_d;}
const FT d() const {return m_d;}
Halfedge_handle he() {return m_he;}
};
//////////////IMPLEMENTATION//////////////////////////
template < class TPoly >
T_Gate<TPoly>::T_Gate( Vertex_handle v, Halfedge_handle he)
: m_he(he)
{
Point_3 p0 = v->point(),
p1 = he->vertex()->point(),
p2 = he->next()->vertex()->point(),
p3 = he->prev()->vertex()->point();
Vector_3 p0p1 = p0 - p1,
p0p2 = p0 - p2,
p0p3 = p0 - p3;
FT d1 = p0p1*p0p1,
d2 = p0p2*p0p2,
d3 = p0p3*p0p3;
m_d = CGAL::sqrt( (std::max)( (std::max)(d1,d2), d3) );
}
//---------------------------------------------------------------------------
// functor for priority queue
// order so than the top element is the smallest in the queue
@ -68,35 +71,34 @@ struct compare_gates
template < class TPoly > class T_PolyhedralSurf_neighbors
{
public:
typedef typename TPoly::Traits::FT FT;
typedef typename TPoly::Traits::Vector_3 Vector_3;
typedef typename TPoly::Traits::Point_3 Point_3;
typedef typename TPoly::Vertex_handle Vertex_handle;
typedef typename TPoly::Halfedge_handle Halfedge_handle;
typedef typename TPoly::Traits::FT FT;
typedef typename TPoly::Traits::Vector_3 Vector_3;
typedef typename TPoly::Traits::Point_3 Point_3;
typedef typename TPoly::Vertex_handle Vertex_handle;
typedef typename TPoly::Halfedge_handle Halfedge_handle;
typedef typename TPoly::Halfedge_around_vertex_circulator
Halfedge_around_vertex_circulator;
typedef typename TPoly::Vertex_iterator Vertex_iterator;
typedef typename TPoly::Vertex_iterator Vertex_iterator;
typedef T_Gate<TPoly> Gate;
public:
T_PolyhedralSurf_neighbors(TPoly& P);
// vertex_neigh stores the vertex v and its 1Ring neighbors contour
// stores halfedges, oriented CW, following the 1Ring disk border
// OneRingSize is the max distance from v to its OneRing
// neighbors. (the tag is_visited is not mofified)
void compute_one_ring(Vertex_handle v,
std::vector<Vertex_handle> &vertex_neigh,
std::list<Halfedge_handle> &contour,
FT &OneRingSize);
std::vector<Vertex_handle> &vertex_neigh,
std::list<Halfedge_handle> &contour,
FT &OneRingSize);
// call compute_one_ring and expand the contour (circle of halfedges
// CW), vertex_neigh are vertices on and inside the contour (there
// tag is_visited is set to true, but reset to false at the end),
// size is such that gates with distance less than size*OneRingSize
// are processed
void compute_neighbors(Vertex_handle v,
std::vector<Vertex_handle> &vertex_neigh,
std::list<Halfedge_handle> &contour,
FT size);
std::vector<Vertex_handle> &vertex_neigh,
std::list<Halfedge_handle> &contour,
FT size);
//vertex tags is_visited are set to false
void reset_is_visited_map(std::vector<Vertex_handle> &vces);
@ -112,7 +114,6 @@ public:
};
//////////////IMPLEMENTATION//////////////////////////
//////////////////////////////////////////////////////
template < class TPoly >
T_PolyhedralSurf_neighbors < TPoly >::
T_PolyhedralSurf_neighbors(TPoly& P)
@ -142,16 +143,6 @@ compute_one_ring(Vertex_handle v,
he_circ++;
} while (he_circ != he_end);
/* //debug check if the contour is closed */
/* list_it itbc = contour.begin(), itec = contour.end(), h_cur, h_next; */
/* for (; itbc != itec; itbc++) */
/* { */
/* h_cur = itbc; */
/* if ( h_cur != (--contour.end()) ) {h_next = ++h_cur; h_cur--;} */
/* else h_next = contour.begin(); */
/* assert( (*h_cur)->vertex() == (*h_next)->opposite()->vertex() ); */
/* } */
//compute OneRingSize = distance(v, 1Ring)
OneRingSize = 0;
typename std::vector<Vertex_handle>::iterator itb = vertex_neigh.begin(),
@ -190,57 +181,39 @@ compute_neighbors(Vertex_handle v,
typename std::list<Halfedge_handle>::iterator itb = contour.begin(),
ite = contour.end();
for (; itb != ite; itb++) {
// cerr << "[" << (*itb) << ' '
if (!( (*itb)->is_border() )) GatePQ.push(Gate(v, *itb));
if (!( (*itb)->is_border() )) GatePQ.push(Gate(v, *itb));
}
// init d_max
// init d_current
Gate firstGate = GatePQ.top();
FT d_current = firstGate.d();
// main loop
// main loop
while ( !GatePQ.empty() && d_current <= d_max ) {
/* //debug check if the contour is closed */
/* typename std::list<Halfedge_handle>::iterator itbc = contour.begin(), */
/* itec = contour.end(), */
/* h_cur, h_next; */
/* for (; itbc != itec; itbc++) */
/* { */
/* h_cur = itbc; */
/* if ( h_cur != (--contour.end()) ) {h_next = ++h_cur; h_cur--;} */
/* else h_next = contour.begin(); */
/* assert( (*h_cur)->vertex() == (*h_next)->opposite()->vertex() ); */
/* //cout << endl << &**itbc ; */
/* } */
/* //debug */
/* //cout << endl; cout << endl; */
Gate gate = GatePQ.top();
GatePQ.pop();
d_current = gate.d();
Halfedge_handle he = gate.he(), he1, he2;
//cerr << '\n' << &(*he) << '\n';
Vertex_handle v1;
// find the gate on the contour
// find the gate on the contour
typename std::list<Halfedge_handle>::iterator pos_he, pos_prev, pos_next, iter;
//pos_he = find(contour.begin(), contour.end(), he);
itb = contour.begin();
ite = contour.end();
for (; itb != ite; itb++)
{
// cout //<< endl << " *he = " << *he << " **itb = " << *(*itb);
// << endl << " he = " << he << " *itb = " << (*itb)
// << endl << "&*he = " << &(*he) << " &**itb = " << &(*(*itb));
if ( &(*(*itb)) == &(*he) )
{//pos_he = itb;
break;}
}
pos_he = itb;
pos_he = find(contour.begin(), contour.end(), he);
iter = pos_he;
// if the gate is not encoutered on the contour (case 3)
/**
there are different cases to expand the contour :
(case 3) he is not on the contour, nothing to do
(case 2) he is on the contour and either the previous or the next
following edge in the triangle is also on the contour, then delete
these 2 he from the contour and add the third one to the contour
and the PQ.
(case1) the vertex opposite to he is not visited, then the he is removed
from the contour, the two others are added to the contour and PQ, the
vertex is set visited.
*/
// if the gate is not encountered on the contour (case 3)
if ( pos_he == contour.end() ) continue;
// simulate a circulator on the contour:
// find the prev and next pos on coutour
// simulate a circulator on the contour:
// find the prev and next pos on coutour
if ( (++iter) != ite ) pos_next = iter;
else pos_next = contour.begin();
iter = pos_he;
@ -269,39 +242,27 @@ compute_neighbors(Vertex_handle v,
if ( !(he1->is_border()) ) GatePQ.push(Gate(v, he1));
continue;
}
v1 = he->next()->vertex();
if ( !is_visited_map.find(v1)->second )
{ // case 1
//vertex
is_visited_map.find(v1)->second = true;
vertex_neigh.push_back(v1);
//contour
he1 = he->prev()->opposite();
he2 = he->next()->opposite();
contour.insert(pos_he, he1);
contour.insert(pos_he, he2);
contour.erase(pos_he);
//GatePQ
if ( !(he1->is_border()) ) GatePQ.push(Gate(v, he1));
if ( !(he2->is_border()) ) GatePQ.push(Gate(v, he2));
continue;
}
//else case non admissible
v1 = he->next()->vertex();
if ( !is_visited_map.find(v1)->second )
{ // case 1
//vertex
is_visited_map.find(v1)->second = true;
vertex_neigh.push_back(v1);
//contour
he1 = he->prev()->opposite();
he2 = he->next()->opposite();
contour.insert(pos_he, he1);
contour.insert(pos_he, he2);
contour.erase(pos_he);
//GatePQ
if ( !(he1->is_border()) ) GatePQ.push(Gate(v, he1));
if ( !(he2->is_border()) ) GatePQ.push(Gate(v, he2));
continue;
}
//else case non admissible
CGAL_postcondition ( "case non admissible" == 0);
}// end while
/* //debug check if the contour is closed */
/* typename std::list<Halfedge_handle>::iterator itbc = contour.begin(), */
/* itec = contour.end(), */
/* h_cur, h_next; */
/* for (; itbc != itec; itbc++) */
/* { */
/* h_cur = itbc; */
/* if ( h_cur != (--contour.end()) ) {h_next = ++h_cur; h_cur--;} */
/* else h_next = contour.begin(); */
/* assert( (*h_cur)->vertex() == (*h_next)->opposite()->vertex() ); */
/* } */
/* //debug */
reset_is_visited_map(vertex_neigh);
}

167
Ridges_3/include/CGAL/Ridges.h
View File

@ -4,7 +4,6 @@
#include <utility>
#include <list>
#include <map>
#include <boost/property_map.hpp>
#include <CGAL/basic.h>
#include <CGAL/Min_sphere_d.h>
@ -22,7 +21,10 @@ enum Ridge_type {NONE=0, BLUE_RIDGE, RED_RIDGE, CREST,
//---------------------------------------------------------------------------
//Ridge_line : a connected sequence of edges of a Poly crossed by a
//ridge (with a barycentric coordinate to compute the crossing point),
//with type and weigths.
//with a Ridge_type and weigths : strength and sharpness. Note
//sharpness is only available (more precisely only meaningful : for
//Tag_3 it keeps its initial value 0) is the Ridge_approximation has
//been computed with the Tag_order Tag_4.
//--------------------------------------------------------------------------
template < class Poly > class Ridge_line
{
@ -33,15 +35,6 @@ public:
typedef typename Poly::Halfedge_handle Halfedge_handle;
typedef std::pair< Halfedge_handle, FT> ridge_he;
protected:
//one of BLUE_ELLIPTIC_RIDGE, BLUE_HYPERBOLIC_RIDGE, BLUE_CREST,
//RED_ELLIPTIC_RIDGE, RED_HYPERBOLIC_RIDGE or RED_CREST
Ridge_type m_line_type;
std::list<ridge_he> m_line;
FT m_strength;
FT m_sharpness;
public:
const Ridge_type line_type() const {return m_line_type;}
Ridge_type& line_type() {return m_line_type;}
@ -57,11 +50,22 @@ public:
//constructor
Ridge_line();
/* The output is :
line_type, strength, sharpness, list of points of the polyline.
/* The output is : line_type, strength, sharpness, list of points of
the polyline. An insert operator << is also available.
*/
void dump_4ogl(std::ostream& out_stream) ;
void dump_4ogl(std::ostream& out_stream) const ;
void dump_verbose(std::ostream& out_stream) const ;
protected:
//one of BLUE_ELLIPTIC_RIDGE, BLUE_HYPERBOLIC_RIDGE, BLUE_CREST,
//RED_ELLIPTIC_RIDGE, RED_HYPERBOLIC_RIDGE or RED_CREST
Ridge_type m_line_type;
std::list<ridge_he> m_line;
FT m_strength;// = integral of ppal curvature along the line
FT m_sharpness;// = (integral of second derivative of curvature
// along the line) divided by the size of the model
// (which is the radius of the smallest enclosing
// ball)
};
//--------------------------------------------------------------------------
@ -76,13 +80,13 @@ Ridge_line() : m_strength(0.), m_sharpness(0.) {}
template < class Poly >
void Ridge_line<Poly>::
dump_4ogl(std::ostream& out_stream)
dump_4ogl(std::ostream& out_stream) const
{
out_stream << line_type() << " "
<< strength() << " "
<< sharpness() << " ";
typename std::list<ridge_he >::iterator
typename std::list<ridge_he >::const_iterator
iter = line()->begin(),
ite = line()->end();
for (;iter!=ite;iter++){
@ -130,7 +134,7 @@ operator<<(std::ostream& out_stream, const Ridge_line<Poly>& ridge_line)
//---------------------------------------------------------------------------
//Differential_quantities
//--------------------------------------------------------------------------
template < class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap>
template <class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap>
class Differential_quantities {
public :
Vertex2FTPropertyMap vertex2k1_pm, vertex2k2_pm,
@ -144,7 +148,6 @@ template < class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap>
//---------------------------------------------------------------------------
//Ridge_approximation
//--------------------------------------------------------------------------
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
class Ridge_approximation
{
@ -161,30 +164,27 @@ template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2
//are ridges tagged as elliptic or hyperbolic using 3rd or 4th order
//differential quantitities?
//tag_3 is not working or badly, to be removed or improved?
//with tag_3 P1 and P2 are not used and the sharpness is not defined.
enum Tag_order {Tag_3 = 3, Tag_4 = 4};
public:
Ridge_approximation(Poly &P,
Vertex2FTPropertyMap vertex2k1_pm, Vertex2FTPropertyMap vertex2k2_pm,
Vertex2FTPropertyMap vertex2b0_pm, Vertex2FTPropertyMap vertex2b3_pm,
Vertex2FTPropertyMap vertex2P1_pm, Vertex2FTPropertyMap vertex2P2_pm,
Vertex2VectorPropertyMap vertex2d1_pm, Vertex2VectorPropertyMap vertex2d2_pm);
OutputIt compute_all_ridges(Poly &P, OutputIt it, Tag_order ord = Tag_3);
OutputIt compute_all_ridges(OutputIt it, Tag_order ord = Tag_3);
//Find BLUE_RIDGE, RED_RIDGE or CREST ridges
//iterate on P facets, find a non-visited, regular, 2Xing triangle,
//follow non-visited, regular, 2Xing triangles in both sens to create
//a Ridge line.
//Each time an edge is added the strength of the current line is updated
// + length(ridge segment in the facet)*|k|
void compute_ridges(Poly &P,
Ridge_type r_type,
//Find BLUE_RIDGE, RED_RIDGE or CREST ridges iterate on P facets,
//find a non-visited, regular (i.e. if there is a coherent
//orientation of ppal dir at the facet vertices), 2Xing triangle,
//follow non-visited, regular, 2Xing triangles in both sens to
//create a Ridge line. Each time an edge is added the strength and
//sharpness(if Tag_4) are updated.
void compute_ridges(Ridge_type r_type,
OutputIt ridge_lines_it,
Tag_order ord = Tag_3);
// void compute_umbilics(Poly &P );//container, class for umbilics?
protected:
protected:
Poly* P;
FT model_size;//radius of the smallest enclosing sphere of the Poly
//used to make the sharpness scale independant and iso indep
@ -209,7 +209,7 @@ protected:
Halfedge_handle& he1,
Halfedge_handle& he2,
Ridge_type r_type,
Tag_order ord = Tag_3);
Tag_order ord);
//is an edge crossed by a BLUE/RED ridge? (color is BLUE_RIDGE or
//RED_RIDGE ). As we only test edges of regular triangles, the ppal
@ -224,7 +224,7 @@ protected:
bool& is_crossed,
Ridge_type color);
//for the computation with tag_order = 3
//for the computation with tag_order == 3 only
//for a ridge segment [r1,r2] in a triangle (v1,v2,v3), let r = r2 -
//r1 and normalize, the projection of a point p on the line (r1,r2)
//is pp=r1+tr, with t=(p-r1)*r then the vector v starting at p is
@ -243,16 +243,19 @@ protected:
Vector_3 r1, Vector_3 r2,
Ridge_type color);
//a ridge line begins with a segment in a triangle
//a ridge line begins with a segment in a triangle given by the 2 he
//crossed
void init_ridge_line(Ridge_line* ridge_line,
Halfedge_handle h1, Halfedge_handle h2,
Ridge_type r_type);
Ridge_type r_type,
Tag_order ord);
//When the line is extended with a he, the bary coord of the
//crossing point is computed, the pair (he,coord) is added and the
//weigths are updated
void addback(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type);
void addfront(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type);
void addback(Ridge_line* ridge_line, Halfedge_handle he,
Ridge_type r_type, Tag_order ord);
void addfront(Ridge_line* ridge_line, Halfedge_handle he,
Ridge_type r_type, Tag_order ord);
//compute the barycentric coordinate of the xing point (blue or red)
//for he: p->q coord is st xing_point = coord*p + (1-coord)*q
@ -273,38 +276,42 @@ Ridge_approximation(Poly &P,
: P(&P), k1(vertex2k1_pm), k2(vertex2k2_pm), b0(vertex2b0_pm), b3(vertex2b3_pm),
P1(vertex2P1_pm), P2(vertex2P2_pm), d1(vertex2d1_pm), d2(vertex2d2_pm)
{
//init the is_visited_map and check that the mesh is a triangular one.
Facet_iterator itb = P.facets_begin(), ite = P.facets_end();
for(;itb!=ite;itb++) {
is_visited_map[itb] = false;
CGAL_precondition( itb->is_triangle() );
}
CGAL::Min_sphere_d<CGAL::Optimisation_d_traits_3<typename Poly::Traits> >
min_sphere(P.points_begin(), P.points_end());
model_size = min_sphere.squared_radius();
//maybe better to use CGAL::Min_sphere_of_spheres_d ?? but need to create spheres?
//init the is_visited_map
Facet_iterator itb = P.facets_begin(), ite = P.facets_end();
for(;itb!=ite;itb++) is_visited_map[itb] = false;
}
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
OutputIt Ridge_approximation<Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap>::
compute_all_ridges(Poly &P, OutputIt it, Tag_order ord)
compute_all_ridges(OutputIt it, Tag_order ord)
{
compute_ridges(P, BLUE_RIDGE, it, ord);
compute_ridges(P, RED_RIDGE, it, ord);
compute_ridges(P, CREST, it, ord);
compute_ridges(BLUE_RIDGE, it, ord);
compute_ridges(RED_RIDGE, it, ord);
compute_ridges(CREST, it, ord);
return it;
}
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
void Ridge_approximation<Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap>::
compute_ridges(Poly &P, Ridge_type r_type,
OutputIt ridge_lines_it, Tag_order ord)
compute_ridges(Ridge_type r_type, OutputIt ridge_lines_it, Tag_order ord)
{
//set all facets non visited
Facet_iterator itb = P.facets_begin(), ite = P.facets_end();
CGAL_precondition( (r_type == BLUE_RIDGE) || (r_type == RED_RIDGE) || (r_type == CREST) );
//reinit the is_visited_map
Facet_iterator itb = P->facets_begin(), ite = P->facets_end();
for(;itb!=ite;itb++) is_visited_map[itb] = false;
itb = P.facets_begin();
itb = P->facets_begin();
for(;itb!=ite;itb++)
{
Facet_handle f = itb;
@ -321,7 +328,7 @@ compute_ridges(Poly &P, Ridge_type r_type,
//a ridge_line is begining and stored
Ridge_line* cur_ridge_line = new Ridge_line();
init_ridge_line(cur_ridge_line, h1, h2, cur_ridge_type);
init_ridge_line(cur_ridge_line, h1, h2, cur_ridge_type, ord);
*ridge_lines_it++ = cur_ridge_line;
//next triangle adjacent to h1 (push_front)
@ -337,7 +344,7 @@ compute_ridges(Poly &P, Ridge_type r_type,
is_visited_map.find(f)->second = true;
if (curhe->opposite() == curhe1) curhe = curhe2;
else curhe = curhe1;//curhe stays at the ridge extremity
addfront(cur_ridge_line, curhe, cur_ridge_type);
addfront(cur_ridge_line, curhe, cur_ridge_type, ord);
if ( !(curhe->is_border_edge()) ) f =
curhe->opposite()->facet();
else break;
@ -361,7 +368,7 @@ compute_ridges(Poly &P, Ridge_type r_type,
is_visited_map.find(f)->second = true;
if (curhe->opposite() == curhe1) curhe = curhe2;
else curhe = curhe1;
addback(cur_ridge_line, curhe, cur_ridge_type);
addback(cur_ridge_line, curhe, cur_ridge_type, ord);
if ( !(curhe->is_border_edge()) ) f =
curhe->opposite()->facet();
else break;
@ -438,9 +445,9 @@ facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle&
he2 = h3;
}
//check there is no other case (just on edge crossed)
assert ( !( (h1_is_crossed && !h2_is_crossed && !h3_is_crossed)
|| (!h1_is_crossed && h2_is_crossed && !h3_is_crossed)
|| (!h1_is_crossed && h2_is_crossed && !h3_is_crossed)) );
CGAL_postcondition ( !( (h1_is_crossed && !h2_is_crossed && !h3_is_crossed)
|| (!h1_is_crossed && h2_is_crossed && !h3_is_crossed)
|| (!h1_is_crossed && h2_is_crossed && !h3_is_crossed)) );
//There is a ridge segment in the triangle, determine its type
Vertex_handle v_p1 = he1->opposite()->vertex(), v_q1 = he1->vertex(),
@ -493,7 +500,7 @@ facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle&
if ( sign_P2 < 0 ) return RED_ELLIPTIC_RIDGE; else return RED_HYPERBOLIC_RIDGE;
}
}
assert(0);//should return before!
CGAL_postcondition ( "should return before!" == 0);//should return before!
}
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
@ -546,9 +553,9 @@ template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2
}
if ( r != CGAL::NULL_VECTOR ) r = r/CGAL::sqrt(r*r);
FT sign1, sign2, sign3;
sign1 = (r1 - (v1->point()-ORIGIN) + (((v1->point()-ORIGIN)-r1)*r)*r )*dv1;
sign2 = (r1 - (v2->point()-ORIGIN) + (((v2->point()-ORIGIN)-r1)*r)*r )*dv2;
sign3 = (r1 - (v3->point()-ORIGIN) + (((v3->point()-ORIGIN)-r1)*r)*r )*dv3;
sign1 = bv1*(r1 - (v1->point()-ORIGIN) + (((v1->point()-ORIGIN)-r1)*r)*r )*dv1;
sign2 = bv2*(r1 - (v2->point()-ORIGIN) + (((v2->point()-ORIGIN)-r1)*r)*r )*dv2;
sign3 = bv3*(r1 - (v3->point()-ORIGIN) + (((v3->point()-ORIGIN)-r1)*r)*r )*dv3;
int compt = 0;
if ( sign1 > 0 ) compt++; else if (sign1 < 0) compt--;
@ -563,16 +570,18 @@ template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2
void Ridge_approximation<Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap>::
init_ridge_line(Ridge_line* ridge_line,
Halfedge_handle h1, Halfedge_handle h2,
Ridge_type r_type)
Ridge_type r_type,
Tag_order ord)
{
ridge_line->line_type() = r_type;
ridge_line->line()->push_back(ridge_he(h1, bary_coord(h1,r_type)));
addback(ridge_line, h2, r_type);
addback(ridge_line, h2, r_type, ord);
}
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
void Ridge_approximation<Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap>::
addback(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type)
addback(Ridge_line* ridge_line, Halfedge_handle he,
Ridge_type r_type, Tag_order ord)
{
Halfedge_handle he_cur = ( --(ridge_line->line()->end()) )->first;
FT coord_cur = ( --(ridge_line->line()->end()) )->second;//bary_coord(he_cur);
@ -583,7 +592,7 @@ addback(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type)
((v_p_cur->point()-ORIGIN)*coord_cur + (v_q_cur->point()-ORIGIN)*(1-coord_cur));
FT k1x, k2x; //abs value of the ppal curvatures at the Xing point on he.
FT k_second; // abs value of the second derivative of the curvature
FT k_second = 0; // abs value of the second derivative of the curvature
// along the line of curvature
k1x = CGAL::abs(k1[v_p]) * coord + CGAL::abs(k1[v_q]) * (1-coord) ;
k2x = CGAL::abs(k2[v_p]) * coord + CGAL::abs(k2[v_q]) * (1-coord) ;
@ -592,15 +601,19 @@ addback(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type)
|| (ridge_line->line_type() == BLUE_HYPERBOLIC_RIDGE)
|| (ridge_line->line_type() == BLUE_CREST) ) {
ridge_line->strength() += k1x * CGAL::sqrt(segment * segment);
k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size;
if (ord == Tag_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; }
}
if ( (ridge_line->line_type() == RED_ELLIPTIC_RIDGE)
|| (ridge_line->line_type() == RED_HYPERBOLIC_RIDGE)
|| (ridge_line->line_type() == RED_CREST) ) {
ridge_line->strength() += k2x * CGAL::sqrt(segment * segment);
k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size;
if (ord == Tag_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; }
}
ridge_line->line()->push_back( ridge_he(he, coord));
}
@ -609,7 +622,8 @@ addback(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type)
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
void Ridge_approximation<Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap>::
addfront(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type)
addfront(Ridge_line* ridge_line, Halfedge_handle he,
Ridge_type r_type, Tag_order ord)
{
Halfedge_handle he_cur = ( ridge_line->line()->begin() )->first;
FT coord_cur = ( ridge_line->line()->begin() )->second;
@ -629,15 +643,19 @@ addfront(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type)
|| (ridge_line->line_type() == BLUE_HYPERBOLIC_RIDGE)
|| (ridge_line->line_type() == BLUE_CREST) ) {
ridge_line->strength() += k1x * CGAL::sqrt(segment * segment);
k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size;
if (ord == Tag_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; }
}
if ( (ridge_line->line_type() == RED_ELLIPTIC_RIDGE)
|| (ridge_line->line_type() == RED_HYPERBOLIC_RIDGE)
|| (ridge_line->line_type() == RED_CREST) ) {
ridge_line->strength() += k2x * CGAL::sqrt(segment * segment);
k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size;
if (ord == Tag_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; }
}
ridge_line->line()->push_front( ridge_he(he, coord));
}
@ -660,6 +678,7 @@ bary_coord(Halfedge_handle he, Ridge_type r_type)
b_p = b3[he->opposite()->vertex()];
b_q = b3[he->vertex()];
}
//denominator cannot be 0 since there is no crossing when both extremalities are 0
return CGAL::abs(b_q) / ( CGAL::abs(b_q) + CGAL::abs(b_p) );
}

137
Ridges_3/include/CGAL/Umbilic.h
View File

@ -1,65 +1,104 @@
#ifndef _UMBILIC_H_
#define _UMBILIC_H_
#include <list>
#include <vector>
#include <math.h>
#include <CGAL/basic.h>
#include <CGAL/PolyhedralSurf_neighbors.h>
//#include <CGAL/number_utils_classes.h>
CGAL_BEGIN_NAMESPACE
enum Umbilic_type { NON_GENERIC = 0, WEDGE, TRISECTOR};
enum Umbilic_type { UMBILIC_NON_GENERIC = 0, UMBILIC_WEDGE, UMBILIC_TRISECTOR};
//-------------------------------------------------------------------
// Umbilic : stores umbilic data
//Umbilic : stores umbilic data, its location given by a vertex, its
//type and a circle of edges bording a disk containing the vertex
//------------------------------------------------------------------
template < class Poly >
class Umbilic
{
public:
typedef typename Poly::Vertex_handle Vertex_handle;
typedef typename Poly::Halfedge_handle Halfedge_handle;
public:
typedef typename Poly::Vertex_handle Vertex_handle;
typedef typename Poly::Halfedge_handle Halfedge_handle;
typedef typename Poly::Traits::Vector_3 Vector_3;
public:
//contructor
Umbilic(Vertex_handle v_init,
std::list<Halfedge_handle> contour_init);
//access fct
const Vertex_handle vertex() const { return v;}
const Umbilic_type umbilic_type() const { return umb_type;}
Umbilic_type& umbilic_type() { return umb_type;}
const std::list<Halfedge_handle> contour_list() const { return contour;}
std::list<Halfedge_handle>& contour_list() { return contour;}
protected:
Vertex_handle v;
Umbilic_type umb_type;
std::list<Halfedge_handle> contour;
//contructor
Umbilic(Vertex_handle v_init,
std::list<Halfedge_handle> contour_init)
: v(v_init), contour(contour_init) {};
};
//contructor
template <class Poly>
Umbilic<Poly>::
Umbilic(Vertex_handle v_init,
std::list<Halfedge_handle> contour_init)
: v(v_init), contour(contour_init) {}
template <class Poly>
std::ostream&
operator<<(std::ostream& out_stream, const Umbilic<Poly>& umbilic)
{
out_stream << "Umbilic at location (" << umbilic.vertex()->point() << ") of type " ;
switch (umbilic.umbilic_type())
{
case CGAL::UMBILIC_NON_GENERIC: out_stream << "non generic" << std::endl; break;
case CGAL::UMBILIC_WEDGE: out_stream << "wedge" << std::endl; break;
case CGAL::UMBILIC_TRISECTOR: out_stream << "trisector" << std::endl; break;
default : out_stream << "Something wrong appends for sure..." << std::endl; break;
}
return out_stream;
}
//---------------------------------------------------------------------------
//Umbilic_approximation
//Umbilic_approximation : enable computation of umbilics of a
//TriangularPolyhedralSurface. It uses the class
//T_PolyhedralSurf_neighbors to compute topological disk patches
//around vertices
//--------------------------------------------------------------------------
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
class Umbilic_approximation
{
public:
typedef typename Poly::Traits::FT FT;
public:
typedef typename Poly::Traits::FT FT;
typedef typename Poly::Traits::Vector_3 Vector_3;
typedef typename Poly::Vertex_handle Vertex_handle;
typedef typename Poly::Halfedge_handle Halfedge_handle;
typedef typename Poly::Facet_handle Facet_handle;
typedef typename Poly::Facet_iterator Facet_iterator;
typedef typename Poly::Vertex_iterator Vertex_iterator;
typedef typename Poly::Vertex_handle Vertex_handle;
typedef typename Poly::Halfedge_handle Halfedge_handle;
typedef typename Poly::Facet_handle Facet_handle;
typedef typename Poly::Facet_iterator Facet_iterator;
typedef typename Poly::Vertex_iterator Vertex_iterator;
typedef Umbilic<Poly> Umbilic;
static FT neigh_size;//the size of neighbourhood for umbilic
// computation is (neigh_size * OneRingSize)
//constructor : sets propertymaps and the poly_neighbors
Umbilic_approximation(Poly &P,
Vertex2FTPropertyMap vertex2k1_pm, Vertex2FTPropertyMap vertex2k2_pm,
Vertex2VectorPropertyMap vertex2d1_pm, Vertex2VectorPropertyMap vertex2d2_pm);
OutputIt compute(Poly &P, OutputIt it, FT size);
//identify umbilics as vertices minimizing the function k1-k2 on
//their patch and for which the index is not 0. We avoid
//potential umbilics whose contours touch the border.
OutputIt compute(OutputIt it, FT size);
protected:
Poly* P;
typedef T_PolyhedralSurf_neighbors<Poly> Poly_neighbors;
Poly_neighbors* poly_neighbors;
CGAL::Abs<FT> cgal_abs;
To_double<FT> To_double;
//Property maps
Vertex2FTPropertyMap k1, k2;
@ -69,9 +108,10 @@ public:
// max dir of an arbitrary starting point, the max dir field is
// oriented on the next point so that the scalar product of the
// consecutive vectors is positive. Adding all the angles between
// consecutive vectors around the contour gives -/+180 for a
// wedge/trisector
void compute_type(Umbilic& umb);
// consecutive vectors around the contour gives ~ -/+180 for a
// wedge/trisector, ~ 0 gives a false umbilic, everything else gives
// a non_generic umbilic.
int compute_type(Umbilic& umb);
};
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
@ -79,24 +119,30 @@ template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2
Umbilic_approximation(Poly &P,
Vertex2FTPropertyMap vertex2k1_pm, Vertex2FTPropertyMap vertex2k2_pm,
Vertex2VectorPropertyMap vertex2d1_pm, Vertex2VectorPropertyMap vertex2d2_pm)
: k1(vertex2k1_pm), k2(vertex2k2_pm),
: P(&P), k1(vertex2k1_pm), k2(vertex2k2_pm),
d1(vertex2d1_pm), d2(vertex2d2_pm)
{
//check that the mesh is a triangular one.
Facet_iterator itb = P.facets_begin(), ite = P.facets_end();
for(;itb!=ite;itb++) CGAL_precondition( itb->is_triangle() );
poly_neighbors = new Poly_neighbors(P);
}
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
OutputIt Umbilic_approximation< Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap >::
compute(Poly &P, OutputIt umbilics_it, FT size)
compute(OutputIt umbilics_it, FT size)
{
CGAL_precondition( size >= 1 );
std::vector<Vertex_handle> vces;
std::list<Halfedge_handle> contour;
double umbilicEstimatorVertex, umbilicEstimatorNeigh;
FT umbilicEstimatorVertex, umbilicEstimatorNeigh;
bool is_umbilic = true;
//MAIN loop on P vertices
Vertex_iterator itb = P.vertices_begin(), ite = P.vertices_end();
Vertex_iterator itb = P->vertices_begin(), ite = P->vertices_end();
for (;itb != ite; itb++) {
Vertex_handle vh = itb;
umbilicEstimatorVertex = cgal_abs(k1[vh]-k2[vh]);
@ -104,10 +150,12 @@ compute(Poly &P, OutputIt umbilics_it, FT size)
vces.clear();
contour.clear();
is_umbilic = true;
//the size of neighbourhood is (size * OneRingSize)
poly_neighbors->compute_neighbors(vh, vces, contour, size);
// OPTIONAL: avoid umbilics whose contours touch the border
// avoid umbilics whose contours touch the border (Note may not be
// desirable?)
typename std::list<Halfedge_handle>::iterator itb_cont = contour.begin(),
ite_cont = contour.end();
for (; itb_cont != ite_cont; itb_cont++)
@ -119,7 +167,6 @@ compute(Poly &P, OutputIt umbilics_it, FT size)
// neigh vertex has a lower umbilicEstimator value
typename std::vector<Vertex_handle>::iterator itbv = vces.begin(),
itev = vces.end();
assert(*itbv == vh);
itbv++;
for (; itbv != itev; itbv++)
{ umbilicEstimatorNeigh = cgal_abs( k1[*itbv] - k2[*itbv] );
@ -128,24 +175,24 @@ compute(Poly &P, OutputIt umbilics_it, FT size)
}
if (is_umbilic == false) continue;
//v is an umbilic, compute the index
//v is an umbilic (wrt the min of k1-k2), compute the index. If
//the index is not 0 then we have actually an umbilic which is output
Umbilic* cur_umbilic = new Umbilic(vh, contour);
compute_type(*cur_umbilic);
*umbilics_it++ = cur_umbilic;
if (compute_type(*cur_umbilic) != 0) *umbilics_it++ = cur_umbilic;
}
return umbilics_it;
}
template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap >
void Umbilic_approximation< Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap >::
int Umbilic_approximation< Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap >::
compute_type(Umbilic& umb)
{
Vector_3 dir, dirnext, normal;
double cosinus, angle=0, angleSum=0;
const double pi=3.141592653589793;
Vertex_handle v;
typename std::list<Halfedge_handle>::iterator itb = umb.contour.begin(),
itlast = --umb.contour.end();
typename std::list<Halfedge_handle>::iterator itb = umb.contour_list().begin(),
itlast = --umb.contour_list().end();
v = (*itb)->vertex();
dir = d1[v];
@ -156,7 +203,7 @@ compute_type(Umbilic& umb)
itb++;
v=(*itb)->vertex();
dirnext = d1[v];
cosinus = dir*dirnext;
cosinus = To_double(dir*dirnext);
if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;}
if (cosinus>1) cosinus = 1;
//orientation of (dir, dirnext, normal)
@ -164,23 +211,25 @@ compute_type(Umbilic& umb)
else angle = -acos(cosinus);
angleSum += angle;
dir = dirnext;
normal = CGAL::cross_product(d1[v], d2[v]);
normal = CGAL::cross_product(d1[v], d2[v]);
}
while (itb != (itlast));
//angle (v_last, v_0)
v=(*umb.contour.begin())->vertex();
v=(*umb.contour_list().begin())->vertex();
dirnext = d1[v];
cosinus = dir*dirnext;
cosinus = To_double(dir*dirnext);
if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;}
if (cosinus>1) cosinus = 1;
if ( (dir * CGAL::cross_product(dirnext, normal)) > 0) angle = acos(cosinus);
else angle = -acos(cosinus);
angleSum += angle;
if ((angleSum > (pi/2)) && (angleSum < (3*pi/2))) umb.umb_type = TRISECTOR ;
else if ((angleSum < (-pi/2)) && (angleSum > (-3*pi/2))) umb.umb_type = WEDGE;
else umb.umb_type = NON_GENERIC;
if ((angleSum > (pi/2)) && (angleSum < (3*pi/2))) umb.umbilic_type() = UMBILIC_TRISECTOR ;
else if ((angleSum < (-pi/2)) && (angleSum > (-3*pi/2))) umb.umbilic_type() = UMBILIC_WEDGE;
else if ((angleSum <= (pi/2)) && (angleSum >= (-pi/2))) return 0;//is not considered as an umbilic
else umb.umbilic_type() = UMBILIC_NON_GENERIC;
return 1;
}
CGAL_END_NAMESPACE