NorMIT-nav  2023.01.05-dev+develop.0da12
An IGT application
cxWirePhantomWidget.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 "cxWirePhantomWidget.h"
13 
14 #include <QPushButton>
15 #include <QTextEdit>
16 #include <vtkPolyData.h>
17 #include <vtkPoints.h>
18 #include "cxMesh.h"
19 #include "cxDataLocations.h"
20 #include "vesselReg/SeansVesselReg.hxx"
21 #include "cxPointMetric.h"
22 #include "cxDistanceMetric.h"
23 #include "cxTool.h"
24 #include "cxTypeConversions.h"
25 #include "cxPipelineWidget.h"
26 #include "cxSmoothingImageFilter.h"
27 #include "cxBinaryThinningImageFilter3DFilter.h"
28 #include "cxBinaryThresholdImageFilter.h"
29 #include "cxSpaceProvider.h"
30 #include "cxReporter.h"
31 #include "cxRegistrationService.h"
32 #include "cxRegistrationTransform.h"
33 #include "cxPatientModelService.h"
34 #include "cxViewGroupData.h"
35 
36 #include "cxProfile.h"
37 #include "cxViewService.h"
38 #include "cxCompositeTimedAlgorithm.h"
39 #include "cxUSReconstructInputData.h"
40 #include "cxTrackingService.h"
41 #include "cxUsReconstructionServiceProxy.h"
42 
43 namespace cx
44 {
45 
46 WirePhantomWidget::WirePhantomWidget(ctkPluginContext *pluginContext, QWidget* parent) :
47  RegistrationBaseWidget(RegServices::create(pluginContext), parent, "WirePhantomWidget", "Wire Phantom Test"),
48  mUsReconstructionService(new UsReconstructionServiceProxy(pluginContext))
49 {
50  mLastRegistration = Transform3D::Identity();
51 
52  // fill the pipeline with filters:
53  mPipeline.reset(new Pipeline(mServices->patient()));
54  XmlOptionFile options = profile()->getXmlSettings().descend("registration").descend("WirePhantomWidget");
55  FilterGroupPtr filters(new FilterGroup(options));
56  filters->append(FilterPtr(new SmoothingImageFilter(mServices)));
57  filters->append(FilterPtr(new BinaryThresholdImageFilter(mServices)));
58  filters->append(FilterPtr(new BinaryThinningImageFilter3DFilter(mServices)));
59  mPipeline->initialize(filters);
60 
61  mPipeline->getNodes()[0]->setValueName("US Image:");
62  mPipeline->getNodes()[0]->setHelp("Select an US volume acquired from the wire phantom.");
63  mPipeline->setOption("Color", QVariant(QColor("red")));
64 
65  mLayout = new QVBoxLayout(this);
66 
67  mPipelineWidget = new PipelineWidget(mServices->view(), mServices->patient(), NULL, mPipeline);
68 
69  QPushButton* showCrossButton = new QPushButton("Show Cross");
70  showCrossButton->setToolTip("Load the centerline description of the wire cross");
71  connect(showCrossButton, SIGNAL(clicked()), this, SLOT(loadNominalCross()));
72 
73  mMeasureButton = new QPushButton("Execute");
74  mMeasureButton->setToolTip("Measure deviation of input volume from nominal wire cross.");
75  connect(mMeasureButton, SIGNAL(clicked()), this, SLOT(measureSlot()));
76 
77  mCalibrationButton = new QPushButton("Probe calib");
78  mCalibrationButton->setToolTip(""
79  "Estimate probe calibration sMt\n"
80  "based on last accuracy test and\n"
81  "the current probe orientation");
82  connect(mCalibrationButton, SIGNAL(clicked()), this, SLOT(generate_sMt()));
83 
84  QLayout* buttonsLayout = new QHBoxLayout;
85  buttonsLayout->addWidget(mMeasureButton);
86  buttonsLayout->addWidget(mCalibrationButton);
87 
88  mResults = new QTextEdit;
89 
90  mLayout->addWidget(showCrossButton);
91  mLayout->addWidget(mPipelineWidget);
92  mLayout->addLayout(buttonsLayout);
93  mLayout->addWidget(mResults, 1);
94  mLayout->addStretch();
95 }
96 
97 WirePhantomWidget::~WirePhantomWidget()
98 {
99 }
100 
101 QString WirePhantomWidget::defaultWhatsThis() const
102 {
103  return "<html>"
104  "<h3>Measure US accuracy using the wire phantom.</h3>"
105  "<p>"
106  "Select a US recording of the wire phantom, then press"
107  "the register button to find deviation from the nominal"
108  "data. The output is given as a distance from measured"
109  "to nominal."
110  "</p>"
111  "</html>";
112 }
113 
114 MeshPtr WirePhantomWidget::loadNominalCross()
115 {
116  QString nominalCrossFilename = DataLocations::getRootConfigPath()+"/models/wire_phantom_cross_pts.vtk";
117  MeshPtr retval;
118 
119  std::map<QString, MeshPtr> meshes = mServices->patient()->getDataOfType<Mesh>();
120  for (std::map<QString, MeshPtr>::iterator iter=meshes.begin(); iter!=meshes.end(); ++iter)
121  if (iter->first.contains("wire_phantom_cross_pts"))
122  retval = iter->second;
123 
124  if (!retval)
125  {
126  QString infoText;
127  retval = boost::dynamic_pointer_cast<Mesh>(mServices->patient()->importData(nominalCrossFilename, infoText));
128  }
129 
130  if (!retval)
131  {
132  reportError(QString("failed to load %s.").arg(nominalCrossFilename));
133  }
134 
135  retval->setColor(QColor("green"));
136  this->showData(retval);
137 
138  return retval;
139 }
140 
141 void WirePhantomWidget::showData(DataPtr data)
142 {
143  mServices->view()->getGroup(0)->addData(data->getUid());
144 }
145 
146 void WirePhantomWidget::measureSlot()
147 {
148  if (mPipeline->getPipelineTimedAlgorithm()->isRunning())
149  return;
150  connect(mPipeline->getPipelineTimedAlgorithm().get(), SIGNAL(finished()), this, SLOT(registration()));
151  mPipeline->execute();
152 }
153 
156 Vector3D WirePhantomWidget::findCentroid(MeshPtr mesh)
157 {
158  vtkPolyDataPtr poly = mesh->getVtkPolyData();
159  vtkPointsPtr points = poly->GetPoints();
160  int N = points->GetNumberOfPoints();
161  if (N==0)
162  return Vector3D(0,0,0);
163 
164  Vector3D acc(0,0,0);
165  for (int i = 0; i < N; ++i)
166  acc += Vector3D(points->GetPoint(i));
167  return acc/N;
168 }
169 
170 void WirePhantomWidget::registration()
171 {
172 // // disconnect the connection that started this registration
173  disconnect(mPipeline->getPipelineTimedAlgorithm().get(), SIGNAL(finished()), this, SLOT(registration()));
174 
175  // Verify that a centerline is available:
176  // - if no centerline, run centerline algo on segmentation,
177  // - if no segmentation, run segmentation
178  //
179  // Set centerline as Registration: moving data.
180  // Load nominal cross
181  // find center in nominal
182  // Execute V2V algo
183  // Find transform matrix fMm
184  // apply fMm to nominal, return diff
185  //
186 
187  DataPtr measuredCross = mPipeline->getNodes().back()->getData();
188  MeshPtr nominalCross = this->loadNominalCross();
189  if (!nominalCross || !measuredCross)
190  {
191  reportError("Missing fixed/moving data. WirePhantom measurement failed.");
192  return;
193  }
194 
195  mServices->registration()->setFixedData(nominalCross);
196  mServices->registration()->setMovingData(measuredCross);
197 
198  this->showData(nominalCross);
199  this->showData(measuredCross);
200 
201  SeansVesselReg vesselReg;
202  vesselReg.mt_auto_lts = true;
203  vesselReg.mt_ltsRatio = 80;
204  vesselReg.mt_doOnlyLinear = true;
205  QString logPath = mServices->patient()->getActivePatientFolder() + "/Logs/";
206 // vesselReg.setDebugOutput(true);
207 
208  bool success = vesselReg.initialize(mServices->registration()->getMovingData(), mServices->registration()->getFixedData(), logPath);
209  success = success && vesselReg.execute();
210  if (!success)
211  {
212  reportWarning("Vessel registration failed.");
213  return;
214  }
215 
216  Transform3D linearTransform = vesselReg.getLinearResult();
217  std::cout << "v2v linear result:\n" << linearTransform << std::endl;
218 
219  vesselReg.checkQuality(linearTransform);
220  mLastRegistration = linearTransform;
221 
222  // The registration is performed in space r. Thus, given an old data position rMd, we find the
223  // new one as rM'd = Q * rMd, where Q is the inverted registration output.
224  // Delta is thus equal to Q:
225  Transform3D delta = linearTransform.inv();
226 // mRegistrationService->restart();
227  mServices->registration()->setLastRegistrationTime(QDateTime::currentDateTime());//Instead of restart
228  mServices->registration()->addImage2ImageRegistration(delta, "Wire Phantom Measurement");
229 
230 
231  Vector3D t_delta = linearTransform.matrix().block<3, 1>(0, 3);
232  Eigen::AngleAxisd angleAxis = Eigen::AngleAxisd(linearTransform.matrix().block<3, 3>(0, 0));
233  double angle = angleAxis.angle();
234 
235  // Compute the centroid of the wire phantom.
236  // This should be the wire centre, given that the
237  // model is symmetrical.
238  Vector3D cross_r = this->findCentroid(nominalCross);
239  // Vector3D cross(134.25, 134.25, 99.50); // old hardcoded value: obsole after new measurements
240  // should be (CA20121022): <134.212 134.338 100.14>
241 // std::cout << "cross2 " << cross2 << std::endl;
242 
243 
244  // find transform to probe space t_us, i.e. the middle position from the us acquisition
245  std::pair<QString, Transform3D> probePos = this->getLastProbePosition();
246  Transform3D rMt_us = probePos.second;
247 
248  Vector3D cross_moving_r = linearTransform.coord(cross_r);
249  Vector3D diff_r = (cross_r - cross_moving_r);
250  Vector3D diff_tus = rMt_us.inv().vector(diff_r);
251 
252  struct Fmt
253  {
254  QString operator()(double val)
255  {
256  return QString("%1").arg(val, 2, 'f', 2);
257  }
258  };
259  Fmt fmt;
260 
261  QString result;
262  result += QString("Results for: %1\n").arg(mServices->registration()->getMovingData()->getName());
263  result += QString("Shift vector (r): \t%1\t%2\t%3\n").arg(fmt(diff_r[0])).arg(fmt(diff_r[1])).arg(fmt(diff_r[2]));
264  if (!probePos.first.isEmpty())
265  result += QString("Shift vector (probe): \t%1\t%2\t%3\t(used tracking data from %4)\n").arg(fmt(diff_tus[0])).arg(fmt(diff_tus[1])).arg(fmt(diff_tus[2])).arg(probePos.first);
266  result += QString("Accuracy |v|: \t%1\tmm\n").arg(fmt(diff_r.length()));
267  result += QString("Angle: \t%1\t*\n").arg(fmt(angle / M_PI * 180.0));
268 
269  mResults->append(result);
270  report("Wire Phantom Test Results:\n"+result);
271 
272  this->showDataMetrics(cross_r);
273 }
274 
280 void WirePhantomWidget::showDataMetrics(Vector3D cross_r)
281 {
282  // add metrics displaying the distance from cross in the nominal and us spaces:
283  Transform3D usMnom = mServices->spaceProvider()->get_toMfrom(
284  mServices->spaceProvider()->getD(mServices->registration()->getFixedData()),
285  mServices->spaceProvider()->getD(mServices->registration()->getMovingData()));
286  Vector3D cross_us = usMnom.coord(cross_r);
287 
288  PointMetricPtr p1 = boost::dynamic_pointer_cast<PointMetric>(mServices->patient()->getData("cross_nominal"));
289  if (!p1)
290  p1 = mServices->patient()->createSpecificData<PointMetric>("cross_nominal");
291 // p1 = PointMetric::create("cross_nominal", "cross_nominal");
292  p1->get_rMd_History()->setParentSpace(mServices->registration()->getFixedData()->getUid());
293  p1->setSpace(mServices->spaceProvider()->getD(mServices->registration()->getFixedData()));
294  p1->setCoordinate(cross_r);
295 // dataManager()->loadData(p1);
296  mServices->patient()->insertData(p1);
297  //this->showData(p1);
298 
299  PointMetricPtr p2 = mServices->patient()->getData<PointMetric>("cross_us");
300  if (!p2)
301  p2 = mServices->patient()->createSpecificData<PointMetric>("cross_us");
302 // p2 = PointMetric::create("cross_us", "cross_us");
303  p2->get_rMd_History()->setParentSpace(mServices->registration()->getMovingData()->getUid());
304  p2->setSpace(mServices->spaceProvider()->getD(mServices->registration()->getMovingData()));
305  p2->setCoordinate(cross_us);
306 // dataManager()->loadData(p2);
307  mServices->patient()->insertData(p2);
308  //this->showData(p2);
309 
310  DistanceMetricPtr d0 = mServices->patient()->getData<DistanceMetric>("accuracy");
311  if (!d0)
312  d0 = mServices->patient()->createSpecificData<DistanceMetric>("accuracy");
313 // d0.reset(new DistanceMetric("accuracy", "accuracy"));
314  d0->get_rMd_History()->setParentSpace("reference");
315  d0->getArguments()->set(0, p1);
316  d0->getArguments()->set(1, p2);
317 // dataManager()->loadData(d0);
318  mServices->patient()->insertData(d0);
319  this->showData(d0);
320 }
321 
328 std::pair<QString, Transform3D> WirePhantomWidget::getLastProbePosition()
329 {
330  // find transform to probe space t_us, i.e. the middle position from the us acquisition
331  USReconstructInputData usData = mUsReconstructionService->getSelectedFileData();
332  Transform3D prMt_us = Transform3D::Identity();
333  if (usData.mPositions.empty())
334  return std::make_pair("", Transform3D::Identity());
335  prMt_us = usData.mPositions[usData.mPositions.size()/2].mPos;
336  return std::make_pair(usData.mFilename, prMt_us);
337 }
338 
339 void WirePhantomWidget::generate_sMt()
340 {
341  bool translateOnly = true;
342 
343 
344  std::pair<QString, Transform3D> probePos = this->getLastProbePosition();
345  Transform3D rMt_us = probePos.second;
346  if (probePos.first.isEmpty())
347  {
348  reportWarning("Cannot find probe position from last recording, aborting calibration test.");
349  return;
350  }
351 
352  USReconstructInputData usData = mUsReconstructionService->getSelectedFileData();
353  ToolPtr probe = mServices->tracking()->getTool(usData.mProbeUid);
354  if (!probe || !probe->hasType(Tool::TOOL_US_PROBE))
355  {
356  reportWarning("Cannot find probe, aborting calibration test.");
357  return;
358  }
359 
360 // ToolPtr probe = toolManager()->getDominantTool();
361 // if (!probe || !probe->getVisible() || !probe->hasType(Tool::TOOL_US_PROBE))
362 // {
363 // reportWarning("Cannot find visible probe, aborting calibration test.");
364 // return;
365 // }
366  if (!mServices->registration()->getMovingData())
367  {
368  reportWarning("Cannot find moving data, aborting calibration test.");
369  return;
370  }
371 
372 
373  if (translateOnly)
374  {
375  Vector3D cross_r = this->findCentroid(this->loadNominalCross());
376  // Vector3D cross(134.25, 134.25, 99.50); // old hardcoded value: obsole after new measurements
377  // should be (CA20121022): <134.212 134.338 100.14>
378  // std::cout << "cross2 " << cross2 << std::endl;
379 
380 
381  // find transform to probe space t_us, i.e. the middle position from the us acquisition
382 // std::pair<QString, Transform3D> probePos = this->getLastProbePosition();
383 // Transform3D rMt_us = probePos.second;
384 
385  Vector3D cross_moving_r = mLastRegistration.coord(cross_r);
386  Vector3D diff_r = (cross_r - cross_moving_r);
387  Vector3D diff_tus = rMt_us.inv().vector(diff_r);
388 
389  Transform3D sMt = probe->getCalibration_sMt();
390  Transform3D sQt; // Q is the new calibration matrix.
391 
392  sQt = sMt * createTransformTranslate(diff_tus);
393 
394  report(QString(""
395  "Calculated new calibration matrix\n"
396  "adding only translation "
397  "from last accuracy test\n"
398  "and raw data %1:\n"
399  "%2").arg(probePos.first).arg(qstring_cast(sQt)));
400  }
401  else
402  {
403  Transform3D prMt = rMt_us;
404  Transform3D sMt = probe->getCalibration_sMt();
405  Transform3D prMs = prMt * sMt.inv();
406  Transform3D usMpr = mServices->registration()->getMovingData()->get_rMd().inv() * mServices->patient()->get_rMpr();
407  Transform3D nomMus = mLastRegistration.inv();
408 
409  Transform3D sQt; // Q is the new calibration matrix.
410  Transform3D usMs = usMpr*prMs;
411 
412  // start with: nomMus * usMpr * prMs * sMt
413  // move usMpr*prMs to the left and collect a new sMt from the remains:
414  // start with: usMpr * prMs * (usMpr * prMs).inv() * nomMus * usMpr * prMs * sMt
415 
416  sQt = usMs.inv() * nomMus * usMs * sMt;
417 
418  report(QString(""
419  "Calculated new calibration matrix\n"
420  "from last accuracy test\n"
421  "and raw data %1.\n"
422  "Old calibration matrix sMt:\n"
423  "%2\n"
424  "New calibration matrix sQt:\n"
425  "%3\n").arg(probePos.first).arg(qstring_cast(sMt)).arg(qstring_cast(sQt)));
426  }
427 
428 }
429 
430 
431 //------------------------------------------------------------------------------
432 }//namespace cx
cx::SmoothingImageFilter
Definition: cxSmoothingImageFilter.h:26
cx::XmlOptionFile
Helper class for xml files used to store ssc/cx data.
Definition: cxXmlOptionItem.h:78
cx::BinaryThinningImageFilter3DFilter
Definition: cxBinaryThinningImageFilter3DFilter.h:27
cxPipelineWidget.h
qstring_cast
QString qstring_cast(const T &val)
Definition: cxTypeConversions.h:46
cx::BinaryThresholdImageFilter
Definition: cxBinaryThresholdImageFilter.h:31
cx
Namespace for all CustusX production code.
Definition: cx_dev_group_definitions.h:13
cx::report
void report(QString msg)
Definition: cxLogger.cpp:69
cxBinaryThresholdImageFilter.h
cxViewGroupData.h
cxUsReconstructionServiceProxy.h
cx::DataLocations::getRootConfigPath
static QString getRootConfigPath()
return path to root config folder. May be replaced with getExistingConfigPath()
Definition: cxDataLocations.cpp:169
cx::RegistrationBaseWidget
Definition: cxRegistrationBaseWidget.h:28
cxReporter.h
cx::FilterGroup
Definition: cxFilterGroup.h:30
cxUSReconstructInputData.h
cx::PipelineWidget
Definition: cxPipelineWidget.h:59
SeansVesselReg.hxx
cxCompositeTimedAlgorithm.h
vtkPolyDataPtr
vtkSmartPointer< class vtkPolyData > vtkPolyDataPtr
Definition: cxAirwaysFilterService.h:35
cxWirePhantomWidget.h
cx::RegistrationBaseWidget::mServices
RegServicesPtr mServices
Definition: cxRegistrationBaseWidget.h:35
cxDistanceMetric.h
cx::MeshPtr
boost::shared_ptr< class Mesh > MeshPtr
Definition: cxForwardDeclarations.h:48
cxBinaryThinningImageFilter3DFilter.h
cxProfile.h
cx::Transform3D
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
Definition: cxLandmarkPatientRegistrationWidget.h:33
cxViewService.h
cx::WirePhantomWidget::~WirePhantomWidget
virtual ~WirePhantomWidget()
Definition: cxWirePhantomWidget.cpp:97
cxTool.h
cxTypeConversions.h
cx::DataPtr
boost::shared_ptr< class Data > DataPtr
Definition: cxRegistrationApplicator.h:22
cx::Tool::TOOL_US_PROBE
@ TOOL_US_PROBE
Ultrasond probe. The tool has a Probe subinterface with a sector and a video stream.
Definition: cxTool.h:87
cxPatientModelService.h
cxSmoothingImageFilter.h
cx::FilterPtr
boost::shared_ptr< class Filter > FilterPtr
Definition: cxForwardDeclarations.h:155
cx::WirePhantomWidget::WirePhantomWidget
WirePhantomWidget(ctkPluginContext *pluginContext, QWidget *parent=0)
Definition: cxWirePhantomWidget.cpp:46
vtkPointsPtr
vtkSmartPointer< class vtkPoints > vtkPointsPtr
Definition: vtkForwardDeclarations.h:110
cxRegistrationTransform.h
cx::profile
cxResource_EXPORT ProfilePtr profile()
Definition: cxProfile.cpp:160
cxPointMetric.h
cx::FilterGroupPtr
boost::shared_ptr< FilterGroup > FilterGroupPtr
Definition: cxFilterGroup.h:65
M_PI
#define M_PI
Definition: cxReconstructParams.cpp:25
cxDataLocations.h
cx::ToolPtr
boost::shared_ptr< class Tool > ToolPtr
Definition: cxVideoConnectionWidget.h:43
cx::Mesh
A mesh data set.
Definition: cxMesh.h:45
cx::UsReconstructionServiceProxy
Definition: cxUsReconstructionServiceProxy.h:27
cxRegistrationService.h
cx::createTransformTranslate
Transform3D createTransformTranslate(const Vector3D &translation)
Definition: cxTransform3D.cpp:164
cx::DistanceMetricPtr
boost::shared_ptr< class DistanceMetric > DistanceMetricPtr
Definition: cxForwardDeclarations.h:67
cx::Pipeline
Definition: cxPipeline.h:78
cx::WirePhantomWidget::mLayout
QVBoxLayout * mLayout
Definition: cxWirePhantomWidget.h:49
cx::PointMetricPtr
boost::shared_ptr< class PointMetric > PointMetricPtr
Definition: cxForwardDeclarations.h:84
cx::reportError
void reportError(QString msg)
Definition: cxLogger.cpp:71
cxMesh.h
cx::Vector3D
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
cxTrackingService.h
cx::WirePhantomWidget::defaultWhatsThis
virtual QString defaultWhatsThis() const
Definition: cxWirePhantomWidget.cpp:101
cx::RegServices
Definition: cxRegServices.h:33
cx::reportWarning
void reportWarning(QString msg)
Definition: cxLogger.cpp:70
cxSpaceProvider.h