src/math/ConvexHull.cpp

/* Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 *
 *  ConvexHull.cpp
 *
 *  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.6 2008-05-14 14:31:48 chuckv Exp $
 *
 */

/* Standard includes independent of library */
#include <iostream>
#include <fstream>
#include <list>
#include <algorithm>
#include <iterator>
#include "math/ConvexHull.hpp"
#include "utils/simError.h"
using namespace oopse;

/* CGAL version of convex hull first then QHULL */
#ifdef HAVE_CGAL

#include <CGAL/basic.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Convex_hull_traits_3.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/double.h>

#include <CGAL/Quotient.h>
#include <CGAL/MP_Float.h>
#include <CGAL/Lazy_exact_nt.h>

typedef double RT;
typedef CGAL::Simple_cartesian<RT>                K;
typedef CGAL::Convex_hull_traits_3<K>             Traits;
typedef Traits::Polyhedron_3                      Polyhedron_3;
typedef K::Point_3                                Point;


ConvexHull::ConvexHull(){}

bool ConvexHull::genHull(std::vector<Vector3d> pos)
{
  
  std::vector<Point> 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()<<std::endl;
  }
  
  
}


#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<Vector3d> pos)
{
  
  
  int numpoints = pos.size();
  coordT* pt_array;
  coordT* surfpt_array;
  std::list<int> 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<char *>(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<char *>(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;
  //    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<int>::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<char *>(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)
{

  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


Revision:	3392
Committed:	Wed May 14 14:31:48 2008 UTC (16 years, 11 months ago) by chuckv
File size:	8524 byte(s)
Log Message:	Changes....
#	User	Rev	Content
1	chuckv	3083	/* Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
2			*
3			* The University of Notre Dame grants you ("Licensee") a
4			* non-exclusive, royalty free, license to use, modify and
5			* redistribute this software in source and binary code form, provided
6			* that the following conditions are met:
7			*
8			* 1. Acknowledgement of the program authors must be made in any
9			* publication of scientific results based in part on use of the
10			* program. An acceptable form of acknowledgement is citation of
11			* the article in which the program was described (Matthew
12			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
13			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
14			* Parallel Simulation Engine for Molecular Dynamics,"
15			* J. Comput. Chem. 26, pp. 252-271 (2005))
16			*
17			* 2. Redistributions of source code must retain the above copyright
18			* notice, this list of conditions and the following disclaimer.
19			*
20			* 3. Redistributions in binary form must reproduce the above copyright
21			* notice, this list of conditions and the following disclaimer in the
22			* documentation and/or other materials provided with the
23			* distribution.
24			*
25			* This software is provided "AS IS," without a warranty of any
26			* kind. All express or implied conditions, representations and
27			* warranties, including any implied warranty of merchantability,
28			* fitness for a particular purpose or non-infringement, are hereby
29			* excluded. The University of Notre Dame and its licensors shall not
30			* be liable for any damages suffered by licensee as a result of
31			* using, modifying or distributing the software or its
32			* derivatives. In no event will the University of Notre Dame or its
33			* licensors be liable for any lost revenue, profit or data, or for
34			* direct, indirect, special, consequential, incidental or punitive
35			* damages, however caused and regardless of the theory of liability,
36			* arising out of the use of or inability to use software, even if the
37			* University of Notre Dame has been advised of the possibility of
38			* such damages.
39			*
40			*
41			* ConvexHull.cpp
42			*
43	chuckv	3137	* Purpose: To calculate convexhull, hull volume libqhull.
44	chuckv	3083	*
45			* Created by Charles F. Vardeman II on 11 Dec 2006.
46			* @author Charles F. Vardeman II
47	chuckv	3392	* @version $Id: ConvexHull.cpp,v 1.6 2008-05-14 14:31:48 chuckv Exp $
48	chuckv	3083	*
49			*/
50
51	chuckv	3392	/* Standard includes independent of library */
52	chuckv	3083	#include <iostream>
53			#include <fstream>
54	chuckv	3275	#include <list>
55			#include <algorithm>
56			#include <iterator>
57	chuckv	3141	#include "math/ConvexHull.hpp"
58	chuckv	3275	#include "utils/simError.h"
59	chuckv	3083	using namespace oopse;
60
61	chuckv	3392	/* CGAL version of convex hull first then QHULL */
62			#ifdef HAVE_CGAL
63	chuckv	3275
64	chuckv	3392	#include <CGAL/basic.h>
65			#include <CGAL/Simple_cartesian.h>
66			#include <CGAL/Convex_hull_traits_3.h>
67			#include <CGAL/convex_hull_3.h>
68			#include <CGAL/double.h>
69	chuckv	3275
70	chuckv	3392	#include <CGAL/Quotient.h>
71			#include <CGAL/MP_Float.h>
72			#include <CGAL/Lazy_exact_nt.h>
73	chuckv	3275
74	chuckv	3392	typedef double RT;
75			typedef CGAL::Simple_cartesian<RT> K;
76			typedef CGAL::Convex_hull_traits_3<K> Traits;
77			typedef Traits::Polyhedron_3 Polyhedron_3;
78			typedef K::Point_3 Point;
79	chuckv	3275
80
81	chuckv	3392	ConvexHull::ConvexHull(){}
82	chuckv	3275
83	chuckv	3392	bool ConvexHull::genHull(std::vector<Vector3d> pos)
84			{
85
86			std::vector<Point> points;
87
88
89			// Copy the positon vector into a points vector for cgal.
90			for (int i = 0; i < pos.size(); ++i)
91			{
92			Point pt(pos[i][0],pos[i][1],pos[i][2]);
93			points.push_back(pt);
94			}
95
96			// define object to hold convex hull
97			Polyhedron_3 ch_object_;
98			// compute convex hull
99			CGAL::convex_hull_3(points.begin(), points.end(), ch_object_);
100
101			for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin(); ch_object_.vertices_end(); ++v){
102			std::cout<< v.point()<<std::endl;
103			}
104
105
106			}
107	chuckv	3275
108
109
110	chuckv	3392	#else
111			#ifdef HAVE_QHULL
112			/* Old options Qt Qu Qg QG0 FA */
113			ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt Qci Tcv Pp")
114			//ConvexHull::ConvexHull() : dim_(3), options_("qhull d Qbb Qt i")
115			{}
116	chuckv	3275
117	chuckv	3392	bool ConvexHull::genHull(std::vector<Vector3d> pos)
118			{
119
120
121			int numpoints = pos.size();
122			coordT* pt_array;
123			coordT* surfpt_array;
124			std::list<int> surface_atoms;
125			FILE *outdummy = NULL;
126			FILE *errdummy = NULL;
127
128			pt_array = (coordT) malloc(sizeof(coordT) (numpoints * dim_));
129
130
131			for (int i = 0; i < numpoints; i++) {
132			pt_array[dim_ * i] = pos[i][0];
133			pt_array[dim_ * i + 1] = pos[i][1];
134			pt_array[dim_ * i + 2] = pos[i][2];
135			}
136
137
138
139
140			/*
141			qh_initflags(const_cast<char *>(options_.c_str()));
142			qh_init_B(pospoints, numpoints, dim_, ismalloc);
143			qh_qhull();
144			qh_check_output();
145	chuckv	3275
146	chuckv	3392	qh_produce_output();
147			*/
148			boolT ismalloc = False;
149
150			if (!qh_new_qhull(dim_, numpoints, pt_array, ismalloc,
151			const_cast<char *>(options_.c_str()), NULL, stderr)) {
152
153			vertexT vertex, *vertexp;
154			facetT *facet;
155			setT *vertices;
156			unsigned int nf = qh num_facets;
157
158			//Matrix idx(nf, dim);
159			/*
160			int j, i = 0, id = 0;
161
162			int id2 = 0;
163			coordT point,*pointp;
164			realT dist;
165			FORALLfacets {
166			j = 0;
167
168			if (!facet->simplicial){
169			// should never happen with Qt
170			sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
171			painCave.isFatal = 0;
172			simError();
173			}
174	chuckv	3275
175	chuckv	3392	vertices = qh_facet3vertex(facet);
176			FOREACHvertex_(vertices){
177			id = qh_pointid(vertex->point);
178			surface_atoms.push_back(id);
179			//std::cout << "Ag " << pos[id][0] << " " << pos[id][1] << " " << pos[id][2]<< std::endl;
180			}
181			qh_settempfree(&vertices);
182
183			FOREACHpoint_(facet->coplanarset){
184			vertex= qh_nearvertex (facet, point, &dist);
185			//id= qh_pointid (vertex->point);
186			id2= qh_pointid (point);
187			surface_atoms.push_back(id2);
188			//std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl;
189			//printf ("%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
190			//std::cout << "Neighbors are: %d $d %d\n" << id << id2 << facet->id;
191	chuckv	3275
192	chuckv	3392	}
193
194			}
195	chuckv	3275
196			*/
197
198
199	chuckv	3392	}
200	chuckv	3275
201
202
203
204	chuckv	3392	qh_getarea(qh facet_list);
205			volume_ = qh totvol;
206			area_ = qh totarea;
207			// FILE *newGeomFile;
208
209
210			/*
211			FORALLfacets {
212			for (int k=0; k < qh hull_dim; k++)
213			printf ("%6.2g ", facet->normal[k]);
214			printf ("\n");
215			}
216			*/
217
218			int curlong,totlong;
219			// geomviewHull("junk.off");
220			qh_freeqhull(!qh_ALL);
221			qh_memfreeshort(&curlong, &totlong);
222			if (curlong \|\| totlong)
223			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
224			<< totlong << curlong << std::endl;
225			free(pt_array);
226			/*
227			int j = 0;
228			surface_atoms.sort();
229			surface_atoms.unique();
230			surfpt_array = (coordT) malloc(sizeof(coordT) (surface_atoms.size() * dim_));
231			for(std::list<int>::iterator list_iter = surface_atoms.begin();
232			list_iter != surface_atoms.end(); list_iter++)
233			{
234			int i = *list_iter;
235			//surfpt_array[dim_ * j] = pos[i][0];
236			//surfpt_array[dim_ * j + 1] = pos[i][1];
237			//surfpt_array[dim_ * j + 2] = pos[i][2];
238			std::cout << "Ag " << pos[i][0] << " " << pos[i][1] << " "<< pos[i][2] << std::endl;
239			j++;
240			}
241			*/
242
243			/*
244	chuckv	3275	std::string deloptions_ = "qhull d Qt";
245			facetT facet, neighbor;
246			ridgeT ridge, *ridgep;
247	chuckv	3392
248	chuckv	3275	if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc,
249			const_cast<char *>(deloptions_.c_str()), NULL, stderr)) {
250	chuckv	3392
251			qh visit_id++;
252			FORALLfacets {
253			if (!facet->upperdelaunay) {
254			facet->visitid= qh visit_id;
255			qh_makeridges(facet);
256			FOREACHridge_(facet->ridges) {
257			neighbor= otherfacet_(ridge, facet);
258			if (neighbor->visitid != qh visit_id) {
259
260			FOREACHvertex_(ridge->vertices)
261			int id2 = qh_pointid (vertex->point);
262			std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl;
263			}
264			}
265			}
266
267	chuckv	3275
268
269	chuckv	3392
270			}
271	chuckv	3275
272	chuckv	3392	qh_freeqhull(!qh_ALL);
273			qh_memfreeshort(&curlong, &totlong);
274			if (curlong \|\| totlong)
275			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
276			<< totlong << curlong << std::endl;
277			free(surfpt_array);
278			*/
279			return true;
280	chuckv	3141	}
281	chuckv	3083
282	chuckv	3143
283	chuckv	3275
284			RealType ConvexHull::getVolume()
285			{
286	chuckv	3392	return volume_;
287	chuckv	3083	}
288
289	chuckv	3275	void ConvexHull::geomviewHull(const std::string& geomFileName)
290			{
291	chuckv	3137
292	chuckv	3392	FILE *newGeomFile;
293
294			//create new .md file based on old .md file
295			newGeomFile = fopen(geomFileName.c_str(), "w");
296			qh_findgood_all(qh facet_list);
297			for (int i = 0; i < qh_PRINTEND; i++)
298			qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
299
300			fclose(newGeomFile);
301	chuckv	3083	}
302	chuckv	3392	#endif //QHULL
303			#endif //CGAL
304	chuckv	3083
305	chuckv	3275
306