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