vtkDataSetTriangleFilter.cxx 11.7 KB
Newer Older
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkDataSetTriangleFilter.cxx

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

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

=========================================================================*/
#include "vtkDataSetTriangleFilter.h"
16

17 18 19 20
#include "vtkCellData.h"
#include "vtkCellType.h"
#include "vtkGenericCell.h"
#include "vtkIdList.h"
21
#include "vtkImageData.h"
22 23
#include "vtkInformation.h"
#include "vtkInformationVector.h"
24
#include "vtkObjectFactory.h"
25
#include "vtkOrderedTriangulator.h"
26
#include "vtkPointData.h"
27 28 29
#include "vtkStructuredGrid.h"
#include "vtkStructuredPoints.h"
#include "vtkUnstructuredGrid.h"
30
#include "vtkRectilinearGrid.h"
31
#include "vtkUnsignedCharArray.h"
32

Brad King's avatar
Brad King committed
33
vtkStandardNewMacro(vtkDataSetTriangleFilter);
34

35 36 37 38 39
vtkDataSetTriangleFilter::vtkDataSetTriangleFilter()
{
  this->Triangulator = vtkOrderedTriangulator::New();
  this->Triangulator->PreSortedOff();
  this->Triangulator->UseTemplatesOn();
40
  this->TetrahedraOnly = 0;
41 42
}

43 44
vtkDataSetTriangleFilter::~vtkDataSetTriangleFilter()
{
45
  this->Triangulator->Delete();
46
  this->Triangulator = nullptr;
47 48
}

49 50 51 52
int vtkDataSetTriangleFilter::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
53
{
54 55 56 57
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

58
  // get the input and output
59 60 61 62 63
  vtkDataSet *input = vtkDataSet::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

64
  if (input->IsA("vtkStructuredPoints") ||
65
      input->IsA("vtkStructuredGrid") ||
66 67
      input->IsA("vtkImageData") ||
      input->IsA("vtkRectilinearGrid"))
68
  {
69
    this->StructuredExecute(input, output);
70
  }
71
  else
72
  {
73
    this->UnstructuredExecute(input, output);
74
  }
75 76

  vtkDebugMacro(<<"Produced " << this->GetOutput()->GetNumberOfCells() << " cells");
77 78

  return 1;
79 80
}

81 82
void vtkDataSetTriangleFilter::StructuredExecute(vtkDataSet *input,
                                                 vtkUnstructuredGrid *output)
83
{
Amy Squillacote's avatar
Amy Squillacote committed
84 85
  int dimensions[3], i, j, k, l, m;
  vtkIdType newCellId, inId;
86 87 88 89 90
  vtkCellData *inCD = input->GetCellData();
  vtkCellData *outCD = output->GetCellData();
  vtkPoints *cellPts = vtkPoints::New();
  vtkPoints *newPoints = vtkPoints::New();
  vtkIdList *cellPtIds = vtkIdList::New();
Amy Squillacote's avatar
Amy Squillacote committed
91 92
  int numSimplices, numPts, dim, type;
  vtkIdType pts[4], num;
93

94 95
  // Create an array of points. This does an explicit creation
  // of each point.
96
  num = input->GetNumberOfPoints();
97
  newPoints->SetNumberOfPoints(num);
98
  for (i = 0; i < num; ++i)
99
  {
100
    newPoints->SetPoint(i,input->GetPoint(i));
101
  }
102 103 104

  outCD->CopyAllocate(inCD,input->GetNumberOfCells()*5);
  output->Allocate(input->GetNumberOfCells()*5);
105

106
  if (input->IsA("vtkStructuredPoints"))
107
  {
108
    static_cast<vtkStructuredPoints*>(input)->GetDimensions(dimensions);
109
  }
110
  else if (input->IsA("vtkStructuredGrid"))
111
  {
112
    static_cast<vtkStructuredGrid*>(input)->GetDimensions(dimensions);
113
  }
114
  else if (input->IsA("vtkImageData"))
115
  {
116
    static_cast<vtkImageData*>(input)->GetDimensions(dimensions);
117
  }
118
  else if (input->IsA("vtkRectilinearGrid"))
119
  {
120
    static_cast<vtkRectilinearGrid*>(input)->GetDimensions(dimensions);
121
  }
122
  else
123
  {
124 125 126 127 128 129 130
    // Every kind of structured data is listed above, this should never happen.
    // Report an error and produce no output.
    vtkErrorMacro("Unrecognized data set " << input->GetClassName());
    // Dimensions of 1x1x1 means a single point, i.e. dimensionality of zero.
    dimensions[0] = 1;
    dimensions[1] = 1;
    dimensions[2] = 1;
131
  }
132

133 134 135
  dimensions[0] = dimensions[0] - 1;
  dimensions[1] = dimensions[1] - 1;
  dimensions[2] = dimensions[2] - 1;
136

137
  vtkIdType numSlices = ( dimensions[2] > 0 ? dimensions[2] : 1 );
138
  int abort=0;
139
  for (k = 0; k < numSlices && !abort; k++)
140
  {
141
    this->UpdateProgress(static_cast<double>(k) / numSlices);
142 143
    abort = this->GetAbortExecute();

144
    for (j = 0; j < dimensions[1]; j++)
145
    {
146
      for (i = 0; i < dimensions[0]; i++)
147
      {
148
        inId = i+(j+(k*dimensions[1]))*dimensions[0];
149
        vtkCell *cell = input->GetCell(i, j, k);
150
        if ((i+j+k)%2 == 0)
151
        {
152
          cell->Triangulate(0, cellPtIds, cellPts);
153
        }
154
        else
155
        {
156
          cell->Triangulate(1, cellPtIds, cellPts);
157
        }
158

159
        dim = cell->GetCellDimension() + 1;
160

161 162 163 164
        numPts = cellPtIds->GetNumberOfIds();
        numSimplices = numPts / dim;
        type = 0;
        switch (dim)
165
        {
166 167 168 169 170 171 172 173
          case 1:
            type = VTK_VERTEX;    break;
          case 2:
            type = VTK_LINE;      break;
          case 3:
            type = VTK_TRIANGLE;  break;
          case 4:
            type = VTK_TETRA;     break;
174
        }
175
        if (!this->TetrahedraOnly || type == VTK_TETRA)
176
        {
177
          for (l = 0; l < numSimplices; l++ )
178
          {
179
            for (m = 0; m < dim; m++)
180
            {
181
              pts[m] = cellPtIds->GetId(dim*l+m);
182
            }
183 184 185
            // copy cell data
            newCellId = output->InsertNextCell(type, dim, pts);
            outCD->CopyData(inCD, inId, newCellId);
186 187 188 189 190
          }//for all simplices
        }
      }//i dimension
    }//j dimension
  }//k dimension
191

192 193 194 195
  // Update output
  output->SetPoints(newPoints);
  output->GetPointData()->PassData(input->GetPointData());
  output->Squeeze();
196

197 198 199 200 201
  newPoints->Delete();
  cellPts->Delete();
  cellPtIds->Delete();
}

