CustusX  22.04-rc5
An IGT application
cxSliceComputer.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 This file is part of CustusX, an Image Guided Therapy Application.
3 
4 Copyright (c) SINTEF Department of Medical Technology.
5 All rights reserved.
6 
7 CustusX is released under a BSD 3-Clause license.
8 
9 See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
10 =========================================================================*/
11 
12 #include "cxSliceComputer.h"
13 #include "cxDefinitions.h"
14 #include <math.h>
15 
16 namespace cx
17 {
18 
19 SlicePlane::SlicePlane(const Vector3D& c_, const Vector3D& i_, const Vector3D& j_) :
20  c(c_), i(i_), j(j_)
21 {
22 }
23 
24 std::ostream& operator<<(std::ostream& s, const SlicePlane& val)
25 {
26  s << "center : " << val.c << std::endl;
27  s << "i_vector : " << val.i << std::endl;
28  s << "j_vector : " << val.j << std::endl;
29  return s;
30 }
31 
32 bool similar(const SlicePlane& a, const SlicePlane& b)
33 {
34  return similar(a.c, b.c) && similar(a.i, b.i) && similar(a.j, b.j);
35 }
36 
40  mClinicalApplication(mdNEUROLOGICAL),
41  mOrientType(otORTHOGONAL),
42  mPlaneType(ptAXIAL),
43  mFollowType(ftFIXED_CENTER),
44  mFixedCenter(Vector3D(0,0,0)),
45  m_rMt(Transform3D::Identity()),
46  mToolOffset(0.0),
47  mUseGravity(false),
48  mGravityDirection(Vector3D(0,0,-1)) ,
49  mUseViewOffset(false),
50  mViewportHeight(1),
51  mViewOffset(0.5)
52 {
53 }
54 
56 {
57 }
58 
61 void SliceComputer::initializeFromPlane(PLANE_TYPE plane, bool useGravity, const Vector3D& gravityDir, bool useViewOffset, double viewportHeight, double toolViewOffset, CLINICAL_VIEW application)
62 {
63  setPlaneType(plane);
64  mClinicalApplication = application;
65 
66  if (plane == ptSAGITTAL || plane == ptCORONAL || plane == ptAXIAL )
67  {
70  }
71  else if (plane == ptANYPLANE || plane == ptINVERSEANYPLANE || plane==ptRADIALPLANE || plane==ptSIDEPLANE)
72  {
75 
76  setGravity(useGravity, gravityDir);
77  setToolViewOffset(useViewOffset, viewportHeight, toolViewOffset); // TODO finish this one
78  }
79  else if (plane==ptTOOLSIDEPLANE)
80  {
83  setGravity(useGravity, gravityDir);
84  }
85 }
86 
87 void SliceComputer::setClinicalApplication(CLINICAL_VIEW application)
88 {
89  mClinicalApplication = application;
90 }
91 
92 ORIENTATION_TYPE SliceComputer::getOrientationType() const
93 {
94  return mOrientType;
95 }
96 
97 PLANE_TYPE SliceComputer::getPlaneType() const
98 {
99  return mPlaneType;
100 }
101 
102 FOLLOW_TYPE SliceComputer::getFollowType() const
103 {
104  return mFollowType;
105 }
106 
108 {
109  return m_rMt;
110 }
111 
116 {
117  m_rMt = rMt;
118 }
119 
123 void SliceComputer::setOrientationType(ORIENTATION_TYPE val)
124 {
125  mOrientType = val;
126 }
127 
130 void SliceComputer::setPlaneType(PLANE_TYPE val)
131 {
132  mPlaneType = val;
133 }
134 
139 {
140  mFixedCenter = center;
141 }
142 
147 void SliceComputer::setFollowType(FOLLOW_TYPE val)
148 {
149  mFollowType = val;
150 }
151 
155 void SliceComputer::setGravity(bool use, const Vector3D& dir)
156 {
157  mUseGravity = use;
158  mGravityDirection = dir;
159 }
160 
164 {
165  mToolOffset = val;
166 }
167 
172 void SliceComputer::setToolViewOffset(bool use, double viewportHeight, double viewOffset)
173 {
174  mUseViewOffset = use;
175  mViewportHeight = viewportHeight;
176  mViewOffset = viewOffset;
177 }
178 
181 void SliceComputer::setToolViewportHeight(double viewportHeight)
182 {
183  mViewportHeight = viewportHeight;
184 }
185 
189 {
190  std::pair<Vector3D,Vector3D> basis = generateBasisVectors();
191  SlicePlane plane;
192  plane.i = basis.first;
193  plane.j = basis.second;
194  plane.c = Vector3D(0,0,mToolOffset);
195 
196  // transform position from tool to reference space
197  plane.c = m_rMt.coord(plane.c);
198  // transform orientation from tool to reference space for the oblique case only
199  if (mOrientType==otOBLIQUE)
200  {
201  plane.i = m_rMt.vector(plane.i);
202  plane.j = m_rMt.vector(plane.j);
203  }
204 
205  if (mPlaneType == ptTOOLSIDEPLANE)
206  {
207  plane = this->orientToGravityAroundToolZAxisAndAlongTheOperatingTable(plane);
208  }
209  else
210  {
211  plane = orientToGravity(plane);
212  }
213 
214  // try to to this also for oblique views, IF the ftFIXED_CENTER is set.
215  // use special acs centermod algo
216  plane.c = generateFixedIJCenter(mFixedCenter, plane.c, plane.i, plane.j);
217 
218  // set center so that it is a fixed distance from viewport top
219  plane = applyViewOffset(plane);
220 
221  return plane;
222 }
223 
233 SlicePlane SliceComputer::applyViewOffset(const SlicePlane& base) const
234 {
235  if (!mUseViewOffset)
236  {
237  return base;
238  }
239 
240  double centerOffset = this->getViewOffsetAbsoluteFromCenter();
241 
242  SlicePlane retval = base;
243 
244  // limit by tooloffset
245  Vector3D toolOffsetCenter = m_rMt.coord(Vector3D(0,0,mToolOffset));
246  Vector3D newCenter = toolOffsetCenter + centerOffset * base.j;
247  double toolOffsetDistance = dot(newCenter - base.c, base.j);
248 
249  // limit by tooltip
250  Vector3D toolCenter = m_rMt.coord(Vector3D(0,0,0));
251  newCenter = toolCenter - centerOffset * base.j;
252  double toolDistance = dot(newCenter - base.c, base.j);
253 
254  // select a dist and apply
255  double usedDistance = std::min(toolOffsetDistance, toolDistance);
256  retval.c = base.c + usedDistance * base.j; // extract j-component of newCenter
257 
258  return retval;
259 }
260 
261 double SliceComputer::getViewOffsetAbsoluteFromCenter() const
262 {
263  if (mPlaneType==ptRADIALPLANE)
264  return 0; // position in the center
265 
266  return mViewportHeight*(0.5-mViewOffset);
267 }
268 
275 std::pair<Vector3D,Vector3D> SliceComputer::generateBasisVectors() const
276 {
277  switch (mClinicalApplication)
278  {
279  case mdRADIOLOGICAL:
280  return this->generateBasisVectorsRadiology();
281  case mdNEUROLOGICAL:
282  default:
283  return this->generateBasisVectorsNeurology();
284  }
285 }
286 
290 std::pair<Vector3D,Vector3D> SliceComputer::generateBasisVectorsNeurology() const
291 {
292  switch (mPlaneType)
293  {
294  // use left-right ordering for axial/coronal
295  case ptAXIAL: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0,-1, 0));
296  case ptCORONAL: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0, 1));
297  case ptSAGITTAL: return std::make_pair(Vector3D( 0, 1, 0), Vector3D( 0, 0, 1));
298 
299  // use planes corresponding to the cx Tool definitions
300  case ptANYPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D( 0, 0,-1));
301  case ptINVERSEANYPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D( 0, 0, 1));
302  case ptSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1));
303  case ptRADIALPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D(-1, 0, 0));
304  case ptTOOLSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1)); //SIDE
305  default:
306  throw std::exception();
307  }
308 }
309 
313 std::pair<Vector3D,Vector3D> SliceComputer::generateBasisVectorsRadiology() const
314 {
315  switch (mPlaneType)
316  {
317  // use right-left ordering for axial/coronal
318  case ptAXIAL: return std::make_pair(Vector3D( 1, 0, 0), Vector3D( 0,-1, 0));
319  case ptCORONAL: return std::make_pair(Vector3D( 1, 0, 0), Vector3D( 0, 0, 1));
320  case ptSAGITTAL: return std::make_pair(Vector3D( 0, 1, 0), Vector3D( 0, 0, 1));
321 
322  // use planes corresponding to the cx Tool definitions
323  case ptANYPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D( 0, 0,-1));
324  case ptINVERSEANYPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D( 0, 0, 1));
325  case ptSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1));
326  case ptRADIALPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D(-1, 0, 0));
327  case ptTOOLSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1)); //SIDE
328  default:
329  throw std::exception();
330  }
331 }
332 
340 Vector3D SliceComputer::generateFixedIJCenter(const Vector3D& center_r, const Vector3D& cross_r, const Vector3D& i, const Vector3D& j) const
341 {
342  if (mFollowType==ftFIXED_CENTER)
343  {
344  // r is REF, s is SLICE
345  Transform3D M_rs = createTransformIJC(i, j, Vector3D(0,0,0)); // transform from data to slice, zero center.
346  Transform3D M_sr = M_rs.inv();
347  Vector3D center_s = M_sr.coord(center_r);
348  Vector3D cross_s = M_sr.coord(cross_r);
349  // in SLICE space, use {xy} values from center and {z} value from cross.
350  Vector3D q_s(center_s[0], center_s[1], cross_s[2]);
351  Vector3D q_r = M_rs.coord(q_s);
352  return q_r;
353  }
354  return cross_r;
355 }
356 
368 SlicePlane SliceComputer::orientToGravity(const SlicePlane& base) const
369 {
370  if (!mUseGravity)
371  {
372  return base;
373  }
374 
375  SlicePlane retval = base;
376  const Vector3D k = cross(base.i, base.j); // plane normal. Constant
377  Vector3D up;
378  up = -mGravityDirection; // normal case
379 
380  // weight of nongravity, 0=<w=<1, 1 means dont use gravity
381  double w_n = dot(up, k);
382  w_n = w_n*w_n; // square to keep stability near normal use.
383 
384  Vector3D i_g = cross(up, k); // |i_g| varies from 0 to 1 depending on 1-w_n
385  Vector3D i_n = base.i; // |i_n|==1 //It seems to me that i_n might need a change to e.g. i_n = -base.j, to be good in the singularity situation. Look into that if using this method later. jone, 20160712
386 
387  // set i vector to a weighted mean of the two definitions
388  // can also experiment with a tanh function or simply a linear interpolation
389  //
390  // Note: i_g is already small here if w_n is small, this will increase
391  // the effect of i_n. Investigate.
392  //
393  retval.i = i_g*(1.0-w_n) + i_n*w_n;
394  retval.i = retval.i.normal(); // |i|==1
395  retval.j = cross(k, retval.i);
396 
397  return retval;
398 }
399 
423 SlicePlane SliceComputer::orientToGravityAroundToolZAxisAndAlongTheOperatingTable(const SlicePlane &base) const
424 {
425  if (!mUseGravity)
426  {
427  return base;
428  }
429 
430  SlicePlane retval = base;
431  Vector3D up = -mGravityDirection;
432 
433  Vector3D k_neg = -cross(base.i, base.j);
434  Vector3D toolVector = m_rMt.vector(Vector3D(0, 0, 1));
435  Vector3D i_perpendicular = cross(up, toolVector).normal();
436 
437  double w_n = this->getWeightForAngularDifference(up, toolVector);
438 
439  i_perpendicular = i_perpendicular*(1.0-w_n) + k_neg*w_n;
440 
441  Vector3D i_mark = cross(i_perpendicular, up).normal();
442  Vector3D j_mark = up;
443 
444  retval.i = i_mark;
445  retval.j = j_mark;
446 
447  return retval;
448 }
449 
455 double SliceComputer::getWeightForAngularDifference(Vector3D a, Vector3D b) const
456 {
457  double w_n = dot(a, b);
458  w_n = fabs(w_n);
459  // w_n = 0 : normal case
460  // w_n = 1 : singularity
461 
462  double cutoff = sqrt(3.0)/2.0; // cutoff at 30*, i.e. use only toolvector up to that angle between up and tool
463  if (w_n<cutoff)
464  w_n = 0;
465  else
466  w_n = (w_n-cutoff)/(1.0-cutoff);
467  return w_n;
468 }
469 
470 
471 } // namespace cx
Vector3D j
defines the second axis of the plane. unit vector
void setToolViewportHeight(double viewportHeight)
ptCORONAL
a slice seen from the front of the patient
Definition: cxDefinitions.h:39
void setToolViewOffset(bool use, double viewportHeight, double viewOffset)
A 2D slice plane in 3D. i,j are perpendicular unit vectors.
ftFOLLOW_TOOL
center follows tool
Definition: cxDefinitions.h:52
otOBLIQUE
orient planes relative to the tool space
Definition: cxDefinitions.h:33
mdRADIOLOGICAL
Definition: cxDefinitions.h:60
ftFIXED_CENTER
center is set.
Definition: cxDefinitions.h:52
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
void setClinicalApplication(CLINICAL_VIEW application)
ptAXIAL
a slice seen from the top of the patient
Definition: cxDefinitions.h:39
PLANE_TYPE getPlaneType() const
void setOrientationType(ORIENTATION_TYPE val)
Vector3D cross(const Vector3D &a, const Vector3D &b)
compute cross product of a and b.
Definition: cxVector3D.cpp:41
void setToolOffset(double val)
FOLLOW_TYPE getFollowType() const
ptSAGITTAL
a slice seen from the side of the patient
Definition: cxDefinitions.h:39
Transform3D createTransformIJC(const Vector3D &ivec, const Vector3D &jvec, const Vector3D &center)
Vector3D c
defines the center of the plane
void setPlaneType(PLANE_TYPE val)
otORTHOGONAL
orient planes relative to the image/reference space.
Definition: cxDefinitions.h:33
void setFixedCenter(const Vector3D &center)
ptTOOLSIDEPLANE
z-rotated 90* relative to anyplane like side plane, but always kept oriented like the plane defined b...
Definition: cxDefinitions.h:39
void setToolPosition(const Transform3D &rMt)
double dot(const Vector3D &a, const Vector3D &b)
compute inner product (or dot product) of a and b.
Definition: cxVector3D.cpp:46
ptINVERSEANYPLANE
a plane aligned with the tool base plane, inverse of tool direction
Definition: cxDefinitions.h:39
std::ostream & operator<<(std::ostream &s, const IntBoundingBox3D &data)
void initializeFromPlane(PLANE_TYPE plane, bool useGravity, const Vector3D &gravityDir, bool useViewOffset, double viewportHeight, double toolViewOffset, CLINICAL_VIEW application)
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
ptRADIALPLANE
y-rotated 90* relative to anyplane (bird&#39;s view)
Definition: cxDefinitions.h:39
Vector3D i
defines the first axis of the plane. unit vector
bool similar(const CameraInfo &lhs, const CameraInfo &rhs, double tol)
Transform3D getToolPosition() const
mdNEUROLOGICAL
Definition: cxDefinitions.h:60
ptANYPLANE
a plane aligned with the tool base plane
Definition: cxDefinitions.h:39
void setFollowType(FOLLOW_TYPE val)
void setGravity(bool use, const Vector3D &dir)
ORIENTATION_TYPE getOrientationType() const
SlicePlane getPlane() const
ptSIDEPLANE
z-rotated 90* relative to anyplane (dual anyplane)
Definition: cxDefinitions.h:39
Namespace for all CustusX production code.