# | Line 40 | Line 40 | |
---|---|---|
40 | * | |
41 | * ConvexHull.cpp | |
42 | * | |
43 | < | * Purpose: To calculate convexhull, hull volume and radius |
44 | < | * using the CGAL library. |
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.1 2006-12-14 19:32:32 chuckv Exp $ |
47 | > | * @version $Id: ConvexHull.cpp,v 1.7 2008-06-18 17:03:30 chuckv Exp $ |
48 | * | |
49 | */ | |
50 | ||
51 | < | #include "math/ConvexHull.hpp" |
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 | < | using namespace oopse; |
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 | < | points.reserve(pos.size()); |
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 | < | CGAL::convex_hull_3(points.begin(), points.end(), ch_object); |
102 | < | // Make sure hull is a polyhedron |
103 | < | if ( CGAL::assign(ch_polyhedron, ch_object) ) |
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 | < | return true; |
108 | > | std::cerr << (*p_it).x() << std::endl; |
109 | } | |
110 | < | else |
111 | < | { |
112 | < | return false; |
113 | < | } |
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 | ||
87 | – | RealType ConvexHull::getVolume() |
88 | – | { |
89 | – | /* |
90 | – | std::list<Point> L; |
91 | – | L.push_front(Point(0,0,0)); |
92 | – | L.push_front(Point(1,0,0)); |
93 | – | L.push_front(Point(0,1,0)); |
119 | ||
95 | – | Triangulation T(L.begin(), L.end()); |
120 | ||
121 | < | int n = T.number_of_vertices(); |
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 | < | // insertion from a vector : |
129 | < | std::vector<Point> V(3); |
130 | < | V[0] = Point(0,0,1); |
131 | < | V[1] = Point(1,1,1); |
132 | < | V[2] = Point(2,2,2); |
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 | < | n = n + T.insert(V.begin(), V.end()); |
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 | ||
107 | – | assert( n == 6 ); // 6 points have been inserted |
108 | – | assert( T.is_valid() ); // checking validity of T |
212 | ||
213 | < | double sum_v = 0; |
214 | < | for(Triangulation::Cell_iterator c_it = T.cells_begin(); c_it != T.cells_end(); ++c_it) |
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 | < | sum_v += T.tetrahedron(c_it).volume(); |
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 | < | std::cout << "sum_v " << sum_v << std::endl; |
253 | < | */ |
254 | < | return 0.0; |
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 | < | std::ofstream newGeomFile; |
304 | < | |
303 | > | FILE *newGeomFile; |
304 | > | |
305 | //create new .md file based on old .md file | |
306 | < | newGeomFile.open(geomFileName.c_str()); |
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 | ||
128 | – | // Write polyhedron in Object File Format (OFF). |
129 | – | CGAL::set_ascii_mode( std::cout); |
130 | – | newGeomFile << "OFF" << std::endl << ch_polyhedron.size_of_vertices() << ' ' |
131 | – | << ch_polyhedron.size_of_facets() << " 0" << std::endl; |
132 | – | std::copy( ch_polyhedron.points_begin(), ch_polyhedron.points_end(), |
133 | – | std::ostream_iterator<Point_3>( newGeomFile, "\n")); |
134 | – | for ( Facet_iterator i = ch_polyhedron.facets_begin(); i != ch_polyhedron.facets_end(); ++i) |
135 | – | { |
136 | – | Halfedge_facet_circulator j = i->facet_begin(); |
137 | – | // Facets in polyhedral surfaces are at least triangles. |
138 | – | CGAL_assertion( CGAL::circulator_size(j) >= 3); |
139 | – | newGeomFile << CGAL::circulator_size(j) << ' '; |
140 | – | do |
141 | – | { |
142 | – | newGeomFile << ' ' << std::distance(ch_polyhedron.vertices_begin(), j->vertex()); |
143 | – | } |
144 | – | while ( ++j != i->facet_begin()); |
145 | – | newGeomFile << std::endl; |
146 | – | } |
316 | ||
148 | – | newGeomFile.close(); |
317 | ||
150 | – | |
151 | – | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |