vtkRectilinearGridToTetrahedra.cxx 17.5 KB
Newer Older
1
2
3
4
5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkRectilinearGridToTetrahedra.cxx

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7
8
9
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

10
11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
13
14
15
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkRectilinearGridToTetrahedra.h"
16
17
18

#include "vtkCellArray.h"
#include "vtkCellData.h"
Clinton Stimpson's avatar
   
Clinton Stimpson committed
19
#include "vtkSignedCharArray.h"
20
#include "vtkExecutive.h"
21
22
23
#include "vtkFloatArray.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
24
#include "vtkInformationVector.h"
25
#include "vtkIntArray.h"
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include "vtkObjectFactory.h"
#include "vtkRectilinearGrid.h"
#include "vtkUnstructuredGrid.h"
#include "vtkVoxel.h"

vtkStandardNewMacro(vtkRectilinearGridToTetrahedra);

// ways to convert to a voxel to tetrahedra.
// Note that the values 0 and 1 and -1 and 2 are important in
// DetermineGridDivisionTypes()
#define VTK_TETRAHEDRALIZE_5                 0
#define VTK_TETRAHEDRALIZE_5_FLIP            1
#define VTK_TETRAHEDRALIZE_6                 6
#define VTK_TETRAHEDRALIZE_12_CONFORM        -1
#define VTK_TETRAHEDRALIZE_12_CONFORM_FLIP   2
#define VTK_TETRAHEDRALIZE_12                10

//-------------------------------------------------------------------------

vtkRectilinearGridToTetrahedra::vtkRectilinearGridToTetrahedra()
{
47
48
  this->TetraPerCell  = VTK_VOXEL_TO_5_TET;
  this->RememberVoxelId = 0;
49
50
51
52
}

//----------------------------------------------------------------------------

Ken Martin's avatar
Ken Martin committed
53
54
55
56
57
58
59
void vtkRectilinearGridToTetrahedra::SetInput(const double ExtentX, 
                                              const double ExtentY,
                                              const double ExtentZ, 
                                              const double SpacingX,
                                              const double SpacingY,
                                              const double SpacingZ,
                                              const double tol)
60
{
Ken Martin's avatar
Ken Martin committed
61
62
  double Extent[3];
  double Spacing[3];
63
64
  Extent[0]  = ExtentX;    Extent[1] = ExtentY;    Extent[2] = ExtentZ;
  Spacing[0] = SpacingX;  Spacing[1] = SpacingY;  Spacing[2] = SpacingZ;
65
  this->SetInput(Extent,Spacing,tol);
66
67
68
69
70
71
}


//----------------------------------------------------------------------------

// Create an input for the filter
Ken Martin's avatar
Ken Martin committed
72
73
74
void vtkRectilinearGridToTetrahedra::SetInput(const double Extent[3], 
                                              const double Spacing[3], 
                                              const double tol)
75
76
77
78
79
80
81
{
  //
  // Determine the number of points in each direction, and the positions
  // The last voxel may have a different spacing to fit inside
  // the selected region
  // 

82
  int i, j;
83
84
  int NumPointsInDir[3];
  vtkFloatArray *Coord[3];
85
  for(i = 0;i<3;i++)
86
    {
Ken Martin's avatar
Ken Martin committed
87
    double NumRegion = Extent[i]/Spacing[i];
88
89
90
91
92
93
94
95
96
97
98
99

    // If we are really close to an integer number of elements, use the 
    // integer number
    if (fabs(NumRegion - floor(NumRegion+0.5)) < tol*Spacing[i])
      NumPointsInDir[i] = ((int) floor(NumRegion+0.5)) + 1;
    else 
      NumPointsInDir[i] = (int) ceil(Extent[i]/Spacing[i])+1;
    Coord[i] = vtkFloatArray::New();
    Coord[i]->SetNumberOfValues(NumPointsInDir[i]+1);

     // The last data point inserted is exactly the Extent
    // Thus avoiding a bit of numerical error.
100
    for(j=0;j<NumPointsInDir[i]-1;j++)
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
      {
      Coord[i]->SetValue(j,Spacing[i]*j);
      }
    Coord[i]->SetValue(NumPointsInDir[i]-1,Extent[i]);
    }

  //
  // Form the grid
  //

  vtkRectilinearGrid *RectGrid = vtkRectilinearGrid::New();
  RectGrid->SetDimensions(NumPointsInDir);
  RectGrid->SetXCoordinates(Coord[0]);
  RectGrid->SetYCoordinates(Coord[1]);
  RectGrid->SetZCoordinates(Coord[2]);

  Coord[0]->Delete();
  Coord[1]->Delete();
  Coord[2]->Delete();

  // Get the reference counting right.
122
  this->Superclass::SetInput(RectGrid);
123
124
125
126
127
128
129
130
  RectGrid->Delete();
}

