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.7 2008-06-18 17:03:30 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_3;


ConvexHull::ConvexHull(){}

bool ConvexHull::genHull(std::vector<Vector3d> pos)
{
  
  std::vector<Point_3> 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<Point_3>::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()<<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:	1261
Committed:	Wed Jun 18 17:03:30 2008 UTC (16 years, 10 months ago) by chuckv
File size:	8826 byte(s)
Log Message:	Added files to module section of Makefile.
#	User	Rev	Content
1	chuckv	1097	/* 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	1137	* Purpose: To calculate convexhull, hull volume libqhull.
44	chuckv	1097	*
45			* Created by Charles F. Vardeman II on 11 Dec 2006.
46			* @author Charles F. Vardeman II
47	chuckv	1261	* @version $Id: ConvexHull.cpp,v 1.7 2008-06-18 17:03:30 chuckv Exp $
48	chuckv	1097	*
49			*/
50
51	chuckv	1242	/* Standard includes independent of library */
52	chuckv	1097	#include <iostream>
53			#include <fstream>
54	chuckv	1188	#include <list>
55			#include <algorithm>
56			#include <iterator>
57	chuckv	1141	#include "math/ConvexHull.hpp"
58	chuckv	1188	#include "utils/simError.h"
59	chuckv	1097	using namespace oopse;
60
61	chuckv	1242	/* CGAL version of convex hull first then QHULL */
62			#ifdef HAVE_CGAL
63	chuckv	1188
64	chuckv	1242	#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	1188
70	chuckv	1242	#include <CGAL/Quotient.h>
71			#include <CGAL/MP_Float.h>
72			#include <CGAL/Lazy_exact_nt.h>
73	chuckv	1188
74	chuckv	1242	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	chuckv	1261	typedef K::Point_3 Point_3;
79	chuckv	1188
80
81	chuckv	1242	ConvexHull::ConvexHull(){}
82	chuckv	1188
83	chuckv	1242	bool ConvexHull::genHull(std::vector<Vector3d> pos)
84			{
85
86	chuckv	1261	std::vector<Point_3> points;
87	chuckv	1242
88
89			// Copy the positon vector into a points vector for cgal.
90			for (int i = 0; i < pos.size(); ++i)
91			{
92	chuckv	1261	Point_3 pt(pos[i][0],pos[i][1],pos[i][2]);
93	chuckv	1242	points.push_back(pt);
94			}
95
96			// define object to hold convex hull
97	chuckv	1261	CGAL::Object ch_object_;
98			Polyhedron_3 polyhedron;
99
100	chuckv	1242	// compute convex hull
101	chuckv	1261	std::cerr << "Creating hull" << std::endl;
102	chuckv	1242	CGAL::convex_hull_3(points.begin(), points.end(), ch_object_);
103	chuckv	1261	std::cerr << "Done Creating hull" << std::endl;
104			std::vector<Point_3>::const_iterator p_it;
105
106			for (p_it = points.begin(); p_it != points.end(); p_it++)
107			{
108			std::cerr << (*p_it).x() << std::endl;
109			}
110
111			/*
112			for (Polyhedron_3::Vertex_iterator v = ch_object_.vertices_begin();
113			ch_object_.vertices_end(); ++v){
114	chuckv	1242	std::cout<< v.point()<<std::endl;
115			}
116	chuckv	1261	*/
117	chuckv	1242	}
118	chuckv	1188
119
120
121	chuckv	1242	#else
122			#ifdef HAVE_QHULL
123			/* Old options Qt Qu Qg QG0 FA */
124			ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt Qci Tcv Pp")
125			//ConvexHull::ConvexHull() : dim_(3), options_("qhull d Qbb Qt i")
126			{}
127	chuckv	1188
128	chuckv	1242	bool ConvexHull::genHull(std::vector<Vector3d> pos)
129			{
130
131
132			int numpoints = pos.size();
133			coordT* pt_array;
134			coordT* surfpt_array;
135			std::list<int> surface_atoms;
136			FILE *outdummy = NULL;
137			FILE *errdummy = NULL;
138
139			pt_array = (coordT) malloc(sizeof(coordT) (numpoints * dim_));
140
141
142			for (int i = 0; i < numpoints; i++) {
143			pt_array[dim_ * i] = pos[i][0];
144			pt_array[dim_ * i + 1] = pos[i][1];
145			pt_array[dim_ * i + 2] = pos[i][2];
146			}
147
148
149
150
151			/*
152			qh_initflags(const_cast<char *>(options_.c_str()));
153			qh_init_B(pospoints, numpoints, dim_, ismalloc);
154			qh_qhull();
155			qh_check_output();
156	chuckv	1188
157	chuckv	1242	qh_produce_output();
158			*/
159			boolT ismalloc = False;
160
161			if (!qh_new_qhull(dim_, numpoints, pt_array, ismalloc,
162			const_cast<char *>(options_.c_str()), NULL, stderr)) {
163
164			vertexT vertex, *vertexp;
165			facetT *facet;
166			setT *vertices;
167			unsigned int nf = qh num_facets;
168
169			//Matrix idx(nf, dim);
170			/*
171			int j, i = 0, id = 0;
172
173			int id2 = 0;
174			coordT point,*pointp;
175			realT dist;
176			FORALLfacets {
177			j = 0;
178
179			if (!facet->simplicial){
180			// should never happen with Qt
181			sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected");
182			painCave.isFatal = 0;
183			simError();
184			}
185	chuckv	1188
186	chuckv	1242	vertices = qh_facet3vertex(facet);
187			FOREACHvertex_(vertices){
188			id = qh_pointid(vertex->point);
189			surface_atoms.push_back(id);
190			//std::cout << "Ag " << pos[id][0] << " " << pos[id][1] << " " << pos[id][2]<< std::endl;
191			}
192			qh_settempfree(&vertices);
193
194			FOREACHpoint_(facet->coplanarset){
195			vertex= qh_nearvertex (facet, point, &dist);
196			//id= qh_pointid (vertex->point);
197			id2= qh_pointid (point);
198			surface_atoms.push_back(id2);
199			//std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl;
200			//printf ("%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
201			//std::cout << "Neighbors are: %d $d %d\n" << id << id2 << facet->id;
202	chuckv	1188
203	chuckv	1242	}
204
205			}
206	chuckv	1188
207			*/
208
209
210	chuckv	1242	}
211	chuckv	1188
212
213
214
215	chuckv	1242	qh_getarea(qh facet_list);
216			volume_ = qh totvol;
217			area_ = qh totarea;
218			// FILE *newGeomFile;
219
220
221			/*
222			FORALLfacets {
223			for (int k=0; k < qh hull_dim; k++)
224			printf ("%6.2g ", facet->normal[k]);
225			printf ("\n");
226			}
227			*/
228
229			int curlong,totlong;
230			// geomviewHull("junk.off");
231			qh_freeqhull(!qh_ALL);
232			qh_memfreeshort(&curlong, &totlong);
233			if (curlong \|\| totlong)
234			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
235			<< totlong << curlong << std::endl;
236			free(pt_array);
237			/*
238			int j = 0;
239			surface_atoms.sort();
240			surface_atoms.unique();
241			surfpt_array = (coordT) malloc(sizeof(coordT) (surface_atoms.size() * dim_));
242			for(std::list<int>::iterator list_iter = surface_atoms.begin();
243			list_iter != surface_atoms.end(); list_iter++)
244			{
245			int i = *list_iter;
246			//surfpt_array[dim_ * j] = pos[i][0];
247			//surfpt_array[dim_ * j + 1] = pos[i][1];
248			//surfpt_array[dim_ * j + 2] = pos[i][2];
249			std::cout << "Ag " << pos[i][0] << " " << pos[i][1] << " "<< pos[i][2] << std::endl;
250			j++;
251			}
252			*/
253
254			/*
255	chuckv	1188	std::string deloptions_ = "qhull d Qt";
256			facetT facet, neighbor;
257			ridgeT ridge, *ridgep;
258	chuckv	1242
259	chuckv	1188	if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc,
260			const_cast<char *>(deloptions_.c_str()), NULL, stderr)) {
261	chuckv	1242
262			qh visit_id++;
263			FORALLfacets {
264			if (!facet->upperdelaunay) {
265			facet->visitid= qh visit_id;
266			qh_makeridges(facet);
267			FOREACHridge_(facet->ridges) {
268			neighbor= otherfacet_(ridge, facet);
269			if (neighbor->visitid != qh visit_id) {
270
271			FOREACHvertex_(ridge->vertices)
272			int id2 = qh_pointid (vertex->point);
273			std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl;
274			}
275			}
276			}
277
278	chuckv	1188
279
280	chuckv	1242
281			}
282	chuckv	1188
283	chuckv	1242	qh_freeqhull(!qh_ALL);
284			qh_memfreeshort(&curlong, &totlong);
285			if (curlong \|\| totlong)
286			std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) "
287			<< totlong << curlong << std::endl;
288			free(surfpt_array);
289			*/
290			return true;
291	chuckv	1141	}
292	chuckv	1097
293	chuckv	1143
294	chuckv	1188
295			RealType ConvexHull::getVolume()
296			{
297	chuckv	1242	return volume_;
298	chuckv	1097	}
299
300	chuckv	1188	void ConvexHull::geomviewHull(const std::string& geomFileName)
301			{
302	chuckv	1137
303	chuckv	1242	FILE *newGeomFile;
304
305			//create new .md file based on old .md file
306			newGeomFile = fopen(geomFileName.c_str(), "w");
307			qh_findgood_all(qh facet_list);
308			for (int i = 0; i < qh_PRINTEND; i++)
309			qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
310
311			fclose(newGeomFile);
312	chuckv	1097	}
313	chuckv	1242	#endif //QHULL
314			#endif //CGAL
315	chuckv	1097
316	chuckv	1188
317