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 (17 years, 4 months ago) by chuckv
File size:	8826 byte(s)
Log Message:	Added files to module section of Makefile.
#	Content
1	/* 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	* Purpose: To calculate convexhull, hull volume libqhull.
44	*
45	* Created by Charles F. Vardeman II on 11 Dec 2006.
46	* @author Charles F. Vardeman II
47	* @version $Id: ConvexHull.cpp,v 1.7 2008-06-18 17:03:30 chuckv Exp $
48	*
49	*/
50
51	/* Standard includes independent of library */
52	#include <iostream>
53	#include <fstream>
54	#include <list>
55	#include <algorithm>
56	#include <iterator>
57	#include "math/ConvexHull.hpp"
58	#include "utils/simError.h"
59	using namespace oopse;
60
61	/* CGAL version of convex hull first then QHULL */
62	#ifdef HAVE_CGAL
63
64	#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
70	#include <CGAL/Quotient.h>
71	#include <CGAL/MP_Float.h>
72	#include <CGAL/Lazy_exact_nt.h>
73
74	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_3;
79
80
81	ConvexHull::ConvexHull(){}
82
83	bool ConvexHull::genHull(std::vector<Vector3d> pos)
84	{
85
86	std::vector<Point_3> 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_3 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	CGAL::Object ch_object_;
98	Polyhedron_3 polyhedron;
99
100	// compute convex hull
101	std::cerr << "Creating hull" << std::endl;
102	CGAL::convex_hull_3(points.begin(), points.end(), ch_object_);
103	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	std::cout<< v.point()<<std::endl;
115	}
116	*/
117	}
118
119
120
121	#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
128	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
157	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
186	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
203	}
204
205	}
206
207	*/
208
209
210	}
211
212
213
214
215	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	std::string deloptions_ = "qhull d Qt";
256	facetT facet, neighbor;
257	ridgeT ridge, *ridgep;
258
259	if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc,
260	const_cast<char *>(deloptions_.c_str()), NULL, stderr)) {
261
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
279
280
281	}
282
283	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	}
292
293
294
295	RealType ConvexHull::getVolume()
296	{
297	return volume_;
298	}
299
300	void ConvexHull::geomviewHull(const std::string& geomFileName)
301	{
302
303	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	}
313	#endif //QHULL
314	#endif //CGAL
315
316
317