//----------------------------------------------------------------------------

// Determine how to Divide each voxel in the vtkRectilinearGrid
void vtkRectilinearGridToTetrahedra::DetermineGridDivisionTypes(
                                             vtkRectilinearGrid *RectGrid,
Clinton Stimpson's avatar
   
Clinton Stimpson committed
131
                                             vtkSignedCharArray *VoxelSubdivisionType,
132
133
134
135
136
137
138
139
140
141
142
                                             const int &TetraPerCell)
{
  int numRec = RectGrid->GetNumberOfCells();
  int NumPointsInDir[3];
  RectGrid->GetDimensions(NumPointsInDir);

  // How to break into Tetrahedra.
  // For division into 5's, we need to flip from one orientation to
  // the next

  int Rec[3];
143
  int flip;
144
  int i;
145
146
147
148

  switch (TetraPerCell)
    {
    case (VTK_VOXEL_TO_12_TET):
149
      for(i=0;i<numRec;i++)
150
151
152
153
154
        {
        VoxelSubdivisionType->SetValue(i,VTK_TETRAHEDRALIZE_12);
        }
    break;
    case (VTK_VOXEL_TO_6_TET):
155
      for(i=0;i<numRec;i++)
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
        {
        VoxelSubdivisionType->SetValue(i,VTK_TETRAHEDRALIZE_6);
        }
    break;
    case (VTK_VOXEL_TO_5_TET):
      for(Rec[0] = 0; Rec[0]<NumPointsInDir[0]-1; Rec[0]++)
        {
        for(Rec[1] = 0; Rec[1]<NumPointsInDir[1]-1; Rec[1]++)
          {
          flip = ( Rec[1] + Rec[0] ) % 2;
          for(Rec[2] = 0; Rec[2]<NumPointsInDir[2]-1; Rec[2]++)
            {
            VoxelSubdivisionType->SetValue(RectGrid->ComputeCellId(Rec),flip);
            flip = 1 - flip;
            }
          }
        }
    break;
    case (VTK_VOXEL_TO_5_AND_12_TET):
      for(Rec[0] = 0; Rec[0]<NumPointsInDir[0]-1; Rec[0]++)
        {
        for(Rec[1] = 0; Rec[1]<NumPointsInDir[1]-1; Rec[1]++)
          {
          flip = ( Rec[1] + Rec[0] ) % 2;
          for(Rec[2] = 0; Rec[2]<NumPointsInDir[2]-1; Rec[2]++)
            {
            int CellId = RectGrid->ComputeCellId(Rec);
            if (VoxelSubdivisionType->GetValue(CellId) == 12)
              VoxelSubdivisionType->SetValue(CellId,3*flip-1);
            else VoxelSubdivisionType->SetValue(CellId,flip);
            flip = 1 - flip;
            }
          }
        }
    break;
    }
}

//----------------------------------------------------------------------------

