Fraxinus  18.10
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==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 ptSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1));
302  case ptRADIALPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D(-1, 0, 0));
303  case ptTOOLSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1)); //SIDE
304  default:
305  throw std::exception();
306  }
307 }
308 
312 std::pair<Vector3D,Vector3D> SliceComputer::generateBasisVectorsRadiology() const
313 {
314  switch (mPlaneType)
315  {
316  // use right-left ordering for axial/coronal
317  case ptAXIAL: return std::make_pair(Vector3D( 1, 0, 0), Vector3D( 0,-1, 0));
318  case ptCORONAL: return std::make_pair(Vector3D( 1, 0, 0), Vector3D( 0, 0, 1));
319  case ptSAGITTAL: return std::make_pair(Vector3D( 0, 1, 0), Vector3D( 0, 0, 1));
320 
321  // use planes corresponding to the cx Tool definitions
322  case ptANYPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D( 0, 0,-1));
323  case ptSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1));
324  case ptRADIALPLANE: return std::make_pair(Vector3D( 0,-1, 0), Vector3D(-1, 0, 0));
325  case ptTOOLSIDEPLANE: return std::make_pair(Vector3D(-1, 0, 0), Vector3D( 0, 0,-1)); //SIDE
326  default:
327  throw std::exception();
328  }
329 }
330 
338 Vector3D SliceComputer::generateFixedIJCenter(const Vector3D& center_r, const Vector3D& cross_r, const Vector3D& i, const Vector3D& j) const
339 {
340  if (mFollowType==ftFIXED_CENTER)
341  {
342  // r is REF, s is SLICE
343  Transform3D M_rs = createTransformIJC(i, j, Vector3D(0,0,0)); // transform from data to slice, zero center.
344  Transform3D M_sr = M_rs.inv();
345  Vector3D center_s = M_sr.coord(center_r);
346  Vector3D cross_s = M_sr.coord(cross_r);
347  // in SLICE space, use {xy} values from center and {z} value from cross.
348  Vector3D q_s(center_s[0], center_s[1], cross_s[2]);
349  Vector3D q_r = M_rs.coord(q_s);
350  return q_r;
351  }
352  return cross_r;
353 }
354 
366 SlicePlane SliceComputer::orientToGravity(const SlicePlane& base) const
367 {
368  if (!mUseGravity)
369  {
370  return base;
371  }
372 
373  SlicePlane retval = base;
374  const Vector3D k = cross(base.i, base.j); // plane normal. Constant
375  Vector3D up;
376  up = -mGravityDirection; // normal case
377 
378  // weight of nongravity, 0=<w=<1, 1 means dont use gravity
379  double w_n = dot(up, k);
380  w_n = w_n*w_n; // square to keep stability near normal use.
381 
382  Vector3D i_g = cross(up, k); // |i_g| varies from 0 to 1 depending on 1-w_n
383  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
384 
385  // set i vector to a weighted mean of the two definitions
386  // can also experiment with a tanh function or simply a linear interpolation
387  //
388  // Note: i_g is already small here if w_n is small, this will increase
389  // the effect of i_n. Investigate.
390  //
391  retval.i = i_g*(1.0-w_n) + i_n*w_n;
392  retval.i = retval.i.normal(); // |i|==1
393  retval.j = cross(k, retval.i);
394 
395  return retval;
396 }
397 
421 SlicePlane SliceComputer::orientToGravityAroundToolZAxisAndAlongTheOperatingTable(const SlicePlane &base) const
422 {
423  if (!mUseGravity)
424  {
425  return base;
426  }
427 
428  SlicePlane retval = base;
429  Vector3D up = -mGravityDirection;
430 
431  Vector3D k_neg = -cross(base.i, base.j);
432  Vector3D toolVector = m_rMt.vector(Vector3D(0, 0, 1));
433  Vector3D i_perpendicular = cross(up, toolVector).normal();
434 
435  double w_n = this->getWeightForAngularDifference(up, toolVector);
436 
437  i_perpendicular = i_perpendicular*(1.0-w_n) + k_neg*w_n;
438 
439  Vector3D i_mark = cross(i_perpendicular, up).normal();
440  Vector3D j_mark = up;
441 
442  retval.i = i_mark;
443  retval.j = j_mark;
444 
445  return retval;
446 }
447 
453 double SliceComputer::getWeightForAngularDifference(Vector3D a, Vector3D b) const
454 {
455  double w_n = dot(a, b);
456  w_n = fabs(w_n);
457  // w_n = 0 : normal case
458  // w_n = 1 : singularity
459 
460  double cutoff = sqrt(3.0)/2.0; // cutoff at 30*, i.e. use only toolvector up to that angle between up and tool
461  if (w_n<cutoff)
462  w_n = 0;
463  else
464  w_n = (w_n-cutoff)/(1.0-cutoff);
465  return w_n;
466 }
467 
468 
469 } // 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:38
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:50
otOBLIQUE
orient planes relative to the tool space
Definition: cxDefinitions.h:32
mdRADIOLOGICAL
Definition: cxDefinitions.h:58
ftFIXED_CENTER
center is set.
Definition: cxDefinitions.h:50
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:38
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:38
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:32
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:38
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
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:38
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:58
ptANYPLANE
a plane aligned with the tool base plane
Definition: cxDefinitions.h:38
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:38
Namespace for all CustusX production code.