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