// Take the grid and make it into a tetrahedral mesh.
void vtkRectilinearGridToTetrahedra::GridToTetMesh(vtkRectilinearGrid *RectGrid,
Clinton Stimpson's avatar
   
Clinton Stimpson committed
198
                                vtkSignedCharArray *VoxelSubdivisionType,
199
200
201
202
                                const int &TetraPerCell,
                                const int &RememberVoxelId,
                                vtkUnstructuredGrid *TetMesh)
{
203
  int i, j;
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
  int numPts = RectGrid->GetNumberOfPoints();
  int numRec = RectGrid->GetNumberOfCells();

  // We need a point list and a cell list
  vtkPoints *NodePoints = vtkPoints::New();
  vtkCellArray *TetList = vtkCellArray::New();

  // Guess number of points and cells!! 
  // For mixture of 5 and 12 tet per cell, 
  // it is easier to way overguess to avoid re-allocation
  // slowness and range checking during insertion.

  switch (TetraPerCell)
    {
    case (VTK_VOXEL_TO_5_TET):
      NodePoints->Allocate(numPts);
      TetList->Allocate(numPts*5*5,numPts);
    break;
    case (VTK_VOXEL_TO_5_AND_12_TET):
    case (VTK_VOXEL_TO_12_TET):
      NodePoints->Allocate(numPts*2);
      TetList->Allocate(numPts*5*12,numPts);
    break;
    }

  // Start by copying over the points
230
  for(i=0;i<numPts;i++)
231
    {
232
    NodePoints->InsertNextPoint(RectGrid->GetPoint(i));
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
    }

  // If they want, we can add Scalar Data
  // to the Tets indicating the Voxel Id the tet
  // came from.
  vtkIntArray *TetOriginalVoxel = NULL;
  if (RememberVoxelId) {
    TetOriginalVoxel = vtkIntArray::New();
    TetOriginalVoxel->Allocate(12*numRec);
  }

  // 9 ids, 8 corners and a possible center to be added later
  //        during the tet creation
  vtkIdList *VoxelCorners = vtkIdList::New();
  VoxelCorners->SetNumberOfIds(9);

  int NumTetFromVoxel;
250
  for(i=0;i<numRec;i++)
251
252
253
254
255
256
    {
    RectGrid->GetCellPoints(i,VoxelCorners);
    NumTetFromVoxel = TetrahedralizeVoxel(VoxelCorners,
                                          (int)VoxelSubdivisionType->GetValue(i),
                                          NodePoints,TetList);
    if (RememberVoxelId)
Michalel Halle's avatar
Michalel Halle committed
257
      {
258
      for(j=0;j<NumTetFromVoxel;j++)
259
260
261
        {
        TetOriginalVoxel->InsertNextValue(i);
        }
Michalel Halle's avatar
Michalel Halle committed
262
      }
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
    }
  
  //
  // It may be there are extra points at the end of the PointList.
  //

  NodePoints->Squeeze();

  //
  // Form the Mesh
  //

  // Need to tell the tet mesh that every cell  is a Tetrahedron
  int numTet = TetList->GetNumberOfCells();
  int *CellTypes = new int[numTet];
278
  for(i=0;i<numTet;i++)
279
    {
280
    CellTypes[i] = VTK_TETRA;
281
282
283
284
285
286
287
288
289
290
291
292
    }

  TetMesh->SetPoints(NodePoints);
  TetMesh->SetCells(CellTypes,TetList);

  //
  // Add Scalar Types if wanted
  //

  if(RememberVoxelId)
    {
    TetOriginalVoxel->Squeeze();
293
294
    int idx = TetMesh->GetCellData()->AddArray(TetOriginalVoxel);
    TetMesh->GetCellData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
    TetOriginalVoxel->Delete();
    }

  //
  // Clean Up
  // 
  delete[] CellTypes;
  NodePoints->Delete();
  TetList->Delete();
  VoxelCorners->Delete();

  TetMesh->Squeeze();
}

//----------------------------------------------------------------------------
// Helper Function for Tetrahedralize Voxel
inline void vtkRectilinearGridToTetrahedra::TetrahedralizeAddCenterPoint(
                                                    vtkIdList *VoxelCorners,
                                                    vtkPoints *NodeList)
{
  // Need to add a center point
Ken Martin's avatar
Ken Martin committed
316
  double c1[3], c2[3];
317
318
  NodeList->GetPoint(VoxelCorners->GetId(0), c2);
  NodeList->GetPoint(VoxelCorners->GetId(7), c1);
Ken Martin's avatar
Ken Martin committed
319
  double center[3];
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
  center[0] = (c1[0] + c2[0])/2.0;
  center[1] = (c1[1] + c2[1])/2.0;
  center[2] = (c1[2] + c2[2])/2.0;
  
  VoxelCorners->InsertId(8,NodeList->InsertNextPoint(center));
}

//----------------------------------------------------------------------------

// Split A Cube into Tetrahedrons
// According to the DivisionType
// There had better be 0..8 voxel corners, though only 0..7 maybe needed.
// Why? This function may add id 8 to VoxelCorners.
// If a point needs to be inserted into the nodelist, itselt
// it at NextPointId. Assume there is space in the nodelist.
// Return the number of Tets Added.


