--- trunk/src/math/ConvexHull.cpp 2007/05/30 18:47:04 1141 +++ trunk/src/math/ConvexHull.cpp 2008/06/18 17:03:30 1261 @@ -44,74 +44,274 @@ * * Created by Charles F. Vardeman II on 11 Dec 2006. * @author Charles F. Vardeman II - * @version $Id: ConvexHull.cpp,v 1.3 2007-05-30 18:47:03 chuckv Exp $ + * @version $Id: ConvexHull.cpp,v 1.7 2008-06-18 17:03:30 chuckv Exp $ * */ +/* Standard includes independent of library */ #include #include +#include +#include +#include #include "math/ConvexHull.hpp" - +#include "utils/simError.h" using namespace oopse; -ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt FA") { +/* CGAL version of convex hull first then QHULL */ +#ifdef HAVE_CGAL + +#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_3; + + +ConvexHull::ConvexHull(){} + +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_3 pt(pos[i][0],pos[i][1],pos[i][2]); + points.push_back(pt); + } + + // define object to hold convex hull + CGAL::Object ch_object_; + Polyhedron_3 polyhedron; + + // compute convex hull + std::cerr << "Creating hull" << std::endl; + CGAL::convex_hull_3(points.begin(), points.end(), ch_object_); + std::cerr << "Done Creating hull" << std::endl; + std::vector::const_iterator p_it; + + for (p_it = points.begin(); p_it != points.end(); p_it++) + { + std::cerr << (*p_it).x() << std::endl; + } + + /* + for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin(); + ch_object_.vertices_end(); ++v){ + std::cout<< v.point()< pos) { - FILE *outfile = stdout; - FILE *errfile = stderr; - facetT *facet; - int exitcode; - boolT ismalloc = False; - int curlong,totlong; + + +#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") +{} + +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; - coordT points[numpoints][dim_]; + pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_)); - for (int i=0; i(options_.c_str())); - qh_init_B (points[0], numpoints, dim_, ismalloc); - qh_qhull(); - qh_check_output(); + /* + 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; + + 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; + + } + + } + +*/ + + + } + + + + qh_getarea(qh facet_list); - volume_ = qh totvol; - area_ = qh totarea; + 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); + qh_memfreeshort(&curlong, &totlong); if (curlong || totlong) - fprintf (errfile, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", - totlong, curlong); + 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_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() { + +RealType ConvexHull::getVolume() +{ return volume_; } -void ConvexHull::geomviewHull(const std::string& geomFileName) { +void ConvexHull::geomviewHull(const std::string& geomFileName) +{ + + FILE *newGeomFile; - std::ofstream 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); - - newGeomFile.close(); + fclose(newGeomFile); } +#endif //QHULL +#endif //CGAL + + +