Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkGantryTiltInformation.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,
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 //#define MBILOG_ENABLE_DEBUG
18 
20 
21 #include "mitkDICOMTag.h"
22 
23 #include "mitkLogMacros.h"
24 
26 : m_ShiftUp(0.0)
27 , m_ShiftRight(0.0)
28 , m_ShiftNormal(0.0)
29 , m_ITKAssumedSliceSpacing(0.0)
30 , m_NumberOfSlicesApart(0)
31 {
32 }
33 
34 
35 #define doublepoint(x) \
36  Point3Dd x; \
37  x[0] = x ## f[0]; \
38  x[1] = x ## f[1]; \
39  x[2] = x ## f[2];
40 
41 
42 #define doublevector(x) \
43  Vector3Dd x; \
44  x[0] = x ## f[0]; \
45  x[1] = x ## f[1]; \
46  x[2] = x ## f[2];
47 
49  const Point3D& origin1f, const Point3D& origin2f,
50  const Vector3D& rightf, const Vector3D& upf,
51  unsigned int numberOfSlicesApart)
52 : m_ShiftUp(0.0)
53 , m_ShiftRight(0.0)
54 , m_ShiftNormal(0.0)
55 , m_NumberOfSlicesApart(numberOfSlicesApart)
56 {
57  assert(numberOfSlicesApart);
58 
59  doublepoint(origin1);
60  doublepoint(origin2);
61  doublevector(right);
62  doublevector(up);
63 
64  // determine if slice 1 (imagePosition1 and imageOrientation1) and slice 2 can be in one orthogonal slice stack:
65  // calculate a line from origin 1, directed along the normal of slice (calculated as the cross product of orientation 1)
66  // check if this line passes through origin 2
67 
68  /*
69  Determine if line (imagePosition2 + l * normal) contains imagePosition1.
70  Done by calculating the distance of imagePosition1 from line (imagePosition2 + l *normal)
71 
72  E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
73 
74  squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2
75  /
76  |pointAlongNormal - origin2| ^ 2
77 
78  ( x meaning the cross product )
79  */
80 
81  Vector3Dd normal = itk::CrossProduct(right, up);
82  Point3Dd pointAlongNormal = origin2 + normal;
83 
84  double numerator = itk::CrossProduct( pointAlongNormal - origin2 , origin2 - origin1 ).GetSquaredNorm();
85  double denominator = (pointAlongNormal - origin2).GetSquaredNorm();
86 
87  double distance = sqrt(numerator / denominator);
88 
89  if ( distance > 0.001 ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt
90  {
91  MITK_DEBUG << " Series seems to contain a tilted (or sheared) geometry";
92  MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance;
93  MITK_DEBUG << " ==> storing this shift for later analysis:";
94  MITK_DEBUG << " v origin1: " << origin1;
95  MITK_DEBUG << " v origin2: " << origin2;
96  MITK_DEBUG << " v right: " << right;
97  MITK_DEBUG << " v up: " << up;
98  MITK_DEBUG << " v normal: " << normal;
99 
100  Point3Dd projectionRight = projectPointOnLine( origin1, origin2, right );
101  Point3Dd projectionNormal = projectPointOnLine( origin1, origin2, normal );
102 
103  m_ShiftRight = (projectionRight - origin2).GetNorm();
104  m_ShiftNormal = (projectionNormal - origin2).GetNorm();
105 
106  /*
107  now also check to which side the image is shifted.
108 
109  Calculation e.g. from
110  http://mathworld.wolfram.com/Point-PlaneDistance.html
111  */
112 
113  Point3Dd testPoint = origin1;
114  Vector3Dd planeNormal = up;
115 
116  double signedDistance = (
117  planeNormal[0] * testPoint[0]
118  + planeNormal[1] * testPoint[1]
119  + planeNormal[2] * testPoint[2]
120  - (
121  planeNormal[0] * origin2[0]
122  + planeNormal[1] * origin2[1]
123  + planeNormal[2] * origin2[2]
124  )
125  )
126  /
127  sqrt( planeNormal[0] * planeNormal[0]
128  + planeNormal[1] * planeNormal[1]
129  + planeNormal[2] * planeNormal[2]
130  );
131 
132  m_ShiftUp = signedDistance;
133 
134  m_ITKAssumedSliceSpacing = (origin2 - origin1).GetNorm();
135  // How do we now this is assumed? See header documentation for ITK code references
136  //double itkAssumedSliceSpacing = sqrt( m_ShiftUp * m_ShiftUp + m_ShiftNormal * m_ShiftNormal );
137 
138  MITK_DEBUG << " calculated from slices " << m_NumberOfSlicesApart << " slices apart";
139  MITK_DEBUG << " shift normal: " << m_ShiftNormal;
140  MITK_DEBUG << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing;
141  MITK_DEBUG << " shift up: " << m_ShiftUp;
142  MITK_DEBUG << " shift right: " << m_ShiftRight;
143 
144  MITK_DEBUG << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535;
145  }
146 }
147 
151  const std::string& origin1String,
152  const std::string& origin2String,
153  const std::string orientationString,
154  unsigned int numberOfSlicesApart)
155 {
156  Vector3D right; right.Fill(0.0);
157  Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point
158  bool orientationConversion(false);
159  DICOMStringToOrientationVectors( orientationString, right, up, orientationConversion );
160 
161  if (orientationConversion
162  && !origin1String.empty() && !origin2String.empty()
163  )
164  {
165  bool firstOriginConversion(false);
166  bool lastOriginConversion(false);
167 
168  Point3D firstOrigin = DICOMStringToPoint3D( origin1String, firstOriginConversion );
169  Point3D lastOrigin = DICOMStringToPoint3D( origin2String, lastOriginConversion );
170 
171  if (firstOriginConversion && lastOriginConversion)
172  {
173  return GantryTiltInformation( firstOrigin, lastOrigin, right, up, numberOfSlicesApart );
174  }
175  }
176 
177  std::stringstream ss;
178  ss << "Invalid tag values when constructing tilt information from origin1 '" << origin1String
179  << "', origin2 '" << origin2String
180  << "', and orientation '" << orientationString << "'";
181 
182  throw std::invalid_argument(ss.str());
183 }
184 
185 void
187 ::Print(std::ostream& os) const
188 {
189  os << " calculated from slices " << m_NumberOfSlicesApart << " slices apart" << std::endl;
190  os << " shift normal: " << m_ShiftNormal << std::endl;
191  os << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing << std::endl;
192  os << " shift up: " << m_ShiftUp << std::endl;
193  os << " shift right: " << m_ShiftRight << std::endl;
194 
195  os << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535 << std::endl;
196 }
197 
199 mitk::GantryTiltInformation::projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection )
200 {
207  Vector3Dd lineOriginToP = p - lineOrigin;
208  double innerProduct = lineOriginToP * lineDirection;
209 
210  double factor = innerProduct / lineDirection.GetSquaredNorm();
211  Point3Dd projection = lineOrigin + factor * lineDirection;
212 
213  return projection;
214 }
215 
216 double
218 {
219  return fabs(m_ShiftUp / static_cast<double>(m_NumberOfSlicesApart) * static_cast<double>(imageSizeZ-1));
220 }
221 
222 double
224 {
225  return atan( fabs(m_ShiftUp) / m_ShiftNormal ) * 180.0 / 3.1415926535;
226 }
227 
228 double
230 {
231  // so many mm need to be shifted per slice!
232  return m_ShiftUp / static_cast<double>(m_NumberOfSlicesApart);
233 }
234 
235 double
237 {
238  return m_ShiftNormal / static_cast<double>(m_NumberOfSlicesApart);
239 }
240 
241 
242 bool
244 {
245  return m_NumberOfSlicesApart &&
246  ( fabs(m_ShiftRight) > 0.001
247  || fabs(m_ShiftUp) > 0.001);
248 }
249 
250 
251 bool
253 {
254  return m_NumberOfSlicesApart &&
255  ( fabs(m_ShiftRight) < 0.001
256  && fabs(m_ShiftUp) > 0.001);
257 }
258 
bool IsRegularGantryTilt() const
Whether the shearing is a gantry tilt or more complicated.
double GetMatrixCoefficientForCorrectionInWorldCoordinates() const
The offset distance in Y direction for each slice in mm (describes the tilt result).
void DICOMStringToOrientationVectors(const std::string &s, Vector3D &right, Vector3D &up, bool &successful)
Convert DICOM string describing a point two Vector3D.
#define MITK_DEBUG
Definition: mitkLogMacros.h:26
#define doublevector(x)
void Print(std::ostream &os) const
#define doublepoint(x)
double GetRealZSpacing() const
The z / inter-slice spacing. Needed to correct ImageSeriesReader's result.
bool IsSheared() const
Whether the slices were sheared.
GantryTiltInformation()
Just so we can create empty instances for assigning results later.
static GantryTiltInformation MakeFromTagValues(const std::string &origin1String, const std::string &origin2String, const std::string orientationString, unsigned int numberOfSlicesApart)
Factory method to create GantryTiltInformation from tag values (strings).
Gantry tilt analysis result.
Point3D DICOMStringToPoint3D(const std::string &s, bool &successful)
Convert DICOM string describing a point to Point3D.
double GetTiltCorrectedAdditionalSize(unsigned int imageSizeZ) const
The shift between first and last slice in mm.
itk::Vector< double, 3 > Vector3Dd
double GetTiltAngleInDegrees() const
Calculated tilt angle in degrees.