int vtkRectilinearGridToTetrahedra::TetrahedralizeVoxel(vtkIdList *VoxelCorners,
                                     const int &DivisionType,
                                     vtkPoints *NodeList,
                                     vtkCellArray *TetList)
{

// See vtkVoxel::Triangulate
/*  Looking at the rect: Corner labeling 

     0  1
     2  3

   Directly behind them:
    4   5
    6   7

and 8 is in the middle of the cube if used

Want right handed Tetrahedra...
*/

  // Split voxel in 2 along diagonal, 3 tets on either side
 static int tet6[6][4] =
   { 
     {1,6,2,3}, {1,6,7,5}, {1,6,3,7},
     {1,6,0,2}, {1,6,5,4}, {1,6,4,0},
   };

 static int tet5[5][4]      = 
   { {0,1,4,2},{1,4,7,5},{1,4,2,7},{1,2,3,7},{2,7,4,6} };
 static int tet5flip[5][4]  = 
   { {3,1,0,5}, {0,3,6,2}, {3,5,6,7}, {0,6,5,4}, {0,3,5,6}};

  // 12 tet to confirm to tet5
  static int tet12_conform[12][4] = {
  /* Left side */
    {8,2,4,0},
    {8,4,2,6},
  /* Back side */
    {8,7,4,6},
    {8,4,7,5},
  /* Bottom side */
    {8,7,2,3},
    {8,2,7,6},
  /* Right side */
    {8,7,1,5},
    {8,1,7,3},
  /* Front side */
    {8,1,2,0},
    {8,2,1,3},
  /* Top side */
    {8,4,1,0}, 
    {8,1,4,5}};

  // 12 tet to confirm to tet5flip
  static int tet12_conform_flip[12][4] = {
  /* Left side */ 
    {8,0,6,4},
    {8,6,0,2},
  /* Back side */ 
    {8,5,6,7},
    {8,6,5,4},
  /* Bottom side */ 
    {8,3,6,2},
    {8,6,3,7},
  /* Right side */
    {8,3,5,7},
    {8,5,3,1},
  /* Front side */
    {8,3,0,1},
    {8,0,3,2},
  /* Top side */
    {8,5,0,4},
    {8,0,5,1}};

  // 12 tet chosen to have the least number of edges per node
  static int tet12[12][4] = {
  /* Left side */
    {8,2,4,0},
    {8,4,2,6},
  /* Back side */                                                     
    {8,7,4,6},
    {8,4,7,5},
  /* Right side */                                                    
    {8,3,5,7},
    {8,5,3,1},
  /* Front side */                                                   
    {8,3,0,1},
    {8,0,3,2},
  /* Top side */                                                     
    {8,5,0,4},
    {8,0,5,1},
  /* Bottom side */                                                  
    {8,7,2,3},
    {8,2,7,6}};


435
  int i, j;
436
437
  // Get the point Ids
  int numTet = 0; // =0 removes warning messages
438
  vtkIdType TetPts[4];
439
440
441
442
443

  switch (DivisionType)
    {
    case (VTK_TETRAHEDRALIZE_6) : 
      numTet = 6;
444
      for(i=0;i<numTet;i++)
445
        {
446
        for(j=0;j<4;j++)
447
448
449
          {
          TetPts[j] = VoxelCorners->GetId(tet6[i][j]);
          }
450
        TetList->InsertNextCell((vtkIdType)4,TetPts);
451
452
453
454
        }
    break; 
    case (VTK_TETRAHEDRALIZE_5) : 
      numTet = 5;
455
      for(i=0;i<numTet;i++)
456
        {
457
        for(j=0;j<4;j++)
458
459
460
          {
          TetPts[j] = VoxelCorners->GetId(tet5[i][j]);
          }
461
        TetList->InsertNextCell((vtkIdType)4,TetPts);
462
463
464
465
        }
    break; 
    case (VTK_TETRAHEDRALIZE_5_FLIP) : 
      numTet = 5;
466
      for(i=0;i<numTet;i++)
467
        {
468
          for(j=0;j<4;j++)
469
470
471
            {
            TetPts[j] = VoxelCorners->GetId(tet5flip[i][j]);
            }
472
          TetList->InsertNextCell((vtkIdType)4,TetPts);
473
474
475
476
477
        }
    break; 
    case (VTK_TETRAHEDRALIZE_12) :
      numTet = 12;
      TetrahedralizeAddCenterPoint(VoxelCorners,NodeList);
478
      for(i=0;i<numTet;i++)
479
        {
480
        for(j=0;j<4;j++)
481
482
483
          {
          TetPts[j] = VoxelCorners->GetId(tet12[i][j]);
          }
484
        TetList->InsertNextCell((vtkIdType)4,TetPts);
485
486
487
488
489
        }
    break;
    case (VTK_TETRAHEDRALIZE_12_CONFORM) :
      numTet = 12;
      TetrahedralizeAddCenterPoint(VoxelCorners,NodeList);
490
      for(i=0;i<numTet;i++)
491
        {
492
        for(j=0;j<4;j++)
493
494
495
          {
          TetPts[j] = VoxelCorners->GetId(tet12_conform[i][j]);
          }
496
        TetList->InsertNextCell((vtkIdType)4,TetPts);
497
498
499
500
501
        }
    break;
    case (VTK_TETRAHEDRALIZE_12_CONFORM_FLIP) :
      numTet = 12;
      TetrahedralizeAddCenterPoint(VoxelCorners,NodeList);
502
      for(i=0;i<numTet;i++)
503
        {
504
        for(j=0;j<4;j++)
505
506
507
          {
          TetPts[j] = VoxelCorners->GetId(tet12_conform_flip[i][j]);
          }
508
        TetList->InsertNextCell((vtkIdType)4,TetPts);
509
510
511
512
513
514
515
516
        }
    break;
    }
  return numTet;
}

