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