Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkMeshUtil.h
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #ifndef MITKMESHUTIL_H_INCLUDED
14 #define MITKMESHUTIL_H_INCLUDED
15 
16 #if (_MSC_VER == 1200)
17 #error MeshUtils currently not supported for MS Visual C++ 6.0. Sorry.
18 #endif
19 
20 //#include <itkMesh.h>
21 #include <itkCellInterface.h>
22 #include <itkLineCell.h>
23 #include <itkPolygonCell.h>
24 #include <itkQuadrilateralCell.h>
25 #include <itkTriangleCell.h>
26 //#include <itkDefaultDynamicMeshTraits.h>
27 #include <itkSphereMeshSource.h>
28 //#include <itkTransformMeshFilter.h>
29 //#include <itkTranslationTransform.h>
30 //#include <itkMinimumMaximumImageCalculator.h>
31 #include <itkAutomaticTopologyMeshSource.h>
32 #include <itkRegularSphereMeshSource.h>
33 #include <vnl/vnl_cross.h>
34 
35 #include <vtkActor.h>
36 #include <vtkCellArray.h>
37 #include <vtkCellData.h>
38 #include <vtkFloatArray.h>
39 #include <vtkPointData.h>
40 #include <vtkPoints.h>
41 #include <vtkPolyData.h>
42 #include <vtkProperty.h>
43 #include <vtkUnstructuredGrid.h>
44 
45 #include <mitkBaseGeometry.h>
46 #include <mitkSurface.h>
47 
48 template <typename MeshType>
50 {
51 public:
52  static inline double GetPointScalar(typename MeshType::PointDataContainer * /*pointData*/,
53  typename MeshType::PointIdentifier /*idx*/,
54  MeshType * /*mesh*/ = nullptr,
55  unsigned int /*type*/ = 0)
56  {
57  return (double)0.0;
58  };
59 
60  static inline double GetCellScalar(typename MeshType::CellDataContainer * /*cellData*/,
61  typename MeshType::CellIdentifier /*idx*/,
62  MeshType * /*mesh*/ = nullptr,
63  unsigned int /*type*/ = 0)
64  {
65  return (double)0.0;
66  };
67 };
68 
69 template <typename MeshType>
71 {
72 public:
73  static inline double GetPointScalar(typename MeshType::PointDataContainer *pointData,
74  typename MeshType::PointIdentifier idx,
75  MeshType * /*mesh*/ = nullptr,
76  unsigned int /*type*/ = 0)
77  {
78  return (double)pointData->GetElement(idx);
79  };
80 
81  static inline double GetCellScalar(typename MeshType::CellDataContainer *cellData,
82  typename MeshType::CellIdentifier idx,
83  MeshType * /*mesh*/ = nullptr,
84  unsigned int /*type*/ = 0)
85  {
86  return (double)cellData->GetElement(idx);
87  };
88 };
89 
90 template <typename MeshType>
91 class MeanCurvatureAccessor : public NullScalarAccessor<MeshType>
92 {
93 public:
94  static inline double GetPointScalar(typename MeshType::PointDataContainer * /*point*/,
95  typename MeshType::PointIdentifier idx,
96  MeshType *mesh,
97  unsigned int /*type*/ = 0)
98  {
99  typename MeshType::PixelType dis = 0;
100  mesh->GetPointData(idx, &dis);
101  return (double)dis;
102  };
103 };
104 
105 template <typename MeshType>
106 class SimplexMeshAccessor : public NullScalarAccessor<MeshType>
107 {
108 public:
109  static inline double GetPointScalar(typename MeshType::PointDataContainer * /*point*/,
110  typename MeshType::PointIdentifier idx,
111  MeshType *mesh,
112  unsigned int type = 0)
113  {
114  typename MeshType::GeometryMapPointer geometryData = mesh->GetGeometryData();
115 
116  if (type == 0)
117  {
118  double val = mesh->GetMeanCurvature(idx);
119  mesh->SetPointData(idx, val);
120  return val;
121  }
122  else if (type == 1)
123  {
124  double val = geometryData->GetElement(idx)->meanTension;
125  mesh->SetPointData(idx, val);
126  return val;
127  }
128  else if (type == 2)
129  {
130  double val = geometryData->GetElement(idx)->externalForce.GetNorm();
131  mesh->SetPointData(idx, val);
132  return val;
133  }
134  else if (type == 3)
135  return geometryData->GetElement(idx)->internalForce.GetNorm();
136  else if (type == 4)
137  return geometryData->GetElement(idx)->externalForce.GetNorm() * mesh->GetDistance(idx);
138  else if (type == 5)
139  {
140  typename MeshType::PixelType dis = 0;
141  mesh->GetPointData(idx, &dis);
142  return (double)dis;
143  }
144  else if (type == 6)
145  {
146  return (double)((geometryData->GetElement(idx))->allowSplitting);
147  }
148  else
149  return (double)0;
150  };
151 };
152 
160 template <typename MeshType, class ScalarAccessor = NullScalarAccessor<MeshType>>
161 class MeshUtil
162 {
186  template <class InsertImplementation>
187  class SwitchByCellType : public InsertImplementation
188  {
189  // typedef the itk cells we are interested in
190  typedef typename itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>
191  CellInterfaceType;
192  typedef itk::LineCell<CellInterfaceType> floatLineCell;
193  typedef itk::TriangleCell<CellInterfaceType> floatTriangleCell;
194  typedef itk::PolygonCell<CellInterfaceType> floatPolygonCell;
195  typedef itk::QuadrilateralCell<CellInterfaceType> floatQuadrilateralCell;
196  typedef itk::TetrahedronCell<CellInterfaceType> floatTetrahedronCell;
197  typedef itk::HexahedronCell<CellInterfaceType> floatHexahedronCell;
198  typedef typename CellInterfaceType::PointIdConstIterator PointIdIterator;
199 
200  public:
204  void Visit(unsigned long cellId, floatLineCell *t)
205  {
206  vtkIdType pts[2];
207  int i = 0;
208  unsigned long num = t->GetNumberOfVertices();
209  vtkIdType vtkCellId = -1;
210  if (num == 2)
211  { // useless because itk::LineCell always returns 2
212  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
213  pts[i++] = *it;
214  vtkCellId = this->InsertLine((vtkIdType *)pts);
215  }
216 
217  if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
218  {
219  this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
220  }
221  }
222 
226  void Visit(unsigned long cellId, floatPolygonCell *t)
227  {
228  vtkIdType pts[4096];
229  int i = 0;
230  unsigned long num = t->GetNumberOfVertices();
231  vtkIdType vtkCellId = -1;
232  if (num > 4096)
233  {
234  MITK_ERROR << "Problem in mitkMeshUtil: Polygon with more than maximum number of vertices encountered."
235  << std::endl;
236  }
237  else if (num > 3)
238  {
239  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
240  pts[i++] = *it;
241  vtkCellId = this->InsertPolygon(num, (vtkIdType *)pts);
242  }
243  else if (num == 3)
244  {
245  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
246  pts[i++] = *it;
247  vtkCellId = this->InsertTriangle((vtkIdType *)pts);
248  }
249  else if (num == 2)
250  {
251  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
252  pts[i++] = *it;
253  vtkCellId = this->InsertLine((vtkIdType *)pts);
254  }
255 
256  if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
257  {
258  this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
259  }
260  }
261 
265  void Visit(unsigned long cellId, floatTriangleCell *t)
266  {
267  vtkIdType pts[3];
268  int i = 0;
269  unsigned long num = t->GetNumberOfVertices();
270  vtkIdType vtkCellId = -1;
271  if (num == 3)
272  {
273  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
274  pts[i++] = *it;
275  vtkCellId = this->InsertTriangle((vtkIdType *)pts);
276  }
277  else if (num == 2)
278  {
279  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
280  pts[i++] = *it;
281  vtkCellId = this->InsertLine((vtkIdType *)pts);
282  }
283 
284  if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
285  {
286  this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
287  }
288  }
289 
293  void Visit(unsigned long cellId, floatQuadrilateralCell *t)
294  {
295  vtkIdType pts[4];
296  int i = 0;
297  unsigned long num = t->GetNumberOfVertices();
298  vtkIdType vtkCellId = -1;
299  if (num == 4)
300  {
301  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
302  {
303  if (i == 2)
304  pts[3] = *it;
305  else if (i == 3)
306  pts[2] = *it;
307  else
308  pts[i] = *it;
309  i++;
310  // pts[i++] = *it;
311  }
312  vtkCellId = this->InsertQuad((vtkIdType *)pts);
313  }
314  else if (num == 3)
315  {
316  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
317  pts[i++] = *it;
318  vtkCellId = this->InsertTriangle((vtkIdType *)pts);
319  }
320  else if (num == 2)
321  {
322  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
323  pts[i++] = *it;
324  vtkCellId = this->InsertLine((vtkIdType *)pts);
325  }
326 
327  if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
328  {
329  this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
330  }
331  }
332 
336  void Visit(unsigned long cellId, floatTetrahedronCell *t)
337  {
338  vtkIdType pts[4];
339  int i = 0;
340  unsigned long num = t->GetNumberOfVertices();
341  vtkIdType vtkCellId = -1;
342  if (num == 4)
343  {
344  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
345  pts[i++] = *it;
346  vtkCellId = this->InsertTetra((vtkIdType *)pts);
347  }
348  else if (num == 3)
349  {
350  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
351  pts[i++] = *it;
352  vtkCellId = this->InsertTriangle((vtkIdType *)pts);
353  }
354  else if (num == 2)
355  {
356  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
357  pts[i++] = *it;
358  vtkCellId = this->InsertLine((vtkIdType *)pts);
359  }
360 
361  if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
362  {
363  this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
364  }
365  }
366 
370  void Visit(unsigned long cellId, floatHexahedronCell *t)
371  {
372  vtkIdType pts[8];
373  int i = 0;
374  unsigned long num = t->GetNumberOfVertices();
375  vtkIdType vtkCellId = -1;
376  if (num == 8)
377  {
378  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
379  {
380  if (i == 2)
381  pts[i++] = *(it + 1);
382  else if (i == 3)
383  pts[i++] = *(it - 1);
384  else if (i == 6)
385  pts[i++] = *(it + 1);
386  else if (i == 7)
387  pts[i++] = *(it - 1);
388  else
389  pts[i++] = *it;
390  }
391  vtkCellId = this->InsertHexahedron((vtkIdType *)pts);
392  }
393  else if (num == 4)
394  {
395  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
396  pts[i++] = *it;
397  vtkCellId = this->InsertQuad((vtkIdType *)pts);
398  }
399  else if (num == 3)
400  {
401  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
402  pts[i++] = *it;
403  vtkCellId = this->InsertTriangle((vtkIdType *)pts);
404  }
405  else if (num == 2)
406  {
407  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
408  pts[i++] = *it;
409  vtkCellId = this->InsertLine((vtkIdType *)pts);
410  }
411 
412  if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
413  {
414  this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
415  }
416  }
417  };
418 
431  template <class InsertImplementation>
432  class ExactSwitchByCellType : public InsertImplementation
433  {
434  // typedef the itk cells we are interested in
435  typedef typename itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>
436  CellInterfaceType;
437  typedef itk::LineCell<CellInterfaceType> floatLineCell;
438  typedef itk::TriangleCell<CellInterfaceType> floatTriangleCell;
439  typedef itk::PolygonCell<CellInterfaceType> floatPolygonCell;
440  typedef itk::QuadrilateralCell<CellInterfaceType> floatQuadrilateralCell;
441  typedef itk::TetrahedronCell<CellInterfaceType> floatTetrahedronCell;
442  typedef itk::HexahedronCell<CellInterfaceType> floatHexahedronCell;
443  typedef typename CellInterfaceType::PointIdConstIterator PointIdIterator;
444 
445  public:
449  void Visit(unsigned long, floatLineCell *t)
450  {
451  unsigned long num = t->GetNumberOfVertices();
452  vtkIdType pts[2];
453  int i = 0;
454 
455  if (num == 2)
456  {
457  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
458  pts[i++] = *it;
459  this->InsertLine(pts);
460  }
461  }
462 
466  void Visit(unsigned long, floatPolygonCell *t)
467  {
468  vtkIdType pts[4096];
469  unsigned long num = t->GetNumberOfVertices();
470  if (num > 4096)
471  {
472  MITK_ERROR << "Problem in mitkMeshUtil: Polygon with more than maximum number of vertices encountered."
473  << std::endl;
474  }
475  int i = 0;
476 
477  if (num > 3)
478  {
479  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
480  pts[i++] = *it;
481  this->InsertPolygon(num, pts);
482  }
483  }
484 
488  void Visit(unsigned long, floatTriangleCell *t)
489  {
490  unsigned long num = t->GetNumberOfVertices();
491  vtkIdType pts[3];
492  int i = 0;
493 
494  if (num == 3)
495  {
496  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
497  pts[i++] = *it;
498  this->InsertTriangle(pts);
499  }
500  }
501 
505  void Visit(unsigned long, floatQuadrilateralCell *t)
506  {
507  unsigned long num = t->GetNumberOfVertices();
508  vtkIdType pts[4];
509  int i = 0;
510 
511  if (num == 4)
512  {
513  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
514  pts[i++] = *it;
515 
516  vtkIdType tmpId = pts[2];
517  pts[2] = pts[3];
518  pts[3] = tmpId;
519  this->InsertQuad(pts);
520  }
521  }
522 
526  void Visit(unsigned long, floatTetrahedronCell *t)
527  {
528  unsigned long num = t->GetNumberOfVertices();
529 
530  vtkIdType pts[4];
531  int i = 0;
532 
533  if (num == 4)
534  {
535  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
536  pts[i++] = *it;
537  this->InsertTetra(pts);
538  }
539  }
540 
544  void Visit(unsigned long, floatHexahedronCell *t)
545  {
546  unsigned long num = t->GetNumberOfVertices();
547  vtkIdType pts[8];
548  int i = 0;
549 
550  if (num == 8)
551  {
552  for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
553  pts[i++] = *it;
554 
555  vtkIdType tmp[8];
556  for (unsigned int i = 0; i < 8; i++)
557  tmp[i] = pts[i];
558  pts[2] = tmp[3];
559  pts[3] = tmp[2];
560  pts[6] = tmp[7];
561  pts[7] = tmp[6];
562  this->InsertHexahedron(pts);
563  }
564  }
565  };
566 
573  class SingleCellArrayInsertImplementation
574  {
575  vtkCellArray *m_Cells;
576  int *m_TypeArray;
577  // vtkIdType cellId;
578 
579  protected:
580  bool m_UseCellScalarAccessor;
581  vtkFloatArray *m_CellScalars;
582  typename MeshType::CellDataContainer::Pointer m_CellData;
583 
584  public:
585  SingleCellArrayInsertImplementation() : m_UseCellScalarAccessor(false) {}
588  void SetCellArray(vtkCellArray *cells) { m_Cells = cells; }
592  void SetTypeArray(int *i) { m_TypeArray = i; }
593  void SetUseCellScalarAccessor(bool flag) { m_UseCellScalarAccessor = flag; }
594  void SetCellScalars(vtkFloatArray *scalars) { m_CellScalars = scalars; }
595  vtkFloatArray *GetCellScalars() { return m_CellScalars; }
596  void SetMeshCellData(typename MeshType::CellDataContainer *data) { m_CellData = data; }
597  vtkIdType InsertLine(vtkIdType *pts)
598  {
599  vtkIdType cellId = m_Cells->InsertNextCell(2, pts);
600  m_TypeArray[cellId] = VTK_LINE;
601  return cellId;
602  }
603 
604  vtkIdType InsertTriangle(vtkIdType *pts)
605  {
606  vtkIdType cellId = m_Cells->InsertNextCell(3, pts);
607  m_TypeArray[cellId] = VTK_TRIANGLE;
608  return cellId;
609  }
610 
611  vtkIdType InsertPolygon(vtkIdType npts, vtkIdType *pts)
612  {
613  vtkIdType cellId = m_Cells->InsertNextCell(npts, pts);
614  m_TypeArray[cellId] = VTK_POLYGON;
615  return cellId;
616  }
617 
618  vtkIdType InsertQuad(vtkIdType *pts)
619  {
620  vtkIdType cellId = m_Cells->InsertNextCell(4, pts);
621  m_TypeArray[cellId] = VTK_QUAD;
622  return cellId;
623  }
624 
625  vtkIdType InsertTetra(vtkIdType *pts)
626  {
627  vtkIdType cellId = m_Cells->InsertNextCell(4, pts);
628  m_TypeArray[cellId] = VTK_TETRA;
629  return cellId;
630  }
631 
632  vtkIdType InsertHexahedron(vtkIdType *pts)
633  {
634  vtkIdType cellId = m_Cells->InsertNextCell(8, pts);
635  m_TypeArray[cellId] = VTK_HEXAHEDRON;
636  return cellId;
637  }
638  };
639 
645  class DistributeInsertImplementation
646  {
647  vtkCellArray *m_LineCells;
648  vtkCellArray *m_TriangleCells;
649  vtkCellArray *m_PolygonCells;
650  vtkCellArray *m_QuadCells;
651 
652  protected:
653  bool m_UseCellScalarAccessor;
654  vtkFloatArray *m_CellScalars;
655  typename MeshType::CellDataContainer::Pointer m_CellData;
656 
657  public:
658  DistributeInsertImplementation() : m_UseCellScalarAccessor(false) {}
661  void SetCellArrays(vtkCellArray *lines, vtkCellArray *triangles, vtkCellArray *polygons, vtkCellArray *quads)
662  {
663  m_LineCells = lines;
664  m_TriangleCells = triangles;
665  m_PolygonCells = polygons;
666  m_QuadCells = quads;
667  }
668 
669  vtkIdType InsertLine(vtkIdType *pts) { return m_LineCells->InsertNextCell(2, pts); }
670  vtkIdType InsertTriangle(vtkIdType *pts) { return m_TriangleCells->InsertNextCell(3, pts); }
671  vtkIdType InsertPolygon(vtkIdType npts, vtkIdType *pts) { return m_PolygonCells->InsertNextCell(npts, pts); }
672  vtkIdType InsertQuad(vtkIdType *pts) { return m_QuadCells->InsertNextCell(4, pts); }
673  vtkIdType InsertTetra(vtkIdType * /*pts*/) { return -1; } // ignored
674  vtkIdType InsertHexahedron(vtkIdType * /*pts*/) { return -1; } // ignored
675  };
676 
677  // typedef typename MeshType::CellType CellType;
678  // typedef typename itk::LineCell< CellType > LineType;
679  // typedef typename itk::PolygonCell< CellType > PolygonType;
680  // typedef typename itk::TriangleCell< CellType > TriangleType;
681 
682  typedef SwitchByCellType<SingleCellArrayInsertImplementation> SingleCellArrayUserVisitorType;
683  typedef SwitchByCellType<DistributeInsertImplementation> DistributeUserVisitorType;
684  typedef ExactSwitchByCellType<DistributeInsertImplementation> ExactUserVisitorType;
685 
686 public:
687  typedef itk::MatrixOffsetTransformBase<typename MeshType::CoordRepType, 3, 3> ITKTransformType;
688  typedef itk::MatrixOffsetTransformBase<mitk::ScalarType, 3, 3> MITKTransformType;
689 
694  static void ConvertTransformToItk(const MITKTransformType *mitkTransform, ITKTransformType *itkTransform)
695  {
696  typename MITKTransformType::MatrixType mitkM = mitkTransform->GetMatrix();
697  typename ITKTransformType::MatrixType itkM;
698 
699  typename MITKTransformType::OffsetType mitkO = mitkTransform->GetOffset();
700  typename ITKTransformType::OffsetType itkO;
701 
702  for (short i = 0; i < 3; ++i)
703  {
704  for (short j = 0; j < 3; ++j)
705  {
706  itkM[i][j] = (double)mitkM[i][j];
707  }
708  itkO[i] = (double)mitkO[i];
709  }
710 
711  itkTransform->SetMatrix(itkM);
712  itkTransform->SetOffset(itkO);
713  }
714 
718  static typename MeshType::Pointer MeshFromPolyData(vtkPolyData *poly,
719  mitk::BaseGeometry *geometryFrame = nullptr,
720  mitk::BaseGeometry *polyDataGeometryFrame = nullptr)
721  {
722  // Create a new mesh
723  typename MeshType::Pointer output = MeshType::New();
724  output->SetCellsAllocationMethod(MeshType::CellsAllocatedDynamicallyCellByCell);
725 
726  typedef typename MeshType::CellDataContainer MeshCellDataContainerType;
727 
728  output->SetCellData(MeshCellDataContainerType::New());
729 
730  // Get the points from vtk
731  vtkPoints *vtkpoints = poly->GetPoints();
732  const unsigned int numPoints = poly->GetNumberOfPoints();
733 
734  // Create a compatible point container for the mesh
735  // the mesh is created with a null points container
736  // MeshType::PointsContainer::Pointer points =
737  // MeshType::PointsContainer::New();
738  // // Resize the point container to be able to fit the vtk points
739  // points->Reserve(numPoints);
740  // // Set the point container on the mesh
741  // output->SetPoints(points);
742  double vtkpoint[3];
743  typename MeshType::PointType itkPhysicalPoint;
744  if (geometryFrame == nullptr)
745  {
746  if (polyDataGeometryFrame == nullptr)
747  {
748  for (unsigned int i = 0; i < numPoints; ++i)
749  {
750  vtkpoints->GetPoint(i, vtkpoint);
751  // MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl;
752  // typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i);
753  mitk::vtk2itk(vtkpoint, itkPhysicalPoint);
754  output->SetPoint(i, itkPhysicalPoint);
755  }
756  }
757  else
758  {
759  for (unsigned int i = 0; i < numPoints; ++i)
760  {
761  vtkpoints->GetPoint(i, vtkpoint);
762  // MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl;
763  // typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i);
764  mitk::Point3D mitkWorldPoint;
765  mitk::vtk2itk(vtkpoint, mitkWorldPoint);
766  polyDataGeometryFrame->IndexToWorld(mitkWorldPoint, mitkWorldPoint);
767  mitk::vtk2itk(mitkWorldPoint, itkPhysicalPoint);
768  output->SetPoint(i, itkPhysicalPoint);
769  }
770  }
771  }
772  else
773  {
774  mitk::Point3D mitkWorldPoint;
775  if (polyDataGeometryFrame == nullptr)
776  {
777  for (unsigned int i = 0; i < numPoints; ++i)
778  {
779  vtkpoints->GetPoint(i, vtkpoint);
780  // MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl;
781  // typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i);
782  mitk::vtk2itk(vtkpoint, mitkWorldPoint);
783  geometryFrame->WorldToItkPhysicalPoint(mitkWorldPoint, itkPhysicalPoint);
784  output->SetPoint(i, itkPhysicalPoint);
785  }
786  }
787  else
788  {
789  for (unsigned int i = 0; i < numPoints; ++i)
790  {
791  vtkpoints->GetPoint(i, vtkpoint);
792  // MITK_INFO << "next point: " << test[0]<< "," << test[1] << "," << test[2] << std::endl;
793  // typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i);
794  mitk::vtk2itk(vtkpoint, mitkWorldPoint);
795  polyDataGeometryFrame->IndexToWorld(mitkWorldPoint, mitkWorldPoint);
796  geometryFrame->WorldToItkPhysicalPoint(mitkWorldPoint, itkPhysicalPoint);
797  output->SetPoint(i, itkPhysicalPoint);
798  }
799  }
800  }
801 
802  vtkCellArray *vtkcells = poly->GetPolys();
803  // vtkCellArray* vtkcells = poly->GetStrips();
804  // MeshType::CellsContainerPointer cells = MeshType::CellsContainer::New();
805  // output->SetCells(cells);
806  // extract the cell id's from the vtkUnstructuredGrid
807  int numcells = vtkcells->GetNumberOfCells();
808  int *vtkCellTypes = new int[numcells];
809  int cellId = 0;
810  // poly ids start after verts and lines!
811  int cellIdOfs = poly->GetNumberOfVerts() + poly->GetNumberOfLines();
812  for (; cellId < numcells; cellId++)
813  {
814  vtkCellTypes[cellId] = poly->GetCellType(cellId + cellIdOfs);
815  }
816 
817  // cells->Reserve(numcells);
818  vtkIdType npts;
819  vtkIdType *pts;
820  cellId = 0;
821 
822  typedef typename MeshType::MeshTraits OMeshTraits;
823  typedef typename OMeshTraits::PixelType OPixelType;
824  typedef typename MeshType::CellTraits CellTraits;
825  typedef typename itk::CellInterface<OPixelType, CellTraits> CellInterfaceType;
826  typedef typename itk::TriangleCell<CellInterfaceType> TriCellType;
827  typedef typename TriCellType::CellAutoPointer TriCellPointer;
828 
829  TriCellPointer newCell;
830  output->GetCells()->Reserve(poly->GetNumberOfPolys() + poly->GetNumberOfStrips());
831  output->GetCellData()->Reserve(poly->GetNumberOfPolys() + poly->GetNumberOfStrips());
832 
833  for (vtkcells->InitTraversal(); vtkcells->GetNextCell(npts, pts); cellId++)
834  {
835  switch (vtkCellTypes[cellId])
836  {
837  case VTK_TRIANGLE:
838  {
839  if (npts != 3)
840  continue; // skip non-triangles;
841  itk::IdentifierType pointIds[3];
842  pointIds[0] = (unsigned long)pts[0];
843  pointIds[1] = (unsigned long)pts[1];
844  pointIds[2] = (unsigned long)pts[2];
845 
846  newCell.TakeOwnership(new TriCellType);
847  newCell->SetPointIds(pointIds); //(unsigned long*)pts);
848  output->SetCell(cellId, newCell);
849  output->SetCellData(cellId, (typename MeshType::PixelType)3);
850  break;
851  }
852 
853  case VTK_QUAD:
854  {
855  if (npts != 4)
856  continue; // skip non-quadrilateral
857  itk::IdentifierType pointIds[3];
858 
859  pointIds[0] = (unsigned long)pts[0];
860  pointIds[1] = (unsigned long)pts[1];
861  pointIds[2] = (unsigned long)pts[2];
862  newCell.TakeOwnership(new TriCellType);
863  newCell->SetPointIds(pointIds);
864  output->SetCell(cellId, newCell);
865  output->SetCellData(cellId, (typename MeshType::PixelType)3);
866  cellId++;
867 
868  pointIds[0] = (unsigned long)pts[2];
869  pointIds[1] = (unsigned long)pts[3];
870  pointIds[2] = (unsigned long)pts[0];
871  newCell.TakeOwnership(new TriCellType);
872  newCell->SetPointIds(pointIds);
873  output->SetCell(cellId, newCell);
874  output->SetCellData(cellId, (typename MeshType::PixelType)3);
875  break;
876  }
877 
878  case VTK_EMPTY_CELL:
879  {
880  if (npts != 3)
881  {
882  MITK_ERROR << "Only empty triangle cell supported by now..." << std::endl; // skip non-triangle empty cells;
883  continue;
884  }
885  itk::IdentifierType pointIds[3];
886  pointIds[0] = (unsigned long)pts[0];
887  pointIds[1] = (unsigned long)pts[1];
888  pointIds[2] = (unsigned long)pts[2];
889 
890  newCell.TakeOwnership(new TriCellType);
891  newCell->SetPointIds(pointIds);
892  output->SetCell(cellId, newCell);
893  output->SetCellData(cellId, (typename MeshType::PixelType)3);
894  break;
895  }
896 
897  // case VTK_VERTEX: // If need to implement use
898  // case VTK_POLY_VERTEX: // the poly->GetVerts() and
899  // case VTK_LINE: // poly->GetLines() routines
900  // case VTK_POLY_LINE: // outside of the switch..case.
901  case VTK_POLYGON:
902  case VTK_PIXEL:
903  {
904  if (npts != 4)
905  continue; // skip non-quadrilateral
906  itk::IdentifierType pointIds[3];
907  for (unsigned int idx = 0; idx <= 1; idx++)
908  {
909  pointIds[0] = (unsigned long)pts[idx];
910  pointIds[1] = (unsigned long)pts[idx + 1];
911  pointIds[2] = (unsigned long)pts[idx + 2];
912  newCell.TakeOwnership(new TriCellType);
913  newCell->SetPointIds(pointIds);
914  output->SetCell(cellId + idx, newCell);
915  output->SetCellData(cellId + idx, (typename MeshType::PixelType)3);
916  }
917  cellId++;
918  break;
919  }
920 
921  case VTK_TETRA:
922  case VTK_VOXEL:
923  case VTK_HEXAHEDRON:
924  case VTK_WEDGE:
925  case VTK_PYRAMID:
926  case VTK_PARAMETRIC_CURVE:
927  case VTK_PARAMETRIC_SURFACE:
928  default:
929  MITK_WARN << "Warning, unhandled cell type " << vtkCellTypes[cellId] << std::endl;
930  }
931  }
932 
933  if (poly->GetNumberOfStrips() != 0)
934  {
935  vtkcells = poly->GetStrips();
936  numcells = vtkcells->GetNumberOfCells();
937  vtkCellTypes = new int[numcells];
938  int stripId = 0;
939  // strip ids start after verts, lines and polys!
940  int stripIdOfs = poly->GetNumberOfVerts() + poly->GetNumberOfLines() + poly->GetNumberOfPolys();
941  for (; stripId < numcells; stripId++)
942  {
943  vtkCellTypes[stripId] = poly->GetCellType(stripId + stripIdOfs);
944  }
945  stripId = 0;
946 
947  vtkcells->InitTraversal();
948  while (vtkcells->GetNextCell(npts, pts))
949  {
950  if (vtkCellTypes[stripId] != VTK_TRIANGLE_STRIP)
951  {
952  MITK_ERROR << "Only triangle strips supported!" << std::endl;
953  continue;
954  }
955  stripId++;
956 
957  unsigned int numberOfTrianglesInStrip = npts - 2;
958  itk::IdentifierType pointIds[3];
959  pointIds[0] = (unsigned long)pts[0];
960  pointIds[1] = (unsigned long)pts[1];
961  pointIds[2] = (unsigned long)pts[2];
962 
963  for (unsigned int t = 0; t < numberOfTrianglesInStrip; t++)
964  {
965  newCell.TakeOwnership(new TriCellType);
966  newCell->SetPointIds(pointIds);
967  output->SetCell(cellId, newCell);
968  output->SetCellData(cellId, (typename MeshType::PixelType)3);
969  cellId++;
970  pointIds[0] = pointIds[1];
971  pointIds[1] = pointIds[2];
972  pointIds[2] = pts[t + 3];
973  }
974  }
975  }
976  // output->Print(std::cout);
977  output->BuildCellLinks();
978  delete[] vtkCellTypes;
979  return output;
980  }
981 
985  static typename MeshType::Pointer MeshFromSurface(mitk::Surface *surface, mitk::BaseGeometry *geometryFrame = nullptr)
986  {
987  if (surface == nullptr)
988  return nullptr;
989  return MeshFromPolyData(surface->GetVtkPolyData(), geometryFrame, surface->GetGeometry());
990  }
991 
995  static vtkUnstructuredGrid *MeshToUnstructuredGrid(MeshType *mesh,
996  bool usePointScalarAccessor = false,
997  bool useCellScalarAccessor = false,
998  unsigned int pointDataType = 0,
999  mitk::BaseGeometry *geometryFrame = nullptr)
1000  {
1004  typedef typename itk::CellInterfaceVisitorImplementation<typename MeshType::CellPixelType,
1005  typename MeshType::CellTraits,
1006  itk::LineCell<typename MeshType::CellType>,
1007  SingleCellArrayUserVisitorType>
1008  SingleCellArrayLineVisitor;
1009 
1013  typedef typename itk::CellInterfaceVisitorImplementation<typename MeshType::CellPixelType,
1014  typename MeshType::CellTraits,
1015  itk::PolygonCell<typename MeshType::CellType>,
1016  SingleCellArrayUserVisitorType>
1017  SingleCellArrayPolygonVisitor;
1018 
1022  typedef typename itk::CellInterfaceVisitorImplementation<
1023  typename MeshType::CellPixelType,
1024  typename MeshType::CellTraits,
1025  itk::TriangleCell<itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>>,
1026  SingleCellArrayUserVisitorType>
1027  SingleCellArrayTriangleVisitor;
1028 
1032  typedef typename itk::CellInterfaceVisitorImplementation<
1033  typename MeshType::CellPixelType,
1034  typename MeshType::CellTraits,
1035  itk::QuadrilateralCell<itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>>,
1036  SingleCellArrayUserVisitorType>
1037  SingleCellArrayQuadrilateralVisitor;
1038 
1042  typedef typename itk::CellInterfaceVisitorImplementation<
1043  typename MeshType::CellPixelType,
1044  typename MeshType::CellTraits,
1045  itk::TetrahedronCell<itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>>,
1046  SingleCellArrayUserVisitorType>
1047  SingleCellArrayTetrahedronVisitor;
1048 
1052  typedef typename itk::CellInterfaceVisitorImplementation<
1053  typename MeshType::CellPixelType,
1054  typename MeshType::CellTraits,
1055  itk::HexahedronCell<itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>>,
1056  SingleCellArrayUserVisitorType>
1057  SingleCellArrayHexahedronVisitor;
1058 
1059  // Get the number of points in the mesh
1060  int numPoints = mesh->GetNumberOfPoints();
1061  if (numPoints == 0)
1062  {
1063  // mesh->Print(std::cerr);
1064  MITK_FATAL << "no points in Grid " << std::endl;
1065  exit(-1);
1066  }
1067  // Create a vtkUnstructuredGrid
1068  vtkUnstructuredGrid *vgrid = vtkUnstructuredGrid::New();
1069  // Create the vtkPoints object and set the number of points
1070  vtkPoints *vpoints = vtkPoints::New(VTK_DOUBLE);
1071 
1072  vtkFloatArray *pointScalars = vtkFloatArray::New();
1073  vtkFloatArray *cellScalars = vtkFloatArray::New();
1074  pointScalars->SetNumberOfComponents(1);
1075  cellScalars->SetNumberOfComponents(1);
1076 
1077  typename MeshType::PointsContainer::Pointer points = mesh->GetPoints();
1078  typename MeshType::PointsContainer::Iterator i;
1079 
1080  // iterate over all the points in the itk mesh to find
1081  // the maximal index
1082  unsigned int maxIndex = 0;
1083  for (i = points->Begin(); i != points->End(); ++i)
1084  {
1085  if (maxIndex < i->Index())
1086  maxIndex = i->Index();
1087  }
1088 
1089  // initialize vtk-classes for points and scalars
1090  vpoints->SetNumberOfPoints(maxIndex + 1);
1091  pointScalars->SetNumberOfTuples(maxIndex + 1);
1092  cellScalars->SetNumberOfTuples(mesh->GetNumberOfCells());
1093 
1094  double vtkpoint[3];
1095  typename MeshType::PointType itkPhysicalPoint;
1096  if (geometryFrame == nullptr)
1097  {
1098  for (i = points->Begin(); i != points->End(); ++i)
1099  {
1100  // Get the point index from the point container iterator
1101  int idx = i->Index();
1102 
1103  itkPhysicalPoint = i->Value();
1104  mitk::itk2vtk(itkPhysicalPoint, vtkpoint);
1105  // Set the vtk point at the index with the the coord array from itk
1106  vpoints->SetPoint(idx, vtkpoint);
1107 
1108  if (usePointScalarAccessor)
1109  {
1110  pointScalars->InsertTuple1(
1111  idx, ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1112  }
1113  }
1114  }
1115  else
1116  {
1117  mitk::Point3D mitkWorldPoint;
1118  for (i = points->Begin(); i != points->End(); ++i)
1119  {
1120  // Get the point index from the point container iterator
1121  int idx = i->Index();
1122 
1123  itkPhysicalPoint = i->Value();
1124  geometryFrame->ItkPhysicalPointToWorld(itkPhysicalPoint, mitkWorldPoint);
1125  mitk::itk2vtk(mitkWorldPoint, vtkpoint);
1126  // Set the vtk point at the index with the the coord array from itk
1127  vpoints->SetPoint(idx, vtkpoint);
1128 
1129  if (usePointScalarAccessor)
1130  {
1131  pointScalars->InsertTuple1(
1132  idx, ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1133  }
1134  }
1135  }
1136  // Set the points on the vtk grid
1137  vgrid->SetPoints(vpoints);
1138  if (usePointScalarAccessor)
1139  vgrid->GetPointData()->SetScalars(pointScalars);
1140 
1141  // Now create the cells using the MultiVisitor
1142  // 1. Create a MultiVisitor
1143  typename MeshType::CellType::MultiVisitor::Pointer mv = MeshType::CellType::MultiVisitor::New();
1144  // 2. Create visitors
1145  typename SingleCellArrayLineVisitor::Pointer lv = SingleCellArrayLineVisitor::New();
1146  typename SingleCellArrayPolygonVisitor::Pointer pv = SingleCellArrayPolygonVisitor::New();
1147  typename SingleCellArrayTriangleVisitor::Pointer tv = SingleCellArrayTriangleVisitor::New();
1148  typename SingleCellArrayQuadrilateralVisitor::Pointer qv = SingleCellArrayQuadrilateralVisitor::New();
1149  typename SingleCellArrayTetrahedronVisitor::Pointer tetv = SingleCellArrayTetrahedronVisitor::New();
1150  typename SingleCellArrayHexahedronVisitor::Pointer hv = SingleCellArrayHexahedronVisitor::New();
1151  // 3. Set up the visitors
1152  // int vtkCellCount = 0; // running counter for current cell being inserted into vtk
1153  int numCells = mesh->GetNumberOfCells();
1154  int *types = new int[numCells]; // type array for vtk
1155  // create vtk cells and estimate the size
1156  vtkCellArray *cells = vtkCellArray::New();
1157  cells->Allocate(numCells);
1158  // Set the TypeArray CellCount and CellArray for the visitors
1159  lv->SetTypeArray(types);
1160  lv->SetCellArray(cells);
1161  pv->SetTypeArray(types);
1162  pv->SetCellArray(cells);
1163  tv->SetTypeArray(types);
1164  // tv->SetCellCounter(&vtkCellCount);
1165  tv->SetCellArray(cells);
1166  qv->SetTypeArray(types);
1167  // qv->SetCellCounter(&vtkCellCount);
1168  qv->SetCellArray(cells);
1169  tetv->SetTypeArray(types);
1170  tetv->SetCellArray(cells);
1171  hv->SetTypeArray(types);
1172  hv->SetCellArray(cells);
1173 
1174  if (useCellScalarAccessor)
1175  {
1176  lv->SetUseCellScalarAccessor(true);
1177  lv->SetCellScalars(cellScalars);
1178  lv->SetMeshCellData(mesh->GetCellData());
1179 
1180  pv->SetUseCellScalarAccessor(true);
1181  pv->SetCellScalars(cellScalars);
1182  pv->SetMeshCellData(mesh->GetCellData());
1183 
1184  tv->SetUseCellScalarAccessor(true);
1185  tv->SetCellScalars(cellScalars);
1186  tv->SetMeshCellData(mesh->GetCellData());
1187 
1188  qv->SetUseCellScalarAccessor(true);
1189  qv->SetCellScalars(cellScalars);
1190  qv->SetMeshCellData(mesh->GetCellData());
1191 
1192  tetv->SetUseCellScalarAccessor(true);
1193  tetv->SetCellScalars(cellScalars);
1194  tetv->SetMeshCellData(mesh->GetCellData());
1195 
1196  hv->SetUseCellScalarAccessor(true);
1197  hv->SetCellScalars(cellScalars);
1198  hv->SetMeshCellData(mesh->GetCellData());
1199  }
1200 
1201  // add the visitors to the multivisitor
1202  mv->AddVisitor(lv);
1203  mv->AddVisitor(pv);
1204  mv->AddVisitor(tv);
1205  mv->AddVisitor(qv);
1206  mv->AddVisitor(tetv);
1207  mv->AddVisitor(hv);
1208  // Now ask the mesh to accept the multivisitor which
1209  // will Call Visit for each cell in the mesh that matches the
1210  // cell types of the visitors added to the MultiVisitor
1211  mesh->Accept(mv);
1212  // Now set the cells on the vtk grid with the type array and cell array
1213 
1214  vgrid->SetCells(types, cells);
1215  vgrid->GetCellData()->SetScalars(cellScalars);
1216 
1217  // Clean up vtk objects (no vtkSmartPointer ... )
1218  cells->Delete();
1219  vpoints->Delete();
1220  delete[] types;
1221 
1222  pointScalars->Delete();
1223  cellScalars->Delete();
1224  // MITK_INFO << "meshToUnstructuredGrid end" << std::endl;
1225  return vgrid;
1226  }
1227 
1231  static vtkPolyData *MeshToPolyData(MeshType *mesh,
1232  bool onlyTriangles = false,
1233  bool useScalarAccessor = false,
1234  unsigned int pointDataType = 0,
1235  mitk::BaseGeometry *geometryFrame = nullptr,
1236  vtkPolyData *polydata = nullptr)
1237  {
1241  typedef typename itk::CellInterfaceVisitorImplementation<typename MeshType::CellPixelType,
1242  typename MeshType::CellTraits,
1243  itk::LineCell<typename MeshType::CellType>,
1244  DistributeUserVisitorType>
1245  DistributeLineVisitor;
1246 
1250  typedef typename itk::CellInterfaceVisitorImplementation<typename MeshType::CellPixelType,
1251  typename MeshType::CellTraits,
1252  itk::PolygonCell<typename MeshType::CellType>,
1253  DistributeUserVisitorType>
1254  DistributePolygonVisitor;
1255 
1259  typedef typename itk::CellInterfaceVisitorImplementation<
1260  typename MeshType::CellPixelType,
1261  typename MeshType::CellTraits,
1262  itk::TriangleCell<itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>>,
1263  DistributeUserVisitorType>
1264  DistributeTriangleVisitor;
1265 
1269  typedef typename itk::CellInterfaceVisitorImplementation<
1270  typename MeshType::CellPixelType,
1271  typename MeshType::CellTraits,
1272  itk::QuadrilateralCell<itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>>,
1273  DistributeUserVisitorType>
1274  DistributeQuadrilateralVisitor;
1275 
1279  typedef typename itk::CellInterfaceVisitorImplementation<
1280  typename MeshType::CellPixelType,
1281  typename MeshType::CellTraits,
1282  itk::TriangleCell<itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>>,
1283  ExactUserVisitorType>
1284  ExactTriangleVisitor;
1285 
1286  // Get the number of points in the mesh
1287  int numPoints = mesh->GetNumberOfPoints();
1288  if (numPoints == 0)
1289  {
1290  // mesh->Print(std::cerr);
1291  MITK_ERROR << "no points in Grid " << std::endl;
1292  }
1293  // Create a vtkPolyData
1294  if (polydata == nullptr)
1295  polydata = vtkPolyData::New();
1296  else
1297  polydata->Initialize();
1298 
1299  // Create the vtkPoints object and set the number of points
1300  vtkPoints *vpoints = vtkPoints::New(VTK_DOUBLE);
1301 
1302  vtkFloatArray *scalars = vtkFloatArray::New();
1303  scalars->SetNumberOfComponents(1);
1304 
1305  typename MeshType::PointsContainer::Pointer points = mesh->GetPoints();
1306  typename MeshType::PointsContainer::Iterator i;
1307 
1308  // iterate over all the points in the itk mesh to find
1309  // the maximal index
1310  unsigned int maxIndex = 0;
1311  for (i = points->Begin(); i != points->End(); ++i)
1312  {
1313  if (maxIndex < i->Index())
1314  maxIndex = i->Index();
1315  }
1316 
1317  // initialize vtk-classes for points and scalars
1318  vpoints->SetNumberOfPoints(maxIndex + 1);
1319  scalars->SetNumberOfTuples(maxIndex + 1);
1320 
1321  // iterate over all the points in the itk mesh filling in
1322  // the vtkPoints object as we go
1323 
1324  double vtkpoint[3];
1325  typename MeshType::PointType itkPhysicalPoint;
1326  if (geometryFrame == nullptr)
1327  {
1328  for (i = points->Begin(); i != points->End(); ++i)
1329  {
1330  // Get the point index from the point container iterator
1331  int idx = i->Index();
1332 
1333  itkPhysicalPoint = i->Value();
1334  mitk::itk2vtk(itkPhysicalPoint, vtkpoint);
1335  // Set the vtk point at the index with the the coord array from itk
1336  // itk returns a const pointer, but vtk is not const correct, so
1337  // we have to use a const cast to get rid of the const
1338  // vpoints->SetPoint(idx, const_cast<DATATYPE*>(i->Value().GetDataPointer()));
1339  vpoints->SetPoint(idx, vtkpoint);
1340 
1341  if (useScalarAccessor)
1342  {
1343  scalars->InsertTuple1(idx,
1344  ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1345  }
1346  }
1347  }
1348  else
1349  {
1350  mitk::Point3D mitkWorldPoint;
1351  for (i = points->Begin(); i != points->End(); ++i)
1352  {
1353  // Get the point index from the point container iterator
1354  int idx = i->Index();
1355 
1356  itkPhysicalPoint = i->Value();
1357  geometryFrame->ItkPhysicalPointToWorld(itkPhysicalPoint, mitkWorldPoint);
1358  mitk::itk2vtk(mitkWorldPoint, vtkpoint);
1359  // Set the vtk point at the index with the the coord array from itk
1360  // itk returns a const pointer, but vtk is not const correct, so
1361  // we have to use a const cast to get rid of the const
1362  // vpoints->SetPoint(idx, const_cast<DATATYPE*>(i->Value().GetDataPointer()));
1363  vpoints->SetPoint(idx, vtkpoint);
1364 
1365  if (useScalarAccessor)
1366  {
1367  scalars->InsertTuple1(idx,
1368  ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1369  }
1370  }
1371  }
1372 
1373  // Set the points on the vtk grid
1374  polydata->SetPoints(vpoints);
1375  if (useScalarAccessor)
1376  polydata->GetPointData()->SetScalars(scalars);
1377  polydata->GetPointData()->CopyAllOn();
1378 
1379  // Now create the cells using the MulitVisitor
1380  // 1. Create a MultiVisitor
1381  typedef typename MeshType::CellType::MultiVisitor MeshMV;
1382  typename MeshMV::Pointer mv = MeshMV::New();
1383 
1384  int numCells = mesh->GetNumberOfCells();
1385 
1386  if (onlyTriangles)
1387  {
1388  // create vtk cells and allocate
1389  vtkCellArray *trianglecells = vtkCellArray::New();
1390  trianglecells->Allocate(numCells);
1391 
1392  // 2. Create a triangle visitor and add it to the multivisitor
1393  typename ExactTriangleVisitor::Pointer tv = ExactTriangleVisitor::New();
1394  tv->SetCellArrays(nullptr, trianglecells, nullptr, nullptr);
1395  mv->AddVisitor(tv);
1396  // 3. Now ask the mesh to accept the multivisitor which
1397  // will Call Visit for each cell in the mesh that matches the
1398  // cell types of the visitors added to the MultiVisitor
1399  mesh->Accept(mv);
1400 
1401  // 4. Set the result into our vtkPolyData
1402  if (trianglecells->GetNumberOfCells() > 0)
1403  polydata->SetStrips(trianglecells);
1404 
1405  // 5. Clean up vtk objects (no vtkSmartPointer ... )
1406  trianglecells->Delete();
1407  }
1408  else
1409  {
1410  // create vtk cells and allocate
1411  vtkCellArray *linecells = vtkCellArray::New();
1412  vtkCellArray *trianglecells = vtkCellArray::New();
1413  vtkCellArray *polygoncells = vtkCellArray::New();
1414  linecells->Allocate(numCells);
1415  trianglecells->Allocate(numCells);
1416  polygoncells->Allocate(numCells);
1417 
1418  // 2. Create visitors
1419  typename DistributeLineVisitor::Pointer lv = DistributeLineVisitor::New();
1420  typename DistributePolygonVisitor::Pointer pv = DistributePolygonVisitor::New();
1421  typename DistributeTriangleVisitor::Pointer tv = DistributeTriangleVisitor::New();
1422  typename DistributeQuadrilateralVisitor::Pointer qv = DistributeQuadrilateralVisitor::New();
1423 
1424  lv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells);
1425  pv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells);
1426  tv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells);
1427  qv->SetCellArrays(linecells, trianglecells, polygoncells, polygoncells);
1428 
1429  // add the visitors to the multivisitor
1430  mv->AddVisitor(tv);
1431  mv->AddVisitor(lv);
1432  mv->AddVisitor(pv);
1433  mv->AddVisitor(qv);
1434  // 3. Now ask the mesh to accept the multivisitor which
1435  // will Call Visit for each cell in the mesh that matches the
1436  // cell types of the visitors added to the MultiVisitor
1437  mesh->Accept(mv);
1438 
1439  // 4. Set the result into our vtkPolyData
1440  if (linecells->GetNumberOfCells() > 0)
1441  polydata->SetLines(linecells);
1442  if (trianglecells->GetNumberOfCells() > 0)
1443  polydata->SetStrips(trianglecells);
1444  if (polygoncells->GetNumberOfCells() > 0)
1445  polydata->SetPolys(polygoncells);
1446 
1447  // 5. Clean up vtk objects (no vtkSmartPointer ... )
1448  linecells->Delete();
1449  trianglecells->Delete();
1450  polygoncells->Delete();
1451  }
1452  vpoints->Delete();
1453  scalars->Delete();
1454 
1455  // MITK_INFO << "meshToPolyData end" << std::endl;
1456  return polydata;
1457  }
1458 
1459  static typename MeshType::Pointer CreateRegularSphereMesh(typename MeshType::PointType center,
1460  typename MeshType::PointType::VectorType scale,
1461  int resolution)
1462  {
1463  typedef itk::RegularSphereMeshSource<MeshType> SphereSourceType;
1464  typename SphereSourceType::Pointer mySphereSource = SphereSourceType::New();
1465 
1466  mySphereSource->SetCenter(center);
1467  mySphereSource->SetScale(scale);
1468  mySphereSource->SetResolution(resolution);
1469  mySphereSource->Update();
1470 
1471  typename MeshType::Pointer resultMesh = mySphereSource->GetOutput();
1472  resultMesh->Register(); // necessary ????
1473  return resultMesh;
1474  }
1475 
1476  static typename MeshType::Pointer CreateSphereMesh(typename MeshType::PointType center,
1477  typename MeshType::PointType scale,
1478  int *resolution)
1479  {
1480  typedef typename itk::SphereMeshSource<MeshType> SphereSource;
1481 
1482  typename SphereSource::Pointer mySphereSource = SphereSource::New();
1483 
1484  mySphereSource->SetCenter(center);
1485  mySphereSource->SetScale(scale);
1486  mySphereSource->SetResolutionX(resolution[0]);
1487  mySphereSource->SetResolutionY(resolution[1]);
1488  mySphereSource->SetSquareness1(1);
1489  mySphereSource->SetSquareness2(1);
1490  mySphereSource->Update();
1491  mySphereSource->GetOutput();
1492 
1493  typename MeshType::Pointer resultMesh = mySphereSource->GetOutput();
1494  resultMesh->Register();
1495 
1496  return resultMesh;
1497  }
1498 
1499  // static typename MeshType::Pointer TranslateMesh(typename MeshType::PointType vec, MeshType* input)
1500  // {
1501  //
1502  // typename MeshType::Pointer output = MeshType::New();
1503  // {
1504  // output->SetPoints(input->GetPoints());
1505  // output->SetPointData(input->GetPointData());
1506  // output->SetCells(input->GetCells());
1507  // output->SetLastCellId( input->GetLastCellId() );
1508  // typename MeshType::GeometryMapIterator pointDataIterator = input->GetGeometryData()->Begin();
1509  // typename MeshType::GeometryMapIterator pointDataEnd = input->GetGeometryData()->End();
1510  //
1511  // typename MeshType::PointType inputPoint,outputPoint;
1512  //
1513  // while (pointDataIterator != pointDataEnd)
1514  // {
1515  // unsigned long pointId = pointDataIterator->Index();
1516  // itk::SimplexMeshGeometry* newGeometry = new itk::SimplexMeshGeometry();
1517  // itk::SimplexMeshGeometry* refGeometry = pointDataIterator->Value();
1518  //
1519  // input->GetPoint(pointId, &inputPoint );
1520  // outputPoint[0] = inputPoint[0] + vec[0];
1521  // outputPoint[1] = inputPoint[1] + vec[1];
1522  // outputPoint[2] = inputPoint[2] + vec[2];
1523  // output->SetPoint( pointId, outputPoint );
1524  //
1525  //
1526  // newGeometry->pos = outputPoint;
1527  // newGeometry->neighborIndices = refGeometry->neighborIndices;
1528  // newGeometry->meanCurvature = refGeometry->meanCurvature;
1529  // newGeometry->neighbors = refGeometry->neighbors;
1530  // newGeometry->oldPos = refGeometry->oldPos;
1531  // newGeometry->eps = refGeometry->eps;
1532  // newGeometry->referenceMetrics = refGeometry->referenceMetrics;
1533  // newGeometry->neighborSet = refGeometry->neighborSet;
1534  // newGeometry->distance = refGeometry->distance;
1535  // newGeometry->externalForce = refGeometry->externalForce;
1536  // newGeometry->internalForce = refGeometry->internalForce;
1537  // output->SetGeometryData(pointId, newGeometry);
1538  // pointDataIterator++;
1539  // }
1540  // }
1542  // return output;
1543  // }
1544 
1545  static typename MeshType::Pointer CreateRegularSphereMesh2(typename MeshType::PointType center,
1546  typename MeshType::PointType scale,
1547  int resolution)
1548  {
1549  typedef typename itk::AutomaticTopologyMeshSource<MeshType> MeshSourceType;
1550  typename MeshSourceType::Pointer mySphereSource = MeshSourceType::New();
1551 
1552  typename MeshType::PointType pnt0, pnt1, pnt2, pnt3, pnt4, pnt5, pnt6, pnt7, pnt8, pnt9, pnt10, pnt11;
1553  double c1 = 0.5 * (1.0 + sqrt(5.0));
1554  double c2 = 1.0;
1555  double len = sqrt(c1 * c1 + c2 * c2);
1556  c1 /= len;
1557  c2 /= len;
1558 
1559  pnt0[0] = center[0] - c1 * scale[0];
1560  pnt0[1] = center[1];
1561  pnt0[2] = center[2] + c2 * scale[2];
1562  pnt1[0] = center[0];
1563  pnt1[1] = center[1] + c2 * scale[1];
1564  pnt1[2] = center[2] - c1 * scale[2];
1565  pnt2[0] = center[0];
1566  pnt2[1] = center[1] + c2 * scale[1];
1567  pnt2[2] = center[2] + c1 * scale[2];
1568  pnt3[0] = center[0] + c1 * scale[0];
1569  pnt3[1] = center[1];
1570  pnt3[2] = center[2] - c2 * scale[2];
1571  pnt4[0] = center[0] - c2 * scale[0];
1572  pnt4[1] = center[1] - c1 * scale[1];
1573  pnt4[2] = center[2];
1574  pnt5[0] = center[0] - c2 * scale[0];
1575  pnt5[1] = center[1] + c1 * scale[1];
1576  pnt5[2] = center[2];
1577  pnt6[0] = center[0];
1578  pnt6[1] = center[1] - c2 * scale[1];
1579  pnt6[2] = center[2] + c1 * scale[2];
1580  pnt7[0] = center[0] + c2 * scale[0];
1581  pnt7[1] = center[1] + c1 * scale[1];
1582  pnt7[2] = center[2];
1583  pnt8[0] = center[0];
1584  pnt8[1] = center[1] - c2 * scale[1];
1585  pnt8[2] = center[2] - c1 * scale[2];
1586  pnt9[0] = center[0] + c1 * scale[0];
1587  pnt9[1] = center[1];
1588  pnt9[2] = center[2] + c2 * scale[2];
1589  pnt10[0] = center[0] + c2 * scale[0];
1590  pnt10[1] = center[1] - c1 * scale[1];
1591  pnt10[2] = center[2];
1592  pnt11[0] = center[0] - c1 * scale[0];
1593  pnt11[1] = center[1];
1594  pnt11[2] = center[2] - c2 * scale[2];
1595 
1596  addTriangle(mySphereSource, scale, pnt9, pnt2, pnt6, resolution);
1597  addTriangle(mySphereSource, scale, pnt1, pnt11, pnt5, resolution);
1598  addTriangle(mySphereSource, scale, pnt11, pnt1, pnt8, resolution);
1599  addTriangle(mySphereSource, scale, pnt0, pnt11, pnt4, resolution);
1600  addTriangle(mySphereSource, scale, pnt3, pnt1, pnt7, resolution);
1601  addTriangle(mySphereSource, scale, pnt3, pnt8, pnt1, resolution);
1602  addTriangle(mySphereSource, scale, pnt9, pnt3, pnt7, resolution);
1603  addTriangle(mySphereSource, scale, pnt0, pnt6, pnt2, resolution);
1604  addTriangle(mySphereSource, scale, pnt4, pnt10, pnt6, resolution);
1605  addTriangle(mySphereSource, scale, pnt1, pnt5, pnt7, resolution);
1606  addTriangle(mySphereSource, scale, pnt7, pnt5, pnt2, resolution);
1607  addTriangle(mySphereSource, scale, pnt8, pnt3, pnt10, resolution);
1608  addTriangle(mySphereSource, scale, pnt4, pnt11, pnt8, resolution);
1609  addTriangle(mySphereSource, scale, pnt9, pnt7, pnt2, resolution);
1610  addTriangle(mySphereSource, scale, pnt10, pnt9, pnt6, resolution);
1611  addTriangle(mySphereSource, scale, pnt0, pnt5, pnt11, resolution);
1612  addTriangle(mySphereSource, scale, pnt0, pnt2, pnt5, resolution);
1613  addTriangle(mySphereSource, scale, pnt8, pnt10, pnt4, resolution);
1614  addTriangle(mySphereSource, scale, pnt3, pnt9, pnt10, resolution);
1615  addTriangle(mySphereSource, scale, pnt6, pnt0, pnt4, resolution);
1616 
1617  return mySphereSource->GetOutput();
1618  }
1619 
1620 private:
1621  static void addTriangle(typename itk::AutomaticTopologyMeshSource<MeshType>::Pointer meshSource,
1622  typename MeshType::PointType scale,
1623  typename MeshType::PointType pnt0,
1624  typename MeshType::PointType pnt1,
1625  typename MeshType::PointType pnt2,
1626  int resolution)
1627  {
1628  if (resolution == 0)
1629  {
1630  // add triangle
1631  meshSource->AddTriangle(meshSource->AddPoint(pnt0), meshSource->AddPoint(pnt1), meshSource->AddPoint(pnt2));
1632  }
1633  else
1634  {
1635  vnl_vector_fixed<typename MeshType::CoordRepType, 3> v1, v2, res, pv;
1636  v1 = (pnt1 - pnt0).Get_vnl_vector();
1637  v2 = (pnt2 - pnt0).Get_vnl_vector();
1638  res = vnl_cross_3d(v1, v2);
1639  pv = pnt0.GetVectorFromOrigin().Get_vnl_vector();
1640  // double d = res[0]*pv[0] + res[1]*pv[1] + res[2]*pv[2];
1641 
1642  // subdivision
1643  typename MeshType::PointType pnt01, pnt12, pnt20;
1644  for (int d = 0; d < 3; d++)
1645  {
1646  pnt01[d] = (pnt0[d] + pnt1[d]) / 2.0;
1647  pnt12[d] = (pnt1[d] + pnt2[d]) / 2.0;
1648  pnt20[d] = (pnt2[d] + pnt0[d]) / 2.0;
1649  }
1650  // map new points to sphere
1651  double lenPnt01 = 0;
1652  for (int d = 0; d < 3; d++)
1653  lenPnt01 += pnt01[d] * pnt01[d];
1654  lenPnt01 = sqrt(lenPnt01);
1655  double lenPnt12 = 0;
1656  for (int d = 0; d < 3; d++)
1657  lenPnt12 += pnt12[d] * pnt12[d];
1658  lenPnt12 = sqrt(lenPnt12);
1659  double lenPnt20 = 0;
1660  for (int d = 0; d < 3; d++)
1661  lenPnt20 += pnt20[d] * pnt20[d];
1662  lenPnt20 = sqrt(lenPnt20);
1663  for (int d = 0; d < 3; d++)
1664  {
1665  pnt01[d] *= scale[d] / lenPnt01;
1666  pnt12[d] *= scale[d] / lenPnt12;
1667  pnt20[d] *= scale[d] / lenPnt20;
1668  }
1669  addTriangle(meshSource, scale, pnt0, pnt01, pnt20, resolution - 1);
1670  addTriangle(meshSource, scale, pnt01, pnt1, pnt12, resolution - 1);
1671  addTriangle(meshSource, scale, pnt20, pnt12, pnt2, resolution - 1);
1672  addTriangle(meshSource, scale, pnt01, pnt12, pnt20, resolution - 1);
1673  }
1674  }
1675 };
1676 
1677 #endif // MITKMESHUTIL_H_INCLUDED
#define MITK_FATAL
Definition: mitkLogMacros.h:21
Class for storing surfaces (vtkPolyData).
Definition: mitkSurface.h:28
virtual vtkPolyData * GetVtkPolyData(unsigned int t=0) const
#define MITK_ERROR
Definition: mitkLogMacros.h:20
static double GetCellScalar(typename MeshType::CellDataContainer *cellData, typename MeshType::CellIdentifier idx, MeshType *=nullptr, unsigned int=0)
Definition: mitkMeshUtil.h:81
static MeshType::Pointer CreateRegularSphereMesh(typename MeshType::PointType center, typename MeshType::PointType::VectorType scale, int resolution)
static MeshType::Pointer CreateSphereMesh(typename MeshType::PointType center, typename MeshType::PointType scale, int *resolution)
The class provides mehtods for ITK - VTK mesh conversion.
Definition: mitkMeshUtil.h:161
#define MITK_WARN
Definition: mitkLogMacros.h:19
static double GetPointScalar(typename MeshType::PointDataContainer *, typename MeshType::PointIdentifier idx, MeshType *mesh, unsigned int=0)
Definition: mitkMeshUtil.h:94
static vtkPolyData * MeshToPolyData(MeshType *mesh, bool onlyTriangles=false, bool useScalarAccessor=false, unsigned int pointDataType=0, mitk::BaseGeometry *geometryFrame=nullptr, vtkPolyData *polydata=nullptr)
static MeshType::Pointer CreateRegularSphereMesh2(typename MeshType::PointType center, typename MeshType::PointType scale, int resolution)
void vtk2itk(const Tin &in, Tout &out)
static double GetPointScalar(typename MeshType::PointDataContainer *pointData, typename MeshType::PointIdentifier idx, MeshType *=nullptr, unsigned int=0)
Definition: mitkMeshUtil.h:73
static void ConvertTransformToItk(const MITKTransformType *mitkTransform, ITKTransformType *itkTransform)
Definition: mitkMeshUtil.h:694
itk::MatrixOffsetTransformBase< typename MeshType::CoordRepType, 3, 3 > ITKTransformType
Definition: mitkMeshUtil.h:687
static double GetPointScalar(typename MeshType::PointDataContainer *, typename MeshType::PointIdentifier, MeshType *=nullptr, unsigned int=0)
Definition: mitkMeshUtil.h:52
void itk2vtk(const Tin &in, Tout &out)
static vtkUnstructuredGrid * MeshToUnstructuredGrid(MeshType *mesh, bool usePointScalarAccessor=false, bool useCellScalarAccessor=false, unsigned int pointDataType=0, mitk::BaseGeometry *geometryFrame=nullptr)
Definition: mitkMeshUtil.h:995
static MeshType::Pointer MeshFromSurface(mitk::Surface *surface, mitk::BaseGeometry *geometryFrame=nullptr)
Definition: mitkMeshUtil.h:985
itk::MatrixOffsetTransformBase< mitk::ScalarType, 3, 3 > MITKTransformType
Definition: mitkMeshUtil.h:688
static double GetPointScalar(typename MeshType::PointDataContainer *, typename MeshType::PointIdentifier idx, MeshType *mesh, unsigned int type=0)
Definition: mitkMeshUtil.h:109
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:138
static double GetCellScalar(typename MeshType::CellDataContainer *, typename MeshType::CellIdentifier, MeshType *=nullptr, unsigned int=0)
Definition: mitkMeshUtil.h:60
BaseGeometry Describes the geometry of a data object.
static MeshType::Pointer MeshFromPolyData(vtkPolyData *poly, mitk::BaseGeometry *geometryFrame=nullptr, mitk::BaseGeometry *polyDataGeometryFrame=nullptr)
Definition: mitkMeshUtil.h:718