Logo ROOT   master
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 /*
23 Implementation of "marching cubes" algortihm for GL module. Used by
24 TGLTF3Painter and TGLIsoPainter.
25 Good and clear algorithm explanation can be found here:
26 http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
27 */
28 
29 class TF3;
30 class TKDEFGT;
31 
32 namespace Rgl {
33 namespace Mc {
34 
35 /*
36 Some routines, values and tables for marching cube method.
37 */
38 extern const UInt_t eInt[256];
39 extern const Float_t vOff[8][3];
40 extern const UChar_t eConn[12][2];
41 extern const Float_t eDir[12][3];
42 extern const Int_t conTbl[256][16];
43 
44 /*
45 "T" prefix in class names is only for code-style checker.
46 */
47 
48 /*
49 TCell is a cube from marching cubes algorithm.
50 It has "type" - defines, which vertices
51 are under iso level, which are above.
52 
53 Vertices 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 
72 TCell's "type" is 8-bit number, one bit per vertex.
73 So, if vertex 1 and 2 are under iso-surface, type
74 will 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 
79 type == 6.
80 
81 Edges 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 
100 There are 12 edges, any of them can be intersected by iso-surface
101 (not all 12 simultaneously). Edge's intersection is a vertex in
102 iso-mesh's vertices array, cell holds index of this vertex in
103 fIds array.
104 fVals holds "scalar field" or density values in vertices [0, 7].
105 
106 "V" parameter is the type to hold such values.
107 */
108 
109 template<class V>
110 class TCell {
111 public:
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 /*
126 TSlice of marching cubes' grid. Has W * H cells.
127 If you have TH3 hist, GetNbinsX() is W and GetNbinsY() is H.
128 */
129 template<class V>
130 class TSlice {
131 public:
133  {
134  }
135 
137  {
138  fCells.resize(w * h);
139  }
140 
141  std::vector<TCell<V> > fCells;
142 private:
143  TSlice(const TSlice &rhs);
144  TSlice & operator = (const TSlice &rhs);
145 };
146 
147 /*
148 Mesh builder requires generic "data source": it can
149 be a wrapped TH3 object, a wrapped TF3 object or some
150 "density estimator" object.
151 Mesh builder inherits this data source type.
152 
153 TH3Adapter is one of such data sources.
154 It has _direct_ access to TH3 internal data.
155 GetBinContent(i, j, k) is a virtual function
156 and it calls two other virtual functions - this
157 is very expensive if you call GetBinContent
158 several 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 
163 For example, H == TH3C, E == Char_t.
164 */
165 
166 template<class H, class E>
167 class TH3Adapter {
168 protected:
169 
170  typedef E ElementType_t;
171 
173  : fSrc(0), fW(0), fH(0), fD(0), fSliceSize(0)
174  {
175  }
176 
177  UInt_t GetW()const
178  {
179  return fW - 2;
180  }
181 
182  UInt_t GetH()const
183  {
184  return fH - 2;
185  }
186 
187  UInt_t GetD()const
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 /*
219 TF3Adapter. Lets TMeshBuilder to use TF3 as a 3d array.
220 TF3Adapter, TF3EdgeSplitter (see below) and TMeshBuilder<TF3, ValueType>
221 need TGridGeometry<ValueType>, so TGridGeometry is a virtual base.
222 */
223 
224 class TF3Adapter : protected virtual TGridGeometry<Double_t> {
225 protected:
227 
228  TF3Adapter() : fTF3(0), fW(0), fH(0), fD(0)
229  {
230  }
231 
232  UInt_t GetW()const
233  {
234  return fW;
235  }
236 
237  UInt_t GetH()const
238  {
239  return fH;
240  }
241 
242  UInt_t GetD()const
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 /*
261 TSourceAdapterSelector is aux. class used by TMeshBuilder to
262 select "data-source" base depending on data-source type.
263 */
264 template<class> class TSourceAdapterSelector;
265 
266 template<>
268 public:
270 };
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 /*
309 Edge splitter is the second base class for TMeshBuilder.
310 Its task is to split cell's edge by adding new vertex
311 into mesh.
312 Default splitter is used by TH3 and KDE.
313 */
314 
315 template<class E, class V>
316 V 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 
324 template<class H, class E, class V>
325 class TDefaultSplitter : protected virtual TGridGeometry<V> {
326 protected:
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 /*
347 TF3's edge splitter. Calculates new vertex and surface normal
348 in this vertex using TF3.
349 */
350 
351 class TF3EdgeSplitter : protected virtual TGridGeometry<Double_t> {
352 protected:
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 /*
369 TSplitterSelector is aux. class to select "edge-splitter" base
370 for TMeshBuilder.
371 */
372 
373 template<class, class> class TSplitterSelector;
374 
375 template<class V>
377 public:
379 };
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 /*
418 Mesh builder. Polygonizes scalar field - TH3, TF3 or
419 something else (some density estimator as data-source).
420 
421 ValueType is Float_t or Double_t - the type of vertex'
422 x,y,z components.
423 */
424 
425 template<class DataSource, class ValueType>
426 class TMeshBuilder : public TSourceAdapterSelector<DataSource>::Type_t,
427  public TSplitterSelector<DataSource, ValueType>::Type_t
428 {
429 private:
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 
447 public:
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 
456 private:
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 
482  TMeshBuilder(const TMeshBuilder &rhs);
483  TMeshBuilder & operator = (const TMeshBuilder &rhs);
484 };
485 
486 }//namespace Mc
487 }//namespace Rgl
488 
489 #endif
TH3Adapter< TH3D, Double_t > Type_t
TSlice< ElementType_t > SliceType_t
void ResizeSlice(UInt_t w, UInt_t h)
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:267
void SplitEdge(TCell< E > &cell, TIsoMesh< V > *mesh, UInt_t i, V x, V y, V z, V iso) const
TH3Adapter< TH3S, Short_t > Type_t
void SetDataSource(const TF3 *f)
Double_t GetData(UInt_t i, UInt_t j, UInt_t k) const
#define H(x, y, z)
void SetNormalEvaluator(const H *)
TDefaultSplitter< TH3F, Float_t, V > Type_t
#define f(i)
Definition: RSha256.hxx:104
int Int_t
Definition: RtypesCore.h:41
TMeshBuilder(Bool_t averagedNormals, ValueType eps=1e-7)
const Float_t eDir[12][3]
TDefaultSplitter< TH3D, Double_t, V > Type_t
V GetOffset(E val1, E val2, V iso)
void FetchDensities() const
UInt_t AddVertex(const V *v)
Definition: TGLIsoMesh.h:40
Double_t x[n]
Definition: legend1.C:17
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.
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 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
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:304
A 3-Dim function with parameters.
Definition: TF3.h:28
void BuildRow(SliceType_t *slice) const
The first row (along x) in the first slice: ny == 0, nz == 0, nx : [1, W - 1].
unsigned int UInt_t
Definition: RtypesCore.h:42
const UInt_t eInt[256]
TSplitterSelector< DataSource, ValueType >::Type_t SplitterBase_t
void SetDataSource(const H *hist)
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
#define h(i)
Definition: RSha256.hxx:106
TSlice & operator=(const TSlice &rhs)
void BuildCol(SliceType_t *slice) const
"Col" (column) consists of cubes along y axis on the first slice (nx == 0, nz == 0).
double Double_t
Definition: RtypesCore.h:55
std::vector< TCell< V > > fCells
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
const Float_t vOff[8][3]
3-D histogram with a byte per channel (see TH1 documentation)
Definition: TH3.h:158
void NextStep(UInt_t depth, const SliceType_t *prevSlice, SliceType_t *curr) const
Fill slice with vertices and triangles.
TCell< ElementType_t > CellType_t
const UChar_t eConn[12][2]
void BuildSlice(SliceType_t *slice) const
Slice with nz == 0.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
void BuildFirstCube(SliceType_t *slice) const
The first cube in a grid. nx == 0, ny == 0, nz ==0.
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
DataSourceBase_t::ElementType_t ElementType_t
ElementType_t GetData(UInt_t i, UInt_t j, UInt_t k) const
void FetchDensities() const
void BuildNormals() const
Build averaged normals using vertices and trinagles.