--- trunk/src/math/ConvexHull.cpp 2006/12/14 19:32:32 1097 +++ trunk/src/math/ConvexHull.cpp 2008/05/14 14:31:48 1242 @@ -40,112 +40,267 @@ * * ConvexHull.cpp * - * Purpose: To calculate convexhull, hull volume and radius - * using the CGAL library. + * Purpose: To calculate convexhull, hull volume libqhull. * * Created by Charles F. Vardeman II on 11 Dec 2006. * @author Charles F. Vardeman II - * @version $Id: ConvexHull.cpp,v 1.1 2006-12-14 19:32:32 chuckv Exp $ + * @version $Id: ConvexHull.cpp,v 1.6 2008-05-14 14:31:48 chuckv Exp $ * */ -#include "math/ConvexHull.hpp" +/* Standard includes independent of library */ #include #include +#include +#include +#include +#include "math/ConvexHull.hpp" +#include "utils/simError.h" +using namespace oopse; +/* CGAL version of convex hull first then QHULL */ +#ifdef HAVE_CGAL -using namespace oopse; +#include +#include +#include +#include +#include +#include +#include +#include +typedef double RT; +typedef CGAL::Simple_cartesian K; +typedef CGAL::Convex_hull_traits_3 Traits; +typedef Traits::Polyhedron_3 Polyhedron_3; +typedef K::Point_3 Point; +ConvexHull::ConvexHull(){} + bool ConvexHull::genHull(std::vector pos) { - - std::vector points; - points.reserve(pos.size()); + + std::vector points; + + // Copy the positon vector into a points vector for cgal. for (int i = 0; i < pos.size(); ++i) { - Point_3 pt(pos[i][0],pos[i][1],pos[i][2]); + Point pt(pos[i][0],pos[i][1],pos[i][2]); points.push_back(pt); } - + + // define object to hold convex hull + Polyhedron_3 ch_object_; // compute convex hull - CGAL::convex_hull_3(points.begin(), points.end(), ch_object); - // Make sure hull is a polyhedron - if ( CGAL::assign(ch_polyhedron, ch_object) ) - { - return true; - } - else - { - return false; - } + CGAL::convex_hull_3(points.begin(), points.end(), ch_object_); + + for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin(); ch_object_.vertices_end(); ++v){ + std::cout<< v.point()< L; - L.push_front(Point(0,0,0)); - L.push_front(Point(1,0,0)); - L.push_front(Point(0,1,0)); - Triangulation T(L.begin(), L.end()); - int n = T.number_of_vertices(); +#else +#ifdef HAVE_QHULL +/* Old options Qt Qu Qg QG0 FA */ +ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt Qci Tcv Pp") + //ConvexHull::ConvexHull() : dim_(3), options_("qhull d Qbb Qt i") +{} - // insertion from a vector : - std::vector V(3); - V[0] = Point(0,0,1); - V[1] = Point(1,1,1); - V[2] = Point(2,2,2); +bool ConvexHull::genHull(std::vector pos) +{ + + + int numpoints = pos.size(); + coordT* pt_array; + coordT* surfpt_array; + std::list surface_atoms; + FILE *outdummy = NULL; + FILE *errdummy = NULL; + + pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_)); + + + for (int i = 0; i < numpoints; i++) { + pt_array[dim_ * i] = pos[i][0]; + pt_array[dim_ * i + 1] = pos[i][1]; + pt_array[dim_ * i + 2] = pos[i][2]; + } + + + + + /* + qh_initflags(const_cast(options_.c_str())); + qh_init_B(pospoints, numpoints, dim_, ismalloc); + qh_qhull(); + qh_check_output(); - n = n + T.insert(V.begin(), V.end()); + qh_produce_output(); + */ + boolT ismalloc = False; + + if (!qh_new_qhull(dim_, numpoints, pt_array, ismalloc, + const_cast(options_.c_str()), NULL, stderr)) { + + vertexT *vertex, **vertexp; + facetT *facet; + setT *vertices; + unsigned int nf = qh num_facets; + + //Matrix idx(nf, dim); + /* + int j, i = 0, id = 0; + + int id2 = 0; + coordT *point,**pointp; + realT dist; + FORALLfacets { + j = 0; + + if (!facet->simplicial){ + // should never happen with Qt + sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); + painCave.isFatal = 0; + simError(); + } + + vertices = qh_facet3vertex(facet); + FOREACHvertex_(vertices){ + id = qh_pointid(vertex->point); + surface_atoms.push_back(id); + //std::cout << "Ag " << pos[id][0] << " " << pos[id][1] << " " << pos[id][2]<< std::endl; + } + qh_settempfree(&vertices); + + FOREACHpoint_(facet->coplanarset){ + vertex= qh_nearvertex (facet, point, &dist); + //id= qh_pointid (vertex->point); + id2= qh_pointid (point); + surface_atoms.push_back(id2); + //std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl; + //printf ("%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist); + //std::cout << "Neighbors are: %d $d %d\n" << id << id2 << facet->id; + + } + + } + +*/ + + + } - assert( n == 6 ); // 6 points have been inserted - assert( T.is_valid() ); // checking validity of T - double sum_v = 0; - for(Triangulation::Cell_iterator c_it = T.cells_begin(); c_it != T.cells_end(); ++c_it) + + + qh_getarea(qh facet_list); + volume_ = qh totvol; + area_ = qh totarea; + // FILE *newGeomFile; + + + /* + FORALLfacets { + for (int k=0; k < qh hull_dim; k++) + printf ("%6.2g ", facet->normal[k]); + printf ("\n"); + } + */ + + int curlong,totlong; + // geomviewHull("junk.off"); + qh_freeqhull(!qh_ALL); + qh_memfreeshort(&curlong, &totlong); + if (curlong || totlong) + std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " + << totlong << curlong << std::endl; + free(pt_array); + /* + int j = 0; + surface_atoms.sort(); + surface_atoms.unique(); + surfpt_array = (coordT*) malloc(sizeof(coordT) * (surface_atoms.size() * dim_)); + for(std::list::iterator list_iter = surface_atoms.begin(); + list_iter != surface_atoms.end(); list_iter++) { - sum_v += T.tetrahedron(c_it).volume(); + int i = *list_iter; + //surfpt_array[dim_ * j] = pos[i][0]; + //surfpt_array[dim_ * j + 1] = pos[i][1]; + //surfpt_array[dim_ * j + 2] = pos[i][2]; + std::cout << "Ag " << pos[i][0] << " " << pos[i][1] << " "<< pos[i][2] << std::endl; + j++; } - std::cout << "sum_v " << sum_v << std::endl; -*/ - return 0.0; + */ + + /* + std::string deloptions_ = "qhull d Qt"; + facetT *facet, *neighbor; + ridgeT *ridge, **ridgep; + + if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc, + const_cast(deloptions_.c_str()), NULL, stderr)) { + + qh visit_id++; + FORALLfacets { + if (!facet->upperdelaunay) { + facet->visitid= qh visit_id; + qh_makeridges(facet); + FOREACHridge_(facet->ridges) { + neighbor= otherfacet_(ridge, facet); + if (neighbor->visitid != qh visit_id) { + + FOREACHvertex_(ridge->vertices) + int id2 = qh_pointid (vertex->point); + std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl; + } + } + } + + + + + } + + qh_freeqhull(!qh_ALL); + qh_memfreeshort(&curlong, &totlong); + if (curlong || totlong) + std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " + << totlong << curlong << std::endl; + free(surfpt_array); + */ + return true; } + + +RealType ConvexHull::getVolume() +{ + return volume_; +} + void ConvexHull::geomviewHull(const std::string& geomFileName) { - std::ofstream newGeomFile; - + FILE *newGeomFile; + //create new .md file based on old .md file - newGeomFile.open(geomFileName.c_str()); + newGeomFile = fopen(geomFileName.c_str(), "w"); + qh_findgood_all(qh facet_list); + for (int i = 0; i < qh_PRINTEND; i++) + qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); + + fclose(newGeomFile); +} +#endif //QHULL +#endif //CGAL - // Write polyhedron in Object File Format (OFF). - CGAL::set_ascii_mode( std::cout); - newGeomFile << "OFF" << std::endl << ch_polyhedron.size_of_vertices() << ' ' - << ch_polyhedron.size_of_facets() << " 0" << std::endl; - std::copy( ch_polyhedron.points_begin(), ch_polyhedron.points_end(), - std::ostream_iterator( newGeomFile, "\n")); - for ( Facet_iterator i = ch_polyhedron.facets_begin(); i != ch_polyhedron.facets_end(); ++i) - { - Halfedge_facet_circulator j = i->facet_begin(); - // Facets in polyhedral surfaces are at least triangles. - CGAL_assertion( CGAL::circulator_size(j) >= 3); - newGeomFile << CGAL::circulator_size(j) << ' '; - do - { - newGeomFile << ' ' << std::distance(ch_polyhedron.vertices_begin(), j->vertex()); - } - while ( ++j != i->facet_begin()); - newGeomFile << std::endl; - } - newGeomFile.close(); - -}