versor  3.0
C++11 library for Geometric algebra
vsr_cubicLattice.h
1 //
2 // vsr_lattice.h
3 // Versor
4 //
5 
6 // Should not assume a metric . . .
7 //
8 // Created by Pablo Colapinto on 3/7/13.
9 // Copyright (c) 2013 __MyCompanyName__. All rights reserved.
10 //
11 
12 #ifndef Versor_vsr_cubiclattice_h_included
13 #define Versor_vsr_cubiclattice_h_included
14 
15 
16 #include <vector>
17 
18 #include "vsr.h"
19 
20 #include "form/vsr_interp.h"
21 #include "util/vsr_patch.h"
22 
23 using namespace std;
24 
25 namespace vsr {
26 
27  #define ITER \
28  for (int i = 0; i < mWidth; ++i){\
29  for (int j = 0; j < mHeight; ++j){\
30  for (int k = 0; k < mDepth; ++k){\
31  int tidx = idx(i,j,k);
32 
33  #define ITEND }}}
34 
35  #define ITERV \
36  for (int i = 0; i < mWidth; ++i){\
37  double u = 1.0 * i/mWidth;\
38  for (int j = 0; j < mHeight; ++j){\
39  double v = 1.0 * j/mHeight;\
40  for (int k = 0; k < mDepth; ++k){\
41  double w = 1.0 * k/mDepth; \
42  int tidx = idx(i,j,k);
43 
44  #define BOUNDITER0\
45  for (int i = 0; i < mWidth-1; ++i){ \
46  for (int j = 0; j < mHeight-1; ++j){ \
47  for (int k = 0; k < mDepth-1; ++k){
48 
49  #define BOUNDITER\
50  for (int i = 1; i < mWidth-1; ++i){ \
51  for (int j = 1; j < mHeight-1; ++j){ \
52  for (int k = 1; k < mDepth-1; ++k){
53 
54  #define BOUNDEND \
55  }}}
56 
57 
58 
60  template<class LPnt> // Template on Lattice Point Type
61  class CubicLattice {
62 
63  public:
64 
65 
66  CubicLattice(int _w = 1, int _h = 1, int _d = 1, double _sx = 1.0,double _sy = 1.0,double _sz = 1.0)
67  :mWidth(_w), mHeight(_h), mDepth(_d), mSpacingX(_sx),mSpacingY(_sy), mSpacingZ(_sz),
68  mNum( mWidth * mHeight * mDepth),
69  mNumVxl( (mWidth-1) * (mHeight-1) * (mDepth-1) ),
70  mPoint(NULL), mVxl(NULL), mNbr(NULL), mNbrVxl(NULL)
71  {
72  alloc(); init();
73  }
74 
75  ~CubicLattice() {
76  onDestroy();
77  }
78 
79  void onDestroy(){
80 
81  if (mPoint) delete[] mPoint;
82  if (mVxl) delete[] mVxl;
83  if (mNbr) delete[] mNbr;
84  if (mNbrVxl) delete[] mNbrVxl;
85 
86  if (!mFace.empty() ) mFace.clear();
87  if (!mEdge.empty() ) mEdge.clear();
88  if (!mCorner.empty()) mCorner.clear();
89  if (!mFaceVxl.empty()) mFaceVxl.clear();
90  if (!mEdgeVxl.empty()) mEdgeVxl.clear();
91  if (!mCornerVxl.empty()) mCornerVxl.clear();
92 
93  mWidth = mHeight = mDepth = mNum = mNumVxl = 0;
94  }
95 
96  CubicLattice& resize( int _w, int _h, int _d, double _sx = 1.0,double _sy=1.0,double _sz=1.0){
97  onDestroy();
98  mWidth = _w; mHeight = _h; mDepth = _d; mSpacingX = _sx; mSpacingY = _sy; mSpacingZ = _sz;
99 
100  mNum = _w * _h * _d;
101  mNumVxl = (_w-1) * (_h-1) * (_d -1);
102  alloc(); init();
103  return *this;
104  }
105 
106  int w() const { return mWidth; }
107  int h() const { return mHeight; }
108  int d() const { return mDepth; }
109 
110  int num() const { return mNum; }
111  int numVxl() const { return mNumVxl; }
113  int idx(int i, int j, int k=0) const {
114  if (i > w() - 1) i = w() -1;
115  if (j > h() - 1) j = h() -1;
116  if (k > d() - 1) k = d() -1;
117 
118  if (i < 0 ) i = 0;
119  if (j < 0 ) j = 0;
120  if (k < 0 ) k = 0;
121 
122  return i * mHeight * mDepth + j * mDepth + k;
123  }
124 
125  //VOXEL IDX
126  Vxl vxl(int ix) const { return mVxl[ix]; }
127  Vxl& vxl(int ix) { return mVxl[ix]; }
128 
129  Nbr nbr(int ix) const { return mNbr[ix]; }
130  Nbr& nbr(int ix) { return mNbr[ix]; }
131  Nbr nbrVxl(int ix) const { return mNbrVxl[ix]; }
132  Nbr& nbrVxl(int ix) { return mNbrVxl[ix]; }
133 
134  /* Totals and Offsets From Center */
136  double tw() const { return (mWidth-1) * mSpacingX; }
138  double ow() const { return tw() / 2.0 ; }
140  double th() const { return (mHeight-1) * mSpacingY; }
142  double oh() const { return th() / 2.0 ; }
144  double td() const { return (mDepth-1) * mSpacingZ; }
146  double od() const { return td() / 2.0 ; }//0;}//
148  double px(int i) const { return -ow() + (mSpacingX * i); }
150  double py(int j) const { return -oh() + (mSpacingY * j); }
152  double pz(int k) const { return od() - (mSpacingZ * k); }
153 
154  void alloc(){
155  mPoint = new LPnt[mNum];
156  mNbr = new Nbr[mNum];
157  mVxl = new Vxl[mNumVxl];
158  mNbrVxl = new Nbr[mNumVxl];
159  }
160 
161  //specialize this if desired
162  void initPoints(){
163  ITER
164  mPoint[ tidx ] = LPnt( px(i), py(j), pz(k) );
165  ITEND
166  }
167 
168  void init(){
169 
170  initPoints();
171 
172  ITER
173  int type = 0;
174  type |= (k==0) ? FRONT : 0;
175  type |= (k==mDepth-1) ? BACK: 0;
176  type |= (j==0) ? BOTTOM : 0;
177  type |= (j==mHeight-1) ? TOP: 0;
178  type |= (i==0) ? LEFT : 0;
179  type |= (i==mWidth-1) ? RIGHT: 0;
180 
181  mNbr[tidx] = Nbr(tidx, mWidth, mHeight, mDepth, type);
182 
183  ITEND
184 
185  //VXLS
186  BOUNDITER0
187 
188  int ix = i * (mHeight-1) * (mDepth-1) + j * (mDepth-1) + k;
189  NEVec<3> v( px(i) + mSpacingX/2.0, py(j) + mSpacingY/2.0, pz(k) - mSpacingZ/2.0);
190  mVxl[ix] = vxlAt(v);
191 
192  //assign face information
193  int type = 0;
194  type |= (k==0) ? FRONT : 0;
195  type |= (k==mDepth-2) ? BACK: 0;
196  type |= (j==0) ? BOTTOM : 0;
197  type |= (j==mHeight-2) ? TOP: 0;
198  type |= (i==0) ? LEFT : 0;
199  type |= (i==mWidth-2) ? RIGHT: 0;
200 
201  mVxl[ix].type = type;
202  mNbrVxl[ix] = Nbr(ix, mWidth, mHeight, mDepth, type);
203 
204  BOUNDEND
205 
206  //Assign Edges and Faces idx
207  for (int i = 0; i < mNum; ++i){
208  FE( mNbr[i] );
209  }
210  //Assign EdgeVxl and FaceVxl idx
211  for (int i = 0; i < mNumVxl; ++i){
212  vxlFE( mNbrVxl[i] );
213  }
214 
215  }
216 
218  template<class V>
219  Vxl vxlAt( const V& tv ) const {
220 
221  Vxl vxl;
222 
223  V v = bound(tv);
224 
225  //Bottom - Left corner Voxel
226  int lw = floor( ( v[0] + ow() ) / mSpacingX );
227  int lh = floor((v[1] + oh() ) / mSpacingY );
228  int ld = floor((-v[2] + od() ) / mSpacingZ );
229 
230  if (lw > w()-2) lw = w() -2;
231  if (lh > h()-2) lh = h() -2;
232  if (ld > d() -2) ld = d() -2;
233 
234  vxl.a = idx(lw,lh,ld);
235  vxl.b = idx(lw+1,lh,ld);
236  vxl.c = idx(lw+1,lh+1,ld);
237  vxl.d = idx(lw,lh+1,ld);
238 
239  vxl.e = idx(lw,lh,ld+1);
240  vxl.f = idx(lw+1,lh,ld+1);
241  vxl.g = idx(lw+1,lh+1,ld+1);
242  vxl.h = idx(lw,lh+1,ld+1);
243 
244  return vxl;
245  }
246 
247  template<class V>
248  V bound ( const V& p) const {
249  V v = p;
250 
251  if ( v[0] < px(0) ) v[0] = px(0);// + f.spacing()/2.0;
252  else if ( v[0] > px( w()-1) ) v[0] = px( w()-1);// - f.spacing()/2.0;
253 
254  if ( v[1] < py(0) ) v[1] = py(0);// + f.spacing()/2.0;
255  else if ( v[1] > py( h()-1) ) v[1] = py( h()-1);// - f.spacing()/2.0;
256 
257  if ( v[2] > pz(0) ) v[2] = pz(0);// - f.spacing()/2.0;
258  else if ( v[2] < pz( d()-1) ) v[2] = pz( d()-1);// + f.spacing()/2.0;
259 
260  return v;
261  }
262 
263  //bound and modify to range [0,1]
264  template<class V>
265  V range2D( const V& v) const{
266  V t = v;
267 
268  double minx = px(0);
269  double maxx = px(mWidth-1);
270  double miny = py(0);
271  double maxy = py(mHeight-1);
272 
273  if (t[0] < minx) t[0] = minx;
274  else if (t[0] > maxx) t[0] = maxx;
275  if (t[1] < miny) t[1] = miny;
276  else if (t[1] > maxy) t[1] = maxy;
277 
278  double dx = (t[0] - minx)/tw();
279  double dy = (t[1] - miny)/th();
280 
281  return V(dx,dy);
282  }
283 
284  //bound and modify to range [0,1]
285  template<class V>
286  V range( const V& v) const{
287  V t = v;
288 
289  double minx = px(0);
290  double maxx = px(mWidth-1);
291  double miny = py(0);
292  double maxy = py(mHeight-1);
293  double maxz = pz(0);
294  double minz = pz(mDepth-1);
295 
296  if (t[0] < minx) t[0] = minx;
297  else if (t[0] > maxx) t[0] = maxx;
298  if (t[1] < miny) t[1] = miny;
299  else if (t[1] > maxy) t[1] = maxy;
300  if (t[2] < minz) t[2] = minz;
301  else if (t[2] > maxz) t[2] = maxz;
302 
303  double dx = (t[0] - minx)/tw();
304  double dy = (t[1] - miny)/th();
305  double dz = -(t[2] - maxz)/td();
306 
307  return V(dx,dy,dz);
308  }
309 
310 
312  void FE( Nbr nb){
313 
314  //one undefined neighbor is a face
315  int i = 0;
316  while (i < 7){
317  if (nb[i] == -1) {
318  mFace.push_back( nb.idx );
319  break;
320  }
321  i++;
322  }
323 
324  //two undefined neighbors is an edge
325  i = 0; int n = 0;
326  while ( i < 7 ){
327  if (nb[i] == -1) {
328  n++;
329  if (n==2){
330  mEdge.push_back( nb.idx );
331  break;
332  }
333  }
334  i++;
335  }
336 
337  //two undefined neighbors is a corner
338  i = 0; n = 0;
339  while ( i < 7 ){
340  if (nb[i] == -1) {
341  n++;
342  if (n==3){
343  mCorner.push_back( nb.idx );
344  break;
345  }
346  }
347  i++;
348  }
349 
350  }
351 
352  //see above
353  void vxlFE( Nbr nb){
354  int i = 0;
355  while (i < 7){
356  if (nb[i] == -1) {
357  mFaceVxl.push_back( nb.idx );
358  break;
359  }
360  i++;
361  }
362 
363  i = 0; int n = 0;
364  while ( i < 7 ){
365  if (nb[i] == -1) {
366  n++;
367  if (n==2){
368  mEdgeVxl.push_back( nb.idx );
369  break;
370  }
371  }
372  i++;
373  }
374 
375  i = 0; n = 0;
376  while ( i < 7 ){
377  if (nb[i] == -1) {
378  n++;
379  if (n==3){
380  mCornerVxl.push_back( nb.idx );
381  break;
382  }
383  }
384  i++;
385  }
386  }
387 
389  Patch surfIdx(double u, double v){
390 
391  double fw = u * (mWidth - 1);
392  double fh = v * (mHeight - 1);
393 
394  int iw = floor ( fw );
395  int ih = floor ( fh );
396 
397  double rw = fw - iw;
398  double rh = fh - ih;
399 
400  if (iw == mWidth -1) { iw = mWidth -2; rw = 1.0; }
401  if (ih == mHeight -1) { ih = mHeight -2; rh = 1.0; }
402 
403  int a= ( idx ( iw, ih, 0 ) );
404  int b= ( idx ( iw + 1, ih, 0 ) );
405  int c= ( idx ( iw + 1, ih + 1, 0 ) );
406  int d= ( idx ( iw, ih + 1, 0 ) );
407 
408  return Patch( a, b, c, d, rw, rh);
409  }
410 
412  VPatch vidx(double u, double v, double w) const{
413  double fw = u * (mWidth - 1);
414  double fh = v * (mHeight - 1);
415  double fd = w * (mDepth-1);
416 
417  int iw = floor ( fw );
418  int ih = floor ( fh );
419  int id = floor ( fd );
420 
421  double rw = fw - iw;
422  double rh = fh - ih;
423  double rd = fd - id;
424 
425  if (iw == mWidth -1) { iw = mWidth -2; rw = 1.0; }
426  if (ih == mHeight -1) { ih = mHeight -2; rh = 1.0; }
427  if (id == mDepth -1) { id = mDepth -2; rd = 1.0; }
428 
429  int a= ( idx ( iw, ih, id ) );
430  int b= ( idx ( iw + 1, ih, id ) );
431  int c= ( idx ( iw + 1, ih + 1, id ) );
432  int d= ( idx ( iw, ih + 1, id ) );
433  int e= ( idx ( iw, ih, id +1) );
434  int f= ( idx ( iw + 1, ih, id +1) );
435  int g= ( idx ( iw + 1, ih + 1, id +1) );
436  int h= ( idx ( iw, ih + 1, id +1) );
437 
438  return VPatch( a, b, c, d, e, f, g, h, rw, rh, rd);
439  }
440 
442  Patch idxU(double t){
443  double fw = t * (mWidth - 1);
444 
445  int iw = floor ( fw );
446  double rw = fw - iw;
447  if (iw == mWidth -1) { iw = mWidth -2; rw = 1.0; }
448 
449  int a = idx(iw, 0, 0);
450  int b = idx(iw+1, 0, 0);
451  int c = idx(iw+2, 0, 0);
452 
453  return Patch(a, b, 0, 0, rw, 0);
454  }
455 
456 
457  //backwards compatibilty
458 
460  LPnt& gridAt(int w = 0, int h = 0, int d = 0) { return mPoint[ idx(w, h, d) ]; }
462  LPnt gridAt(int w = 0, int h = 0, int d = 0) const { return mPoint[ idx(w, h, d) ]; }
464  LPnt& grid(int i) { return mPoint[i]; }
466  LPnt grid(int i) const { return mPoint[i]; }
467 
468  LPnt surf(double u, double v){
469 
470  Patch p = surfIdx(u,v);
471 
472  LPnt a = mPoint[ p.a ];//gridAt ( iw, ih, 0 );
473  LPnt b = mPoint[ p.b ];//gridAt ( iw + 1, ih, 0 );
474  LPnt c = mPoint[ p.c ];//gridAt ( iw + 1, ih + 1, 0 );
475  LPnt d = mPoint[ p.d ];//gridAt ( iw, ih + 1, 0 );
476 
477  return nga::Round::null( Interp::surface<LPnt>( a,b,c,d, p.rw, p.rh) ) ;
478  }
479 
480  LPnt surfGrid(double u, double v) { return surfPnt(u,v); }
481  LPnt surfPnt(double u, double v){
482 
483  double pw = 1.0 / ( mWidth-1);
484  double ph = 1.0/ ( mHeight-1);
485  //double pd = 1.0 / mDepth;
486 
487  double fw = u / pw;
488  double fh = v / ph;
489  //double fd = w / pd;
490 
491  int iw = floor ( fw );
492  int ih = floor ( fh );
493  //int id = floor ( fd );
494  // cout << iw << " " << ih << endl;
495 
496  double rw = fw - iw;
497  double rh = fh - ih;
498  //double rd = fd - id;
499  // cout << rw << " " << rh << endl;
500 
501  LPnt a = grid( idx ( iw, ih, 0 ) );
502  LPnt d =grid( idx ( iw, ih + 1, 0 ));
503  LPnt b = grid( idx ( iw + 1, ih, 0 ));
504  LPnt c = grid( idx ( iw + 1, ih + 1, 0 ));
505 
506  return nga::Round::null( Interp::surface<LPnt> (a,b,c,d, rw, rh) );
507  }
508 
509  vector<int>& face() { return mFace; }
510  int face(int ix) const { return mFace[ix]; }
511 
512 
513  vector<int>& faceVxl() { return mFaceVxl; }
514  Vxl faceVxl(int ix) { return mVxl[ mFaceVxl[ix] ]; }
515 
516  LPnt * gridPtr() const { return mPoint; }
517  void gridPtr(LPnt * lp) { mPoint = lp; }
518 
519  protected:
520 
521  int mWidth, mHeight, mDepth;
522  int mNum, mNumVxl;
523  double mSpacingX,mSpacingY,mSpacingZ;
524 
525  //Points in Space
526  LPnt * mPoint;
527 
528  //vxl access
529  Vxl *mVxl;
530 
531  //neighbor access (i.e. mNbr[idx] returns list of neighbor indices)
532  Nbr *mNbr, *mNbrVxl;
533 
534  //INDICES
535  vector <int> mFace;
536  vector <int> mEdge;
537  vector <int> mCorner;
538  vector <int> mFaceVxl;
539  vector <int> mEdgeVxl;
540  vector <int> mCornerVxl;
541 
542 
543  };
544 
545 // extern template class CubicLattice< NPnt<5> >;
546 // extern template void CubicLattice< NPnt<5> >::initPoints();
547  template<> void CubicLattice< NPnt<5> > :: initPoints();
548 
549  #undef ITER
550  #undef ITERV
551  #undef ITEND
552  #undef BOUNDITER
553  #undef BOUNDITER0
554  #undef BOUNDEND
555 }
556 #endif
vsr::CubicLattice::py
double py(int j) const
Spatial Positions of jth element in y direction.
Definition: vsr_cubicLattice.h:150
vsr::CubicLattice::oh
double oh() const
Offset Height.
Definition: vsr_cubicLattice.h:142
vsr::CubicLattice::vidx
VPatch vidx(double u, double v, double w) const
Volume Index at u,v,w [0-1].
Definition: vsr_cubicLattice.h:412
vsr::CubicLattice::px
double px(int i) const
Spatial Positions of ith element in x direction.
Definition: vsr_cubicLattice.h:148
vsr::Nbr
Data Structure of Neighbors in a cartesian volume (left, right, bottom, top, front, back)
Definition: vsr_patch.h:65
vsr::CubicLattice::ow
double ow() const
Offset Width.
Definition: vsr_cubicLattice.h:138
vsr::CubicLattice::od
double od() const
Offset Depth.
Definition: vsr_cubicLattice.h:146
vsr::CubicLattice::FE
void FE(Nbr nb)
Routines to Find Face and Edge Boundary.
Definition: vsr_cubicLattice.h:312
vsr::CubicLattice::pz
double pz(int k) const
Spatial Positions of kth element in z direction.
Definition: vsr_cubicLattice.h:152
vsr::VPatch
Volume Patch Info Container for Euler integration of a 3d Field.
Definition: vsr_patch.h:54
vsr::CubicLattice::tw
double tw() const
Total Width.
Definition: vsr_cubicLattice.h:136
vsr::CubicLattice::td
double td() const
Total Depth.
Definition: vsr_cubicLattice.h:144
vsr::Patch
Info Container for Euler integration of a 2d Field.
Definition: vsr_patch.h:44
vsr::CubicLattice::grid
LPnt grid(int i) const
Get Grid (position) Data.
Definition: vsr_cubicLattice.h:466
vsr::CubicLattice::gridAt
LPnt & gridAt(int w=0, int h=0, int d=0)
Set grid data by Coordinate.
Definition: vsr_cubicLattice.h:460
vsr::CubicLattice
Discretized Volume Indexing (Isometric Cubic Lattice w/o data)
Definition: vsr_cubicLattice.h:61
vsr::CubicLattice::idx
int idx(int i, int j, int k=0) const
Get Index In Array.
Definition: vsr_cubicLattice.h:113
vsr
the versor library namespace
Definition: vsr_algebra.h:29
vsr::CubicLattice::idxU
Patch idxU(double t)
Indices of Line at T.
Definition: vsr_cubicLattice.h:442
vsr::CubicLattice::surfIdx
Patch surfIdx(double u, double v)
Indicex of surface at u, v [0, 1].
Definition: vsr_cubicLattice.h:389
vsr::CubicLattice::gridAt
LPnt gridAt(int w=0, int h=0, int d=0) const
Get grid data by Coordinate.
Definition: vsr_cubicLattice.h:462
vsr::CubicLattice::th
double th() const
Total Height.
Definition: vsr_cubicLattice.h:140
vsr::Vxl
Definition: vsr_patch.h:107
vsr::CubicLattice::grid
LPnt & grid(int i)
Set Grid (position) Data.
Definition: vsr_cubicLattice.h:464
vsr::CubicLattice::vxlAt
Vxl vxlAt(const V &tv) const
Voxel of Vector v.
Definition: vsr_cubicLattice.h:219
Generated on Tue Jul 28 2015 16:26:05 for versor by   doxygen 1.8.10