include/root/TGLMarchingCubes.h

0001 // @(#)root/gl:$Id$
0002 // Author:  Timur Pocheptsov  06/01/2009
0003
0004 /*************************************************************************
0005  * Copyright (C) 1995-2009, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011
0012 #ifndef ROOT_TGLMarchingCubes
0013 #define ROOT_TGLMarchingCubes
0014
0015 #include <vector>
0016
0017 #include "TH3.h"
0018
0019 #include "TGLIsoMesh.h"
0020 #include "TKDEAdapter.h"
0021
0022 /*
0023 Implementation of "marching cubes" algortihm for GL module. Used by
0024 TGLTF3Painter and TGLIsoPainter.
0025 Good and clear algorithm explanation can be found here:
0026 http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
0027 */
0028
0029 class TF3;
0030 class TKDEFGT;
0031
0032 namespace Rgl {
0033 namespace Mc {
0034
0035 /*
0036 Some routines, values and tables for marching cube method.
0037 */
0038 extern const UInt_t  eInt[256];
0039 extern const Float_t vOff[8][3];
0040 extern const UChar_t eConn[12][2];
0041 extern const Float_t eDir[12][3];
0042 extern const Int_t   conTbl[256][16];
0043
0044 /*
0045 "T" prefix in class names is only for code-style checker.
0046 */
0047
0048 /*
0049 TCell is a cube from marching cubes algorithm.
0050 It has "type" - defines, which vertices
0051 are under iso level, which are above.
0052
0053 Vertices numeration:
0054
0055            |z
0056            |
0057            4____________7
0058           /|           /|
0059          / |          / |
0060         /  |         /  |
0061        /   |        /   |
0062       5____|_______6    |
0063       |    0_______|____3______ y
0064       |   /        |   /
0065       |  /         |  /
0066       | /          | /
0067       |/           |/
0068       1____________2
0069      /
0070     /x
0071
0072 TCell's "type" is 8-bit number, one bit per vertex.
0073 So, if vertex 1 and 2 are under iso-surface, type
0074 will be:
0075
0076  7 6 5 4 3 2 1 0 (bit number)
0077 [0 0 0 0 0 1 1 0] bit pattern
0078
0079 type == 6.
0080
0081 Edges numeration:
0082
0083            |z
0084            |
0085            |_____7______
0086           /|           /|
0087          / |          / |
0088        4/  8         6  11
0089        /   |        /   |
0090       /____|5______/    |
0091       |    |_____3_|____|______ y
0092       |   /        |   /
0093       9  /        10  /
0094       | /0         | /2
0095       |/           |/
0096       /____________/
0097      /      1
0098     /x
0099
0100 There are 12 edges, any of them can be intersected by iso-surface
0101 (not all 12 simultaneously). Edge's intersection is a vertex in
0102 iso-mesh's vertices array, cell holds index of this vertex in
0103 fIds array.
0104 fVals holds "scalar field" or density values in vertices [0, 7].
0105
0106 "V" parameter is the type to hold such values.
0107 */
0108
0109 template<class V>
0110 class TCell {
0111 public:
0112    TCell() : fType(), fIds(), fVals()
0113    {
0114       //TCell ctor.
0115       //Such mem-initializer list can produce
0116       //warnings with some versions of MSVC,
0117       //but this list is what I want.
0118    }
0119
0120    UInt_t     fType;
0121    UInt_t     fIds[12];
0122    V          fVals[8];
0123 };
0124
0125 /*
0126 TSlice of marching cubes' grid. Has W * H cells.
0127 If you have TH3 hist, GetNbinsX() is W and GetNbinsY() is H.
0128 */
0129 template<class V>
0130 class TSlice {
0131 public:
0132    TSlice()
0133    {
0134    }
0135
0136    void ResizeSlice(UInt_t w, UInt_t h)
0137    {
0138       fCells.resize(w * h);
0139    }
0140
0141    std::vector<TCell<V> > fCells;
0142 private:
0143    TSlice(const TSlice &rhs) = delete;
0144    TSlice & operator = (const TSlice &rhs) = delete;
0145 };
0146
0147 /*
0148 Mesh builder requires generic "data source": it can
0149 be a wrapped TH3 object, a wrapped TF3 object or some
0150 "density estimator" object.
0151 Mesh builder inherits this data source type.
0152
0153 TH3Adapter is one of such data sources.
0154 It has _direct_ access to TH3 internal data.
0155 GetBinContent(i, j, k) is a virtual function
0156 and it calls two other virtual functions - this
0157 is very expensive if you call GetBinContent
0158 several million times as I do in marching cubes.
0159
0160 "H" parameter is one of TH3 classes,
0161 "E" is the type of internal data.
0162
0163 For example, H == TH3C, E == Char_t.
0164 */
0165
0166 template<class H, class E>
0167 class TH3Adapter {
0168 protected:
0169
0170    typedef E    ElementType_t;
0171
0172    TH3Adapter()
0173       : fSrc(nullptr), fW(0), fH(0), fD(0), fSliceSize(0)
0174    {
0175    }
0176
0177    UInt_t GetW()const
0178    {
0179       return fW - 2;
0180    }
0181
0182    UInt_t GetH()const
0183    {
0184       return fH - 2;
0185    }
0186
0187    UInt_t GetD()const
0188    {
0189       return fD - 2;
0190    }
0191
0192    void SetDataSource(const H *hist)
0193    {
0194       fSrc = hist->GetArray();
0195       fW   = hist->GetNbinsX() + 2;
0196       fH   = hist->GetNbinsY() + 2;
0197       fD   = hist->GetNbinsZ() + 2;
0198       fSliceSize = fW * fH;
0199    }
0200
0201    void FetchDensities()const{}//Do nothing.
0202
0203    ElementType_t GetData(UInt_t i, UInt_t j, UInt_t k)const
0204    {
0205       i += 1;
0206       j += 1;
0207       k += 1;
0208       return fSrc[k * fSliceSize + j * fW + i];
0209    }
0210
0211    const ElementType_t *fSrc;
0212    UInt_t fW;
0213    UInt_t fH;
0214    UInt_t fD;
0215    UInt_t fSliceSize;
0216 };
0217
0218 /*
0219 TF3Adapter. Lets TMeshBuilder to use TF3 as a 3d array.
0220 TF3Adapter, TF3EdgeSplitter (see below) and TMeshBuilder<TF3, ValueType>
0221 need TGridGeometry<ValueType>, so TGridGeometry is a virtual base.
0222 */
0223
0224 class TF3Adapter : protected virtual TGridGeometry<Double_t> {
0225 protected:
0226    typedef Double_t ElementType_t;
0227
0228    TF3Adapter() : fTF3(nullptr), fW(0), fH(0), fD(0)
0229    {
0230    }
0231
0232    UInt_t GetW()const
0233    {
0234       return fW;
0235    }
0236
0237    UInt_t GetH()const
0238    {
0239       return fH;
0240    }
0241
0242    UInt_t GetD()const
0243    {
0244       return fD;
0245    }
0246
0247    void SetDataSource(const TF3 *f);
0248
0249    void FetchDensities()const{}//Do nothing.
0250
0251    Double_t GetData(UInt_t i, UInt_t j, UInt_t k)const;
0252
0253    const TF3 *fTF3;//TF3 data source.
0254    //TF3 grid's dimensions.
0255    UInt_t     fW;
0256    UInt_t     fH;
0257    UInt_t     fD;
0258 };
0259
0260 /*
0261 TSourceAdapterSelector is aux. class used by TMeshBuilder to
0262 select "data-source" base depending on data-source type.
0263 */
0264 template<class> class TSourceAdapterSelector;
0265
0266 template<>
0267 class TSourceAdapterSelector<TH3C> {
0268 public:
0269    typedef TH3Adapter<TH3C, Char_t> Type_t;
0270 };
0271
0272 template<>
0273 class TSourceAdapterSelector<TH3S> {
0274 public:
0275    typedef TH3Adapter<TH3S, Short_t> Type_t;
0276 };
0277
0278 template<>
0279 class TSourceAdapterSelector<TH3I> {
0280 public:
0281    typedef TH3Adapter<TH3I, Int_t> Type_t;
0282 };
0283
0284 template<>
0285 class TSourceAdapterSelector<TH3F> {
0286 public:
0287    typedef TH3Adapter<TH3F, Float_t> Type_t;
0288 };
0289
0290 template<>
0291 class TSourceAdapterSelector<TH3D> {
0292 public:
0293    typedef TH3Adapter<TH3D, Double_t> Type_t;
0294 };
0295
0296 template<>
0297 class TSourceAdapterSelector<TF3> {
0298 public:
0299    typedef TF3Adapter Type_t;
0300 };
0301
0302 template<>
0303 class TSourceAdapterSelector<TKDEFGT> {
0304 public:
0305    typedef Fgt::TKDEAdapter Type_t;
0306 };
0307
0308 /*
0309 Edge splitter is the second base class for TMeshBuilder.
0310 Its task is to split cell's edge by adding new vertex
0311 into mesh.
0312 Default splitter is used by TH3 and KDE.
0313 */
0314
0315 template<class E, class V>
0316 V GetOffset(E val1, E val2, V iso)
0317 {
0318    const V delta = val2 - val1;
0319    if (!delta)
0320       return 0.5f;
0321    return (iso - val1) / delta;
0322 }
0323
0324 template<class H, class E, class V>
0325 class TDefaultSplitter : protected virtual TGridGeometry<V> {
0326 protected:
0327    void SetNormalEvaluator(const H * /*source*/)
0328    {
0329    }
0330    void SplitEdge(TCell<E> & cell, TIsoMesh<V> * mesh, UInt_t i,
0331                   V x, V y, V z, V iso)const
0332 {
0333    V v[3];
0334    const V offset = GetOffset(cell.fVals[eConn[i][0]],
0335                               cell.fVals[eConn[i][1]],
0336                               iso);
0337    v[0] = x + (vOff[eConn[i][0]][0] + offset * eDir[i][0]) * this->fStepX;
0338    v[1] = y + (vOff[eConn[i][0]][1] + offset * eDir[i][1]) * this->fStepY;
0339    v[2] = z + (vOff[eConn[i][0]][2] + offset * eDir[i][2]) * this->fStepZ;
0340    cell.fIds[i] = mesh->AddVertex(v);
0341 }
0342
0343
0344 };
0345
0346 /*
0347 TF3's edge splitter. Calculates new vertex and surface normal
0348 in this vertex using TF3.
0349 */
0350
0351 class TF3EdgeSplitter : protected virtual TGridGeometry<Double_t> {
0352 protected:
0353    TF3EdgeSplitter() : fTF3(nullptr)
0354    {
0355    }
0356
0357    void SetNormalEvaluator(const TF3 *tf3)
0358    {
0359       fTF3 = tf3;
0360    }
0361
0362    void SplitEdge(TCell<Double_t> & cell, TIsoMesh<Double_t> * mesh, UInt_t i,
0363                   Double_t x, Double_t y, Double_t z, Double_t iso)const;
0364
0365    const TF3 *fTF3;
0366 };
0367
0368 /*
0369 TSplitterSelector is aux. class to select "edge-splitter" base
0370 for TMeshBuilder.
0371 */
0372
0373 template<class, class> class TSplitterSelector;
0374
0375 template<class V>
0376 class TSplitterSelector<TH3C, V> {
0377 public:
0378    typedef TDefaultSplitter<TH3C, Char_t, V> Type_t;
0379 };
0380
0381 template<class V>
0382 class TSplitterSelector<TH3S, V> {
0383 public:
0384    typedef TDefaultSplitter<TH3S, Short_t, V> Type_t;
0385 };
0386
0387 template<class V>
0388 class TSplitterSelector<TH3I, V> {
0389 public:
0390    typedef TDefaultSplitter<TH3I, Int_t, V> Type_t;
0391 };
0392
0393 template<class V>
0394 class TSplitterSelector<TH3F, V> {
0395 public:
0396    typedef TDefaultSplitter<TH3F, Float_t, V> Type_t;
0397 };
0398
0399 template<class V>
0400 class TSplitterSelector<TH3D, V> {
0401 public:
0402    typedef TDefaultSplitter<TH3D, Double_t, V> Type_t;
0403 };
0404
0405 template<class V>
0406 class TSplitterSelector<TKDEFGT, V> {
0407 public:
0408    typedef TDefaultSplitter<TKDEFGT, Float_t, Float_t> Type_t;
0409 };
0410
0411 template<class V>
0412 class TSplitterSelector<TF3, V> {
0413 public:
0414    typedef TF3EdgeSplitter Type_t;
0415 };
0416
0417 /*
0418 Mesh builder. Polygonizes scalar field - TH3, TF3 or
0419 something else (some density estimator as data-source).
0420
0421 ValueType is Float_t or Double_t - the type of vertex'
0422 x,y,z components.
0423 */
0424
0425 template<class DataSource, class ValueType>
0426 class TMeshBuilder : public TSourceAdapterSelector<DataSource>::Type_t,
0427                      public TSplitterSelector<DataSource, ValueType>::Type_t
0428 {
0429 private:
0430    //Two base classes.
0431    typedef typename TSourceAdapterSelector<DataSource>::Type_t       DataSourceBase_t;
0432    typedef typename TSplitterSelector<DataSource, ValueType>::Type_t SplitterBase_t;
0433    //Using declarations required, since these are
0434    //type-dependant names in template.
0435    using DataSourceBase_t::GetW;
0436    using DataSourceBase_t::GetH;
0437    using DataSourceBase_t::GetD;
0438    using DataSourceBase_t::GetData;
0439    using SplitterBase_t::SplitEdge;
0440
0441    typedef typename DataSourceBase_t::ElementType_t ElementType_t;
0442
0443    typedef TCell<ElementType_t>  CellType_t;
0444    typedef TSlice<ElementType_t> SliceType_t;
0445    typedef TIsoMesh<ValueType>   MeshType_t;
0446
0447 public:
0448    TMeshBuilder(Bool_t averagedNormals, ValueType eps = 1e-7)
0449       : fAvgNormals(averagedNormals), fMesh(nullptr), fIso(), fEpsilon(eps)
0450    {
0451    }
0452
0453    void BuildMesh(const DataSource *src, const TGridGeometry<ValueType> &geom,
0454                   MeshType_t *mesh, ValueType iso);
0455
0456 private:
0457
0458    Bool_t      fAvgNormals;
0459    SliceType_t fSlices[2];
0460    MeshType_t *fMesh;
0461    ValueType   fIso;
0462    ValueType   fEpsilon;
0463
0464    void NextStep(UInt_t depth, const SliceType_t *prevSlice,
0465                  SliceType_t *curr)const;
0466
0467    void BuildFirstCube(SliceType_t *slice)const;
0468    void BuildRow(SliceType_t *slice)const;
0469    void BuildCol(SliceType_t *slice)const;
0470    void BuildSlice(SliceType_t *slice)const;
0471    void BuildFirstCube(UInt_t depth, const SliceType_t *prevSlice,
0472                        SliceType_t *slice)const;
0473    void BuildRow(UInt_t depth, const SliceType_t *prevSlice,
0474                  SliceType_t *slice)const;
0475    void BuildCol(UInt_t depth, const SliceType_t *prevSlice,
0476                  SliceType_t *slice)const;
0477    void BuildSlice(UInt_t depth, const SliceType_t *prevSlice,
0478                    SliceType_t *slice)const;
0479
0480    void BuildNormals()const;
0481
0482    TMeshBuilder(const TMeshBuilder &rhs);
0483    TMeshBuilder & operator = (const TMeshBuilder &rhs);
0484 };
0485
0486 }//namespace Mc
0487 }//namespace Rgl
0488
0489 #endif