--- trunk/src/math/ConvexHull.cpp 2008/04/25 15:14:47 1241 +++ trunk/src/math/ConvexHull.cpp 2008/05/14 14:31:48 1242 @@ -44,10 +44,11 @@ * * Created by Charles F. Vardeman II on 11 Dec 2006. * @author Charles F. Vardeman II - * @version $Id: ConvexHull.cpp,v 1.5 2007-11-22 16:39:45 chuckv Exp $ + * @version $Id: ConvexHull.cpp,v 1.6 2008-05-14 14:31:48 chuckv Exp $ * */ +/* Standard includes independent of library */ #include #include #include @@ -57,197 +58,249 @@ using namespace oopse; #include "utils/simError.h" using namespace oopse; -/* 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") -{} +/* CGAL version of convex hull first then QHULL */ +#ifdef HAVE_CGAL -bool ConvexHull::genHull(std::vector pos) -{ +#include +#include +#include +#include +#include +#include +#include +#include - int numpoints = pos.size(); - coordT* pt_array; - coordT* surfpt_array; - std::list surface_atoms; - FILE *outdummy = NULL; - FILE *errdummy = NULL; +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; - pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_)); +ConvexHull::ConvexHull(){} - 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]; - } +bool ConvexHull::genHull(std::vector pos) +{ + + std::vector points; + + + // Copy the positon vector into a points vector for cgal. + for (int i = 0; i < pos.size(); ++i) + { + 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_); + + for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin(); ch_object_.vertices_end(); ++v){ + std::cout<< v.point()<(options_.c_str())); - qh_init_B(pospoints, numpoints, dim_, ismalloc); - qh_qhull(); - qh_check_output(); +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(); - 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; + 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(); + } - 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; + 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; - } - - } + } + + } */ - } + } - 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++) - { - 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++; - } -*/ - -/* + 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++) + { + 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::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 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; } - 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_; + return volume_; } void ConvexHull::geomviewHull(const std::string& geomFileName) { - FILE *newGeomFile; - - //create new .md file based on old .md file - 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); + FILE *newGeomFile; + + //create new .md file based on old .md file + 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 -