# | Line 44 | Line 44 | |
---|---|---|
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.4 2007-05-30 19:51:07 chuckv Exp $ |
47 | > | * @version $Id: ConvexHull.cpp,v 1.6 2008-05-14 14:31:48 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 | < | |
58 | > | #include "utils/simError.h" |
59 | using namespace oopse; | |
60 | ||
61 | < | ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt FA") { |
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; |
79 | > | |
80 | > | |
81 | > | ConvexHull::ConvexHull(){} |
82 | > | |
83 | > | 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 | ||
108 | < | bool ConvexHull::genHull(std::vector<Vector3d> pos) { |
109 | < | FILE *outfile = stdout; |
110 | < | FILE *errfile = stderr; |
111 | < | facetT *facet; |
112 | < | int exitcode; |
113 | < | boolT ismalloc = False; |
114 | < | int curlong,totlong; |
108 | > | |
109 | > | |
110 | > | #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 | > | |
117 | > | bool ConvexHull::genHull(std::vector<Vector3d> pos) |
118 | > | { |
119 | ||
120 | + | |
121 | int numpoints = pos.size(); | |
122 | < | |
123 | < | coordT* points; |
124 | < | points = (coordT*) malloc(sizeof(coordT) * (numpoints*dim_)); |
125 | < | |
126 | < | for (int i=0; i<numpoints; i++) { |
127 | < | points[dim_ * i] = pos[i][0]; |
128 | < | points[dim_ * i + 1] = pos[i][1]; |
129 | < | points[dim_ * i + 2] = pos[i][2]; |
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 | ||
80 | – | qh_initflags (const_cast<char *>(options_.c_str())); |
81 | – | qh_init_B (points, numpoints, dim_, ismalloc); |
82 | – | qh_qhull(); |
83 | – | qh_check_output(); |
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 | + | |
146 | + | 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 | + | |
175 | + | 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 | + | |
192 | + | } |
193 | + | |
194 | + | } |
195 | + | |
196 | + | */ |
197 | + | |
198 | + | |
199 | + | } |
200 | + | |
201 | + | |
202 | + | |
203 | + | |
204 | qh_getarea(qh facet_list); | |
205 | < | volume_ = qh totvol; |
206 | < | area_ = qh totarea; |
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); |
221 | > | qh_memfreeshort(&curlong, &totlong); |
222 | if (curlong || totlong) | |
223 | < | fprintf (errfile, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", |
224 | < | totlong, curlong); |
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 | < | |
243 | > | /* |
244 | > | std::string deloptions_ = "qhull d Qt"; |
245 | > | facetT *facet, *neighbor; |
246 | > | ridgeT *ridge, **ridgep; |
247 | > | |
248 | > | if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc, |
249 | > | const_cast<char *>(deloptions_.c_str()), NULL, stderr)) { |
250 | > | |
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 | > | |
268 | > | |
269 | > | |
270 | > | } |
271 | > | |
272 | > | 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 | } | |
281 | ||
282 | ||
283 | < | RealType ConvexHull::getVolume() { |
283 | > | |
284 | > | RealType ConvexHull::getVolume() |
285 | > | { |
286 | return volume_; | |
287 | } | |
288 | ||
289 | < | void ConvexHull::geomviewHull(const std::string& geomFileName) { |
289 | > | void ConvexHull::geomviewHull(const std::string& geomFileName) |
290 | > | { |
291 | > | |
292 | > | FILE *newGeomFile; |
293 | ||
110 | – | std::ofstream newGeomFile; |
111 | – | |
294 | //create new .md file based on old .md file | |
295 | < | newGeomFile.open(geomFileName.c_str()); |
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 | < | |
116 | < | newGeomFile.close(); |
300 | > | fclose(newGeomFile); |
301 | } | |
302 | + | #endif //QHULL |
303 | + | #endif //CGAL |
304 | + | |
305 | + | |
306 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |