Medical Imaging Interaction Toolkit  2018.4.99-c4b6bb11
Medical Imaging Interaction Toolkit
mitkACVD.cpp
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 #include "mitkACVD.h"
14 #include <mitkExceptionMacro.h>
15 #include <vtkIdList.h>
16 #include <vtkIntArray.h>
17 #include <vtkIsotropicDiscreteRemeshing.h>
18 #include <vtkMultiThreader.h>
19 #include <vtkPolyData.h>
20 #include <vtkPolyDataNormals.h>
21 #include <vtkSmartPointer.h>
22 #include <vtkSurface.h>
23 
24 struct ClustersQuadrics
25 {
26  explicit ClustersQuadrics(int size) : Elements(new double *[size]), Size(size)
27  {
28  for (int i = 0; i < size; ++i)
29  {
30  Elements[i] = new double[9];
31 
32  for (int j = 0; j < 9; ++j)
33  Elements[i][j] = 0.0;
34  }
35  }
36 
37  ~ClustersQuadrics()
38  {
39  for (int i = 0; i < Size; ++i)
40  delete[] Elements[i];
41 
42  delete Elements;
43  }
44 
45  double **Elements;
46  int Size;
47 
48 private:
49  ClustersQuadrics(const ClustersQuadrics &);
50  ClustersQuadrics &operator=(const ClustersQuadrics &);
51 };
52 
53 static void ValidateSurface(mitk::Surface::ConstPointer surface, unsigned int t)
54 {
55  if (surface.IsNull())
56  mitkThrow() << "Input surface is nullptr!";
57 
58  if (t >= surface->GetSizeOfPolyDataSeries())
59  mitkThrow() << "Input surface doesn't have data at time step " << t << "!";
60 
61  vtkPolyData *polyData = const_cast<mitk::Surface *>(surface.GetPointer())->GetVtkPolyData(t);
62 
63  if (polyData == nullptr)
64  mitkThrow() << "PolyData of input surface at time step " << t << " is nullptr!";
65 
66  if (polyData->GetNumberOfPolys() == 0)
67  mitkThrow() << "Input surface has no polygons at time step " << t << "!";
68 }
69 
71  unsigned int t,
72  int numVertices,
73  double gradation,
74  int subsampling,
75  double edgeSplitting,
76  int optimizationLevel,
77  bool forceManifold,
78  bool boundaryFixing)
79 {
80  ValidateSurface(surface, t);
81 
82  MITK_INFO << "Start remeshing...";
83 
84  vtkSmartPointer<vtkPolyData> surfacePolyData = vtkSmartPointer<vtkPolyData>::New();
85  surfacePolyData->DeepCopy(const_cast<Surface *>(surface.GetPointer())->GetVtkPolyData(t));
86 
87  vtkSmartPointer<vtkSurface> mesh = vtkSmartPointer<vtkSurface>::New();
88 
89  mesh->CreateFromPolyData(surfacePolyData);
90  mesh->GetCellData()->Initialize();
91  mesh->GetPointData()->Initialize();
92 
93  mesh->DisplayMeshProperties();
94 
95  if (numVertices == 0)
96  numVertices = surfacePolyData->GetNumberOfPoints();
97 
98  if (edgeSplitting != 0.0)
99  mesh->SplitLongEdges(edgeSplitting);
100 
101  vtkSmartPointer<vtkIsotropicDiscreteRemeshing> remesher = vtkSmartPointer<vtkIsotropicDiscreteRemeshing>::New();
102 
103  remesher->GetMetric()->SetGradation(gradation);
104  remesher->SetBoundaryFixing(boundaryFixing);
105  remesher->SetConsoleOutput(1);
106  remesher->SetForceManifold(forceManifold);
107  remesher->SetInput(mesh);
108  remesher->SetNumberOfClusters(numVertices);
109  remesher->SetNumberOfThreads(vtkMultiThreader::GetGlobalDefaultNumberOfThreads());
110  remesher->SetSubsamplingThreshold(subsampling);
111 
112  remesher->Remesh();
113 
114  // Optimization: Minimize distance between input surface and remeshed surface
115  if (optimizationLevel != 0)
116  {
117  ClustersQuadrics clustersQuadrics(numVertices);
118 
119  vtkSmartPointer<vtkIdList> faceList = vtkSmartPointer<vtkIdList>::New();
120  vtkSmartPointer<vtkIntArray> clustering = remesher->GetClustering();
121  vtkSmartPointer<vtkSurface> remesherInput = remesher->GetInput();
122  int clusteringType = remesher->GetClusteringType();
123  int numItems = remesher->GetNumberOfItems();
124  int numMisclassifiedItems = 0;
125 
126  for (int i = 0; i < numItems; ++i)
127  {
128  int cluster = clustering->GetValue(i);
129 
130  if (cluster >= 0 && cluster < numVertices)
131  {
132  if (clusteringType != 0)
133  {
134  remesherInput->GetVertexNeighbourFaces(i, faceList);
135  int numIds = static_cast<int>(faceList->GetNumberOfIds());
136 
137  for (int j = 0; j < numIds; ++j)
138  vtkQuadricTools::AddTriangleQuadric(
139  clustersQuadrics.Elements[cluster], remesherInput, faceList->GetId(j), false);
140  }
141  else
142  {
143  vtkQuadricTools::AddTriangleQuadric(clustersQuadrics.Elements[cluster], remesherInput, i, false);
144  }
145  }
146  else
147  {
148  ++numMisclassifiedItems;
149  }
150  }
151 
152  if (numMisclassifiedItems != 0)
153  std::cout << numMisclassifiedItems << " items with wrong cluster association" << std::endl;
154 
155  vtkSmartPointer<vtkSurface> remesherOutput = remesher->GetOutput();
156  double point[3];
157 
158  for (int i = 0; i < numVertices; ++i)
159  {
160  remesherOutput->GetPoint(i, point);
161  vtkQuadricTools::ComputeRepresentativePoint(clustersQuadrics.Elements[i], point, optimizationLevel);
162  remesherOutput->SetPointCoordinates(i, point);
163  }
164 
165  std::cout << "After quadrics post-processing:" << std::endl;
166  remesherOutput->DisplayMeshProperties();
167  }
168 
169  vtkSmartPointer<vtkPolyDataNormals> normals = vtkSmartPointer<vtkPolyDataNormals>::New();
170 
171  normals->SetInputData(remesher->GetOutput());
172  normals->AutoOrientNormalsOn();
173  normals->ComputeCellNormalsOff();
174  normals->ComputePointNormalsOn();
175  normals->ConsistencyOff();
176  normals->FlipNormalsOff();
177  normals->NonManifoldTraversalOff();
178  normals->SplittingOff();
179 
180  normals->Update();
181 
182  Surface::Pointer remeshedSurface = Surface::New();
183  remeshedSurface->SetVtkPolyData(normals->GetOutput());
184 
185  MITK_INFO << "Finished remeshing";
186 
187  return remeshedSurface;
188 }
189 
190 mitk::ACVD::RemeshFilter::RemeshFilter()
191  : m_TimeStep(0),
192  m_NumVertices(0),
193  m_Gradation(1.0),
194  m_Subsampling(10),
195  m_EdgeSplitting(0.0),
196  m_OptimizationLevel(1),
197  m_ForceManifold(false),
198  m_BoundaryFixing(false)
199 {
200  Surface::Pointer output = Surface::New();
201  this->SetNthOutput(0, output);
202 }
203 
204 mitk::ACVD::RemeshFilter::~RemeshFilter()
205 {
206 }
207 
209 {
210  Surface::Pointer output = Remesh(this->GetInput(),
211  m_TimeStep,
212  m_NumVertices,
213  m_Gradation,
214  m_Subsampling,
215  m_EdgeSplitting,
216  m_OptimizationLevel,
217  m_ForceManifold,
218  m_BoundaryFixing);
219  this->SetNthOutput(0, output);
220 }
MITKREMESHING_EXPORT Surface::Pointer Remesh(Surface::ConstPointer surface, unsigned int t, int numVertices, double gradation, int subsampling=10, double edgeSplitting=0.0, int optimizationLevel=1, bool forceManifold=false, bool boundaryFixing=false)
Remesh a surface and store the result in a new surface.
Definition: mitkACVD.cpp:70
Class for storing surfaces (vtkPolyData).
Definition: mitkSurface.h:28
#define MITK_INFO
Definition: mitkLogMacros.h:18
static void ValidateSurface(mitk::Surface::ConstPointer surface, unsigned int t)
Definition: mitkACVD.cpp:53
#define mitkThrow()
void GenerateData() override
Definition: mitkACVD.cpp:208