Fraxinus  17.12
An IGT application
cxTransform3D.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) 2008-2014, SINTEF Department of Medical Technology
5 All rights reserved.
6 
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions are met:
9 
10 1. Redistributions of source code must retain the above copyright notice,
11  this list of conditions and the following disclaimer.
12 
13 2. Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16 
17 3. Neither the name of the copyright holder nor the names of its contributors
18  may be used to endorse or promote products derived from this software
19  without specific prior written permission.
20 
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
25 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
29 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 =========================================================================*/
32 
33 #include "cxTransform3D.h"
34 
35 #include <sstream>
36 #include <vector>
37 #include <vtkMatrix4x4.h>
38 #include <vtkTransform.h>
39 #include "cxTypeConversions.h"
40 #include "cxBoundingBox3D.h"
41 #include "vtkForwardDeclarations.h"
42 
43 // --------------------------------------------------------
45 {
46 
47 boost::array<double, 16> flatten(const Eigen::Affine3d* self)
48 {
49  boost::array<double, 16> retval;
50  boost::array<double, 16>::iterator raw = retval.begin();
51 
52  for (int r = 0; r < 4; ++r)
53  for (int c = 0; c < 4; ++c)
54  *raw++ = (*self)(r, c);
55 
56  return retval;
57 }
58 
59 void fill(Eigen::Affine3d* self, vtkMatrix4x4Ptr m)
60 {
61  if (!m)
62  return;
63  for (int r = 0; r < 4; ++r)
64  for (int c = 0; c < 4; ++c)
65  (*self)(r, c) = m->GetElement(r, c);
66 }
67 
68 void fill(Eigen::Affine3d* self, float m[4][4])
69 {
70  if (!m)
71  return;
72  for (int r = 0; r < 4; ++r)
73  for (int c = 0; c < 4; ++c)
74  (*self)(r, c) = m[r][c];
75 }
76 
77 
81 void fill(Eigen::Affine3d* self, const double* raw)
82 {
83  for (int r = 0; r < 4; ++r)
84  for (int c = 0; c < 4; ++c)
85  (*self)(r, c) = *raw++;
86 }
87 
88 vtkMatrix4x4Ptr getVtkMatrix(const Eigen::Affine3d* self)
89 {
90  vtkMatrix4x4Ptr m = vtkMatrix4x4Ptr::New();
91 
92  for (int r = 0; r < 4; ++r)
93  for (int c = 0; c < 4; ++c)
94  m->SetElement(r, c, (*self)(r, c));
95  ;
96 
97  return m;
98 }
99 
100 vtkTransformPtr getVtkTransform(const Eigen::Affine3d* self)
101 {
102  vtkTransformPtr retval = vtkTransform::New();
103  retval->SetMatrix(self->getVtkMatrix());
104  retval->Update();
105  return retval;
106 }
107 
108 std::ostream& put(const Eigen::Affine3d* self, std::ostream& s, int indent, char newline)
109 {
110  QString ind(indent, ' ');
111 
112  std::ostringstream ss; // avoid changing state of input stream
113  ss << setprecision(3) << std::fixed;
114 
115  for (unsigned i = 0; i < 4; ++i)
116  {
117  ss << ind;
118  for (unsigned j = 0; j < 4; ++j)
119  {
120  ss << setw(10) << (*self)(i, j) << " ";
121  }
122  if (i != 3)
123  {
124  ss << newline;
125  }
126  }
127 
128  s << ss.str();
129 
130  return s;
131 }
132 
133 Eigen::Affine3d fromString(const QString& text, bool* _ok)
134 {
135  bool okval = false; // if input _ok is null, we still need a flag
136  bool* ok = &okval;
137  if (_ok)
138  ok = _ok;
139 
140  std::vector<double> raw = convertQString2DoubleVector(text, ok);
141  if (raw.size() != 16)
142  *ok = false;
143  if (!*ok)
144  return Eigen::Affine3d();
145 
146  Eigen::Affine3d retval;
147  fill(&retval, &*raw.begin());
148  return retval;
149 }
150 
151 } // namespace cx_transform3D_internal
152 
153 namespace cx
154 {
155 
156 bool similar(const Transform3D& a, const Transform3D& b, double tol)
157 {
158  boost::array<double, 16> m = a.flatten();
159  boost::array<double, 16> n = b.flatten();
160  for (int j = 0; j < 16; ++j)
161  {
162  if (!similar(n[j], m[j], tol))
163  {
164  return false;
165  }
166  }
167  return true;
168 }
169 // --------------------------------------------------------
170 
172 {
173  Vector3D a = m.coord(bb.bottomLeft());
174  Vector3D b = m.coord(bb.topRight());
175  return DoubleBoundingBox3D(a, b);
176 }
177 
179 {
180  Transform3D retval = Transform3D::Identity();
181  retval.scale(scale_);
182  return retval;
183 }
184 
186 {
187  Transform3D retval = Transform3D::Identity();
188  retval.translate(translation);
189  return retval;
190 }
191 
193 {
194  Transform3D retval = Transform3D::Identity();
195  retval.rotate(Eigen::AngleAxisd(angle, Vector3D::UnitX()));
196  return retval;
197 }
198 
200 {
201  Transform3D retval = Transform3D::Identity();
202  retval.rotate(Eigen::AngleAxisd(angle, Vector3D::UnitY()));
203  return retval;
204 }
205 
207 {
208  Transform3D retval = Transform3D::Identity();
209  retval.rotate(Eigen::AngleAxisd(angle, Vector3D::UnitZ()));
210  return retval;
211 }
212 
214 {
215  Vector3D k = cross(from, to);
216 
217 
218  // handle special cases
219  if (similar(k.length(), 0.0))
220  {
221  // dot==1 -> point in the same direction
222  if (similar(dot(from, to), 1.0))
223  {
224  return Transform3D::Identity();
225  }
226 
227  // dot==-1 ->point in opposite directions, cross product will fail.
228  // Find an arbitrary vector perpendicular to from, rotate 180 around that one.
229  if (similar(dot(from, to), -1.0))
230  {
231  Vector3D e_x = Vector3D::UnitX();
232  Vector3D kk = cross(from, e_x);
233  if (similar(kk.length(), 0.0))
234  {
235  Vector3D e_y = Vector3D::UnitY();
236  kk = cross(from, e_y);
237  }
238 
239  Transform3D retval = Transform3D::Identity();
240  retval.rotate(Eigen::AngleAxisd(M_PI, kk));
241  return retval;
242  }
243  }
244 
245  double dotnormal = dot(from, to)/from.length()/to.length();
246  double angle = acos(dotnormal);
247 
248  Transform3D retval = Transform3D::Identity();
249  retval.rotate(Eigen::AngleAxisd(angle, k.normal()));
250  return retval;
251 }
252 
253 
255 {
256  // translate input bottomleft to origin, scale, translate back to output bottomleft.
258  Vector3D inrange = in.range();
259  Vector3D outrange = out.range();
260  Vector3D scale;
261 
262  // check for zero input dimensions
263  for (unsigned i = 0; i < scale.size(); ++i)
264  {
265  if (fabs(inrange[i]) < 1.0E-5)
266  scale[i] = 0;
267  else
268  scale[i] = outrange[i] / inrange[i];
269  }
272  Transform3D M = T1 * S * T0;
273  return M;
274 }
275 
276 Transform3D createTransformIJC(const Vector3D& ivec, const Vector3D& jvec, const Vector3D& center)
277 {
278  Transform3D t = Transform3D::Identity();
279  t.matrix().col(0).head(3) = ivec;
280  t.matrix().col(1).head(3) = jvec;
281  t.matrix().col(2).head(3) = cross(ivec, jvec);
282  t.matrix().col(3).head(3) = center;
283  return t;
284 }
285 
286 
288 {
290 }
291 
292 cxResource_EXPORT Transform3D createTransformFromReferenceToExternal(PATIENT_COORDINATE_SYSTEM external)
293 {
294  if (external==pcsRAS)
295  return createTransformLPS2RAS();
296  else
297  return Transform3D::Identity();
298 }
299 
301 {
302  Eigen::IOFormat singleLineFormatting(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", " | ", "", "", "", "");
303  std::ostringstream stream;
304  stream << transform.matrix().format(singleLineFormatting);
305 
306  return stream.str();
307 }
308 
309 
310 } // namespace cx
311 // --------------------------------------------------------
312 
vtkSmartPointer< class vtkMatrix4x4 > vtkMatrix4x4Ptr
Definition: cxMathBase.h:58
pcsRAS
Right-Anterior-Superior, used by Slicer3D, ITK-Snap, nifti, MINC.
Transform3D createTransformRotateY(const double angle)
bool similar(const Transform3D &a, const Transform3D &b, double tol)
Vector3D topRight() const
std::ostream & put(const Eigen::Affine3d *self, std::ostream &s, int indent, char newline)
Transform3D createTransformScale(const Vector3D &scale_)
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
boost::array< double, 16 > flatten(const Eigen::Affine3d *self)
cxResource_EXPORT Transform3D createTransformFromReferenceToExternal(PATIENT_COORDINATE_SYSTEM external)
Transform3D createTransformRotationBetweenVectors(Vector3D from, Vector3D to)
cxResource_EXPORT Transform3D createTransformLPS2RAS()
std::string matrixAsSingleLineString(cx::Transform3D transform)
Vector3D bottomLeft() const
vtkSmartPointer< class vtkTransform > vtkTransformPtr
Definition: cxMathBase.h:62
Transform3D createTransformNormalize(const DoubleBoundingBox3D &in, const DoubleBoundingBox3D &out)
Vector3D cross(const Vector3D &a, const Vector3D &b)
compute cross product of a and b.
Definition: cxVector3D.cpp:62
vtkMatrix4x4Ptr getVtkMatrix(const Eigen::Affine3d *self)
Transform3D createTransformIJC(const Vector3D &ivec, const Vector3D &jvec, const Vector3D &center)
Eigen::Affine3d fromString(const QString &text, bool *_ok)
Transform3D createTransformTranslate(const Vector3D &translation)
double dot(const Vector3D &a, const Vector3D &b)
compute inner product (or dot product) of a and b.
Definition: cxVector3D.cpp:67
Representation of a floating-point bounding box in 3D. The data are stored as {xmin,xmax,ymin,ymax,zmin,zmax}, in order to simplify communication with vtk.
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
void fill(Eigen::Affine3d *self, vtkMatrix4x4Ptr m)
Transform3D createTransformRotateZ(const double angle)
vtkTransformPtr getVtkTransform(const Eigen::Affine3d *self)
Transform3D createTransformRotateX(const double angle)
std::vector< double > convertQString2DoubleVector(const QString &input, bool *ok)
float4 transform(float16 matrix, float4 voxel)
#define M_PI
Namespace for all CustusX production code.