Fraxinus  16.5.0-fx-rc9
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cxProbeSector.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 
34 #include "cxProbeSector.h"
35 
36 #include "vtkImageData.h"
37 #include <vtkPointData.h>
38 #include <vtkUnsignedCharArray.h>
39 #include <vtkPolyData.h>
40 #include <vtkCellArray.h>
41 #include <vtkFloatArray.h>
42 #include <vtkPolyLine.h>
43 #include <vtkClipPolyData.h>
44 #include <vtkBox.h>
45 #include <vtkPlane.h>
46 #include <vtkPlanes.h>
47 #include <vtkCutter.h>
48 #include <vtkAppendPolyData.h>
49 #include "cxBoundingBox3D.h"
50 #include "cxVolumeHelpers.h"
51 #include "cxUtilHelpers.h"
52 
53 typedef vtkSmartPointer<class vtkPlanes> vtkPlanesPtr;
54 typedef vtkSmartPointer<class vtkPlane> vtkPlanePtr;
55 typedef vtkSmartPointer<class vtkBox> vtkBoxPtr;
56 typedef vtkSmartPointer<class vtkCutter> vtkCutterPtr;
57 typedef vtkSmartPointer<class vtkAppendPolyData> vtkAppendPolyDataPtr;
58 typedef vtkSmartPointer<class vtkFloatArray> vtkFloatArrayPtr;
59 
60 namespace cx
61 {
62 
64 {
65  mPolyData = vtkPolyDataPtr::New();
67 }
68 
70 {
71  mData = data;
72  // this->test();
73 }
74 
79 {
80 public:
82  mData(data), m_vMu(uMv.inv())
83  {
84  mCachedCenter_v = m_vMu.coord(mData.getOrigin_u());
85  mClipRect_v = transform(m_vMu, mData.getClipRect_u());
86  mClipRect_v[4] = -1;
87  mClipRect_v[5] = 1;
88  }
89  bool operator ()(int x, int y) const
90  {
91  Vector3D p_v = multiply_elems(Vector3D(x, y, 0), mData.getSpacing());
92 
93  return this->insideClipRect(p_v) && this->insideSector(p_v);
94  }
95 
96 private:
102  bool insideClipRect(const Vector3D& p_v) const
103  {
104  return mClipRect_v.contains(p_v);
105  }
106 
112  bool insideSector(const Vector3D& p_v) const
113  {
114  Vector3D d = p_v - mCachedCenter_v;
115 
116  if (mData.getType() == ProbeDefinition::tSECTOR)
117  {
118  double angle = atan2(d[1], d[0]);
119  angle -= M_PI_2; // center angle on us probe axis at 90*.
120  if (angle < -M_PI)
121  angle += 2.0 * M_PI;
122 
123  if (fabs(angle) > mData.getWidth() / 2.0)
124  return false;
125  if (d.length() < mData.getDepthStart())
126  return false;
127  if (d.length() > mData.getDepthEnd())
128  return false;
129  return true;
130  }
131  else // tLINEAR
132  {
133  if (fabs(d[0]) > mData.getWidth() / 2.0)
134  return false;
135  if (d[1] < mData.getDepthStart())
136  return false;
137  if (d[1] > mData.getDepthEnd())
138  return false;
139  return true;
140  }
141  }
142 
143  ProbeDefinition mData;
144  Transform3D m_vMu;
145  Vector3D mCachedCenter_v;
146  DoubleBoundingBox3D mClipRect_v;
147 };
148 
153 {
155  return vtkImageDataPtr();
156  InsideMaskFunctor checkInside(mData, this->get_uMv());
157  vtkImageDataPtr retval;
158  retval = generateVtkImageData(Eigen::Array3i(mData.getSize().width(), mData.getSize().height(), 1),
159  mData.getSpacing(), 0);
160 
161  int* dim(retval->GetDimensions());
162  unsigned char* dataPtr = static_cast<unsigned char*> (retval->GetScalarPointer());
163  for (int x = 0; x < dim[0]; x++)
164  for (int y = 0; y < dim[1]; y++)
165  {
166  dataPtr[x + y * dim[0]] = checkInside(x, y) ? 1 : 0;
167  }
168 
169  return retval;
170 }
171 
173 {
174  Transform3D tMu = this->get_tMu();
175  Vector3D e_x(1, 0, 0);
176  Vector3D e_y(0, 1, 0);
177  Vector3D e_z(0, 0, 1);
178 
179  // zero = tMu * mOrigin_u
180  std::cout << "zero = tMu * mOrigin_u, zero: " << tMu.coord(mData.getOrigin_u()) << ", mOrigin_u: "
181  << mData.getOrigin_u() << std::endl;
182 
183  // e_z = tMu * -e_y
184  std::cout << "e_z = tMu * -e_y " << tMu.vector(-e_y) << std::endl;
185 
186  // e_y = tMu * -e_x
187  std::cout << "e_y = tMu * -e_x " << tMu.vector(-e_x) << std::endl;
188 
189  // tMu * e_x
190  std::cout << "tMu * e_x = <0,-1,0>" << tMu.vector(e_x) << std::endl;
191  // tMu * e_y
192  std::cout << "tMu * e_y = <0,0,-1> " << tMu.vector(e_y) << std::endl;
193 
194 }
195 
197 {
200  Transform3D R = (Rx * Rz);
202 
203  Transform3D tMu = R * T;
204  return tMu;
205 }
206 
208 {
209  // use H-1 because we count between pixel centers.
210  double H = (mData.getSize().height() - 1) * mData.getSpacing()[1];
212 }
213 
215 {
216  this->updateSector();
217  return mPolyData;
218 }
219 
224 bool ProbeSector::clipRectIntersectsSector() const
225 {
226  DoubleBoundingBox3D s(mPolyData->GetPoints()->GetBounds());
228 
229  double tol = 1; // assume 1mm tolerance
230  bool outside = ( (c[0] < s[0]) || similar(c[0],s[0], tol) )
231  && ( (s[1] < c[1]) || similar(s[1],c[1], tol) )
232  && ( (c[2] < s[2]) || similar(c[2],s[2], tol) )
233  && ( (s[3] < c[3]) || similar(s[3],c[3], tol) );
234  return !outside;
235 }
236 
238 {
239  this->updateSector();
241  return mPolyData;
242 
243  vtkPolyDataPtr output = vtkPolyDataPtr::New();
244  output->SetPoints(mPolyData->GetPoints());
245  output->SetLines(mPolyData->GetLines());
246 
247  // also display the cliprect
248  vtkAppendPolyDataPtr retval = vtkAppendPolyDataPtr::New();
249  retval->AddInputData(output);
250 
251  if (this->clipRectIntersectsSector())
252  retval->AddInputData(this->getClipRectPolyData());
253 
254  retval->Update();
255  return retval->GetOutput();
256 }
257 
259 {
260  this->updateSector();
262  return mPolyData;
263 
264  vtkPolyDataPtr output = vtkPolyDataPtr::New();
265  output->SetPoints(mPolyData->GetPoints());
266  output->SetLines(mPolyData->GetLines());
267 
268 // output->Update();
269  return output;
270 }
271 
273 {
274  return this->getClipRectPolyData();
275 }
276 
277 vtkPolyDataPtr ProbeSector::getClipRectPolyData()
278 {
279  vtkPointsPtr points = vtkPointsPtr::New();
280  vtkCellArrayPtr sides = vtkCellArrayPtr::New();
281 
282  vtkIdType cells[5] =
283  { 0, 1, 2, 3, 0 };
284  sides->InsertNextCell(5, cells);
285 
286  DoubleBoundingBox3D bb = mData.getClipRect_u();
287  points->InsertNextPoint(bb.corner(0, 0, 0).begin());
288  points->InsertNextPoint(bb.corner(1, 0, 0).begin());
289  points->InsertNextPoint(bb.corner(1, 1, 0).begin());
290  points->InsertNextPoint(bb.corner(0, 1, 0).begin());
291 
292  vtkPolyDataPtr polydata = vtkPolyDataPtr::New();
293  polydata->SetPoints(points);
294  polydata->SetLines(sides);
295 
296  return polydata;
297 }
298 
303 {
304  vtkPointsPtr points = vtkPointsPtr::New();
305  vtkCellArrayPtr sides = vtkCellArrayPtr::New();
306 
307  vtkIdType cells[4] =
308  { 1, 0, 2, 3 };
309  sides->InsertNextCell(4, cells);
310 
311  Vector3D o_u = mData.getOrigin_u();
312  double length = (mData.getDepthStart() - mData.getDepthEnd())/15;
313  length = constrainValue(length, 2, 10);
314  Vector3D tip = o_u + Vector3D(0, -length, 0);
315  Vector3D left = o_u + Vector3D(-length/3, 0, 0);
316  Vector3D right = o_u + Vector3D(length/3, 0, 0);
317 
318  points->InsertNextPoint(o_u.begin());
319  points->InsertNextPoint(tip.begin());
320  points->InsertNextPoint(left.begin());
321  points->InsertNextPoint(right.begin());
322 
323  vtkPolyDataPtr polydata = vtkPolyDataPtr::New();
324  polydata->SetPoints(points);
325  polydata->SetLines(sides);
326 
327  return polydata;
328 }
329 
331 {
333  {
334  mPolyData = vtkPolyDataPtr::New();
335  return;
336  }
337 
338  Vector3D bounds = Vector3D(mData.getSize().width() - 1, mData.getSize().height() - 1, 1);
339  bounds = multiply_elems(bounds, mData.getSpacing());
340 
341  vtkFloatArrayPtr newTCoords = vtkFloatArrayPtr::New();
342  newTCoords->SetNumberOfComponents(2);
343 
344  Vector3D p(0, 0, 0); // tool position in local space
345  // first define the shape of the probe in a xy-plane.
346  // after that, transform into yz-plane because thats the tool plane (probe point towards positive z)
347  // then transform to global space.
348  Transform3D tMl = createTransformIJC(Vector3D(0, 1, 0), Vector3D(0, 0, 1), Vector3D(0, 0, 0));
349  Transform3D texMu = createTransformNormalize(DoubleBoundingBox3D(0, bounds[0], 0, bounds[1], 0, 1), DoubleBoundingBox3D(0, 1, 0, 1, 0, 1));
350  Transform3D uMt = this->get_tMu().inv();
351  Transform3D texMl = texMu * uMt * tMl;
352  Transform3D uMl = uMt * tMl;
353 
354  //Transform3D M = tMl;
355  Vector3D e_x = unitVector(0);
356  Vector3D e_y = unitVector(M_PI_2);
357  Vector3D e_z(0, 0, 1);
358 
359  vtkPointsPtr points = vtkPointsPtr::New();
360  vtkCellArrayPtr sides = vtkCellArrayPtr::New();
361  vtkCellArrayPtr strips = vtkCellArrayPtr::New();
362  vtkCellArrayPtr polys = vtkCellArrayPtr::New();
363 
364  DoubleBoundingBox3D bb_u;
365 
367  {
368  Vector3D cr = mData.getDepthStart() * e_y + mData.getWidth() / 2 * e_x;
369  Vector3D cl = mData.getDepthStart() * e_y - mData.getWidth() / 2 * e_x;
370  Vector3D pr = mData.getDepthEnd() * e_y + mData.getWidth() / 2 * e_x;
371  Vector3D pl = mData.getDepthEnd() * e_y - mData.getWidth() / 2 * e_x;
372 
373  points->Allocate(4);
374  points->InsertNextPoint(uMl.coord(cl).begin());
375  points->InsertNextPoint(uMl.coord(cr).begin());
376  points->InsertNextPoint(uMl.coord(pr).begin());
377  points->InsertNextPoint(uMl.coord(pl).begin());
378 
379  newTCoords->Allocate(4);
380  newTCoords->InsertNextTuple(texMl.coord(cl).begin());
381  newTCoords->InsertNextTuple(texMl.coord(cr).begin());
382  newTCoords->InsertNextTuple(texMl.coord(pr).begin());
383  newTCoords->InsertNextTuple(texMl.coord(pl).begin());
384 
385  vtkIdType cells[5] = { 0, 1, 2, 3, 0 };
386  sides->InsertNextCell(5, cells);
387  polys->InsertNextCell(5, cells);
388  vtkIdType s_cells[5] = { 0, 3, 1, 2 };
389  strips->InsertNextCell(4, s_cells);
390  }
391  else if (mData.getType() == ProbeDefinition::tSECTOR)
392  {
393  Vector3D c(0, 0, 0); // arc center point
394  c += mData.getCenterOffset() * e_y; // arc center point
395 
396  int arcRes = 20;//Number of points in arc
397  double angleIncrement = mData.getWidth() / arcRes;
398  double startAngle = M_PI_2 - mData.getWidth() / 2.0;
399  double stopAngle = M_PI_2 + mData.getWidth() / 2.0;
400  int N = 2 * (arcRes + 1); // total number of points
401 
402  points->Allocate(N);
403  newTCoords->Allocate(2 * N);
404 
405  for (int i = 0; i <= arcRes; i++)
406  {
407  double theta = startAngle + i * angleIncrement;
408  Vector3D startTheta;
409  if (mData.getCenterOffset() == 0)
410  {
411  startTheta = c + (mData.getDepthStart()-mData.getCenterOffset()) * unitVector(theta);
412  }
413  else
414  {
415  startTheta = c + (mData.getDepthStart()-mData.getCenterOffset()) * unitVector(theta) / sin(theta);
416  }
417  newTCoords->InsertNextTuple(texMl.coord(startTheta).begin());
418  points->InsertNextPoint(uMl.coord(startTheta).begin());
419  }
420 
421  for (int i = 0; i <= arcRes; i++)
422  {
423  double theta = stopAngle - i * angleIncrement;
424  Vector3D endTheta;
425  double offset = 0;
426 
427  if (mData.getCenterOffset() != 0)
428  {
429  offset = (mData.getDepthStart()-mData.getCenterOffset()) / sin(theta) - (mData.getDepthStart()-mData.getCenterOffset());
430  }
431  endTheta = c + (mData.getDepthEnd()-mData.getCenterOffset()+offset) * unitVector(theta);
432  newTCoords->InsertNextTuple(texMl.coord(endTheta).begin());
433  points->InsertNextPoint(uMl.coord(endTheta).begin());
434  }
435 
436  sides->InsertNextCell(N + 1);
437  for (int i = 0; i < N; i++)
438  sides->InsertCellPoint(i);
439  sides->InsertCellPoint(0);
440 
441  polys->InsertNextCell(N + 1);
442  for (int i = 0; i < arcRes * 2 + 2; i++)
443  polys->InsertCellPoint(i);
444  polys->InsertCellPoint(0);
445 
446  strips->InsertNextCell(N);
447  for (int i = 0; i <= arcRes; ++i)
448  {
449  strips->InsertCellPoint(i);
450  strips->InsertCellPoint(N - 1 - i);
451  }
452  }
453 
454  vtkPolyDataPtr polydata = vtkPolyDataPtr::New();
455  polydata->SetPoints(points);
456  polydata->SetStrips(strips);
457  polydata->GetPointData()->SetTCoords(newTCoords);
458  polydata->SetLines(sides);
459  mPolyData = polydata;
460 }
461 
462 }
double getWidth() const
DoubleBoundingBox3D transform(const Transform3D &m, const DoubleBoundingBox3D &bb)
Transform3D createTransformRotateY(const double angle)
vtkSmartPointer< class vtkPlanes > vtkPlanesPtr
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
vtkPolyDataPtr getSector()
get a polydata representation of the us sector
vtkSmartPointer< class vtkCellArray > vtkCellArrayPtr
double getCenterOffset() const
US beam is emitted straight forward.
Vector3D unitVector(double thetaXY, double thetaZ)
compute a unit vector given angles xy in the xy plane and z meaning the elevation from the xy plane...
Definition: cxVector3D.cpp:77
vtkSmartPointer< class vtkAppendPolyData > vtkAppendPolyDataPtr
vtkSmartPointer< class vtkFloatArray > vtkFloatArrayPtr
vtkPolyDataPtr getSectorLinesOnly()
get a polydata representation of the us sector
vtkSmartPointer< class vtkBox > vtkBoxPtr
double getDepthStart() const
Transform3D createTransformNormalize(const DoubleBoundingBox3D &in, const DoubleBoundingBox3D &out)
bool similar(const DoubleBoundingBox3D &a, const DoubleBoundingBox3D &b, double tol)
double getDepthEnd() const
vtkSmartPointer< class vtkPolyData > vtkPolyDataPtr
Transform3D createTransformIJC(const Vector3D &ivec, const Vector3D &jvec, const Vector3D &center)
Transform3D get_uMv() const
get transform from inverted image space v (origin in ul corner) to image space u. ...
InsideMaskFunctor(ProbeDefinition data, Transform3D uMv)
US beam is emitted radially in a flat cone.
vtkPolyDataPtr getSectorSectorOnlyLinesOnly()
get a polydata representation of the us sector
double constrainValue(double val, double min, double max)
DoubleBoundingBox3D getClipRect_u() const
sector clipping rect in image space u. (lower-left corner origin)
Transform3D createTransformTranslate(const Vector3D &translation)
vtkPolyDataPtr getClipRectLinesOnly()
get a polydata representation of the us clip rect
Transform3D get_tMu() const
get transform from image space u to probe tool space t.
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.
vtkImageDataPtr generateVtkImageData(Eigen::Array3i dim, Vector3D spacing, const unsigned char initValue, int components)
vtkSmartPointer< class vtkCutter > vtkCutterPtr
Vector3D getOrigin_u() const
probe origin in image space u. (lower-left corner origin)
ProbeDefinition mData
Definition: cxProbeSector.h:75
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
Definition of characteristics for an Ultrasound Probe Sector.
Vector3D multiply_elems(const Vector3D &a, const Vector3D &b)
perform element-wise multiplication of a and b.
Definition: cxVector3D.cpp:52
Vector3D getSpacing() const
bool operator()(int x, int y) const
vtkImageDataPtr getMask()
RealScalar length() const
Transform3D createTransformRotateX(const double angle)
vtkPolyDataPtr getOriginPolyData()
get a polydata representation of the origin
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
vtkSmartPointer< class vtkPlane > vtkPlanePtr
#define M_PI
vtkSmartPointer< class vtkPoints > vtkPointsPtr
bool contains(const Vector3D &p) const
void setData(ProbeDefinition data)