202 203 204 205
// 3D cells use the ordered triangulator. The ordered triangulator is used
// to create templates on the fly. Once the templates are created then they
// are used to produce the final triangulation.
//
206 207
void vtkDataSetTriangleFilter::UnstructuredExecute(vtkDataSet *dataSetInput,
                                                   vtkUnstructuredGrid *output)
208
{
209
  vtkPointSet *input = static_cast<vtkPointSet*>(dataSetInput); //has to be
Amy Squillacote's avatar
Amy Squillacote committed
210
  vtkIdType numCells = input->GetNumberOfCells();
211
  vtkGenericCell *cell;
212
  vtkIdType newCellId, j;
Amy Squillacote's avatar
Amy Squillacote committed
213
  int k;
214 215
  vtkCellData *inCD=input->GetCellData();
  vtkCellData *outCD=output->GetCellData();
216 217
  vtkPoints *cellPts;
  vtkIdList *cellPtIds;
Amy Squillacote's avatar
Amy Squillacote committed
218
  vtkIdType ptId, numTets, ncells;
219 220
  int numPts, type;
  int numSimplices, dim;
221
  vtkIdType pts[4];
222
  double x[3];
223

224
  if (numCells == 0)
225
  {
226
    return;
227
  }
228

229 230 231
  vtkUnstructuredGrid * inUgrid =
    vtkUnstructuredGrid::SafeDownCast(dataSetInput);
  if (inUgrid)
232
  {
233 234 235
    //avoid doing cell simplification if all cells are already simplices
    vtkUnsignedCharArray* cellTypes = inUgrid->GetCellTypesArray();
    if (cellTypes)
236
    {
237 238
      int allsimplices = 1;
      for (vtkIdType cellId = 0; cellId < cellTypes->GetSize() && allsimplices; cellId++)
239
      {
240
        switch (cellTypes->GetValue(cellId))
241
        {
242 243 244 245
          case VTK_TETRA:
            break;
          case VTK_VERTEX:
          case VTK_LINE:
246
          case VTK_TRIANGLE:
247
            if (this->TetrahedraOnly)
248
            {
249
              allsimplices = 0; //don't shallowcopy need to stip non tets
250
            }
251 252 253 254 255
            break;
          default:
            allsimplices = 0;
            break;
        }
256
      }
257
      if (allsimplices)
258
      {
259 260 261 262
        output->ShallowCopy(input);
        return;
      }
    }
263
  }
264

265 266 267 268
  cell = vtkGenericCell::New();
  cellPts = vtkPoints::New();
  cellPtIds = vtkIdList::New();

269
  // Create an array of points
270 271
  vtkCellData *tempCD = vtkCellData::New();
  tempCD->ShallowCopy(inCD);
272
  tempCD->SetActiveGlobalIds(nullptr);
273 274

  outCD->CopyAllocate(tempCD, input->GetNumberOfCells()*5);
275
  output->Allocate(input->GetNumberOfCells()*5);
276

277 278 279 280
  // Points are passed through
  output->SetPoints(input->GetPoints());
  output->GetPointData()->PassData(input->GetPointData());

281 282 283
  int abort=0;
  vtkIdType updateTime = numCells/20 + 1;  // update roughly every 5%
  for (vtkIdType cellId=0; cellId < numCells && !abort; cellId++)
284
  {
285
    if ( !(cellId % updateTime) )
286
    {
287
      this->UpdateProgress(static_cast<double>(cellId) / numCells);
288
      abort = this->GetAbortExecute();
289
    }
290 291

    input->GetCell(cellId, cell);
292
    dim = cell->GetCellDimension();
293

294
    if (cell->GetCellType() == VTK_POLYHEDRON) //polyhedron
295
    {
296 297 298
      dim = 4;
      cell->Triangulate(0, cellPtIds, cellPts);
      numPts = cellPtIds->GetNumberOfIds();
299

300 301 302 303
      numSimplices = numPts / dim;
      type = VTK_TETRA;

      for ( j=0; j < numSimplices; j++ )
304
      {
305
        for (k=0; k<dim; k++)
306
        {
307
          pts[k] = cellPtIds->GetId(dim*j+k);
308
        }
309 310
        // copy cell data
        newCellId = output->InsertNextCell(type, dim, pts);
311
        outCD->CopyData(tempCD, cellId, newCellId);
312
      }
313
    }
314 315

    else if ( dim == 3 ) //use ordered triangulation
316
    {
317
      numPts = cell->GetNumberOfPoints();
318
      double *p, *pPtr=cell->GetParametricCoords();
319 320
      this->Triangulator->InitTriangulation(0.0,1.0, 0.0,1.0, 0.0,1.0, numPts);
      for (p=pPtr, j=0; j<numPts; j++, p+=3)
321
      {
322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
        // the wedge is "flipped" compared to other cells in that
        // the normal of the first face points out instead of in
        // so we flip the way we pass the points to the triangulator
        const vtkIdType wedgemap[18] = {3, 4, 5, 0, 1, 2, 9, 10, 11, 6, 7, 8, 12, 13, 14, 15, 16, 17};
        type = cell->GetCellType();
        if (type == VTK_WEDGE || type == VTK_QUADRATIC_WEDGE || type == VTK_QUADRATIC_LINEAR_WEDGE ||
            type == VTK_BIQUADRATIC_QUADRATIC_WEDGE)
        {
          ptId = cell->PointIds->GetId(wedgemap[j]);
          cell->Points->GetPoint(wedgemap[j], x);
        }
        else
        {
          ptId = cell->PointIds->GetId(j);
          cell->Points->GetPoint(j, x);
        }
338
        this->Triangulator->InsertPoint(ptId, x, p, 0);
339
      }//for all cell points
340
      if ( cell->IsPrimaryCell() ) //use templates if topology is fixed
341
      {
342 343 344
        int numEdges=cell->GetNumberOfEdges();
        this->Triangulator->TemplateTriangulate(cell->GetCellType(),
                                                numPts,numEdges);
345
      }
346
      else //use ordered triangulator
347
      {
348
        this->Triangulator->Triangulate();
349
      }
350 351 352

      ncells = output->GetNumberOfCells();
      numTets = this->Triangulator->AddTetras(0,output);
353

354
      for (j=0; j < numTets; j++)
355
      {
356
        outCD->CopyData(tempCD, cellId, ncells+j);
357
      }
358
    }
359

360
    else if (!this->TetrahedraOnly) //2D or lower dimension
361
    {
362 363
      dim++;
      cell->Triangulate(0, cellPtIds, cellPts);
364
      numPts = cellPtIds->GetNumberOfIds();
365

366
      numSimplices = numPts / dim;
367 368
      type = 0;
      switch (dim)
369
      {
370 371 372 373 374 375
        case 1:
          type = VTK_VERTEX;    break;
        case 2:
          type = VTK_LINE;      break;
        case 3:
          type = VTK_TRIANGLE;  break;
376
      }
377 378

      for ( j=0; j < numSimplices; j++ )
379
      {
380
        for (k=0; k<dim; k++)
381
        {
382
          pts[k] = cellPtIds->GetId(dim*j+k);
383
        }
384 385
        // copy cell data
        newCellId = output->InsertNextCell(type, dim, pts);
386
        outCD->CopyData(tempCD, cellId, newCellId);
387 388 389
      }
    } //if 2D or less cell
  } //for all cells
390

391 392
  // Update output
  output->Squeeze();
393

394 395
  tempCD->Delete();

396 397 398 399 400
  cellPts->Delete();
  cellPtIds->Delete();
  cell->Delete();
}

401 402 403 404 405 406
int vtkDataSetTriangleFilter::FillInputPortInformation(int, vtkInformation *info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
  return 1;
}

407 408
void vtkDataSetTriangleFilter::PrintSelf(ostream& os, vtkIndent indent)
{
Brad King's avatar
Brad King committed
409
  this->Superclass::PrintSelf(os,indent);
410
  os << indent << "TetrahedraOnly: " << (this->TetrahedraOnly ? "On":"Off") << "\n";
411 412
}