diff --git a/src/fgt/filter_qhull/filter_qhull.cpp b/src/fgt/filter_qhull/filter_qhull.cpp index 14cce0f52..59b535a4a 100644 --- a/src/fgt/filter_qhull/filter_qhull.cpp +++ b/src/fgt/filter_qhull/filter_qhull.cpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -45,6 +46,9 @@ using namespace std; using namespace vcg; +//Internal prototype +bool create_mesh_convex_hull(MeshDocument &md,MeshModel &m,FilterParameterSet & par); + // Constructor usually performs only two simple tasks of filling the two lists // - typeList: with all the possible id of the filtering actions // - actionList with the corresponding actions. If you want to add icons to your filtering actions you can do here by construction the QActions accordingly @@ -87,20 +91,26 @@ const QString QhullPlugin::filterName(FilterIDType filterId) const QString QhullPlugin::filterInfo(FilterIDType filterId) { switch(filterId) { - case FP_QHULL_CONVEX_HULL : return QString("Calculate mesh convex hull with Qhull library.
" - "The convex hull of a set of points is the smallest convex set containing the points."); - case FP_QHULL_DELAUNAY_TRIANGULATION : return QString("Calculate mesh Delaunay triangulations with Qhull library.
" - "The Delaunay triangulation of a set of points in d-dimensional spaces is " - "the projection of the convex hull of the projections of the points onto a " - "(d+1)-dimensional paraboloid.
"); - case FP_QHULL_VORONOI_FILTERING : return QString("Compute mesh Voronoi filtering with Qhull library.
" - "The filter reconstructs the mesh surface starting from its vertices through a double " - "Delaunay triangulation.
The second one takes in input some Voronoi vertices too.
"); - case FP_QHULL_ALPHA_COMPLEX_AND_SHAPE: return QString("Calculate Alpha Shape.
" - "The Alpha Shape is the boundary of the alpha complex, that is a subcomplex of the Delaunay triangulation.
" + case FP_QHULL_CONVEX_HULL : return QString("Calculate the convex hull with Qhull library (http://www.qhull.org/html/qconvex.htm).

" + "The convex hull of a set of points is the boundary of the minimal convex set containing the given non-empty finite set of points."); + case FP_QHULL_DELAUNAY_TRIANGULATION : return QString("Calculate the Delaunay triangulation with Qhull library (http://www.qhull.org/html/qdelaun.htm).

" + "The Delaunay triangulation of a set of points in d-dimensional spaces is a triangulation of the convex hull.
" + "The result is non 2-manifold by definition."); + case FP_QHULL_VORONOI_FILTERING : return QString("Compute a Voronoi filtering (Amenta and Bern 1998) with Qhull library (http://www.qhull.org/).

" + "The filter computes a piecewise-linear approximation of a smooth surface from a finite set of sample points." + "The algorithm uses Voronoi vertices to remove triangles from the Delaunay triangulation.
" + "After computing the Voronoi diagram, foreach sample point it chooses the two farthest opposite Voronoi vertices." + "Then computes a Delaunay triangulation of the sample points and the selected Voronoi vertices, and keep " + "only those triangles in witch all three vertices are sample points."); + case FP_QHULL_ALPHA_COMPLEX_AND_SHAPE: return QString("Calculate the Alpha Shape (Edelsbrunner and P.Mucke 1994) with Qhull library (http://www.qhull.org/).

" + "From a given finite point set in the space it computes 'the shape' of the set." + "The Alpha Shape is the boundary of the alpha complex, that is a subcomplex of the Delaunay triangulation of the given point set.
" "For a given value of 'alpha', the alpha complex includes all the simplices in the Delaunay " "triangulation which have an empty circumsphere with radius equal or smaller than 'alpha'"); - case FP_QHULL_VISIBLE_POINTS: return QString("Select points visible from camera viewpoint"); + case FP_QHULL_VISIBLE_POINTS: return QString("Select the visible points in a point cloud, as viewed from a given viewpoint.
" + "It uses the Qhull library (http://www.qhull.org/

" + "The algorithm used (Katz, Tal and Basri 2007) determines visibility without reconstructing a surface or estimating normals." + "The visible points are those that corresponds to the ones that reside on the convex hull of a transformed point cloud."); default : assert(0); } return QString("Error: Unknown Filter"); @@ -119,7 +129,7 @@ const QhullPlugin::FilterClass QhullPlugin::getClass(QAction *a) case FP_QHULL_ALPHA_COMPLEX_AND_SHAPE: return FilterClass (MeshFilterInterface::Remeshing) ; case FP_QHULL_VISIBLE_POINTS: - return FilterClass (MeshFilterInterface::Selection); + return FilterClass (MeshFilterInterface::Selection + MeshFilterInterface::PointSet); default : assert(0); } return FilterClass(0); @@ -137,6 +147,7 @@ void QhullPlugin::initParameterSet(QAction *action,MeshModel &m, FilterParameter switch(ID(action)) { case FP_QHULL_CONVEX_HULL : { + parlst.addBool("reorient", false,"Re-orient all faces coherentely","Re-orient all faces coherentely"); break; } case FP_QHULL_DELAUNAY_TRIANGULATION : @@ -145,6 +156,7 @@ void QhullPlugin::initParameterSet(QAction *action,MeshModel &m, FilterParameter } case FP_QHULL_VORONOI_FILTERING : { + parlst.addDynamicFloat("threshold",100.0f, 100.0f, 2000.0f,MeshModel::MM_NONE,"threshold ","????."); break; } case FP_QHULL_ALPHA_COMPLEX_AND_SHAPE: @@ -160,7 +172,7 @@ void QhullPlugin::initParameterSet(QAction *action,MeshModel &m, FilterParameter { parlst.addDynamicFloat("radiusThreshold", 0.0f, 0.0f, 7.0f,MeshModel::MM_VERTFLAGSELECT, - "radius threshold ","Bounds the radius oh the sphere used to select visible points.
" + "radius threshold ","Bounds the radius of the sphere used to select visible points.
" "Use a big threshold for dense point clouds, a small one for sparse clouds"); parlst.addBool ("usecamera", @@ -172,7 +184,11 @@ void QhullPlugin::initParameterSet(QAction *action,MeshModel &m, FilterParameter "ViewPoint", "if UseCamera is true, this value is ignored"); - parlst.addBool("convex_hull",false,"Show Partial Convex Hull of flipped points", "Show Partial Convex Hull of the transformed point cloud"); + parlst.addBool("convex_hullFP",false,"Show Partial Convex Hull of flipped points", "Show Partial Convex Hull of the transformed point cloud"); + parlst.addBool("convex_hullNFP",false,"Show the Convex Hull of the given points", "Show the Convex Hull of the given points"); + parlst.addBool("reorient", false,"Re-orient all faces of the CH coherentely","Re-orient all faces of the CH coherentely." + "If no Convex Hulls are selected , this value is ignored"); + break; break; } default : assert(0); @@ -187,79 +203,13 @@ bool QhullPlugin::applyFilter(QAction *filter, MeshDocument &md, FilterParameter case FP_QHULL_CONVEX_HULL : { MeshModel &m=*md.mm(); - MeshModel &pm =*md.addNewMesh("Convex Hull"); + bool result= create_mesh_convex_hull(md,m,par); - if (m.hasDataMask(MeshModel::MM_WEDGTEXCOORD)){ - m.clearDataMask(MeshModel::MM_WEDGTEXCOORD); - } - if (m.hasDataMask(MeshModel::MM_VERTTEXCOORD)){ - m.clearDataMask(MeshModel::MM_VERTTEXCOORD); - } - - int dim= 3; /* dimension of points */ - int numpoints= m.cm.vn; /* number of mesh vertices */ - - //facet_list contains the convex hull - //facet_list contains simplicial (triangulated) facets. - //The convex hull of a set of points is the smallest convex set containing the points. - - facetT *facet_list = compute_convex_hull(dim,numpoints,m); - - if(facet_list!=NULL){ - - int convexNumFaces = qh num_facets; - int convexNumVert = qh_setsize(qh_facetvertices (facet_list, NULL, false)); - assert( qh num_vertices == convexNumVert); - - tri::Allocator::AddVertices(pm.cm,convexNumVert); - tri::Allocator::AddFaces(pm.cm,convexNumFaces); - - - /*ivp length is 'numpoints' because each vertex is accessed through its ID whose range is - 0<=qh_pointid(vertex->point)::VertexPointer> ivp(numpoints); - vertexT *vertex; - int i=0; - FORALLvertices{ - if ((*vertex).point){ - pm.cm.vert[i].P()[0] = (*vertex).point[0]; - pm.cm.vert[i].P()[1] = (*vertex).point[1]; - pm.cm.vert[i].P()[2] = (*vertex).point[2]; - ivp[qh_pointid(vertex->point)] = &pm.cm.vert[i]; - i++; - } - } - - facetT *facet; - i=0; - FORALLfacet_(facet_list){ - vertexT *vertex; - int vertex_n, vertex_i; - FOREACHvertex_i_((*facet).vertices){ - pm.cm.face[i].V(vertex_i)= ivp[qh_pointid(vertex->point)]; - } - i++; - } - - assert( pm.cm.fn == convexNumFaces); - Log(GLLogStream::FILTER,"Successfully created a mesh of %i vert and %i faces",convexNumVert,convexNumFaces); - //m.cm.Clear(); - - vcg::tri::UpdateBounding::Box(pm.cm); - vcg::tri::UpdateNormals::PerVertexNormalizedPerFace(pm.cm); - - int curlong, totlong; /* memory remaining after qh_memfreeshort */ - qh_freeqhull(!qh_ALL); - qh_memfreeshort (&curlong, &totlong); - if (curlong || totlong) - fprintf (stderr, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", - totlong, curlong); - + if(result){ + Log(GLLogStream::FILTER,"Successfully created a mesh of %i vert and %i faces",m.cm.vn,m.cm.fn); return true; - } - else - return false; + else return false; } case FP_QHULL_DELAUNAY_TRIANGULATION: { @@ -356,8 +306,10 @@ bool QhullPlugin::applyFilter(QAction *filter, MeshDocument &md, FilterParameter int dim= 3; /* dimension of points */ int numpoints= m.cm.vn; /* number of mesh vertices */ + + float threshold = par.getDynamicFloat("threshold"); - bool result = compute_voronoi(dim,numpoints,m,pm); + bool result = compute_voronoi(dim,numpoints,m,pm,threshold); if(result){ vcg::tri::UpdateBounding::Box(pm.cm); @@ -433,6 +385,9 @@ bool QhullPlugin::applyFilter(QAction *filter, MeshDocument &md, FilterParameter { m.updateDataMask(MeshModel::MM_VERTFLAGSELECT); } + + //Clear old selection + tri::UpdateSelection::ClearVertex(m.cm); int dim= 3; /* dimension of points */ int numpoints= m.cm.vn; /* number of mesh vertices */ @@ -454,25 +409,43 @@ bool QhullPlugin::applyFilter(QAction *filter, MeshDocument &md, FilterParameter MeshModel &pm =*md.addNewMesh("CH Flipped Points"); - if (m.hasDataMask(MeshModel::MM_WEDGTEXCOORD)){ - m.clearDataMask(MeshModel::MM_WEDGTEXCOORD); + if (pm.hasDataMask(MeshModel::MM_WEDGTEXCOORD)){ + pm.clearDataMask(MeshModel::MM_WEDGTEXCOORD); } - if (m.hasDataMask(MeshModel::MM_VERTTEXCOORD)){ - m.clearDataMask(MeshModel::MM_VERTTEXCOORD); + if (pm.hasDataMask(MeshModel::MM_VERTTEXCOORD)){ + pm.clearDataMask(MeshModel::MM_VERTTEXCOORD); } + bool convex_hullFP = par.getBool("convex_hullFP"); + int result = visible_points(dim,numpoints,m,pm,viewpoint,threshold,convex_hullFP ); - bool convex_hull = par.getBool("convex_hull"); - int result = visible_points(dim,numpoints,m,pm,viewpoint,threshold,convex_hull); - - if(!convex_hull) + if(!convex_hullFP) md.delMesh(&pm); else{ + //Re-orient normals + bool reorient= par.getBool("reorient"); + if (reorient){ + pm.updateDataMask(MeshModel::MM_FACEFACETOPO); + bool oriented,orientable; + + tri::Clean::IsOrientedMesh(pm.cm, oriented,orientable); + vcg::tri::UpdateTopology::FaceFace(pm.cm); + vcg::tri::UpdateTopology::TestFaceFace(pm.cm); + pm.clearDataMask(MeshModel::MM_FACEFACETOPO); + tri::UpdateSelection::ClearFace(pm.cm); + } vcg::tri::UpdateBounding::Box(pm.cm); vcg::tri::UpdateNormals::PerVertexNormalizedPerFace(pm.cm); Log(GLLogStream::FILTER,"Successfully created a mesh of %i vert and %i faces",pm.cm.vn,pm.cm.fn); } - + + bool convex_hullNFP = par.getBool("convex_hullNFP"); + if(convex_hullNFP){ + bool result =create_mesh_convex_hull(md,m,par); + if(result) + Log(GLLogStream::FILTER,"Successfully created a mesh of %i vert and %i faces",(*md.mm()).cm.vn,(*md.mm()).cm.fn); + } + if(result>=0){ Log(GLLogStream::FILTER,"Selected %i visible points", result); return true; @@ -483,4 +456,92 @@ bool QhullPlugin::applyFilter(QAction *filter, MeshDocument &md, FilterParameter } } +bool create_mesh_convex_hull(MeshDocument &md,MeshModel &m,FilterParameterSet & par){ + MeshModel &pm =*md.addNewMesh("Convex Hull"); + + if (m.hasDataMask(MeshModel::MM_WEDGTEXCOORD)){ + m.clearDataMask(MeshModel::MM_WEDGTEXCOORD); + } + if (m.hasDataMask(MeshModel::MM_VERTTEXCOORD)){ + m.clearDataMask(MeshModel::MM_VERTTEXCOORD); + } + + int dim= 3; /* dimension of points */ + int numpoints= m.cm.vn; /* number of mesh vertices */ + + //facet_list contains the convex hull + //facet_list contains simplicial (triangulated) facets. + //The convex hull of a set of points is the smallest convex set containing the points. + + facetT *facet_list = compute_convex_hull(dim,numpoints,m); + + if(facet_list!=NULL){ + + int convexNumFaces = qh num_facets; + int convexNumVert = qh_setsize(qh_facetvertices (facet_list, NULL, false)); + assert( qh num_vertices == convexNumVert); + + tri::Allocator::AddVertices(pm.cm,convexNumVert); + tri::Allocator::AddFaces(pm.cm,convexNumFaces); + + + /*ivp length is 'numpoints' because each vertex is accessed through its ID whose range is + 0<=qh_pointid(vertex->point)::VertexPointer> ivp(numpoints); + vertexT *vertex; + int i=0; + FORALLvertices{ + if ((*vertex).point){ + pm.cm.vert[i].P()[0] = (*vertex).point[0]; + pm.cm.vert[i].P()[1] = (*vertex).point[1]; + pm.cm.vert[i].P()[2] = (*vertex).point[2]; + ivp[qh_pointid(vertex->point)] = &pm.cm.vert[i]; + i++; + } + } + + facetT *facet; + i=0; + FORALLfacet_(facet_list){ + vertexT *vertex; + int vertex_n, vertex_i; + FOREACHvertex_i_((*facet).vertices){ + pm.cm.face[i].V(vertex_i)= ivp[qh_pointid(vertex->point)]; + } + i++; + } + + assert( pm.cm.fn == convexNumFaces); + //m.cm.Clear(); + + //Re-orient normals + bool reorient= par.getBool("reorient"); + if (reorient){ + pm.updateDataMask(MeshModel::MM_FACEFACETOPO); + bool oriented, orientable; + + tri::Clean::IsOrientedMesh(pm.cm, oriented,orientable); + vcg::tri::UpdateTopology::FaceFace(pm.cm); + vcg::tri::UpdateTopology::TestFaceFace(pm.cm); + pm.clearDataMask(MeshModel::MM_FACEFACETOPO); + tri::UpdateSelection::ClearFace(pm.cm); + } + + vcg::tri::UpdateBounding::Box(pm.cm); + vcg::tri::UpdateNormals::PerVertexNormalizedPerFace(pm.cm); + + int curlong, totlong; /* memory remaining after qh_memfreeshort */ + qh_freeqhull(!qh_ALL); + qh_memfreeshort (&curlong, &totlong); + if (curlong || totlong) + fprintf (stderr, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", + totlong, curlong); + + return true; + + } + else + return false; +} + Q_EXPORT_PLUGIN(QhullPlugin) diff --git a/src/fgt/filter_qhull/qhull_tools.cpp b/src/fgt/filter_qhull/qhull_tools.cpp index f2a3d2c42..7c5f8369a 100644 --- a/src/fgt/filter_qhull/qhull_tools.cpp +++ b/src/fgt/filter_qhull/qhull_tools.cpp @@ -69,7 +69,7 @@ double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim); m --> original mesh compute_convex_hull(int dim, int numpoints, MeshModel &m) - build convex hull from a set of vertices of a mesh + build convex hull from a set of vertices of a mesh with Qhull library (http://www.qhull.org/html/qconvex.htm). returns the convex hull as a list of simplicial (triangulated) facets, if there are no errors; @@ -110,10 +110,10 @@ facetT *compute_convex_hull(int dim, int numpoints, MeshModel &m) m --> original mesh compute_delaunay(int dim, int numpoints, MeshModel &m) - build Delauanay triangulation from a set of vertices of a mesh. + build Delauanay triangulation from a set of vertices of a mesh with Qhull library (http://www.qhull.org/html/qdelaun.htm). - The Delaunay triangulation of a set of points in d-dimensional spaces is the projection of the convex - hull of the projections of the points onto a (d+1)-dimensional paraboloid. + The Delaunay triangulation of a set of points in d-dimensional spaces is a triangulation of the convex hull. + The result is non 2-manifold by definition. By default, qdelaunay merges regions with cocircular or cospherical input sites. If you want a simplicial triangulation use triangulated output ('Qt') or joggled input ('QJ'). @@ -159,10 +159,16 @@ facetT *compute_delaunay(int dim, int numpoints, MeshModel &m) numpoints --> number of points m --> original mesh pm --> new mesh + threshold --> ?????????????????????????????????????? compute_voronoi(int dim, int numpoints, MeshModel &m, MeshModel &pm) - Reconstructs the mesh surface starting from its vertices through a double - Delaunay triangulation. The second one takes in input some Voronoi vertices too. + Compute a Voronoi filtering(Amenta and Bern 1998) with Qhull library (http://www.qhull.org/). + The filter computes a piecewise-linear approximation of a smooth surface from a finite set of sample points + The algorithm uses Voronoi vertices to remove triangles from the Delaunay triangulation. + + After computing the Voronoi diagram, foreach sample point it chooses the two farthest opposite Voronoi vertices. + Then computes a Delaunay triangulation of the sample points and the selected Voronoi vertices, and keep + only those triangles in witch all three vertices are sample points. Options 'Qt' (triangulated output) and 'QJ' (joggled input) may produce unexpected results. @@ -177,7 +183,7 @@ facetT *compute_delaunay(int dim, int numpoints, MeshModel &m) false otherwise. */ -bool compute_voronoi(int dim, int numpoints, MeshModel &m, MeshModel &pm) +bool compute_voronoi(int dim, int numpoints, MeshModel &m, MeshModel &pm, float threshold) { coordT *points; /* array of coordinates for each point*/ boolT ismalloc= True; /* True if qhull should free points in qh_freeqhull() or reallocation */ @@ -292,7 +298,7 @@ bool compute_voronoi(int dim, int numpoints, MeshModel &m, MeshModel &pm) for(int i=0;i(500*pm.cm.bbox.Diag())) + if(qh_pointdist(bbCenter,pole,dim+1)>(threshold*pm.cm.bbox.Diag())) discard=true; } if(!discard) @@ -342,7 +348,7 @@ bool compute_voronoi(int dim, int numpoints, MeshModel &m, MeshModel &pm) for(int i=0;i(500*pm.cm.bbox.Diag())) + if(qh_pointdist(bbCenter,pole,dim+1)>(threshold*pm.cm.bbox.Diag())) discard=true; } if(!discard) @@ -463,7 +469,7 @@ bool compute_voronoi(int dim, int numpoints, MeshModel &m, MeshModel &pm) alphashape --> true to calculate alpha shape, false alpha complex compute_alpha_shapes(int dim, int numpoints, MeshModel &m, MeshModel &pm, double alpha, bool alphashape) - build alpha complex or alpha shapes from a set of vertices of a mesh. + build alpha complex or alpha shapes (Edelsbrunner and P.Mucke 1994)from a set of vertices of a mesh with Qhull library (http://www.qhull.org/ Insert the minimum value of alpha (radius of the circumsphere) in attribute quality foreach face. The Alpha Shape is the boundary of the alpha complex, that is a subcomplex of the Delaunay triangulation. @@ -657,21 +663,27 @@ bool compute_alpha_shapes(int dim, int numpoints, MeshModel &m, MeshModel &pm, d dim --> dimension of points numpoints --> number of points m --> original mesh - pm --> new mesh that is the convex hull of the flipped points, if convex_hull is true + pm --> new mesh that is the convex hull of the flipped points, if convex_hullFP is true + pm2 --> new mesh that is the convex hull of the original points, if convex_hullNFP is true viewpoint -> the viewpoint threshold -> bounds the radius oh the sphere used to select visible points - convex_hull --> true if you want to show the partial Convex Hull of the transformed points cloud + convex_hullFP --> true if you want to show the partial Convex Hull of the transformed points cloud + false otherwise + convex_hullNFP --> true if you want to show the Convex Hull of the original points cloud false otherwise bool visible_points(int dim, int numpoints, MeshModel &m, MeshModel &pm) - Select the visibile points from a given viewpoint. + Select the visible points in a point cloud, as viewed from a given viewpoint. + It uses the Qhull library (http://www.qhull.org/. + The algorithm used (Katz, Tal and Basri 2007) determines visibility without reconstructing a surface or estimating normals. + The visible points are those that corresponds to the ones that reside on the convex hull of a transformed point cloud. returns the number of visible points if no errors occurred; -1 otherwise. */ -int visible_points(int dim, int numpoints, MeshModel &m,MeshModel &pm,Point3f viewpointP,float threshold, bool convex_hull){ +int visible_points(int dim, int numpoints, MeshModel &m, MeshModel &pm,vcg::Point3f viewpointP,float threshold,bool convex_hullFP){ boolT ismalloc= True; /* True if qhull should free points in qh_freeqhull() or reallocation */ char flags[]= "qhull Tcv"; /* option flags for qhull, see qh_opt.htm */ FILE *outfile= NULL; /* output from qh_produce_output() @@ -764,7 +776,7 @@ int visible_points(int dim, int numpoints, MeshModel &m,MeshModel &pm,Point3f vi //mark the vertex in the mesh that is correspondent to the flipped one (share the same index) ivp[qh_pointid(vertex->point)]->SetS(); - if(convex_hull){ //Add flipped points to the new mesh + if(convex_hullFP){ //Add flipped points to the new mesh pm.cm.vert[i].P()[0] = (*vertex).point[0]; pm.cm.vert[i].P()[1] = (*vertex).point[1]; pm.cm.vert[i].P()[2] = (*vertex).point[2]; @@ -775,7 +787,7 @@ int visible_points(int dim, int numpoints, MeshModel &m,MeshModel &pm,Point3f vi } vcg::tri::UpdateColor::VertexSelected(m.cm); - if(convex_hull){ + if(convex_hullFP){ //Add to the new mesh only the faces of the convex hull whose vertices aren't the viewpoint facetT *facet; FORALLfacet_(qh facet_list){ @@ -852,90 +864,4 @@ double calculate_circumradius(pointT* p0,pointT* p1,pointT* p2, int dim){ coordT sum =(a + b + c)*0.5; coordT area = sum*(a+b-sum)*(a+c-sum)*(b+c-sum); return (double) (a*b*c)/(4*sqrt(area)); -} - -//-------------------------------------------------------------------------------- - -bool test(facetT* facet, ridgeT* triangle ,double alpha){ - vertexT* vertex = (vertexT *)(facet->vertices->e[0].p); - double* center = facet->center; - double px = (*vertex).point[0]; - double qx = *center; - double py = (*vertex).point[1]; - double qy = *center++; - double pz = (*vertex).point[2]; - double qz = *center++; - double radius = sqrt(pow(px-qx,2)+pow(py-qy,2)+pow(pz-qz,2)); - int ridgesCount=0; - if (radius>alpha){ - //facet->visitid= qh visit_id; - //qh_makeridges(facet); - //ridgeT *ridge, **ridgep; - //FOREACHridge_(facet->ridges) { - // if (ridge!=triangle) { - // //Calcola il raggio della circosfera - // pointT* p0 = ((vertexT*) (ridge->vertices->e[0].p))->point; - // pointT* p1 = ((vertexT*) (ridge->vertices->e[1].p))->point; - // pointT* p2 = ((vertexT*) (ridge->vertices->e[2].p))->point; - - // coordT a = qh_pointdist(p0,p1,3); - // coordT b = qh_pointdist(p1,p2,3); - // coordT c = qh_pointdist(p2,p0,3); - - // coordT sum =(a + b + c)*0.5; - // coordT area = sum*(a+b-sum)*(a+c-sum)*(b+c-sum); - // radius = (double) (a*b*c)/(4*sqrt(area)); - // - // if(radius <=alpha){ - // ridgesCount++; - // } - // } - //if(ridgesCount==3) - // return true; - //else - return false; - } - else return true; -} - -void addDelaunay (coordT *points, int numpoints, int numnew, int dim) { - int j; - coordT *point; - facetT *facet; - realT bestdist; - boolT isoutside; - - for (j= 0; j < numnew ; j++) { - point= points + (numpoints+j)*dim; - if (points == qh first_point) /* in case of 'QRn' */ - qh num_points= numpoints+j+1; - /* qh num_points sets the size of the points array. You may - allocate the point elsewhere. If so, qh_addpoint records - the point's address in qh other_points - */ - - qh_setdelaunay (dim, 1, point); - facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside); - if (isoutside) { - if (!qh_addpoint (point, facet, False)) - break; /* user requested an early exit with 'TVn' or 'TCn' */ - } - - /* qh_produce_output(); */ - } - if (qh DOcheckmax) - qh_check_maxout(); - else if (qh KEEPnearinside) - qh_nearcoplanar(); -} /*.addDelaunay.*/ - -coordT pointdist(pointT *point1, pointT *point2, int dim) { - coordT dist, diff; - dist= 0.0; - - for (int k=0;k