Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
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) 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 
13 #include "cxProbeSector.h"
14 
15 #include "vtkImageData.h"
16 #include <vtkPointData.h>
17 #include <vtkUnsignedCharArray.h>
18 #include <vtkPolyData.h>
19 #include <vtkCellArray.h>
20 #include <vtkFloatArray.h>
21 #include <vtkPolyLine.h>
22 #include <vtkClipPolyData.h>
23 #include <vtkBox.h>
24 #include <vtkPlane.h>
25 #include <vtkPlanes.h>
26 #include <vtkCutter.h>
27 #include <vtkAppendPolyData.h>
28 #include "cxBoundingBox3D.h"
29 #include "cxVolumeHelpers.h"
30 #include "cxUtilHelpers.h"
31 
32 typedef vtkSmartPointer<class vtkPlanes> vtkPlanesPtr;
33 typedef vtkSmartPointer<class vtkPlane> vtkPlanePtr;
34 typedef vtkSmartPointer<class vtkBox> vtkBoxPtr;
35 typedef vtkSmartPointer<class vtkCutter> vtkCutterPtr;
36 typedef vtkSmartPointer<class vtkAppendPolyData> vtkAppendPolyDataPtr;
37 typedef vtkSmartPointer<class vtkFloatArray> vtkFloatArrayPtr;
38 
39 namespace cx
40 {
41 
43 {
44  mPolyData = vtkPolyDataPtr::New();
46 }
47 
49 {
50  mData = data;
51  // this->test();
52 }
53 
58 {
59 public:
61  mData(data), m_vMu(uMv.inv())
62  {
63  mCachedCenter_v = m_vMu.coord(mData.getOrigin_u());
64  mClipRect_v = transform(m_vMu, mData.getClipRect_u());
65  mClipRect_v[4] = -1;
66  mClipRect_v[5] = 1;
67  }
68  bool operator ()(int x, int y) const
69  {
70  Vector3D p_v = multiply_elems(Vector3D(x, y, 0), mData.getSpacing());
71 
72  return this->insideClipRect(p_v) && this->insideSector(p_v);
73  }
74 
75 private:
81  bool insideClipRect(const Vector3D& p_v) const
82  {
83  return mClipRect_v.contains(p_v);
84  }
85 
91  bool insideSector(const Vector3D& p_v) const
92  {
93  Vector3D d = p_v - mCachedCenter_v;
94 
96  {
97  double angle = atan2(d[1], d[0]);
98  angle -= M_PI_2; // center angle on us probe axis at 90*.
99  if (angle < -M_PI)
100  angle += 2.0 * M_PI;
101 
102  if (fabs(angle) > mData.getWidth() / 2.0)
103  return false;
104  if (d.length() < mData.getDepthStart())
105  return false;
106  if (d.length() > mData.getDepthEnd())
107  return false;
108  return true;
109  }
110  else // tLINEAR
111  {
112  if (fabs(d[0]) > mData.getWidth() / 2.0)
113  return false;
114  if (d[1] < mData.getDepthStart())
115  return false;
116  if (d[1] > mData.getDepthEnd())
117  return false;
118  return true;
119  }
120  }
121 
123  Transform3D m_vMu;
124  Vector3D mCachedCenter_v;
125  DoubleBoundingBox3D mClipRect_v;
126 };
127 
132 {
134  return vtkImageDataPtr();
135  InsideMaskFunctor checkInside(mData, this->get_uMv());
136  vtkImageDataPtr retval;
137  retval = generateVtkImageData(Eigen::Array3i(mData.getSize().width(), mData.getSize().height(), 1),
138  mData.getSpacing(), 0);
139 
140  int* dim(retval->GetDimensions());
141  unsigned char* dataPtr = static_cast<unsigned char*> (retval->GetScalarPointer());
142  for (int x = 0; x < dim[0]; x++)
143  for (int y = 0; y < dim[1]; y++)
144  {
145  dataPtr[x + y * dim[0]] = checkInside(x, y) ? 1 : 0;
146  }
147 
148  return retval;
149 }
150 
152 {
153  Transform3D tMu = this->get_tMu();
154  Vector3D e_x(1, 0, 0);
155  Vector3D e_y(0, 1, 0);
156  Vector3D e_z(0, 0, 1);
157 
158  // zero = tMu * mOrigin_u
159  std::cout << "zero = tMu * mOrigin_u, zero: " << tMu.coord(mData.getOrigin_u()) << ", mOrigin_u: "
160  << mData.getOrigin_u() << std::endl;
161 
162  // e_z = tMu * -e_y
163  std::cout << "e_z = tMu * -e_y " << tMu.vector(-e_y) << std::endl;
164 
165  // e_y = tMu * -e_x
166  std::cout << "e_y = tMu * -e_x " << tMu.vector(-e_x) << std::endl;
167 
168  // tMu * e_x
169  std::cout << "tMu * e_x = <0,-1,0>" << tMu.vector(e_x) << std::endl;
170  // tMu * e_y
171  std::cout << "tMu * e_y = <0,0,-1> " << tMu.vector(e_y) << std::endl;
172 
173 }
174 
176 {
179  Transform3D R = (Rx * Rz);
181 
182  Transform3D tMu = R * T;
183  return tMu;
184 }
185 
187 {
188  // use H-1 because we count between pixel centers.
189  double H = (mData.getSize().height() - 1) * mData.getSpacing()[1];
191 }
192 
194 {
195  this->updateSector();
196  return mPolyData;
197 }
198 
203 bool ProbeSector::clipRectIntersectsSector() const
204 {
205  DoubleBoundingBox3D s(mPolyData->GetPoints()->GetBounds());
207 
208  double tol = 1; // assume 1mm tolerance
209  bool outside = ( (c[0] < s[0]) || similar(c[0],s[0], tol) )
210  && ( (s[1] < c[1]) || similar(s[1],c[1], tol) )
211  && ( (c[2] < s[2]) || similar(c[2],s[2], tol) )
212  && ( (s[3] < c[3]) || similar(s[3],c[3], tol) );
213  return !outside;
214 }
215 
217 {
218  this->updateSector();
220  return mPolyData;
221 
222  vtkPolyDataPtr output = vtkPolyDataPtr::New();
223  output->SetPoints(mPolyData->GetPoints());
224  output->SetLines(mPolyData->GetLines());
225 
226  // also display the cliprect
227  vtkAppendPolyDataPtr retval = vtkAppendPolyDataPtr::New();
228  retval->AddInputData(output);
229 
230  if (this->clipRectIntersectsSector())
231  retval->AddInputData(this->getClipRectPolyData());
232 
233  retval->Update();
234  return retval->GetOutput();
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 // output->Update();
248  return output;
249 }
250 
252 {
253  return this->getClipRectPolyData();
254 }
255 
256 vtkPolyDataPtr ProbeSector::getClipRectPolyData()
257 {
258  vtkPointsPtr points = vtkPointsPtr::New();
259  vtkCellArrayPtr sides = vtkCellArrayPtr::New();
260 
261  vtkIdType cells[5] =
262  { 0, 1, 2, 3, 0 };
263  sides->InsertNextCell(5, cells);
264 
266  points->InsertNextPoint(bb.corner(0, 0, 0).begin());
267  points->InsertNextPoint(bb.corner(1, 0, 0).begin());
268  points->InsertNextPoint(bb.corner(1, 1, 0).begin());
269  points->InsertNextPoint(bb.corner(0, 1, 0).begin());
270 
271  vtkPolyDataPtr polydata = vtkPolyDataPtr::New();
272  polydata->SetPoints(points);
273  polydata->SetLines(sides);
274 
275  return polydata;
276 }
277 
282 {
283  vtkPointsPtr points = vtkPointsPtr::New();
284  vtkCellArrayPtr sides = vtkCellArrayPtr::New();
285 
286  vtkIdType cells[4] =
287  { 1, 0, 2, 3 };
288  sides->InsertNextCell(4, cells);
289 
290  Vector3D o_u = mData.getOrigin_u();
291  double length = (mData.getDepthStart() - mData.getDepthEnd())/15;
292  length = constrainValue(length, 2.0, 10.0);
293  Vector3D tip = o_u + Vector3D(0, -length, 0);
294  Vector3D left = o_u + Vector3D(-length/3, 0, 0);
295  Vector3D right = o_u + Vector3D(length/3, 0, 0);
296 
297  points->InsertNextPoint(o_u.begin());
298  points->InsertNextPoint(tip.begin());
299  points->InsertNextPoint(left.begin());
300  points->InsertNextPoint(right.begin());
301 
302  vtkPolyDataPtr polydata = vtkPolyDataPtr::New();
303  polydata->SetPoints(points);
304  polydata->SetLines(sides);
305 
306  return polydata;
307 }
308 
310 {
312  {
313  mPolyData = vtkPolyDataPtr::New();
314  return;
315  }
316 
317  Vector3D bounds = Vector3D(mData.getSize().width() - 1, mData.getSize().height() - 1, 1);
318  bounds = multiply_elems(bounds, mData.getSpacing());
319 
320  vtkFloatArrayPtr newTCoords = vtkFloatArrayPtr::New();
321  newTCoords->SetNumberOfComponents(2);
322 
323  Vector3D p(0, 0, 0); // tool position in local space
324  // first define the shape of the probe in a xy-plane.
325  // after that, transform into yz-plane because thats the tool plane (probe point towards positive z)
326  // then transform to global space.
327  Transform3D tMl = createTransformIJC(Vector3D(0, 1, 0), Vector3D(0, 0, 1), Vector3D(0, 0, 0));
328  Transform3D texMu = createTransformNormalize(DoubleBoundingBox3D(0, bounds[0], 0, bounds[1], 0, 1), DoubleBoundingBox3D(0, 1, 0, 1, 0, 1));
329  Transform3D uMt = this->get_tMu().inv();
330  Transform3D texMl = texMu * uMt * tMl;
331  Transform3D uMl = uMt * tMl;
332 
333  //Transform3D M = tMl;
334  Vector3D e_x = unitVector(0);
335  Vector3D e_y = unitVector(M_PI_2);
336  Vector3D e_z(0, 0, 1);
337 
338  vtkPointsPtr points = vtkPointsPtr::New();
339  vtkCellArrayPtr sides = vtkCellArrayPtr::New();
340  vtkCellArrayPtr strips = vtkCellArrayPtr::New();
341  vtkCellArrayPtr polys = vtkCellArrayPtr::New();
342 
343  DoubleBoundingBox3D bb_u;
344 
346  {
347  Vector3D cr = mData.getDepthStart() * e_y + mData.getWidth() / 2 * e_x;
348  Vector3D cl = mData.getDepthStart() * e_y - mData.getWidth() / 2 * e_x;
349  Vector3D pr = mData.getDepthEnd() * e_y + mData.getWidth() / 2 * e_x;
350  Vector3D pl = mData.getDepthEnd() * e_y - mData.getWidth() / 2 * e_x;
351 
352  points->Allocate(4);
353  points->InsertNextPoint(uMl.coord(cl).begin());
354  points->InsertNextPoint(uMl.coord(cr).begin());
355  points->InsertNextPoint(uMl.coord(pr).begin());
356  points->InsertNextPoint(uMl.coord(pl).begin());
357 
358  newTCoords->Allocate(4);
359  newTCoords->InsertNextTuple(texMl.coord(cl).begin());
360  newTCoords->InsertNextTuple(texMl.coord(cr).begin());
361  newTCoords->InsertNextTuple(texMl.coord(pr).begin());
362  newTCoords->InsertNextTuple(texMl.coord(pl).begin());
363 
364  vtkIdType cells[5] = { 0, 1, 2, 3, 0 };
365  sides->InsertNextCell(5, cells);
366  polys->InsertNextCell(5, cells);
367  vtkIdType s_cells[5] = { 0, 3, 1, 2 };
368  strips->InsertNextCell(4, s_cells);
369  }
370  else if (mData.getType() == ProbeDefinition::tSECTOR)
371  {
372  Vector3D c(0, 0, 0); // arc center point
373  c += mData.getCenterOffset() * e_y; // arc center point
374 
375  int arcRes = 20;//Number of points in arc
376  double angleIncrement = mData.getWidth() / arcRes;
377  double startAngle = M_PI_2 - mData.getWidth() / 2.0;
378  double stopAngle = M_PI_2 + mData.getWidth() / 2.0;
379  int N = 2 * (arcRes + 1); // total number of points
380 
381  points->Allocate(N);
382  newTCoords->Allocate(2 * N);
383 
384  for (int i = 0; i <= arcRes; i++)
385  {
386  double theta = startAngle + i * angleIncrement;
387  Vector3D startTheta;
388  if (mData.getCenterOffset() == 0)
389  {
390  startTheta = c + (mData.getDepthStart()-mData.getCenterOffset()) * unitVector(theta);
391  }
392  else
393  {
394  startTheta = c + (mData.getDepthStart()-mData.getCenterOffset()) * unitVector(theta) / sin(theta);
395  }
396  newTCoords->InsertNextTuple(texMl.coord(startTheta).begin());
397  points->InsertNextPoint(uMl.coord(startTheta).begin());
398  }
399 
400  for (int i = 0; i <= arcRes; i++)
401  {
402  double theta = stopAngle - i * angleIncrement;
403  Vector3D endTheta;
404  double offset = 0;
405 
406  if (mData.getCenterOffset() != 0)
407  {
408  offset = (mData.getDepthStart()-mData.getCenterOffset()) / sin(theta) - (mData.getDepthStart()-mData.getCenterOffset());
409  }
410  endTheta = c + (mData.getDepthEnd()-mData.getCenterOffset()+offset) * unitVector(theta);
411  newTCoords->InsertNextTuple(texMl.coord(endTheta).begin());
412  points->InsertNextPoint(uMl.coord(endTheta).begin());
413  }
414 
415  sides->InsertNextCell(N + 1);
416  for (int i = 0; i < N; i++)
417  sides->InsertCellPoint(i);
418  sides->InsertCellPoint(0);
419 
420  polys->InsertNextCell(N + 1);
421  for (int i = 0; i < arcRes * 2 + 2; i++)
422  polys->InsertCellPoint(i);
423  polys->InsertCellPoint(0);
424 
425  strips->InsertNextCell(N);
426  for (int i = 0; i <= arcRes; ++i)
427  {
428  strips->InsertCellPoint(i);
429  strips->InsertCellPoint(N - 1 - i);
430  }
431  }
432 
433  vtkPolyDataPtr polydata = vtkPolyDataPtr::New();
434  polydata->SetPoints(points);
435  polydata->SetStrips(strips);
436  polydata->GetPointData()->SetTCoords(newTCoords);
437  polydata->SetLines(sides);
438  mPolyData = polydata;
439 }
440 
441 }
double getWidth() const
DoubleBoundingBox3D transform(const Transform3D &m, const DoubleBoundingBox3D &bb)
Transform3D createTransformRotateY(const double angle)
vtkSmartPointer< class vtkPlanes > vtkPlanesPtr
Vector3D corner(int x, int y, int z) const
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:56
vtkSmartPointer< class vtkAppendPolyData > vtkAppendPolyDataPtr
vtkSmartPointer< class vtkFloatArray > vtkFloatArrayPtr
vtkSmartPointer< vtkPoints > vtkPointsPtr
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)
double getDepthEnd() const
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)
vtkSmartPointer< vtkPolyData > vtkPolyDataPtr
ProbeDefinition mData
Definition: cxProbeSector.h:54
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
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:31
Vector3D getSpacing() const
bool similar(const CameraInfo &lhs, const CameraInfo &rhs, double tol)
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
void setData(ProbeDefinition data)
Namespace for all CustusX production code.