Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkClippedSurfaceBoundsCalculator.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 ============================================================================*/
13 #include "mitkLine.h"
14 
15 #define ROUND_P(x) ((x) >= 0 ? (int)((x) + 0.5) : (int)((x)-0.5))
16 
19  : m_PlaneGeometry(nullptr), m_Geometry3D(nullptr), m_Image(nullptr)
20 {
21  this->InitializeOutput();
22 
23  this->SetInput(geometry, image);
24 }
25 
28  : m_PlaneGeometry(nullptr), m_Geometry3D(nullptr), m_Image(nullptr)
29 {
30  this->InitializeOutput();
31 
32  this->SetInput(geometry, image);
33 }
34 
37  : m_PlaneGeometry(nullptr), m_Geometry3D(nullptr), m_Image(image)
38 {
39  this->InitializeOutput();
40 
42 }
43 
45 {
46  // initialize with meaningless slice indices
47  m_MinMaxOutput.clear();
48  m_MinMaxOutput.reserve(3);
49  for (int i = 0; i < 3; i++)
50  {
52  }
53 }
54 
56 {
57 }
58 
60 {
61  if (geometry && image)
62  {
63  this->m_PlaneGeometry = geometry;
64  this->m_Image = image;
65  this->m_Geometry3D = nullptr; // Not possible to set both
67  }
68 }
69 
71 {
72  if (geometry && image)
73  {
74  this->m_Geometry3D = geometry;
75  this->m_Image = image;
76  this->m_PlaneGeometry = nullptr; // Not possible to set both
78  }
79 }
80 
82 {
83  if (!pointlist.empty() && image)
84  {
85  m_Geometry3D = nullptr;
86  m_PlaneGeometry = nullptr;
87  m_Image = image;
89  }
90 }
91 
93 {
94  return this->m_MinMaxOutput[0];
95 }
96 
98 {
99  return this->m_MinMaxOutput[1];
100 }
101 
103 {
104  return this->m_MinMaxOutput[2];
105 }
106 
108 {
109  this->m_MinMaxOutput.clear();
110  m_MinMaxOutput.reserve(3);
111  for (int i = 0; i < 3; i++)
112  {
114  }
115 
116  if (m_PlaneGeometry.IsNotNull())
117  {
119  }
120  else if (m_Geometry3D.IsNotNull())
121  {
122  // go through all slices of the image, ...
123  const auto *slicedGeometry3D =
124  dynamic_cast<const mitk::SlicedGeometry3D *>(m_Geometry3D.GetPointer());
125  int allSlices = slicedGeometry3D->GetSlices();
126  this->CalculateIntersectionPoints(dynamic_cast<mitk::PlaneGeometry *>(slicedGeometry3D->GetPlaneGeometry(0)));
128  dynamic_cast<mitk::PlaneGeometry *>(slicedGeometry3D->GetPlaneGeometry(allSlices - 1)));
129  }
130  else if (!m_ObjectPointsInWorldCoordinates.empty())
131  {
133  this->EnforceImageBounds();
134  }
135 }
136 
138 {
139  // SEE HEADER DOCUMENTATION for explanation
140 
141  const mitk::BaseGeometry::Pointer imageGeometry = m_Image->GetGeometry()->Clone();
142 
143  // the cornerpoint(0) is the corner based Origin, which is original center based
144  Point3D origin = imageGeometry->GetCornerPoint(0); // Left, bottom, front
145 
146  // Get axis vector for the spatial directions
147  const Vector3D xDirection = imageGeometry->GetAxisVector(0);
148  const Vector3D yDirection = imageGeometry->GetAxisVector(1);
149  const Vector3D zDirection = imageGeometry->GetAxisVector(2);
150 
151  const Point3D leftBottomFront = origin;
152  const Point3D leftTopFront = origin + yDirection;
153  const Point3D leftBottomBack = origin + zDirection;
154  const Point3D leftTopBack = origin + yDirection + zDirection;
155  const Point3D rightBottomFront = origin + xDirection;
156  const Point3D rightTopFront = origin + xDirection + yDirection;
157  const Point3D rightBottomBack = origin + xDirection + zDirection;
158  const Point3D rightTopBack = origin + xDirection + yDirection + zDirection;
159 
160  typedef std::vector<std::pair<mitk::Point3D, mitk::Point3D>> EdgesVector;
161  EdgesVector edgesOf3DBox;
162  edgesOf3DBox.reserve(12);
163 
164  edgesOf3DBox.push_back(std::make_pair(leftBottomFront, // x = left=xfront, y=bottom=yfront, z=front=zfront
165  leftTopFront)); // left, top, front
166 
167  edgesOf3DBox.push_back(std::make_pair(leftBottomFront, // left, bottom, front
168  leftBottomBack)); // left, bottom, back
169 
170  edgesOf3DBox.push_back(std::make_pair(leftBottomFront, // left, bottom, front
171  rightBottomFront)); // right, bottom, front
172 
173  edgesOf3DBox.push_back(std::make_pair(leftTopFront, // left, top, front
174  rightTopFront)); // right, top, front
175 
176  edgesOf3DBox.push_back(std::make_pair(leftTopFront, // left, top, front
177  leftTopBack)); // left, top, back
178 
179  edgesOf3DBox.push_back(std::make_pair(rightTopFront, // right, top, front
180  rightTopBack)); // right, top, back
181 
182  edgesOf3DBox.push_back(std::make_pair(rightTopFront, // right, top, front
183  rightBottomFront)); // right, bottom, front
184 
185  edgesOf3DBox.push_back(std::make_pair(rightBottomFront, // right, bottom, front
186  rightBottomBack)); // right, bottom, back
187 
188  edgesOf3DBox.push_back(std::make_pair(rightBottomBack, // right, bottom, back
189  leftBottomBack)); // left, bottom, back
190 
191  edgesOf3DBox.push_back(std::make_pair(rightBottomBack, // right, bottom, back
192  rightTopBack)); // right, top, back
193 
194  edgesOf3DBox.push_back(std::make_pair(rightTopBack, // right, top, back
195  leftTopBack)); // left, top, back
196 
197  edgesOf3DBox.push_back(std::make_pair(leftTopBack, // left, top, back
198  leftBottomBack)); // left, bottom, back
199 
200  for (auto iterator = edgesOf3DBox.cbegin(); iterator != edgesOf3DBox.cend(); ++iterator)
201  {
202  const Point3D startPoint = (*iterator).first; // start point of the line
203  const Point3D endPoint = (*iterator).second; // end point of the line
204  const Vector3D lineDirection = endPoint - startPoint;
205 
206  const mitk::Line3D line(startPoint, lineDirection);
207 
208  // Get intersection point of line and plane geometry
209  Point3D intersectionWorldPoint(std::numeric_limits<int>::min());
210 
211  double t = -1.0;
212  bool doesLineIntersectWithPlane(false);
213 
214  const double norm = line.GetDirection().GetNorm();
215  const double dist = geometry->Distance(line.GetPoint1());
216  if (norm < mitk::eps && dist < mitk::sqrteps)
217  {
218  t = 1.0;
219  doesLineIntersectWithPlane = true;
220  intersectionWorldPoint = line.GetPoint1();
221  }
222  else
223  {
224  geometry->IntersectionPoint(line, intersectionWorldPoint);
225  doesLineIntersectWithPlane = geometry->IntersectionPointParam(line, t);
226  }
227 
228  // Get index point
229  mitk::Point3D intersectionIndexPoint;
230  imageGeometry->WorldToIndex(intersectionWorldPoint, intersectionIndexPoint);
231 
232  const bool lowerBoundGood = (0 - mitk::sqrteps) <= t;
233  const bool upperBoundGood = t <= 1.0 + mitk::sqrteps;
234  if (doesLineIntersectWithPlane && lowerBoundGood && upperBoundGood)
235  {
236  for (int dim = 0; dim < 3; ++dim)
237  {
238  m_MinMaxOutput[dim].first = std::min(m_MinMaxOutput[dim].first, ROUND_P(intersectionIndexPoint[dim]));
239  m_MinMaxOutput[dim].second = std::max(m_MinMaxOutput[dim].second, ROUND_P(intersectionIndexPoint[dim]));
240  }
241  this->EnforceImageBounds();
242  }
243  }
244 }
245 
247 {
248  PointListType::const_iterator pointIterator;
249 
250  const mitk::SlicedGeometry3D::Pointer imageGeometry = m_Image->GetSlicedGeometry();
251  for (pointIterator = pointList.cbegin(); pointIterator != pointList.cend(); ++pointIterator)
252  {
253  mitk::Point3D pntInIndexCoordinates;
254  imageGeometry->WorldToIndex((*pointIterator), pntInIndexCoordinates);
255 
256  m_MinMaxOutput[0].first =
257  pntInIndexCoordinates[0] < m_MinMaxOutput[0].first ? ROUND_P(pntInIndexCoordinates[0]) : m_MinMaxOutput[0].first;
258  m_MinMaxOutput[0].second = pntInIndexCoordinates[0] > m_MinMaxOutput[0].second ? ROUND_P(pntInIndexCoordinates[0]) :
259  m_MinMaxOutput[0].second;
260 
261  m_MinMaxOutput[1].first =
262  pntInIndexCoordinates[1] < m_MinMaxOutput[1].first ? ROUND_P(pntInIndexCoordinates[1]) : m_MinMaxOutput[1].first;
263  m_MinMaxOutput[1].second = pntInIndexCoordinates[1] > m_MinMaxOutput[1].second ? ROUND_P(pntInIndexCoordinates[1]) :
264  m_MinMaxOutput[1].second;
265 
266  m_MinMaxOutput[2].first =
267  pntInIndexCoordinates[2] < m_MinMaxOutput[2].first ? ROUND_P(pntInIndexCoordinates[2]) : m_MinMaxOutput[2].first;
268  m_MinMaxOutput[2].second = pntInIndexCoordinates[2] > m_MinMaxOutput[2].second ? ROUND_P(pntInIndexCoordinates[2]) :
269  m_MinMaxOutput[2].second;
270  }
271 
272  // this->EnforceImageBounds();
273 }
274 
276 {
277  m_MinMaxOutput[0].first = std::max(m_MinMaxOutput[0].first, 0);
278  m_MinMaxOutput[1].first = std::max(m_MinMaxOutput[1].first, 0);
279  m_MinMaxOutput[2].first = std::max(m_MinMaxOutput[2].first, 0);
280 
281  m_MinMaxOutput[0].first = std::min(m_MinMaxOutput[0].first, (int)m_Image->GetDimension(0) - 1);
282  m_MinMaxOutput[1].first = std::min(m_MinMaxOutput[1].first, (int)m_Image->GetDimension(1) - 1);
283  m_MinMaxOutput[2].first = std::min(m_MinMaxOutput[2].first, (int)m_Image->GetDimension(2) - 1);
284 
285  m_MinMaxOutput[0].second = std::min(m_MinMaxOutput[0].second, (int)m_Image->GetDimension(0) - 1);
286  m_MinMaxOutput[1].second = std::min(m_MinMaxOutput[1].second, (int)m_Image->GetDimension(1) - 1);
287  m_MinMaxOutput[2].second = std::min(m_MinMaxOutput[2].second, (int)m_Image->GetDimension(2) - 1);
288 
289  m_MinMaxOutput[0].second = std::max(m_MinMaxOutput[0].second, 0);
290  m_MinMaxOutput[1].second = std::max(m_MinMaxOutput[1].second, 0);
291  m_MinMaxOutput[2].second = std::max(m_MinMaxOutput[2].second, 0);
292 }
OutputType GetMinMaxSpatialDirectionZ()
What Z coordinates (slice indices) are cut/visible in given plane.
void CalculateIntersectionPoints(const mitk::PlaneGeometry *geometry)
Descibes a line.
Definition: mitkLine.h:28
static char * line
Definition: svm.cpp:2870
void SetInput(const mitk::PlaneGeometry *geometry, mitk::Image *image)
bool IntersectionPointParam(const Line3D &line, double &t) const
Calculate line parameter of intersection point between the plane and a line.
virtual unsigned int GetSlices() const
Get the number of slices.
bool IntersectionPoint(const Line3D &line, Point3D &intersectionPoint) const
Calculate intersection point between the plane and a line.
const itk::Point< TCoordRep, NPointDimension > & GetPoint1() const
Get start point of the line.
Definition: mitkLine.h:111
const itk::Vector< TCoordRep, NPointDimension > & GetDirection() const
Get the direction vector of the line.
Definition: mitkLine.h:70
void EnforceImageBounds()
Clips the resulting index-coordinates to make sure they do not exceed the imagebounds.
Image class for storing images.
Definition: mitkImage.h:72
ClippedSurfaceBoundsCalculator(const mitk::PlaneGeometry *geometry=nullptr, mitk::Image::Pointer image=nullptr)
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
Describes the geometry of a data object consisting of slices.
static T min(T x, T y)
Definition: svm.cpp:53
MITKCORE_EXPORT const ScalarType sqrteps
MITKCORE_EXPORT const ScalarType eps
ScalarType Distance(const Point3D &pt3d_mm) const
Distance of the point from the geometry (bounding-box not considered)
OutputType GetMinMaxSpatialDirectionY()
What Y coordinates (slice indices) are cut/visible in given plane.
Describes a two-dimensional, rectangular plane.
OutputType GetMinMaxSpatialDirectionX()
What X coordinates (slice indices) are cut/visible in given plane.
std::pair< int, int > OutputType
Minimum (first) and maximum (second) slice index.
BaseGeometry Describes the geometry of a data object.