# | 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.2 2007-05-29 22:50:14 chuckv Exp $ |
47 | > | * @version $Id: ConvexHull.cpp,v 1.5 2007-11-22 16:39:45 chuckv Exp $ |
48 | * | |
49 | */ | |
50 | ||
51 | – | #include "math/ConvexHull.hpp" |
51 | #include <iostream> | |
52 | #include <fstream> | |
53 | < | |
54 | < | char options[] = "qhull Qt FA"; |
55 | < | int dim_ = 3; |
56 | < | |
53 | > | #include <list> |
54 | > | #include <algorithm> |
55 | > | #include <iterator> |
56 | > | #include "math/ConvexHull.hpp" |
57 | > | #include "utils/simError.h" |
58 | using namespace oopse; | |
59 | ||
60 | < | ConvexHull::ConvexHull(){} |
60 | > | /* Old options Qt Qu Qg QG0 FA */ |
61 | > | ConvexHull::ConvexHull() : dim_(3), options_("qhull Qt Qci Tcv Pp") |
62 | > | //ConvexHull::ConvexHull() : dim_(3), options_("qhull d Qbb Qt i") |
63 | > | {} |
64 | ||
62 | – | |
65 | bool ConvexHull::genHull(std::vector<Vector3d> pos) | |
66 | { | |
67 | < | FILE *outfile = stdout; |
68 | < | FILE *errfile = stderr; |
67 | < | facetT *facet; |
68 | < | int exitcode; |
69 | < | boolT ismalloc = False; |
70 | < | int curlong,totlong; |
71 | < | |
67 | > | |
68 | > | |
69 | int numpoints = pos.size(); | |
70 | < | |
71 | < | coordT points[numpoints][dim_]; |
72 | < | |
73 | < | for (int i=0; i<numpoints; i++) |
74 | < | { |
75 | < | points[i][0] = pos[i][0]; |
76 | < | points[i][1] = pos[i][1]; |
77 | < | points[i][2] = pos[i][2]; |
70 | > | coordT* pt_array; |
71 | > | coordT* surfpt_array; |
72 | > | std::list<int> surface_atoms; |
73 | > | FILE *outdummy = NULL; |
74 | > | FILE *errdummy = NULL; |
75 | > | |
76 | > | pt_array = (coordT*) malloc(sizeof(coordT) * (numpoints * dim_)); |
77 | > | |
78 | > | |
79 | > | for (int i = 0; i < numpoints; i++) { |
80 | > | pt_array[dim_ * i] = pos[i][0]; |
81 | > | pt_array[dim_ * i + 1] = pos[i][1]; |
82 | > | pt_array[dim_ * i + 2] = pos[i][2]; |
83 | } | |
82 | – | |
84 | ||
85 | ||
86 | < | qh_initflags (options); |
87 | < | qh_init_B (points[0], numpoints, dim_, ismalloc); |
88 | < | qh_qhull(); |
89 | < | qh_check_output(); |
90 | < | |
91 | < | |
92 | < | |
93 | < | qh_getarea(qh facet_list); |
94 | < | volume_ = qh totvol; |
95 | < | area_ = qh totarea; |
96 | < | |
97 | < | |
98 | < | |
99 | < | qh_freeqhull(!qh_ALL); |
100 | < | qh_memfreeshort (&curlong, &totlong); |
86 | > | |
87 | > | |
88 | > | /* |
89 | > | qh_initflags(const_cast<char *>(options_.c_str())); |
90 | > | qh_init_B(pospoints, numpoints, dim_, ismalloc); |
91 | > | qh_qhull(); |
92 | > | qh_check_output(); |
93 | > | |
94 | > | qh_produce_output(); |
95 | > | */ |
96 | > | boolT ismalloc = False; |
97 | > | |
98 | > | if (!qh_new_qhull(dim_, numpoints, pt_array, ismalloc, |
99 | > | const_cast<char *>(options_.c_str()), NULL, stderr)) { |
100 | > | |
101 | > | vertexT *vertex, **vertexp; |
102 | > | facetT *facet; |
103 | > | setT *vertices; |
104 | > | unsigned int nf = qh num_facets; |
105 | > | |
106 | > | //Matrix idx(nf, dim); |
107 | > | /* |
108 | > | int j, i = 0, id = 0; |
109 | > | |
110 | > | int id2 = 0; |
111 | > | coordT *point,**pointp; |
112 | > | realT dist; |
113 | > | FORALLfacets { |
114 | > | j = 0; |
115 | > | |
116 | > | if (!facet->simplicial){ |
117 | > | // should never happen with Qt |
118 | > | sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); |
119 | > | painCave.isFatal = 0; |
120 | > | simError(); |
121 | > | } |
122 | > | |
123 | > | vertices = qh_facet3vertex(facet); |
124 | > | FOREACHvertex_(vertices){ |
125 | > | id = qh_pointid(vertex->point); |
126 | > | surface_atoms.push_back(id); |
127 | > | //std::cout << "Ag " << pos[id][0] << " " << pos[id][1] << " " << pos[id][2]<< std::endl; |
128 | > | } |
129 | > | qh_settempfree(&vertices); |
130 | > | |
131 | > | FOREACHpoint_(facet->coplanarset){ |
132 | > | vertex= qh_nearvertex (facet, point, &dist); |
133 | > | //id= qh_pointid (vertex->point); |
134 | > | id2= qh_pointid (point); |
135 | > | surface_atoms.push_back(id2); |
136 | > | //std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl; |
137 | > | //printf ("%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist); |
138 | > | //std::cout << "Neighbors are: %d $d %d\n" << id << id2 << facet->id; |
139 | > | |
140 | > | } |
141 | > | |
142 | > | } |
143 | > | |
144 | > | */ |
145 | > | |
146 | > | |
147 | > | } |
148 | > | |
149 | > | |
150 | > | |
151 | > | |
152 | > | qh_getarea(qh facet_list); |
153 | > | volume_ = qh totvol; |
154 | > | area_ = qh totarea; |
155 | > | // FILE *newGeomFile; |
156 | > | |
157 | > | |
158 | > | /* |
159 | > | FORALLfacets { |
160 | > | for (int k=0; k < qh hull_dim; k++) |
161 | > | printf ("%6.2g ", facet->normal[k]); |
162 | > | printf ("\n"); |
163 | > | } |
164 | > | */ |
165 | > | |
166 | > | int curlong,totlong; |
167 | > | // geomviewHull("junk.off"); |
168 | > | qh_freeqhull(!qh_ALL); |
169 | > | qh_memfreeshort(&curlong, &totlong); |
170 | if (curlong || totlong) | |
171 | < | fprintf (errfile, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", |
172 | < | totlong, curlong); |
171 | > | std::cerr << "qhull internal warning (main): did not free %d bytes of long memory (%d pieces) " |
172 | > | << totlong << curlong << std::endl; |
173 | > | free(pt_array); |
174 | > | /* |
175 | > | int j = 0; |
176 | > | surface_atoms.sort(); |
177 | > | surface_atoms.unique(); |
178 | > | surfpt_array = (coordT*) malloc(sizeof(coordT) * (surface_atoms.size() * dim_)); |
179 | > | for(std::list<int>::iterator list_iter = surface_atoms.begin(); |
180 | > | list_iter != surface_atoms.end(); list_iter++) |
181 | > | { |
182 | > | int i = *list_iter; |
183 | > | //surfpt_array[dim_ * j] = pos[i][0]; |
184 | > | //surfpt_array[dim_ * j + 1] = pos[i][1]; |
185 | > | //surfpt_array[dim_ * j + 2] = pos[i][2]; |
186 | > | std::cout << "Ag " << pos[i][0] << " " << pos[i][1] << " "<< pos[i][2] << std::endl; |
187 | > | j++; |
188 | > | } |
189 | > | */ |
190 | ||
191 | + | /* |
192 | + | std::string deloptions_ = "qhull d Qt"; |
193 | + | facetT *facet, *neighbor; |
194 | + | ridgeT *ridge, **ridgep; |
195 | ||
196 | + | if (!qh_new_qhull(dim_, surface_atoms.size(), surfpt_array, ismalloc, |
197 | + | const_cast<char *>(deloptions_.c_str()), NULL, stderr)) { |
198 | + | |
199 | + | qh visit_id++; |
200 | + | FORALLfacets { |
201 | + | if (!facet->upperdelaunay) { |
202 | + | facet->visitid= qh visit_id; |
203 | + | qh_makeridges(facet); |
204 | + | FOREACHridge_(facet->ridges) { |
205 | + | neighbor= otherfacet_(ridge, facet); |
206 | + | if (neighbor->visitid != qh visit_id) { |
207 | + | |
208 | + | FOREACHvertex_(ridge->vertices) |
209 | + | int id2 = qh_pointid (vertex->point); |
210 | + | std::cout << "Ag " << pos[id2][0] << " " << pos[id2][1] << " " << pos[id2][2]<< std::endl; |
211 | + | } |
212 | + | } |
213 | + | } |
214 | + | |
215 | + | |
216 | + | |
217 | + | |
218 | + | } |
219 | + | |
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(surfpt_array); |
226 | + | */ |
227 | return true; | |
228 | } | |
229 | ||
230 | ||
231 | + | |
232 | RealType ConvexHull::getVolume() | |
233 | { | |
234 | return volume_; | |
# | Line 114 | Line 237 | void ConvexHull::geomviewHull(const std::string& geomF | |
237 | void ConvexHull::geomviewHull(const std::string& geomFileName) | |
238 | { | |
239 | ||
240 | < | std::ofstream newGeomFile; |
240 | > | FILE *newGeomFile; |
241 | ||
242 | < | //create new .md file based on old .md file |
243 | < | newGeomFile.open(geomFileName.c_str()); |
242 | > | //create new .md file based on old .md file |
243 | > | newGeomFile = fopen(geomFileName.c_str(), "w"); |
244 | > | qh_findgood_all(qh facet_list); |
245 | > | for (int i = 0; i < qh_PRINTEND; i++) |
246 | > | qh_printfacets(newGeomFile, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); |
247 | ||
248 | + | fclose(newGeomFile); |
249 | + | } |
250 | ||
123 | – | newGeomFile.close(); |
251 | ||
252 | ||
253 | < | } |
253 | > |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |