From 3b4f578ae61d1a42acd9d979fd3d03f7a376eeed Mon Sep 17 00:00:00 2001 From: Marc Pouget Date: Thu, 24 Aug 2006 12:54:27 +0000 Subject: [PATCH] CGAL_pre and post condition added --- .gitattributes | 2 +- Ridges_3/demo/Ridges_3/README | 2 + Ridges_3/demo/Ridges_3/SketchSample.cpp | 51 +++-- Ridges_3/demo/Ridges_3/SketchSample.h | 7 +- Ridges_3/demo/Ridges_3/makefile | 2 - Ridges_3/demo/Ridges_3/visu.cpp | 77 ++++--- Ridges_3/doc_tex/Ridges_3/Ridges_3_user.tex | 160 ++++++++------ Ridges_3/doc_tex/Ridges_3/main.tex | 2 +- Ridges_3/examples/Ridges_3/README | 5 +- Ridges_3/examples/Ridges_3/blind.cpp | 37 ++-- ...0.062500-off => poly2x^2+y^2-0.062500.off} | 0 Ridges_3/examples/Ridges_3/makefile | 5 +- .../include/CGAL/PolyhedralSurf_neighbors.h | 203 +++++++----------- Ridges_3/include/CGAL/Ridges.h | 169 ++++++++------- Ridges_3/include/CGAL/Umbilic.h | 139 ++++++++---- 15 files changed, 474 insertions(+), 387 deletions(-) rename Ridges_3/examples/Ridges_3/data/{poly2x^2+y^2-0.062500-off => poly2x^2+y^2-0.062500.off} (100%) diff --git a/.gitattributes b/.gitattributes index 3a36878bd39..eaaf3e38502 100644 --- a/.gitattributes +++ b/.gitattributes @@ -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 diff --git a/Ridges_3/demo/Ridges_3/README b/Ridges_3/demo/Ridges_3/README index 1e18cc24b23..0cbdfd1df27 100644 --- a/Ridges_3/demo/Ridges_3/README +++ b/Ridges_3/demo/Ridges_3/README @@ -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 + + diff --git a/Ridges_3/demo/Ridges_3/SketchSample.cpp b/Ridges_3/demo/Ridges_3/SketchSample.cpp index e4f4e9de86d..53e3657f297 100644 --- a/Ridges_3/demo/Ridges_3/SketchSample.cpp +++ b/Ridges_3/demo/Ridges_3/SketchSample.cpp @@ -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::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::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(); -} diff --git a/Ridges_3/demo/Ridges_3/SketchSample.h b/Ridges_3/demo/Ridges_3/SketchSample.h index 801d9aab16e..e4026c136ee 100644 --- a/Ridges_3/demo/Ridges_3/SketchSample.h +++ b/Ridges_3/demo/Ridges_3/SketchSample.h @@ -7,11 +7,16 @@ #include "Sketcher.h" #include +#include +#include +#include //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{ diff --git a/Ridges_3/demo/Ridges_3/makefile b/Ridges_3/demo/Ridges_3/makefile index 511c111c131..4e17ab42994 100644 --- a/Ridges_3/demo/Ridges_3/makefile +++ b/Ridges_3/demo/Ridges_3/makefile @@ -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 diff --git a/Ridges_3/demo/Ridges_3/visu.cpp b/Ridges_3/demo/Ridges_3/visu.cpp index 61510915a42..096419a924c 100644 --- a/Ridges_3/demo/Ridges_3/visu.cpp +++ b/Ridges_3/demo/Ridges_3/visu.cpp @@ -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 < jvw(out_jvx); - jvw.set_title("ridges"); - jvw.write_header(); +// std::ofstream out_jvx("debug.jvx"); +// CGAL::Javaview_writer jvw(out_jvx); +// jvw.set_title("ridges"); +// jvw.write_header(); -// //first the polysurf -// jvw.set_geometry_name("polysurf"); -// jvw.begin_geometry(); -// polyhedron_javaview_writer(jvw, P); -// jvw.end_geometry(); +// // //first the polysurf +// // jvw.set_geometry_name("polysurf"); +// // jvw.begin_geometry(); +// // polyhedron_javaview_writer(jvw, P); +// // jvw.end_geometry(); - int compt = 0; +// int compt = 0; - 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()); +// 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()); - //color - if ((*iter_lines)->ridge_type ==BC) jvw.set_color(CGAL::BLUE); - else jvw.set_color(CGAL::RED); +// //color +// if ((*iter_lines)->ridge_type == CGAL::BLUE_CREST) jvw.set_color(CGAL::BLUE); +// else jvw.set_color(CGAL::RED); - //lines - jvw.begin_geometry(); - polyline_javaview_writer(jvw, (*iter_lines)->ridge_points.begin(), - (*iter_lines)->ridge_points.end()); - jvw.end_geometry(); - } +// //lines +// jvw.begin_geometry(); +// polyline_javaview_writer(jvw, (*iter_lines)->ridge_points.begin(), +// (*iter_lines)->ridge_points.end()); +// jvw.end_geometry(); +// } - jvw.write_footer(); +// jvw.write_footer(); diff --git a/Ridges_3/doc_tex/Ridges_3/Ridges_3_user.tex b/Ridges_3/doc_tex/Ridges_3/Ridges_3_user.tex index 38e439aba00..e8016764fd4 100644 --- a/Ridges_3/doc_tex/Ridges_3/Ridges_3_user.tex +++ b/Ridges_3/doc_tex/Ridges_3/Ridges_3_user.tex @@ -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} +
+
+\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} +
+
+\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} diff --git a/Ridges_3/doc_tex/Ridges_3/main.tex b/Ridges_3/doc_tex/Ridges_3/main.tex index cdee29859c5..a41587c9c38 100644 --- a/Ridges_3/doc_tex/Ridges_3/main.tex +++ b/Ridges_3/doc_tex/Ridges_3/main.tex @@ -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} diff --git a/Ridges_3/examples/Ridges_3/README b/Ridges_3/examples/Ridges_3/README index 802528b79d1..02f0e1cd695 100644 --- a/Ridges_3/examples/Ridges_3/README +++ b/Ridges_3/examples/Ridges_3/README @@ -28,7 +28,8 @@ 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 + ./blind -f data/ellipe0.003.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 diff --git a/Ridges_3/examples/Ridges_3/blind.cpp b/Ridges_3/examples/Ridges_3/blind.cpp index 72946c5d46c..4e15677fd69 100644 --- a/Ridges_3/examples/Ridges_3/blind.cpp +++ b/Ridges_3/examples/Ridges_3/blind.cpp @@ -1,3 +1,4 @@ + #include #include @@ -11,6 +12,7 @@ #include #include + // #include //does not work since not in the release yet // #include #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 > 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::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 umbilics; -// back_insert_iterator > 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 umbilics; + back_insert_iterator > umb_it(umbilics); + umbilic_approximation.compute(umb_it, umb_size); -// std::vector::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::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; } diff --git a/Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500-off b/Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500.off similarity index 100% rename from Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500-off rename to Ridges_3/examples/Ridges_3/data/poly2x^2+y^2-0.062500.off diff --git a/Ridges_3/examples/Ridges_3/makefile b/Ridges_3/examples/Ridges_3/makefile index ac9657316bf..ca5f1190a96 100644 --- a/Ridges_3/examples/Ridges_3/makefile +++ b/Ridges_3/examples/Ridges_3/makefile @@ -7,7 +7,7 @@ # Choose the right include file from the /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 diff --git a/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h b/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h index 2d98405f278..a0f25f89f94 100644 --- a/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h +++ b/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h @@ -2,6 +2,7 @@ #define _POLYHEDRALSURF_NEIGHBORS_H_ #include +#include #include 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; -private: - double m_d; - Halfedge_handle m_he; + 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;} -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;} +private: + FT m_d; + Halfedge_handle m_he; }; +//////////////IMPLEMENTATION////////////////////////// +template < class TPoly > +T_Gate::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 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_neigh, - std::list &contour, - FT &OneRingSize); + std::vector &vertex_neigh, + std::list &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_neigh, - std::list &contour, - FT size); + std::vector &vertex_neigh, + std::list &contour, + FT size); //vertex tags is_visited are set to false void reset_is_visited_map(std::vector &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::iterator itb = vertex_neigh.begin(), @@ -190,57 +181,39 @@ compute_neighbors(Vertex_handle v, typename std::list::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::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::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::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); } diff --git a/Ridges_3/include/CGAL/Ridges.h b/Ridges_3/include/CGAL/Ridges.h index 10403858f92..409f23ea4ff 100644 --- a/Ridges_3/include/CGAL/Ridges.h +++ b/Ridges_3/include/CGAL/Ridges.h @@ -4,7 +4,6 @@ #include #include #include -#include #include #include @@ -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 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 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:: -dump_4ogl(std::ostream& out_stream) +dump_4ogl(std::ostream& out_stream) const { out_stream << line_type() << " " << strength() << " " << sharpness() << " "; - typename std::list::iterator + typename std::list::const_iterator iter = line()->begin(), ite = line()->end(); for (;iter!=ite;iter++){ @@ -130,7 +134,7 @@ operator<<(std::ostream& out_stream, const Ridge_line& ridge_line) //--------------------------------------------------------------------------- //Differential_quantities //-------------------------------------------------------------------------- -template < class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap> +template 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 > 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:: -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:: -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:: 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:: -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:: -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) ); } diff --git a/Ridges_3/include/CGAL/Umbilic.h b/Ridges_3/include/CGAL/Umbilic.h index 2268786d43c..cf1a4ca9d05 100644 --- a/Ridges_3/include/CGAL/Umbilic.h +++ b/Ridges_3/include/CGAL/Umbilic.h @@ -1,65 +1,104 @@ #ifndef _UMBILIC_H_ #define _UMBILIC_H_ +#include +#include #include #include #include +//#include + + 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; + + //contructor + Umbilic(Vertex_handle v_init, + std::list 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 contour_list() const { return contour;} + std::list& contour_list() { return contour;} -public: + protected: Vertex_handle v; Umbilic_type umb_type; std::list contour; - //contructor - Umbilic(Vertex_handle v_init, - std::list contour_init) - : v(v_init), contour(contour_init) {}; }; +//contructor +template +Umbilic:: +Umbilic(Vertex_handle v_init, + std::list contour_init) + : v(v_init), contour(contour_init) {} + +template +std::ostream& +operator<<(std::ostream& out_stream, const Umbilic& 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 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_neighbors; Poly_neighbors* poly_neighbors; CGAL::Abs cgal_abs; + To_double 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 vces; std::list 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::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::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::iterator itb = umb.contour.begin(), - itlast = --umb.contour.end(); + typename std::list::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