17 #ifndef MITKMESHUTIL_H_INCLUDED
18 #define MITKMESHUTIL_H_INCLUDED
20 #if (_MSC_VER == 1200)
21 #error MeshUtils currently not supported for MS Visual C++ 6.0. Sorry.
25 #include <itkCellInterface.h>
26 #include <itkLineCell.h>
27 #include <itkPolygonCell.h>
28 #include <itkQuadrilateralCell.h>
29 #include <itkTriangleCell.h>
31 #include <itkSphereMeshSource.h>
35 #include <itkAutomaticTopologyMeshSource.h>
36 #include <itkRegularSphereMeshSource.h>
37 #include <vnl/vnl_cross.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>
52 template <
typename MeshType>
56 static inline double GetPointScalar(
typename MeshType::PointDataContainer * ,
57 typename MeshType::PointIdentifier ,
64 static inline double GetCellScalar(
typename MeshType::CellDataContainer * ,
65 typename MeshType::CellIdentifier ,
73 template <
typename MeshType>
77 static inline double GetPointScalar(
typename MeshType::PointDataContainer *pointData,
78 typename MeshType::PointIdentifier idx,
82 return (
double)pointData->GetElement(idx);
85 static inline double GetCellScalar(
typename MeshType::CellDataContainer *cellData,
86 typename MeshType::CellIdentifier idx,
90 return (
double)cellData->GetElement(idx);
94 template <
typename MeshType>
98 static inline double GetPointScalar(
typename MeshType::PointDataContainer * ,
99 typename MeshType::PointIdentifier idx,
104 mesh->GetPointData(idx, &dis);
109 template <
typename MeshType>
114 typename MeshType::PointIdentifier idx,
116 unsigned int type = 0)
118 typename MeshType::GeometryMapPointer geometryData = mesh->GetGeometryData();
122 double val = mesh->GetMeanCurvature(idx);
123 mesh->SetPointData(idx, val);
128 double val = geometryData->GetElement(idx)->meanTension;
129 mesh->SetPointData(idx, val);
134 double val = geometryData->GetElement(idx)->externalForce.GetNorm();
135 mesh->SetPointData(idx, val);
139 return geometryData->GetElement(idx)->internalForce.GetNorm();
141 return geometryData->GetElement(idx)->externalForce.GetNorm() * mesh->GetDistance(idx);
145 mesh->GetPointData(idx, &dis);
150 return (
double)((geometryData->GetElement(idx))->allowSplitting);
164 template <
typename MeshType,
class ScalarAccessor = NullScalarAccessor<MeshType>>
190 template <
class InsertImplementation>
191 class SwitchByCellType :
public InsertImplementation
194 typedef typename itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>
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;
208 void Visit(
unsigned long cellId, floatLineCell *t)
212 unsigned long num = t->GetNumberOfVertices();
213 vtkIdType vtkCellId = -1;
216 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
218 vtkCellId = this->InsertLine((vtkIdType *)pts);
221 if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
223 this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
230 void Visit(
unsigned long cellId, floatPolygonCell *t)
234 unsigned long num = t->GetNumberOfVertices();
235 vtkIdType vtkCellId = -1;
238 MITK_ERROR <<
"Problem in mitkMeshUtil: Polygon with more than maximum number of vertices encountered."
243 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
245 vtkCellId = this->InsertPolygon(num, (vtkIdType *)pts);
249 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
251 vtkCellId = this->InsertTriangle((vtkIdType *)pts);
255 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
257 vtkCellId = this->InsertLine((vtkIdType *)pts);
260 if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
262 this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
269 void Visit(
unsigned long cellId, floatTriangleCell *t)
273 unsigned long num = t->GetNumberOfVertices();
274 vtkIdType vtkCellId = -1;
277 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
279 vtkCellId = this->InsertTriangle((vtkIdType *)pts);
283 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
285 vtkCellId = this->InsertLine((vtkIdType *)pts);
288 if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
290 this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
297 void Visit(
unsigned long cellId, floatQuadrilateralCell *t)
301 unsigned long num = t->GetNumberOfVertices();
302 vtkIdType vtkCellId = -1;
305 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
316 vtkCellId = this->InsertQuad((vtkIdType *)pts);
320 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
322 vtkCellId = this->InsertTriangle((vtkIdType *)pts);
326 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
328 vtkCellId = this->InsertLine((vtkIdType *)pts);
331 if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
333 this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
340 void Visit(
unsigned long cellId, floatTetrahedronCell *t)
344 unsigned long num = t->GetNumberOfVertices();
345 vtkIdType vtkCellId = -1;
348 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
350 vtkCellId = this->InsertTetra((vtkIdType *)pts);
354 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
356 vtkCellId = this->InsertTriangle((vtkIdType *)pts);
360 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
362 vtkCellId = this->InsertLine((vtkIdType *)pts);
365 if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
367 this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
374 void Visit(
unsigned long cellId, floatHexahedronCell *t)
378 unsigned long num = t->GetNumberOfVertices();
379 vtkIdType vtkCellId = -1;
382 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
385 pts[i++] = *(it + 1);
387 pts[i++] = *(it - 1);
389 pts[i++] = *(it + 1);
391 pts[i++] = *(it - 1);
395 vtkCellId = this->InsertHexahedron((vtkIdType *)pts);
399 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
401 vtkCellId = this->InsertQuad((vtkIdType *)pts);
405 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
407 vtkCellId = this->InsertTriangle((vtkIdType *)pts);
411 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
413 vtkCellId = this->InsertLine((vtkIdType *)pts);
416 if (this->m_UseCellScalarAccessor && vtkCellId >= 0)
418 this->m_CellScalars->InsertTuple1(vtkCellId, ScalarAccessor::GetCellScalar(this->m_CellData, cellId));
435 template <
class InsertImplementation>
436 class ExactSwitchByCellType :
public InsertImplementation
439 typedef typename itk::CellInterface<typename MeshType::CellPixelType, typename MeshType::CellTraits>
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;
453 void Visit(
unsigned long, floatLineCell *t)
455 unsigned long num = t->GetNumberOfVertices();
461 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
463 this->InsertLine(pts);
470 void Visit(
unsigned long, floatPolygonCell *t)
473 unsigned long num = t->GetNumberOfVertices();
476 MITK_ERROR <<
"Problem in mitkMeshUtil: Polygon with more than maximum number of vertices encountered."
483 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
485 this->InsertPolygon(num, pts);
492 void Visit(
unsigned long, floatTriangleCell *t)
494 unsigned long num = t->GetNumberOfVertices();
500 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
502 this->InsertTriangle(pts);
509 void Visit(
unsigned long, floatQuadrilateralCell *t)
511 unsigned long num = t->GetNumberOfVertices();
517 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
520 vtkIdType tmpId = pts[2];
523 this->InsertQuad(pts);
530 void Visit(
unsigned long, floatTetrahedronCell *t)
532 unsigned long num = t->GetNumberOfVertices();
539 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
541 this->InsertTetra(pts);
548 void Visit(
unsigned long, floatHexahedronCell *t)
550 unsigned long num = t->GetNumberOfVertices();
556 for (PointIdIterator it = t->PointIdsBegin(); it != t->PointIdsEnd(); it++)
560 for (
unsigned int i = 0; i < 8; i++)
566 this->InsertHexahedron(pts);
577 class SingleCellArrayInsertImplementation
579 vtkCellArray *m_Cells;
584 bool m_UseCellScalarAccessor;
585 vtkFloatArray *m_CellScalars;
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)
603 vtkIdType cellId = m_Cells->InsertNextCell(2, pts);
604 m_TypeArray[cellId] = VTK_LINE;
608 vtkIdType InsertTriangle(vtkIdType *pts)
610 vtkIdType cellId = m_Cells->InsertNextCell(3, pts);
611 m_TypeArray[cellId] = VTK_TRIANGLE;
615 vtkIdType InsertPolygon(vtkIdType npts, vtkIdType *pts)
617 vtkIdType cellId = m_Cells->InsertNextCell(npts, pts);
618 m_TypeArray[cellId] = VTK_POLYGON;
622 vtkIdType InsertQuad(vtkIdType *pts)
624 vtkIdType cellId = m_Cells->InsertNextCell(4, pts);
625 m_TypeArray[cellId] = VTK_QUAD;
629 vtkIdType InsertTetra(vtkIdType *pts)
631 vtkIdType cellId = m_Cells->InsertNextCell(4, pts);
632 m_TypeArray[cellId] = VTK_TETRA;
636 vtkIdType InsertHexahedron(vtkIdType *pts)
638 vtkIdType cellId = m_Cells->InsertNextCell(8, pts);
639 m_TypeArray[cellId] = VTK_HEXAHEDRON;
649 class DistributeInsertImplementation
651 vtkCellArray *m_LineCells;
652 vtkCellArray *m_TriangleCells;
653 vtkCellArray *m_PolygonCells;
654 vtkCellArray *m_QuadCells;
657 bool m_UseCellScalarAccessor;
658 vtkFloatArray *m_CellScalars;
662 DistributeInsertImplementation() : m_UseCellScalarAccessor(
false) {}
665 void SetCellArrays(vtkCellArray *lines, vtkCellArray *triangles, vtkCellArray *polygons, vtkCellArray *quads)
668 m_TriangleCells = triangles;
669 m_PolygonCells = polygons;
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 * ) {
return -1; }
678 vtkIdType InsertHexahedron(vtkIdType * ) {
return -1; }
686 typedef SwitchByCellType<SingleCellArrayInsertImplementation> SingleCellArrayUserVisitorType;
687 typedef SwitchByCellType<DistributeInsertImplementation> DistributeUserVisitorType;
688 typedef ExactSwitchByCellType<DistributeInsertImplementation> ExactUserVisitorType;
691 typedef itk::MatrixOffsetTransformBase<typename MeshType::CoordRepType, 3, 3>
ITKTransformType;
700 typename MITKTransformType::MatrixType mitkM = mitkTransform->GetMatrix();
701 typename ITKTransformType::MatrixType itkM;
703 typename MITKTransformType::OffsetType mitkO = mitkTransform->GetOffset();
704 typename ITKTransformType::OffsetType itkO;
706 for (
short i = 0; i < 3; ++i)
708 for (
short j = 0; j < 3; ++j)
710 itkM[i][j] = (double)mitkM[i][j];
712 itkO[i] = (double)mitkO[i];
715 itkTransform->SetMatrix(itkM);
716 itkTransform->SetOffset(itkO);
728 output->SetCellsAllocationMethod(MeshType::CellsAllocatedDynamicallyCellByCell);
730 typedef typename MeshType::CellDataContainer MeshCellDataContainerType;
735 vtkPoints *vtkpoints = poly->GetPoints();
736 const unsigned int numPoints = poly->GetNumberOfPoints();
748 if (geometryFrame == NULL)
750 if (polyDataGeometryFrame == NULL)
752 for (
unsigned int i = 0; i < numPoints; ++i)
754 vtkpoints->GetPoint(i, vtkpoint);
758 output->SetPoint(i, itkPhysicalPoint);
763 for (
unsigned int i = 0; i < numPoints; ++i)
765 vtkpoints->GetPoint(i, vtkpoint);
770 polyDataGeometryFrame->IndexToWorld(mitkWorldPoint, mitkWorldPoint);
772 output->SetPoint(i, itkPhysicalPoint);
779 if (polyDataGeometryFrame == NULL)
781 for (
unsigned int i = 0; i < numPoints; ++i)
783 vtkpoints->GetPoint(i, vtkpoint);
787 geometryFrame->WorldToItkPhysicalPoint(mitkWorldPoint, itkPhysicalPoint);
788 output->SetPoint(i, itkPhysicalPoint);
793 for (
unsigned int i = 0; i < numPoints; ++i)
795 vtkpoints->GetPoint(i, vtkpoint);
799 polyDataGeometryFrame->IndexToWorld(mitkWorldPoint, mitkWorldPoint);
800 geometryFrame->WorldToItkPhysicalPoint(mitkWorldPoint, itkPhysicalPoint);
801 output->SetPoint(i, itkPhysicalPoint);
806 vtkCellArray *vtkcells = poly->GetPolys();
811 int numcells = vtkcells->GetNumberOfCells();
812 int *vtkCellTypes =
new int[numcells];
815 int cellIdOfs = poly->GetNumberOfVerts() + poly->GetNumberOfLines();
816 for (; cellId < numcells; cellId++)
818 vtkCellTypes[cellId] = poly->GetCellType(cellId + cellIdOfs);
826 typedef typename MeshType::MeshTraits OMeshTraits;
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;
833 TriCellPointer newCell;
834 output->GetCells()->Reserve(poly->GetNumberOfPolys() + poly->GetNumberOfStrips());
835 output->GetCellData()->Reserve(poly->GetNumberOfPolys() + poly->GetNumberOfStrips());
837 for (vtkcells->InitTraversal(); vtkcells->GetNextCell(npts, pts); cellId++)
839 switch (vtkCellTypes[cellId])
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];
850 newCell.TakeOwnership(
new TriCellType);
851 newCell->SetPointIds(pointIds);
852 output->SetCell(cellId, newCell);
861 unsigned long pointIds[3];
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);
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);
886 MITK_ERROR <<
"Only empty triangle cell supported by now..." << std::endl;
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];
894 newCell.TakeOwnership(
new TriCellType);
895 newCell->SetPointIds(pointIds);
896 output->SetCell(cellId, newCell);
910 unsigned long pointIds[3];
911 for (
unsigned int idx = 0; idx <= 1; idx++)
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);
930 case VTK_PARAMETRIC_CURVE:
931 case VTK_PARAMETRIC_SURFACE:
933 MITK_WARN <<
"Warning, unhandled cell type " << vtkCellTypes[cellId] << std::endl;
937 if (poly->GetNumberOfStrips() != 0)
939 vtkcells = poly->GetStrips();
940 numcells = vtkcells->GetNumberOfCells();
941 vtkCellTypes =
new int[numcells];
944 int stripIdOfs = poly->GetNumberOfVerts() + poly->GetNumberOfLines() + poly->GetNumberOfPolys();
945 for (; stripId < numcells; stripId++)
947 vtkCellTypes[stripId] = poly->GetCellType(stripId + stripIdOfs);
951 vtkcells->InitTraversal();
952 while (vtkcells->GetNextCell(npts, pts))
954 if (vtkCellTypes[stripId] != VTK_TRIANGLE_STRIP)
956 MITK_ERROR <<
"Only triangle strips supported!" << std::endl;
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];
967 for (
unsigned int t = 0; t < numberOfTrianglesInStrip; t++)
969 newCell.TakeOwnership(
new TriCellType);
970 newCell->SetPointIds(pointIds);
971 output->SetCell(cellId, newCell);
974 pointIds[0] = pointIds[1];
975 pointIds[1] = pointIds[2];
976 pointIds[2] = pts[t + 3];
981 output->BuildCellLinks();
982 delete[] vtkCellTypes;
1000 bool usePointScalarAccessor =
false,
1001 bool useCellScalarAccessor =
false,
1002 unsigned int pointDataType = 0,
1008 typedef typename itk::CellInterfaceVisitorImplementation<
typename MeshType::CellPixelType,
1009 typename MeshType::CellTraits,
1010 itk::LineCell<typename MeshType::CellType>,
1011 SingleCellArrayUserVisitorType>
1012 SingleCellArrayLineVisitor;
1017 typedef typename itk::CellInterfaceVisitorImplementation<
typename MeshType::CellPixelType,
1018 typename MeshType::CellTraits,
1019 itk::PolygonCell<typename MeshType::CellType>,
1020 SingleCellArrayUserVisitorType>
1021 SingleCellArrayPolygonVisitor;
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;
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;
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;
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;
1064 int numPoints = mesh->GetNumberOfPoints();
1068 MITK_FATAL <<
"no points in Grid " << std::endl;
1078 pointScalars->SetNumberOfComponents(1);
1079 cellScalars->SetNumberOfComponents(1);
1082 typename MeshType::PointsContainer::Iterator i;
1086 unsigned int maxIndex = 0;
1087 for (i = points->Begin(); i != points->End(); ++i)
1089 if (maxIndex < i->Index())
1090 maxIndex = i->Index();
1094 vpoints->SetNumberOfPoints(maxIndex + 1);
1095 pointScalars->SetNumberOfTuples(maxIndex + 1);
1096 cellScalars->SetNumberOfTuples(mesh->GetNumberOfCells());
1100 if (geometryFrame == 0)
1102 for (i = points->Begin(); i != points->End(); ++i)
1105 int idx = i->Index();
1107 itkPhysicalPoint = i->Value();
1110 vpoints->SetPoint(idx, vtkpoint);
1112 if (usePointScalarAccessor)
1114 pointScalars->InsertTuple1(
1115 idx, ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1122 for (i = points->Begin(); i != points->End(); ++i)
1125 int idx = i->Index();
1127 itkPhysicalPoint = i->Value();
1128 geometryFrame->ItkPhysicalPointToWorld(itkPhysicalPoint, mitkWorldPoint);
1131 vpoints->SetPoint(idx, vtkpoint);
1133 if (usePointScalarAccessor)
1135 pointScalars->InsertTuple1(
1136 idx, ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1141 vgrid->SetPoints(vpoints);
1142 if (usePointScalarAccessor)
1143 vgrid->GetPointData()->SetScalars(pointScalars);
1157 int numCells = mesh->GetNumberOfCells();
1158 int *types =
new int[numCells];
1161 cells->Allocate(numCells);
1163 lv->SetTypeArray(types);
1164 lv->SetCellArray(cells);
1165 pv->SetTypeArray(types);
1166 pv->SetCellArray(cells);
1167 tv->SetTypeArray(types);
1169 tv->SetCellArray(cells);
1170 qv->SetTypeArray(types);
1172 qv->SetCellArray(cells);
1173 tetv->SetTypeArray(types);
1174 tetv->SetCellArray(cells);
1175 hv->SetTypeArray(types);
1176 hv->SetCellArray(cells);
1178 if (useCellScalarAccessor)
1180 lv->SetUseCellScalarAccessor(
true);
1181 lv->SetCellScalars(cellScalars);
1182 lv->SetMeshCellData(mesh->GetCellData());
1184 pv->SetUseCellScalarAccessor(
true);
1185 pv->SetCellScalars(cellScalars);
1186 pv->SetMeshCellData(mesh->GetCellData());
1188 tv->SetUseCellScalarAccessor(
true);
1189 tv->SetCellScalars(cellScalars);
1190 tv->SetMeshCellData(mesh->GetCellData());
1192 qv->SetUseCellScalarAccessor(
true);
1193 qv->SetCellScalars(cellScalars);
1194 qv->SetMeshCellData(mesh->GetCellData());
1196 tetv->SetUseCellScalarAccessor(
true);
1197 tetv->SetCellScalars(cellScalars);
1198 tetv->SetMeshCellData(mesh->GetCellData());
1200 hv->SetUseCellScalarAccessor(
true);
1201 hv->SetCellScalars(cellScalars);
1202 hv->SetMeshCellData(mesh->GetCellData());
1210 mv->AddVisitor(tetv);
1218 vgrid->SetCells(types, cells);
1219 vgrid->GetCellData()->SetScalars(cellScalars);
1226 pointScalars->Delete();
1227 cellScalars->Delete();
1236 bool onlyTriangles =
false,
1237 bool useScalarAccessor =
false,
1238 unsigned int pointDataType = 0,
1240 vtkPolyData *polydata = NULL)
1245 typedef typename itk::CellInterfaceVisitorImplementation<
typename MeshType::CellPixelType,
1246 typename MeshType::CellTraits,
1247 itk::LineCell<typename MeshType::CellType>,
1248 DistributeUserVisitorType>
1249 DistributeLineVisitor;
1254 typedef typename itk::CellInterfaceVisitorImplementation<
typename MeshType::CellPixelType,
1255 typename MeshType::CellTraits,
1256 itk::PolygonCell<typename MeshType::CellType>,
1257 DistributeUserVisitorType>
1258 DistributePolygonVisitor;
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;
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;
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;
1291 int numPoints = mesh->GetNumberOfPoints();
1295 MITK_ERROR <<
"no points in Grid " << std::endl;
1298 if (polydata == NULL)
1301 polydata->Initialize();
1307 scalars->SetNumberOfComponents(1);
1310 typename MeshType::PointsContainer::Iterator i;
1314 unsigned int maxIndex = 0;
1315 for (i = points->Begin(); i != points->End(); ++i)
1317 if (maxIndex < i->Index())
1318 maxIndex = i->Index();
1322 vpoints->SetNumberOfPoints(maxIndex + 1);
1323 scalars->SetNumberOfTuples(maxIndex + 1);
1330 if (geometryFrame == NULL)
1332 for (i = points->Begin(); i != points->End(); ++i)
1335 int idx = i->Index();
1337 itkPhysicalPoint = i->Value();
1343 vpoints->SetPoint(idx, vtkpoint);
1345 if (useScalarAccessor)
1347 scalars->InsertTuple1(idx,
1348 ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1355 for (i = points->Begin(); i != points->End(); ++i)
1358 int idx = i->Index();
1360 itkPhysicalPoint = i->Value();
1361 geometryFrame->ItkPhysicalPointToWorld(itkPhysicalPoint, mitkWorldPoint);
1367 vpoints->SetPoint(idx, vtkpoint);
1369 if (useScalarAccessor)
1371 scalars->InsertTuple1(idx,
1372 ScalarAccessor::GetPointScalar(mesh->GetPointData(), i->Index(), mesh, pointDataType));
1378 polydata->SetPoints(vpoints);
1379 if (useScalarAccessor)
1380 polydata->GetPointData()->SetScalars(scalars);
1381 polydata->GetPointData()->CopyAllOn();
1385 typedef typename MeshType::CellType::MultiVisitor MeshMV;
1388 int numCells = mesh->GetNumberOfCells();
1394 trianglecells->Allocate(numCells);
1398 tv->SetCellArrays(NULL, trianglecells, NULL, NULL);
1406 if (trianglecells->GetNumberOfCells() > 0)
1407 polydata->SetStrips(trianglecells);
1410 trianglecells->Delete();
1418 linecells->Allocate(numCells);
1419 trianglecells->Allocate(numCells);
1420 polygoncells->Allocate(numCells);
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);
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);
1452 linecells->Delete();
1453 trianglecells->Delete();
1454 polygoncells->Delete();
1467 typedef itk::RegularSphereMeshSource<MeshType> SphereSourceType;
1470 mySphereSource->SetCenter(center);
1471 mySphereSource->SetScale(scale);
1472 mySphereSource->SetResolution(resolution);
1473 mySphereSource->Update();
1476 resultMesh->Register();
1484 typedef typename itk::SphereMeshSource<MeshType> SphereSource;
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();
1498 resultMesh->Register();
1553 typedef typename itk::AutomaticTopologyMeshSource<MeshType> MeshSourceType;
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));
1559 double len = sqrt(c1 * c1 + c2 * c2);
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];
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);
1621 return mySphereSource->GetOutput();
1632 if (resolution == 0)
1635 meshSource->AddTriangle(meshSource->AddPoint(pnt0), meshSource->AddPoint(pnt1), meshSource->AddPoint(pnt2));
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();
1648 for (
int d = 0; d < 3; d++)
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;
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++)
1669 pnt01[d] *= scale[d] / lenPnt01;
1670 pnt12[d] *= scale[d] / lenPnt12;
1671 pnt20[d] *= scale[d] / lenPnt20;
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);
1681 #endif // MITKMESHUTIL_H_INCLUDED
Class for storing surfaces (vtkPolyData).
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)
static double GetPointScalar(typename MeshType::PointDataContainer *pointData, typename MeshType::PointIdentifier idx, MeshType *=NULL, unsigned int=0)
virtual vtkPolyData * GetVtkPolyData(unsigned int t=0) const
static double GetCellScalar(typename MeshType::CellDataContainer *, typename MeshType::CellIdentifier, MeshType *=NULL, unsigned int=0)
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)
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)
The class provides mehtods for ITK - VTK mesh conversion.
itk::Vector< float, 3 > VectorType
static double GetPointScalar(typename MeshType::PointDataContainer *, typename MeshType::PointIdentifier idx, MeshType *mesh, unsigned int=0)
static MeshType::Pointer MeshFromSurface(mitk::Surface *surface, mitk::BaseGeometry *geometryFrame=NULL)
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)
static void ConvertTransformToItk(const MITKTransformType *mitkTransform, ITKTransformType *itkTransform)
itk::MatrixOffsetTransformBase< typename MeshType::CoordRepType, 3, 3 > ITKTransformType
void itk2vtk(const Tin &in, Tout &out)
itk::MatrixOffsetTransformBase< mitk::ScalarType, 3, 3 > MITKTransformType
static double GetPointScalar(typename MeshType::PointDataContainer *, typename MeshType::PointIdentifier idx, MeshType *mesh, unsigned int type=0)
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
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.