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