Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Analysis.AlgorithmBehavior/qhull-2012.1/src/libqhullcpp/QhullFacet.cpp @ 10207

Last change on this file since 10207 was 10207, checked in by ascheibe, 11 years ago

#1886 added a unit test for volume calculation and the qhull library
File size: 15.8 KB
Line 
1/****************************************************************************
2**
3** Copyright (c) 2008-2012 C.B. Barber. All rights reserved.
4** $Id: //main/2011/qhull/src/libqhullcpp/QhullFacet.cpp#7 $$Change: 1464 $
5** $DateTime: 2012/01/25 22:58:41 $$Author: bbarber $
6**
7****************************************************************************/
8
9#//! QhullFacet -- Qhull's facet structure, facetT, as a C++ class
10
11#include "QhullError.h"
12#include "QhullSet.h"
13#include "QhullPoint.h"
14#include "QhullPointSet.h"
15#include "QhullRidge.h"
16#include "QhullFacet.h"
17#include "QhullFacetSet.h"
18#include "QhullVertex.h"
19
20#include <ostream>
21
22using std::endl;
23using std::string;
24using std::vector;
25using std::ostream;
26
27#ifdef _MSC_VER  // Microsoft Visual C++ -- warning level 4
28#pragma warning( disable : 4611)  // interaction between '_setjmp' and C++ object destruction is non-portable
29#pragma warning( disable : 4996)  // function was declared deprecated(strcpy, localtime, etc.)
30#endif
31
32namespace orgQhull {
33
34#//class statics
35facetT QhullFacet::
36        s_empty_facet= {0,0,0,0,{0},
37                0,0,0,0,0,
38                0,0,0,0,0,
39                0,0,0,0,0,
40                0,0,0,0,0,
41                0,0,0,0,0,
42                0,0,0,0,0,
43                0,0,0,0};
44
45#//GetSet
46
47int QhullFacet::
48dimension() const
49{
50    if(qh_facet->ridges){
51        setT *s= qh_facet->ridges;
52        ridgeT *r= reinterpret_cast<ridgeT *>(SETfirst_(s));
53        return r ? QhullSetBase::count(r->vertices)+1 : 0;
54    }else{
55        return QhullSetBase::count(qh_facet->vertices);
56    }
57}//dimension
58
59//! Return voronoi center or facet centrum.  Derived from qh_printcenter [io.c]
60//! printFormat=qh_PRINTtriangles if return centrum of a Delaunay facet
61//! Sets center if needed
62//! Code duplicated for PrintCenter and getCenter
63//! .empty() if none or qh_INFINITE
64QhullPoint QhullFacet::
65getCenter(int qhRunId, qh_PRINT printFormat)
66{
67    UsingLibQhull q(qhRunId);
68
69    if(qh CENTERtype==qh_ASvoronoi){
70        if(!qh_facet->normal || !qh_facet->upperdelaunay || !qh ATinfinity){
71            if(!qh_facet->center){
72                int exitCode = setjmp(qh errexit);
73                if(!exitCode){ // no object creation -- destructors skipped on longjmp()
74                    qh_facet->center= qh_facetcenter(qh_facet->vertices);
75                }
76                q.maybeThrowQhullMessage(exitCode);
77            }
78            return QhullPoint(qh hull_dim-1, qh_facet->center);
79        }
80    }else if(qh CENTERtype==qh_AScentrum){
81        volatile int numCoords= qh hull_dim;
82        if(printFormat==qh_PRINTtriangles && qh DELAUNAY){
83            numCoords--;
84        }
85        if(!qh_facet->center){
86            int exitCode = setjmp(qh errexit);
87            if(!exitCode){ // no object creation -- destructors skipped on longjmp()
88                qh_facet->center= qh_getcentrum(getFacetT());
89            }
90            q.maybeThrowQhullMessage(exitCode);
91        }
92        return QhullPoint(numCoords, qh_facet->center);
93    }
94    return QhullPoint();
95 }//getCenter
96
97//! Return innerplane clearly below the vertices
98//! from io.c[qh_PRINTinner]
99QhullHyperplane QhullFacet::
100innerplane(int qhRunId) const{
101    UsingLibQhull q(qhRunId);
102    realT inner;
103    // Does not error
104    qh_outerinner(const_cast<facetT *>(getFacetT()), NULL, &inner);
105    QhullHyperplane h= hyperplane();
106    h.setOffset(h.offset()-inner); //inner is negative
107    return h;
108}//innerplane
109
110//! Return outerplane clearly above all points
111//! from io.c[qh_PRINTouter]
112QhullHyperplane QhullFacet::
113outerplane(int qhRunId) const{
114    UsingLibQhull q(qhRunId);
115    realT outer;
116    // Does not error
117    qh_outerinner(const_cast<facetT *>(getFacetT()), &outer, NULL);
118    QhullHyperplane h= hyperplane();
119    h.setOffset(h.offset()-outer); //outer is positive
120    return h;
121}//outerplane
122
123//! Set by qh_triangulate for option 'Qt'.
124//! Errors if tricoplanar and facetArea() or qh_getarea() called first.
125QhullFacet QhullFacet::
126tricoplanarOwner() const
127{
128    if(qh_facet->tricoplanar){
129        if(qh_facet->isarea){
130            throw QhullError(10018, "Qhull error: facetArea() or qh_getarea() previously called.  triCoplanarOwner() is not available.");
131        }
132        return qh_facet->f.triowner;
133    }
134    return 0; // FIXUP QH11009 Should false be the NULL facet or empty facet
135}//tricoplanarOwner
136
137QhullPoint QhullFacet::
138voronoiVertex(int qhRunId)
139{
140    if(
141#if qh_QHpointer
142      !qh_qh ||
143#endif
144      qh CENTERtype!=qh_ASvoronoi){
145          throw QhullError(10052, "Error: QhullFacet.voronoiVertex() requires option 'v' (qh_ASvoronoi)");
146    }
147    return getCenter(qhRunId);
148}//voronoiVertex
149
150#//Value
151
152//! Disables tricoplanarOwner()
153double QhullFacet::
154facetArea(int qhRunId)
155{
156    if(!qh_facet->isarea){
157        UsingLibQhull q(qhRunId);
158        int exitCode = setjmp(qh errexit);
159        if(!exitCode){ // no object creation -- destructors skipped on longjmp()
160            qh_facet->f.area= qh_facetarea(qh_facet);
161            qh_facet->isarea= True;
162        }
163        q.maybeThrowQhullMessage(exitCode);
164    }
165    return qh_facet->f.area;
166}//facetArea
167
168#//Foreach
169
170QhullPointSet QhullFacet::
171coplanarPoints() const
172{
173    return QhullPointSet(dimension(), qh_facet->coplanarset);
174}//coplanarPoints
175
176QhullFacetSet QhullFacet::
177neighborFacets() const
178{
179    return QhullFacetSet(qh_facet->neighbors);
180}//neighborFacets
181
182QhullPointSet QhullFacet::
183outsidePoints() const
184{
185    return QhullPointSet(dimension(), qh_facet->outsideset);
186}//outsidePoints
187
188QhullRidgeSet QhullFacet::
189ridges() const
190{
191    return QhullRidgeSet(qh_facet->ridges);
192}//ridges
193
194QhullVertexSet QhullFacet::
195vertices() const
196{
197    return QhullVertexSet(qh_facet->vertices);
198}//vertices
199
200
201}//namespace orgQhull
202
203#//operator<<
204
205using std::ostream;
206
207using orgQhull::QhullFacet;
208using orgQhull::QhullFacetSet;
209using orgQhull::QhullPoint;
210using orgQhull::QhullPointSet;
211using orgQhull::QhullRidge;
212using orgQhull::QhullRidgeSet;
213using orgQhull::QhullSetBase;
214using orgQhull::QhullVertexSet;
215using orgQhull::UsingLibQhull;
216
217ostream &
218operator<<(ostream &os, const QhullFacet::PrintFacet &pr)
219{
220    QhullFacet f= *pr.facet;
221    if(f.getFacetT()==0){ // Special values from set iterator
222        os << " NULLfacet" << endl;
223        return os;
224    }
225    if(f.getFacetT()==qh_MERGEridge){
226        os << " MERGEridge" << endl;
227        return os;
228    }
229    if(f.getFacetT()==qh_DUPLICATEridge){
230        os << " DUPLICATEridge" << endl;
231        return os;
232    }
233    os << f.printHeader(pr.run_id);
234    if(!f.ridges().isEmpty()){
235        os << f.printRidges(pr.run_id);
236    }
237    return os;
238}//operator<< PrintFacet
239
240//! Print Voronoi center or facet centrum to stream.  Same as qh_printcenter [io.c]
241//! Code duplicated for PrintCenter and getCenter
242//! Sets center if needed
243ostream &
244operator<<(ostream &os, const QhullFacet::PrintCenter &pr)
245{
246    facetT *f= pr.facet->getFacetT();
247    if(qh CENTERtype!=qh_ASvoronoi && qh CENTERtype!=qh_AScentrum){
248        return os;
249    }
250    if (pr.message){
251        os << pr.message;
252    }
253    int numCoords;
254    if(qh CENTERtype==qh_ASvoronoi){
255        numCoords= qh hull_dim-1;
256        if(!f->normal || !f->upperdelaunay || !qh ATinfinity){
257            if(!f->center){
258                f->center= qh_facetcenter(f->vertices);
259            }
260            for(int k=0; k<numCoords; k++){
261                os << f->center[k] << " "; // FIXUP QH11010 qh_REAL_1
262            }
263        }else{
264            for(int k=0; k<numCoords; k++){
265                os << qh_INFINITE << " "; // FIXUP QH11010 qh_REAL_1
266            }
267        }
268    }else{ // qh CENTERtype==qh_AScentrum
269        numCoords= qh hull_dim;
270        if(pr.print_format==qh_PRINTtriangles && qh DELAUNAY){
271            numCoords--;
272        }
273        if(!f->center){
274            f->center= qh_getcentrum(f);
275        }
276        for(int k=0; k<numCoords; k++){
277            os << f->center[k] << " "; // FIXUP QH11010 qh_REAL_1
278        }
279    }
280    if(pr.print_format==qh_PRINTgeom && numCoords==2){
281        os << " 0";
282    }
283    os << endl;
284    return os;
285}//operator<< PrintCenter
286
287//! Print flags for facet to stream.  Space prefix.  From qh_printfacetheader [io.c]
288ostream &
289operator<<(ostream &os, const QhullFacet::PrintFlags &p)
290{
291    const facetT *f= p.facet->getFacetT();
292    if(p.message){
293        os << p.message;
294    }
295
296    os << (p.facet->isTopOrient() ? " top" : " bottom");
297    if(p.facet->isSimplicial()){
298        os << " simplicial";
299    }
300    if(p.facet->isTriCoplanar()){
301        os << " tricoplanar";
302    }
303    if(p.facet->isUpperDelaunay()){
304        os << " upperDelaunay";
305    }
306    if(f->visible){
307        os << " visible";
308    }
309    if(f->newfacet){
310        os << " new";
311    }
312    if(f->tested){
313        os << " tested";
314    }
315    if(!f->good){
316        os << " notG";
317    }
318    if(f->seen){
319        os << " seen";
320    }
321    if(f->coplanar){
322        os << " coplanar";
323    }
324    if(f->mergehorizon){
325        os << " mergehorizon";
326    }
327    if(f->keepcentrum){
328        os << " keepcentrum";
329    }
330    if(f->dupridge){
331        os << " dupridge";
332    }
333    if(f->mergeridge && !f->mergeridge2){
334        os << " mergeridge1";
335    }
336    if(f->mergeridge2){
337        os << " mergeridge2";
338    }
339    if(f->newmerge){
340        os << " newmerge";
341    }
342    if(f->flipped){
343        os << " flipped";
344    }
345    if(f->notfurthest){
346        os << " notfurthest";
347    }
348    if(f->degenerate){
349        os << " degenerate";
350    }
351    if(f->redundant){
352        os << " redundant";
353    }
354    os << endl;
355    return os;
356}//operator<< PrintFlags
357
358//! Print header for facet to stream. Space prefix.  From qh_printfacetheader [io.c]
359ostream &
360operator<<(ostream &os, const QhullFacet::PrintHeader &pr)
361{
362    QhullFacet facet= *pr.facet;
363    facetT *f= facet.getFacetT();
364    os << "- f" << facet.id() << endl;
365    os << facet.printFlags("    - flags:");
366    if(f->isarea){
367        os << "    - area: " << f->f.area << endl; //FIXUP QH11010 2.2g
368    }else if(qh NEWfacets && f->visible && f->f.replace){
369        os << "    - replacement: f" << f->f.replace->id << endl;
370    }else if(f->newfacet){
371        if(f->f.samecycle && f->f.samecycle != f){
372            os << "    - shares same visible/horizon as f" << f->f.samecycle->id << endl;
373        }
374    }else if(f->tricoplanar /* !isarea */){
375        if(f->f.triowner){
376            os << "    - owner of normal & centrum is facet f" << f->f.triowner->id << endl;
377        }
378    }else if(f->f.newcycle){
379        os << "    - was horizon to f" << f->f.newcycle->id << endl;
380    }
381    if(f->nummerge){
382        os << "    - merges: " << f->nummerge << endl;
383    }
384    os << facet.hyperplane().print("    - normal: ", "\n    - offset: "); // FIXUP QH11010 %10.7g
385    if(qh CENTERtype==qh_ASvoronoi || f->center){
386        os << facet.printCenter(pr.run_id, qh_PRINTfacets, "    - center: ");
387    }
388#if qh_MAXoutside
389    if(f->maxoutside > qh DISTround){
390        os << "    - maxoutside: " << f->maxoutside << endl; //FIXUP QH11010 %10.7g
391    }
392#endif
393    QhullPointSet ps= facet.outsidePoints();
394    if(!ps.isEmpty()){
395        QhullPoint furthest= ps.last();
396        if (ps.size() < 6) {
397            os << "    - outside set(furthest p" << furthest.id(pr.run_id) << "):" << endl;
398            for(QhullPointSet::iterator i=ps.begin(); i!=ps.end(); ++i){
399                QhullPoint p= *i;
400                os << p.print(pr.run_id, "     ");
401            }
402        }else if(ps.size()<21){
403            os << ps.print(pr.run_id, "    - outside set:");
404        }else{
405            os << "    - outside set:  " << ps.size() << " points.";
406            os << furthest.print(pr.run_id, "  Furthest");
407        }
408#if !qh_COMPUTEfurthest
409        os << "    - furthest distance= " << f->furthestdist << endl; //FIXUP QH11010 %2.2g
410#endif
411    }
412    QhullPointSet cs= facet.coplanarPoints();
413    if(!cs.isEmpty()){
414        QhullPoint furthest= cs.last();
415        if (cs.size() < 6) {
416            os << "    - coplanar set(furthest p" << furthest.id(pr.run_id) << "):" << endl;
417            for(QhullPointSet::iterator i=cs.begin(); i!=cs.end(); ++i){
418                QhullPoint p= *i;
419                os << p.print(pr.run_id, "     ");
420            }
421        }else if(cs.size()<21){
422            os << cs.print(pr.run_id, "    - coplanar set:");
423        }else{
424            os << "    - coplanar set:  " << cs.size() << " points.";
425            os << furthest.print(pr.run_id, "  Furthest");
426        }
427        zinc_(Zdistio);
428        double d= facet.distance(furthest);
429        os << "      furthest distance= " << d << endl; //FIXUP QH11010 %2.2g
430    }
431    QhullVertexSet vs= facet.vertices();
432    if(!vs.isEmpty()){
433        os << vs.print(pr.run_id, "    - vertices:");
434    }
435    QhullFacetSet fs= facet.neighborFacets();
436    fs.selectAll();
437    if(!fs.isEmpty()){
438        os << fs.printIdentifiers("    - neighboring facets:");
439    }
440    return os;
441}//operator<< PrintHeader
442
443
444//! Print ridges of facet to stream.  Same as qh_printfacetridges [io.c]
445//! If qhRunId==UsingLibQhull::NOqhRunId, does not use qh
446ostream &
447operator<<(ostream &os, const QhullFacet::PrintRidges &pr)
448{
449    const QhullFacet facet= *pr.facet;
450    facetT *f= facet.getFacetT();
451    QhullRidgeSet rs= facet.ridges();
452    if(!rs.isEmpty()){
453        if(pr.run_id!=UsingLibQhull::NOqhRunId){
454            UsingLibQhull q(pr.run_id);
455            // No calls to libqhull
456            if(f->visible && qh NEWfacets){
457                os << "    - ridges(ids may be garbage):";
458                for(QhullRidgeSet::iterator i=rs.begin(); i!=rs.end(); ++i){
459                    QhullRidge r= *i;
460                    os << " r" << r.id();
461                }
462                os << endl;
463            }else{
464                os << "    - ridges:" << endl;
465            }
466        }else{
467            os << "    - ridges:" << endl;
468        }
469
470        // Keep track of printed ridges
471        for(QhullRidgeSet::iterator i=rs.begin(); i!=rs.end(); ++i){
472            QhullRidge r= *i;
473            r.getRidgeT()->seen= false;
474        }
475        int ridgeCount= 0;
476        if(facet.dimension()==3){
477            for(QhullRidge r= rs.first(); !r.getRidgeT()->seen; r= r.nextRidge3d(facet)){
478                r.getRidgeT()->seen= true;
479                os << r.print(pr.run_id);
480                ++ridgeCount;
481                if(!r.hasNextRidge3d(facet)){
482                    break;
483                }
484            }
485        }else {
486            QhullFacetSet ns(facet.neighborFacets());
487            for(QhullFacetSet::iterator i=ns.begin(); i!=ns.end(); ++i){
488                QhullFacet neighbor= *i;
489                QhullRidgeSet nrs(neighbor.ridges());
490                for(QhullRidgeSet::iterator j=nrs.begin(); j!=nrs.end(); ++j){
491                    QhullRidge r= *j;
492                    if(r.otherFacet(neighbor)==facet){
493                        r.getRidgeT()->seen= true;
494                        os << r.print(pr.run_id);
495                        ridgeCount++;
496                    }
497                }
498            }
499        }
500        if(ridgeCount!=rs.count()){
501            os << "     - all ridges:";
502            for(QhullRidgeSet::iterator i=rs.begin(); i!=rs.end(); ++i){
503                QhullRidge r= *i;
504                os << " r" << r.id();
505            }
506            os << endl;
507        }
508        for(QhullRidgeSet::iterator i=rs.begin(); i!=rs.end(); ++i){
509            QhullRidge r= *i;
510            if(!r.getRidgeT()->seen){
511                os << r.print(pr.run_id);
512            }
513        }
514    }
515    return os;
516}//operator<< PrintRidges
517
518// "No conversion" error if defined inline
519ostream &
520operator<<(ostream &os, QhullFacet &f)
521{
522    os << f.print(UsingLibQhull::NOqhRunId);
523    return os;
524}//<< QhullFacet
Note: See TracBrowser for help on using the repository browser.