00001
00002
00003
00004
00005
00006
00007
00008
00009 #//! Qhull -- invoke qhull from C++
00010 #//! Compile libqhull and Qhull together due to use of setjmp/longjmp()
00011
00012 #include "QhullError.h"
00013 #include "UsingLibQhull.h"
00014 #include "RboxPoints.h"
00015 #include "QhullQh.h"
00016 #include "QhullFacet.h"
00017 #include "QhullFacetList.h"
00018 #include "Qhull.h"
00019 extern "C" {
00020 #include "libqhull/qhull_a.h"
00021 }
00022
00023 #include <iostream>
00024
00025 using std::cerr;
00026 using std::string;
00027 using std::vector;
00028 using std::ostream;
00029
00030 #ifdef _MSC_VER // Microsoft Visual C++ -- warning level 4
00031 #pragma warning( disable : 4611) // interaction between '_setjmp' and C++ object destruction is non-portable
00032 #pragma warning( disable : 4996) // function was declared deprecated(strcpy, localtime, etc.)
00033 #endif
00034
00035 namespace orgQhull {
00036
00037 #//Global variables
00038
00039 char s_unsupported_options[]=" Fd TI ";
00040 char s_not_output_options[]= " Fd TI A C d E H P Qb QbB Qbb Qc Qf Qg Qi Qm QJ Qr QR Qs Qt Qv Qx Qz Q0 Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 R Tc TC TM TP TR Tv TV TW U v V W ";
00041
00042 #//Constructor, destructor, etc.
00043 Qhull::
00044 Qhull()
00045 : qhull_qh(0)
00046 , qhull_run_id(UsingLibQhull::NOqhRunId)
00047 , origin_point()
00048 , qhull_status(qh_ERRnone)
00049 , qhull_dimension(0)
00050 , run_called(false)
00051 , qh_active(false)
00052 , qhull_message()
00053 , error_stream(0)
00054 , output_stream(0)
00055 , feasiblePoint()
00056 , useOutputStream(false)
00057 {
00058 initializeQhull();
00059 }
00060
00061 Qhull::
00062 Qhull(const RboxPoints &rboxPoints, const char *qhullCommand2)
00063 : qhull_qh(0)
00064 , qhull_run_id(UsingLibQhull::NOqhRunId)
00065 , origin_point()
00066 , qhull_status(qh_ERRnone)
00067 , qhull_dimension(0)
00068 , run_called(false)
00069 , qh_active(false)
00070 , qhull_message()
00071 , error_stream(0)
00072 , output_stream(0)
00073 , feasiblePoint()
00074 , useOutputStream(false)
00075 {
00076 initializeQhull();
00077 runQhull(rboxPoints, qhullCommand2);
00078 }
00079
00080 Qhull::
00081 Qhull(const char *rboxCommand2, int pointDimension, int pointCount, const realT *pointCoordinates, const char *qhullCommand2)
00082 : qhull_qh(0)
00083 , qhull_run_id(UsingLibQhull::NOqhRunId)
00084 , origin_point()
00085 , qhull_status(qh_ERRnone)
00086 , qhull_dimension(0)
00087 , run_called(false)
00088 , qh_active(false)
00089 , qhull_message()
00090 , error_stream(0)
00091 , output_stream(0)
00092 , feasiblePoint()
00093 , useOutputStream(false)
00094 {
00095 initializeQhull();
00096 runQhull(rboxCommand2, pointDimension, pointCount, pointCoordinates, qhullCommand2);
00097 }
00098
00099 Qhull::
00100 Qhull(const Qhull &other)
00101 : qhull_qh(0)
00102 , qhull_run_id(UsingLibQhull::NOqhRunId)
00103 , origin_point()
00104 , qhull_status(qh_ERRnone)
00105 , qhull_dimension(0)
00106 , run_called(false)
00107 , qh_active(false)
00108 , qhull_message(other.qhull_message)
00109 , error_stream(other.error_stream)
00110 , output_stream(other.output_stream)
00111 , feasiblePoint(other.feasiblePoint)
00112 , useOutputStream(other.useOutputStream)
00113 {
00114 if(other.initialized()){
00115 throw QhullError(10069, "Qhull error: can not use Qhull copy constructor if initialized() is true");
00116 }
00117 initializeQhull();
00118 }
00119
00120 Qhull & Qhull::
00121 operator=(const Qhull &other)
00122 {
00123 if(other.initialized() || initialized()){
00124 throw QhullError(10070, "Qhull error: can not use Qhull copy assignment if initialized() is true");
00125 }
00126 qhull_message= other.qhull_message;
00127 error_stream= other.error_stream;
00128 output_stream= other.output_stream;
00129 feasiblePoint= other.feasiblePoint;
00130 useOutputStream= other.useOutputStream;
00131 return *this;
00132 }
00133
00134 void Qhull::
00135 initializeQhull()
00136 {
00137 #if qh_QHpointer
00138 qhull_qh= new QhullQh;
00139 qhull_qh->old_qhstat= qh_qhstat;
00140 qhull_qh->old_tempstack= qhmem.tempstack;
00141 qh_qh= 0;
00142 qh_qhstat= 0;
00143 #else
00144 qhull_qh= new (&qh_qh) QhullQh;
00145 qhull_qh->old_qhstat= &qh_qhstat;
00146 qhull_qh->old_tempstack= qhmem.tempstack;
00147 #endif
00148 qhmem.tempstack= 0;
00149 qhull_run_id= qhull_qh->run_id;
00150 }
00151
00152 Qhull::
00153 ~Qhull() throw()
00154 {
00156 UsingLibQhull q(this, QhullError::NOthrow);
00157 if(q.defined()){
00158 int exitCode = setjmp(qh errexit);
00159 if(!exitCode){
00160 #if qh_QHpointer
00161 delete qhull_qh;
00162
00163 qh_qh= 0;
00164 #else
00165 qhull_qh->~QhullQh();
00166 qhull_qh= 0;
00167 #endif
00168 qhull_run_id= UsingLibQhull::NOqhRunId;
00169
00170 if(hasQhullMessage()){
00171 cerr<< "\nQhull output at end\n";
00172 cerr<<qhullMessage();
00173 clearQhullMessage();
00174 }
00175 }
00176 maybeThrowQhullMessage(exitCode, QhullError::NOthrow);
00177 }
00178 s_qhull_output= 0;
00179 }
00180
00181 #//Messaging
00182
00183 void Qhull::
00184 appendQhullMessage(const string &s)
00185 {
00186 if(output_stream && useOutputStream && qh USEstdout){
00187 *output_stream << s;
00188 }else if(error_stream){
00189 *error_stream << s;
00190 }else{
00191 qhull_message += s;
00192 }
00193 }
00194
00196 void Qhull::
00197 clearQhullMessage()
00198 {
00199 qhull_status= qh_ERRnone;
00200 qhull_message.clear();
00201 RoadError::clearGlobalLog();
00202 }
00203
00205 bool Qhull::
00206 hasQhullMessage() const
00207 {
00208 return (!qhull_message.empty() || qhull_status!=qh_ERRnone);
00209
00210 }
00211
00213 std::string Qhull::
00214 qhullMessage() const
00215 {
00216 if(qhull_message.empty() && qhull_status!=qh_ERRnone){
00217 return "qhull: no message for error. Check cerr or error stream\n";
00218 }else{
00219 return qhull_message;
00220 }
00221 }
00222
00223 int Qhull::
00224 qhullStatus() const
00225 {
00226 return qhull_status;
00227 }
00228
00229 void Qhull::
00230 setErrorStream(ostream *os)
00231 {
00232 error_stream= os;
00233 }
00234
00236 void Qhull::
00237 setOutputStream(ostream *os)
00238 {
00239 output_stream= os;
00240 useOutputStream= (os!=0);
00241 }
00242
00243 #//GetSet
00244
00245 void Qhull::
00246 checkIfQhullInitialized()
00247 {
00248 if(!initialized()){
00249 throw QhullError(10023, "Qhull error: checkIfQhullInitialized failed. Call runQhull() first.");
00250 }
00251 }
00252
00254 int Qhull::
00255 runId()
00256 {
00257 UsingLibQhull u(this);
00258 QHULL_UNUSED(u);
00259
00260 return qhull_run_id;
00261 }
00262
00263
00264 #//GetValue
00265
00266 double Qhull::
00267 area(){
00268 checkIfQhullInitialized();
00269 UsingLibQhull q(this);
00270 if(!qh hasAreaVolume){
00271 int exitCode = setjmp(qh errexit);
00272 if(!exitCode){
00273 qh_getarea(qh facet_list);
00274 }
00275 maybeThrowQhullMessage(exitCode);
00276 }
00277 return qh totarea;
00278 }
00279
00280 double Qhull::
00281 volume(){
00282 checkIfQhullInitialized();
00283 UsingLibQhull q(this);
00284 if(!qh hasAreaVolume){
00285 int exitCode = setjmp(qh errexit);
00286 if(!exitCode){
00287 qh_getarea(qh facet_list);
00288 }
00289 maybeThrowQhullMessage(exitCode);
00290 }
00291 return qh totvol;
00292 }
00293
00294 #//ForEach
00295
00299 void Qhull::
00300 defineVertexNeighborFacets(){
00301 checkIfQhullInitialized();
00302 UsingLibQhull q(this);
00303 if(!qh hasAreaVolume){
00304 int exitCode = setjmp(qh errexit);
00305 if(!exitCode){
00306 qh_vertexneighbors();
00307 }
00308 maybeThrowQhullMessage(exitCode);
00309 }
00310 }
00311
00312 QhullFacetList Qhull::
00313 facetList() const{
00314 return QhullFacetList(beginFacet(), endFacet());
00315 }
00316
00317 QhullPoints Qhull::
00318 points() const
00319 {
00320 return QhullPoints(hullDimension(), qhull_qh->num_points*hullDimension(), qhull_qh->first_point);
00321 }
00322
00323 QhullPointSet Qhull::
00324 otherPoints() const
00325 {
00326 return QhullPointSet(hullDimension(), qhull_qh->other_points);
00327 }
00328
00330 QhullVertexList Qhull::
00331 vertexList() const{
00332 return QhullVertexList(beginVertex(), endVertex());
00333 }
00334
00335 #//Modify
00336
00337 void Qhull::
00338 outputQhull()
00339 {
00340 checkIfQhullInitialized();
00341 UsingLibQhull q(this);
00342 int exitCode = setjmp(qh errexit);
00343 if(!exitCode){
00344 qh_produce_output2();
00345 }
00346 maybeThrowQhullMessage(exitCode);
00347 }
00348
00349 void Qhull::
00350 outputQhull(const char *outputflags)
00351 {
00352 checkIfQhullInitialized();
00353 UsingLibQhull q(this);
00354 string cmd(" ");
00355 cmd += outputflags;
00356 char *command= const_cast<char*>(cmd.c_str());
00357 int exitCode = setjmp(qh errexit);
00358 if(!exitCode){
00359 qh_clear_outputflags();
00360 char *s = qh qhull_command + strlen(qh qhull_command) + 1;
00361 strncat(qh qhull_command, command, sizeof(qh qhull_command)-strlen(qh qhull_command)-1);
00362 qh_checkflags(command, s_not_output_options);
00363 qh_initflags(s);
00364 qh_initqhull_outputflags();
00365 if(qh KEEPminArea < REALmax/2
00366 || (0 != qh KEEParea + qh KEEPmerge + qh GOODvertex
00367 + qh GOODthreshold + qh GOODpoint + qh SPLITthresholds)){
00368 facetT *facet;
00369 qh ONLYgood= False;
00370 FORALLfacet_(qh facet_list) {
00371 facet->good= True;
00372 }
00373 qh_prepare_output();
00374 }
00375 qh_produce_output2();
00376 if(qh VERIFYoutput && !qh STOPpoint && !qh STOPcone){
00377 qh_check_points();
00378 }
00379 }
00380 maybeThrowQhullMessage(exitCode);
00381 }
00382
00383 void Qhull::
00384 runQhull(const RboxPoints &rboxPoints, const char *qhullCommand2)
00385 {
00386 runQhull(rboxPoints.comment().c_str(), rboxPoints.dimension(), rboxPoints.count(), &*rboxPoints.coordinates(), qhullCommand2);
00387 }
00388
00391 void Qhull::
00392 runQhull(const char *rboxCommand2, int pointDimension, int pointCount, const realT *rboxPoints, const char *qhullCommand2)
00393 {
00394 if(run_called){
00395 throw QhullError(10027, "Qhull error: runQhull called twice. Only one call allowed.");
00396 }
00397 run_called= true;
00398 string s("qhull ");
00399 s += qhullCommand2;
00400 char *command= const_cast<char*>(s.c_str());
00401 UsingLibQhull q(this);
00402 int exitCode = setjmp(qh errexit);
00403 if(!exitCode){
00404 qh_checkflags(command, s_unsupported_options);
00405 qh_initflags(command);
00406 *qh rbox_command= '\0';
00407 strncat( qh rbox_command, rboxCommand2, sizeof(qh rbox_command)-1);
00408 if(qh DELAUNAY){
00409 qh PROJECTdelaunay= True;
00410 }
00411 pointT *newPoints= const_cast<pointT*>(rboxPoints);
00412 int newDimension= pointDimension;
00413 int newIsMalloc= False;
00414 if(qh HALFspace){
00415 --newDimension;
00416 initializeFeasiblePoint(newDimension);
00417 newPoints= qh_sethalfspace_all(pointDimension, pointCount, newPoints, qh feasible_point);
00418 newIsMalloc= True;
00419 }
00420 qh_init_B(newPoints, pointCount, newDimension, newIsMalloc);
00421 qhull_dimension= (qh DELAUNAY ? qh hull_dim - 1 : qh hull_dim);
00422 qh_qhull();
00423 qh_check_output();
00424 qh_prepare_output();
00425 if(qh VERIFYoutput && !qh STOPpoint && !qh STOPcone){
00426 qh_check_points();
00427 }
00428 }
00429 for(int k= qhull_dimension; k--; ){
00430 origin_point << 0.0;
00431 }
00432 maybeThrowQhullMessage(exitCode);
00433 }
00434
00435 #//Helpers -- be careful of allocating C++ objects due to setjmp/longjmp() error handling by qh_... routines
00436
00437 void Qhull::
00438 initializeFeasiblePoint(int hulldim)
00439 {
00440 if(qh feasible_string){
00441 qh_setfeasible(hulldim);
00442 }else{
00443 if(feasiblePoint.empty()){
00444 qh_fprintf(qh ferr, 6209, "qhull error: missing feasible point for halfspace intersection. Use option 'Hn,n' or set qh.feasiblePoint\n");
00445 qh_errexit(qh_ERRmem, NULL, NULL);
00446 }
00447 if(feasiblePoint.size()!=(size_t)hulldim){
00448 qh_fprintf(qh ferr, 6210, "qhull error: dimension of feasiblePoint should be %d. It is %u", hulldim, feasiblePoint.size());
00449 qh_errexit(qh_ERRmem, NULL, NULL);
00450 }
00451 if (!(qh feasible_point= (coordT*)qh_malloc(hulldim * sizeof(coordT)))) {
00452 qh_fprintf(qh ferr, 6202, "qhull error: insufficient memory for feasible point\n");
00453 qh_errexit(qh_ERRmem, NULL, NULL);
00454 }
00455 coordT *t= qh feasible_point;
00456
00457 for(Coordinates::ConstIterator p=feasiblePoint.begin(); p<feasiblePoint.end(); p++){
00458 *t++= *p;
00459 }
00460 }
00461 }
00462
00463 void Qhull::
00464 maybeThrowQhullMessage(int exitCode)
00465 {
00466 if(qhull_status==qh_ERRnone){
00467 qhull_status= exitCode;
00468 }
00469 if(qhull_status!=qh_ERRnone){
00470 QhullError e(qhull_status, qhull_message);
00471 clearQhullMessage();
00472 throw e;
00473 }
00474 }
00475
00476 void Qhull::
00477 maybeThrowQhullMessage(int exitCode, int noThrow) throw()
00478 {
00479 QHULL_UNUSED(noThrow);
00480
00481 if(qhull_status==qh_ERRnone){
00482 qhull_status= exitCode;
00483 }
00484 if(qhull_status!=qh_ERRnone){
00485 QhullError e(qhull_status, qhull_message);
00486 e.logError();
00487 }
00488 }
00489
00490 }
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505 extern "C"
00506 void qh_fprintf(FILE *fp, int msgcode, const char *fmt, ... ) {
00507 va_list args;
00508
00509 using namespace orgQhull;
00510
00511 if(!s_qhull_output){
00512 fprintf(stderr, "QH10025 Qhull error: UsingLibQhull not declared prior to calling qh_...(). s_qhull_output==NULL.\n");
00513 qh_exit(10025);
00514 }
00515 Qhull *out= s_qhull_output;
00516 va_start(args, fmt);
00517 if(msgcode<MSG_OUTPUT || fp == qh_FILEstderr){
00518 if(msgcode>=MSG_ERROR && msgcode<MSG_WARNING){
00519 if(out->qhull_status<MSG_ERROR || out->qhull_status>=MSG_WARNING){
00520 out->qhull_status= msgcode;
00521 }
00522 }
00523 char newMessage[MSG_MAXLEN];
00524
00525 vsnprintf(newMessage, sizeof(newMessage), fmt, args);
00526 out->appendQhullMessage(newMessage);
00527 va_end(args);
00528 return;
00529 }
00530 if(out->output_stream && out->useOutputStream){
00531 char newMessage[MSG_MAXLEN];
00532 vsnprintf(newMessage, sizeof(newMessage), fmt, args);
00533 *out->output_stream << newMessage;
00534 va_end(args);
00535 return;
00536 }
00537
00538 char newMessage[MSG_MAXLEN];
00539 vsnprintf(newMessage, sizeof(newMessage), fmt, args);
00540 out->appendQhullMessage(newMessage);
00541 va_end(args);
00542 }
00543