//----------------------------------------------------------------------------

517
518
519
520
int vtkRectilinearGridToTetrahedra::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
521
{
522
523
524
525
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

526
  // get the input and output
527
528
529
530
  vtkRectilinearGrid *RectGrid = vtkRectilinearGrid::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));
531
532
533

  // Create internal version of VoxelSubdivisionType
  // VoxelSubdivisionType indicates how to subdivide each cell
Clinton Stimpson's avatar
   
Clinton Stimpson committed
534
535
  vtkSignedCharArray *VoxelSubdivisionType;
  VoxelSubdivisionType = vtkSignedCharArray::New();
536
537
538
539
540
541
542

  // If we have a mixture of 5 and 12 Tet, we need to get the information from
  // the scalars of the Input. Note that we will modify the array internally
  // so we need to copy it.
  if (this->TetraPerCell == VTK_VOXEL_TO_5_AND_12_TET)
    {
    vtkDataArray *TempVoxelSubdivisionType = RectGrid->GetCellData()->GetScalars();
543
    if(TempVoxelSubdivisionType == NULL)
544
545
      {
      vtkErrorMacro(<< "Scalars to input Should be set!");
546
      return 1;
547
548
      }
    VoxelSubdivisionType->SetNumberOfValues(RectGrid->GetNumberOfCells());
Clinton Stimpson's avatar
   
Clinton Stimpson committed
549
    VoxelSubdivisionType->vtkSignedCharArray::DeepCopy(TempVoxelSubdivisionType);
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
    }
  else
    { // Otherwise, just create the GridDivisionTypes
    VoxelSubdivisionType->SetNumberOfValues(RectGrid->GetNumberOfCells());
    }

  vtkDebugMacro(<<"Number of points: " 
                << RectGrid->GetNumberOfPoints());
  vtkDebugMacro(<< "Number of voxels in input: " 
                << RectGrid->GetNumberOfCells());

  // Determine how each Cell should be subdivided
  DetermineGridDivisionTypes(RectGrid,VoxelSubdivisionType,
                             this->TetraPerCell);

  // Subdivide each cell to a tetrahedron, forming the TetMesh
  GridToTetMesh(RectGrid,VoxelSubdivisionType,this->TetraPerCell,
567
                this->RememberVoxelId,output);
568
569

  vtkDebugMacro(<< "Number of output points: " 
570
                << output->GetNumberOfPoints());
571
  vtkDebugMacro(<< "Number of output tetrahedra: " 
572
                << output->GetNumberOfCells());
573
574
575

  // Clean Up
  VoxelSubdivisionType->Delete();
576
577

  return 1;
578
579
}

580
581
582
583
584
585
586
587
588
589
590
591
592
//----------------------------------------------------------------------------
int
vtkRectilinearGridToTetrahedra
::FillInputPortInformation(int port, vtkInformation* info)
{
  if(!this->Superclass::FillInputPortInformation(port, info))
    {
    return 0;
    }
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkRectilinearGrid");
  return 1;
}

593
594
595
596
//----------------------------------------------------------------------------

void vtkRectilinearGridToTetrahedra::PrintSelf(ostream& os, vtkIndent indent)
{
597
  this->Superclass::PrintSelf(os,indent);
598

599
600
  os << indent << "Mesh Type: " << this->TetraPerCell << "\n";
  os << indent << "RememberVoxel Id: " << this->RememberVoxelId << "\n";
601
}