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