Logo ROOT   6.16/01
Reference Guide
TGLMarchingCubes.h
Go to the documentation of this file.
1// @(#)root/gl:$Id$
2// Author: Timur Pocheptsov 06/01/2009
3
4/*************************************************************************
5 * Copyright (C) 1995-2009, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12#ifndef ROOT_TGLMarchingCubes
13#define ROOT_TGLMarchingCubes
14
15#include <vector>
16
17#include "TH3.h"
18
19#include "TGLIsoMesh.h"
20#include "TKDEAdapter.h"
21
22/*
23Implementation of "marching cubes" algortihm for GL module. Used by
24TGLTF3Painter and TGLIsoPainter.
25Good and clear algorithm explanation can be found here:
26http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
27*/
28
29class TF3;
30class TKDEFGT;
31
32namespace Rgl {
33namespace Mc {
34
35/*
36Some routines, values and tables for marching cube method.
37*/
38extern const UInt_t eInt[256];
39extern const Float_t vOff[8][3];
40extern const UChar_t eConn[12][2];
41extern const Float_t eDir[12][3];
42extern const Int_t conTbl[256][16];
43
44/*
45"T" prefix in class names is only for code-style checker.
46*/
47
48/*
49TCell is a cube from marching cubes algorithm.
50It has "type" - defines, which vertices
51are under iso level, which are above.
52
53Vertices numeration:
54
55 |z
56 |
57 4____________7
58 /| /|
59 / | / |
60 / | / |
61 / | / |
62 5____|_______6 |
63 | 0_______|____3______ y
64 | / | /
65 | / | /
66 | / | /
67 |/ |/
68 1____________2
69 /
70 /x
71
72TCell's "type" is 8-bit number, one bit per vertex.
73So, if vertex 1 and 2 are under iso-surface, type
74will be:
75
76 7 6 5 4 3 2 1 0 (bit number)
77[0 0 0 0 0 1 1 0] bit pattern
78
79type == 6.
80
81Edges numeration:
82
83 |z
84 |
85 |_____7______
86 /| /|
87 / | / |
88 4/ 8 6 11
89 / | / |
90 /____|5______/ |
91 | |_____3_|____|______ y
92 | / | /
93 9 / 10 /
94 | /0 | /2
95 |/ |/
96 /____________/
97 / 1
98 /x
99
100There are 12 edges, any of them can be intersected by iso-surface
101(not all 12 simultaneously). Edge's intersection is a vertex in
102iso-mesh's vertices array, cell holds index of this vertex in
103fIds array.
104fVals holds "scalar field" or density values in vertices [0, 7].
105
106"V" parameter is the type to hold such values.
107*/
108
109template<class V>
110class TCell {
111public:
112 TCell() : fType(), fIds(), fVals()
113 {
114 //TCell ctor.
115 //Such mem-initializer list can produce
116 //warnings with some versions of MSVC,
117 //but this list is what I want.
118 }
119
122 V fVals[8];
123};
124
125/*
126TSlice of marching cubes' grid. Has W * H cells.
127If you have TH3 hist, GetNbinsX() is W and GetNbinsY() is H.
128*/
129template<class V>
130class TSlice {
131public:
133 {
134 }
135
137 {
138 fCells.resize(w * h);
139 }
140
141 std::vector<TCell<V> > fCells;
142private:
143 TSlice(const TSlice &rhs);
144 TSlice & operator = (const TSlice &rhs);
145};
146
147/*
148Mesh builder requires generic "data source": it can
149be a wrapped TH3 object, a wrapped TF3 object or some
150"density estimator" object.
151Mesh builder inherits this data source type.
152
153TH3Adapter is one of such data sources.
154It has _direct_ access to TH3 internal data.
155GetBinContent(i, j, k) is a virtual function
156and it calls two other virtual functions - this
157is very expensive if you call GetBinContent
158several million times as I do in marching cubes.
159
160"H" parameter is one of TH3 classes,
161"E" is the type of internal data.
162
163For example, H == TH3C, E == Char_t.
164*/
165
166template<class H, class E>
168protected:
169
171
173 : fSrc(0), fW(0), fH(0), fD(0), fSliceSize(0)
174 {
175 }
176
178 {
179 return fW - 2;
180 }
181
183 {
184 return fH - 2;
185 }
186
188 {
189 return fD - 2;
190 }
191
192 void SetDataSource(const H *hist)
193 {
194 fSrc = hist->GetArray();
195 fW = hist->GetNbinsX() + 2;
196 fH = hist->GetNbinsY() + 2;
197 fD = hist->GetNbinsZ() + 2;
198 fSliceSize = fW * fH;
199 }
200
201 void FetchDensities()const{}//Do nothing.
202
204 {
205 i += 1;
206 j += 1;
207 k += 1;
208 return fSrc[k * fSliceSize + j * fW + i];
209 }
210
216};
217
218/*
219TF3Adapter. Lets TMeshBuilder to use TF3 as a 3d array.
220TF3Adapter, TF3EdgeSplitter (see below) and TMeshBuilder<TF3, ValueType>
221need TGridGeometry<ValueType>, so TGridGeometry is a virtual base.
222*/
223
224class TF3Adapter : protected virtual TGridGeometry<Double_t> {
225protected:
227
228 TF3Adapter() : fTF3(0), fW(0), fH(0), fD(0)
229 {
230 }
231
233 {
234 return fW;
235 }
236
238 {
239 return fH;
240 }
241
243 {
244 return fD;
245 }
246
247 void SetDataSource(const TF3 *f);
248
249 void FetchDensities()const{}//Do nothing.
250
251 Double_t GetData(UInt_t i, UInt_t j, UInt_t k)const;
252
253 const TF3 *fTF3;//TF3 data source.
254 //TF3 grid's dimensions.
258};
259
260/*
261TSourceAdapterSelector is aux. class used by TMeshBuilder to
262select "data-source" base depending on data-source type.
263*/
264template<class> class TSourceAdapterSelector;
265
266template<>
268public:
270};
271
272template<>
274public:
276};
277
278template<>
280public:
282};
283
284template<>
286public:
288};
289
290template<>
292public:
294};
295
296template<>
298public:
300};
301
302template<>
304public:
306};
307
308/*
309Edge splitter is the second base class for TMeshBuilder.
310Its task is to split cell's edge by adding new vertex
311into mesh.
312Default splitter is used by TH3 and KDE.
313*/
314
315template<class E, class V>
316V GetOffset(E val1, E val2, V iso)
317{
318 const V delta = val2 - val1;
319 if (!delta)
320 return 0.5f;
321 return (iso - val1) / delta;
322}
323
324template<class H, class E, class V>
325class TDefaultSplitter : protected virtual TGridGeometry<V> {
326protected:
327 void SetNormalEvaluator(const H * /*source*/)
328 {
329 }
330 void SplitEdge(TCell<E> & cell, TIsoMesh<V> * mesh, UInt_t i,
331 V x, V y, V z, V iso)const
332{
333 V v[3];
334 const V offset = GetOffset(cell.fVals[eConn[i][0]],
335 cell.fVals[eConn[i][1]],
336 iso);
337 v[0] = x + (vOff[eConn[i][0]][0] + offset * eDir[i][0]) * this->fStepX;
338 v[1] = y + (vOff[eConn[i][0]][1] + offset * eDir[i][1]) * this->fStepY;
339 v[2] = z + (vOff[eConn[i][0]][2] + offset * eDir[i][2]) * this->fStepZ;
340 cell.fIds[i] = mesh->AddVertex(v);
341}
342
343
344};
345
346/*
347TF3's edge splitter. Calculates new vertex and surface normal
348in this vertex using TF3.
349*/
350
351class TF3EdgeSplitter : protected virtual TGridGeometry<Double_t> {
352protected:
354 {
355 }
356
357 void SetNormalEvaluator(const TF3 *tf3)
358 {
359 fTF3 = tf3;
360 }
361
362 void SplitEdge(TCell<Double_t> & cell, TIsoMesh<Double_t> * mesh, UInt_t i,
363 Double_t x, Double_t y, Double_t z, Double_t iso)const;
364
365 const TF3 *fTF3;
366};
367
368/*
369TSplitterSelector is aux. class to select "edge-splitter" base
370for TMeshBuilder.
371*/
372
373template<class, class> class TSplitterSelector;
374
375template<class V>
377public:
379};
380
381template<class V>
383public:
385};
386
387template<class V>
389public:
391};
392
393template<class V>
395public:
397};
398
399template<class V>
401public:
403};
404
405template<class V>
407public:
409};
410
411template<class V>
413public:
415};
416
417/*
418Mesh builder. Polygonizes scalar field - TH3, TF3 or
419something else (some density estimator as data-source).
420
421ValueType is Float_t or Double_t - the type of vertex'
422x,y,z components.
423*/
424
425template<class DataSource, class ValueType>
426class TMeshBuilder : public TSourceAdapterSelector<DataSource>::Type_t,
427 public TSplitterSelector<DataSource, ValueType>::Type_t
428{
429private:
430 //Two base classes.
433 //Using declarations required, since these are
434 //type-dependant names in template.
435 using DataSourceBase_t::GetW;
436 using DataSourceBase_t::GetH;
437 using DataSourceBase_t::GetD;
438 using DataSourceBase_t::GetData;
439 using SplitterBase_t::SplitEdge;
440
441 typedef typename DataSourceBase_t::ElementType_t ElementType_t;
442
446
447public:
448 TMeshBuilder(Bool_t averagedNormals, ValueType eps = 1e-7)
449 : fAvgNormals(averagedNormals), fMesh(0), fIso(), fEpsilon(eps)
450 {
451 }
452
453 void BuildMesh(const DataSource *src, const TGridGeometry<ValueType> &geom,
454 MeshType_t *mesh, ValueType iso);
455
456private:
457
461 ValueType fIso;
462 ValueType fEpsilon;
463
464 void NextStep(UInt_t depth, const SliceType_t *prevSlice,
465 SliceType_t *curr)const;
466
467 void BuildFirstCube(SliceType_t *slice)const;
468 void BuildRow(SliceType_t *slice)const;
469 void BuildCol(SliceType_t *slice)const;
470 void BuildSlice(SliceType_t *slice)const;
471 void BuildFirstCube(UInt_t depth, const SliceType_t *prevSlice,
472 SliceType_t *slice)const;
473 void BuildRow(UInt_t depth, const SliceType_t *prevSlice,
474 SliceType_t *slice)const;
475 void BuildCol(UInt_t depth, const SliceType_t *prevSlice,
476 SliceType_t *slice)const;
477 void BuildSlice(UInt_t depth, const SliceType_t *prevSlice,
478 SliceType_t *slice)const;
479
480 void BuildNormals()const;
481
484};
485
486}//namespace Mc
487}//namespace Rgl
488
489#endif
SVector< double, 2 > v
Definition: Dict.h:5
#define f(i)
Definition: RSha256.hxx:104
#define h(i)
Definition: RSha256.hxx:106
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:41
unsigned char UChar_t
Definition: RtypesCore.h:34
unsigned int UInt_t
Definition: RtypesCore.h:42
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
float Float_t
Definition: RtypesCore.h:53
void SetNormalEvaluator(const H *)
void SplitEdge(TCell< E > &cell, TIsoMesh< V > *mesh, UInt_t i, V x, V y, V z, V iso) const
Double_t GetData(UInt_t i, UInt_t j, UInt_t k) const
void SetDataSource(const TF3 *f)
void FetchDensities() const
void SetNormalEvaluator(const TF3 *tf3)
void SplitEdge(TCell< Double_t > &cell, TIsoMesh< Double_t > *mesh, UInt_t i, Double_t x, Double_t y, Double_t z, Double_t iso) const
Split the edge and find normal in a new vertex.
void SetDataSource(const H *hist)
const ElementType_t * fSrc
void FetchDensities() const
ElementType_t GetData(UInt_t i, UInt_t j, UInt_t k) const
UInt_t AddVertex(const V *v)
Definition: TGLIsoMesh.h:40
void BuildMesh(const DataSource *src, const TGridGeometry< ValueType > &geom, MeshType_t *mesh, ValueType iso)
Build iso-mesh using marching cubes.
TMeshBuilder & operator=(const TMeshBuilder &rhs)
void NextStep(UInt_t depth, const SliceType_t *prevSlice, SliceType_t *curr) const
Fill slice with vertices and triangles.
TMeshBuilder(const TMeshBuilder &rhs)
TCell< ElementType_t > CellType_t
TMeshBuilder(Bool_t averagedNormals, ValueType eps=1e-7)
DataSourceBase_t::ElementType_t ElementType_t
void BuildRow(SliceType_t *slice) const
The first row (along x) in the first slice: ny == 0, nz == 0, nx : [1, W - 1].
void BuildCol(SliceType_t *slice) const
"Col" (column) consists of cubes along y axis on the first slice (nx == 0, nz == 0).
TSplitterSelector< DataSource, ValueType >::Type_t SplitterBase_t
void BuildSlice(SliceType_t *slice) const
Slice with nz == 0.
TSlice< ElementType_t > SliceType_t
TIsoMesh< ValueType > MeshType_t
TSourceAdapterSelector< DataSource >::Type_t DataSourceBase_t
void BuildNormals() const
Build averaged normals using vertices and trinagles.
void BuildFirstCube(SliceType_t *slice) const
The first cube in a grid. nx == 0, ny == 0, nz ==0.
TSlice(const TSlice &rhs)
TSlice & operator=(const TSlice &rhs)
std::vector< TCell< V > > fCells
void ResizeSlice(UInt_t w, UInt_t h)
TH3Adapter< TH3D, Double_t > Type_t
TH3Adapter< TH3F, Float_t > Type_t
TH3Adapter< TH3S, Short_t > Type_t
TDefaultSplitter< TH3C, Char_t, V > Type_t
TDefaultSplitter< TH3D, Double_t, V > Type_t
TDefaultSplitter< TH3F, Float_t, V > Type_t
TDefaultSplitter< TH3I, Int_t, V > Type_t
TDefaultSplitter< TH3S, Short_t, V > Type_t
TDefaultSplitter< TKDEFGT, Float_t, Float_t > Type_t
A 3-Dim function with parameters.
Definition: TF3.h:28
3-D histogram with a byte per channel (see TH1 documentation)
Definition: TH3.h:158
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:304
3-D histogram with a float per channel (see TH1 documentation)}
Definition: TH3.h:267
3-D histogram with an int per channel (see TH1 documentation)}
Definition: TH3.h:230
3-D histogram with a short per channel (see TH1 documentation)
Definition: TH3.h:194
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
#define H(x, y, z)
const UChar_t eConn[12][2]
const Int_t conTbl[256][16]
const Float_t vOff[8][3]
V GetOffset(E val1, E val2, V iso)
const Float_t eDir[12][3]
const UInt_t eInt[256]
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:97