NorMIT-nav  16.5
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cxElastixExecuter.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 "cxElastixExecuter.h"
35 
36 #include <QProcess>
37 #include <QFile>
38 #include <QDir>
39 #include <QtGlobal>
40 
41 #include "cxTypeConversions.h"
42 #include "cxTime.h"
43 #include "cxData.h"
44 #include "cxBoundingBox3D.h"
45 #include "cxTransformFile.h"
46 #include "cxCustomMetaImage.h"
47 
48 #include "cxPatientModelService.h"
49 
50 namespace cx
51 {
52 
54  TimedBaseAlgorithm("ElastiX", 1000),
55  mServices(services)
56 {
57  mProcess = new QProcess(this);
58  connect(mProcess, SIGNAL(stateChanged(QProcess::ProcessState)), this, SLOT(processStateChanged(QProcess::ProcessState)));
59  connect(mProcess, SIGNAL(error(QProcess::ProcessError)), this, SLOT(processError(QProcess::ProcessError)));
60  connect(mProcess, SIGNAL(finished(int, QProcess::ExitStatus)), this, SLOT(processFinished(int, QProcess::ExitStatus)));
61 }
62 
64 {
65  disconnect(mProcess, SIGNAL(error(QProcess::ProcessError)), this, SLOT(processError(QProcess::ProcessError)));
66  mProcess->close();
67 }
68 
70 {
71  disconnect(mProcess, SIGNAL(readyRead()), this, SLOT(processReadyRead()));
72 
73  if (on)
74  {
75  connect(mProcess, SIGNAL(readyRead()), this, SLOT(processReadyRead()));
76  // mProcess->setProcessChannelMode(QProcess::ForwardedChannels);
77  mProcess->setProcessChannelMode(QProcess::MergedChannels);
78  mProcess->setReadChannel(QProcess::StandardOutput);
79  }
80 }
81 
82 
83 bool ElastixExecuter::setInput(QString application,
84  DataPtr fixed,
85  DataPtr moving,
86  QString outdir,
87  QStringList parameterfiles)
88 {
89  mFixed = fixed;
90  mMoving = moving;
91 
92  if (!fixed || !moving)
93  {
94  reportWarning("Failed to start elastiX registration, fixed or missing image missing.");
95  return false;
96  }
97 
98  if (mProcess->state() != QProcess::NotRunning)
99  {
100  reportWarning("Failed to start elastiX registration, process already running");
101  return false;
102  }
103 
104  // create output dir (required by ElastiX)
105  QDir().mkpath(outdir);
106 
107  QString initFilename = this->writeInitTransformToElastixfile(fixed, moving, outdir);
108  this->writeInitTransformToCalfile(fixed, moving, outdir);
109 
110  mLastOutdir = outdir;
111 
112  QStringList cmd;
113  cmd << "\"" + application + "\"";
114  cmd << "-f" << mServices->patient()->getActivePatientFolder()+"/"+fixed->getFilename();
115  cmd << "-m" << mServices->patient()->getActivePatientFolder()+"/"+moving->getFilename();
116  cmd << "-out" << outdir;
117  cmd << "-t0" << initFilename;
118  for (int i=0; i<parameterfiles.size(); ++i)
119  cmd << "-p" << parameterfiles[i];
120 
121  QString commandLine = cmd.join(" ");
122  report(QString("Executing registration with command line: [%1]").arg(commandLine));
123 
124 #ifdef Q_OS_LINUX
125  // hack that inserts . into library path for linux. Solveds issue with elastix lib not being fixed up on linux.
126  QString path = QFileInfo(application).absolutePath();
127  QProcessEnvironment env = QProcessEnvironment::systemEnvironment();
128  env.insert("LD_LIBRARY_PATH", path);
129  mProcess->setProcessEnvironment(env);
130 #endif
131  mProcess->start(commandLine);
132  return true;
133 }
134 
136 {
137  emit aboutToStart(); // if properly set: triggers the setInput()
138 }
139 
141 {
142  std::cout << "TODO: check ElastixExecuter::isFinished()" << std::endl;
143  return mProcess->atEnd();
144 }
145 
147 {
148  return mProcess->state()!=QProcess::NotRunning;
149 }
150 
151 QString ElastixExecuter::writeInitTransformToElastixfile(
152  DataPtr fixed,
153  DataPtr moving,
154  QString outdir)
155 {
156  // elastiX uses the transforms present in the mhd files. If the mhd info is up to date
157  // with the dataManager info, the T0=I, i.e we dont need to tell elastiX anything.
158  // If NOT up to date, then compare file and dataManager, then insert the difference as T0:
159  //
160  // Let f be fixed, m be moving, mm and ff be the intermediate spaces between the file and r:
161  // rMd is the transform stored in ssc.
162  // ddMd is the transform stored in the mhd file.
163  // T0 is the remainder to be sent to elastiX
164  //
165  // mMf = mMr*rMf = mMmm*mmMr*rMff*ffMf = mMmm*T0*ffMf
166  // -->
167  // T0 = mmMm*mMr*rMf*fMff
168  //
169  Transform3D rMf = fixed->get_rMd();
170  Transform3D rMm = moving->get_rMd();
171  Transform3D ffMf = this->getFileTransform_ddMd(mFixed);
172  Transform3D mmMm = this->getFileTransform_ddMd(mMoving);
173 // Transform3D mMf = rMm.inv() * rMf;
174  // -->
175  // The remainder transform, not stored in mhd files, must be sent to elastiX:
176  Transform3D T0 = mmMm*rMm.inv()*rMf*ffMf.inv();
177 
178 // Transform3D mMf = moving->get_rMd().inv() * fixed->get_rMd();
179  ElastixEulerTransform E = ElastixEulerTransform::create(T0, fixed->boundingBox().center());
180 
181  QString elastiXText = QString(""
182  "// Input transform file\n"
183  "// Auto-generated by CustusX on %1\n"
184  "\n"
185  "(Transform \"EulerTransform\")\n"
186  "(NumberOfParameters 6)\n"
187  "(TransformParameters %2 %3)\n"
188  "(InitialTransformParametersFileName \"NoInitialTransform\")\n"
189  "(HowToCombineTransforms \"Compose\")\n"
190  "\n"
191  "// EulerTransform specific\n"
192  "(CenterOfRotationPoint %4)\n"
193  "(ComputeZYX \"false\")\n"
194  "").arg(QDateTime::currentDateTime().toString(timestampSecondsFormatNice()),
198 
199  QFile initTransformFile(outdir+"/t0.txt");
200  if (!initTransformFile.open(QIODevice::WriteOnly | QIODevice::Text))
201  {
202  reportWarning(QString("Failed to open file %1 for writing.").arg(initTransformFile.fileName()));
203  }
204  initTransformFile.write(elastiXText.toLatin1());
205  return initTransformFile.fileName();
206 }
207 
208 QString ElastixExecuter::writeInitTransformToCalfile(
209  DataPtr fixed,
210  DataPtr moving,
211  QString outdir)
212 {
213  Transform3D mMf = moving->get_rMd().inv() * fixed->get_rMd();
214 
215  TransformFile file(outdir+"/moving_M_fixed_initial.cal");
216  file.write(mMf);
217 
218  return file.fileName();
219 }
220 
221 void ElastixExecuter::processReadyRead()
222 {
223  report(QString(mProcess->readAllStandardOutput()));
224 }
225 
226 void ElastixExecuter::processError(QProcess::ProcessError error)
227 {
228  QString msg;
229  msg += "Registration process reported an error: ";
230 
231  switch (error)
232  {
233  case QProcess::FailedToStart:
234  msg += "Failed to start";
235  break;
236  case QProcess::Crashed:
237  msg += "Crashed";
238  break;
239  case QProcess::Timedout:
240  msg += "Timed out";
241  break;
242  case QProcess::WriteError:
243  msg += "Write Error";
244  break;
245  case QProcess::ReadError:
246  msg += "Read Error";
247  break;
248  case QProcess::UnknownError:
249  msg += "Unknown Error";
250  break;
251  default:
252  msg += "Invalid error";
253  }
254 
255  reportError(msg);
256 }
257 
258 void ElastixExecuter::processFinished(int code, QProcess::ExitStatus status)
259 {
260  if (status == QProcess::CrashExit)
261  reportError("Registration process crashed");
262 }
263 
264 
265 void ElastixExecuter::processStateChanged(QProcess::ProcessState newState)
266 {
267 // removed text: to much information
268 
269  QString msg;
270  msg += "Registration process";
271 
272  if (newState == QProcess::Running)
273  {
274 // report(msg + " running.");
275  emit started(0);
276  }
277  if (newState == QProcess::NotRunning)
278  {
279  emit finished();
280 // report(msg + " not running.");
281  }
282  if (newState == QProcess::Starting)
283  {
284 // report(msg + " starting.");
285  }
286 }
287 
288 QString ElastixExecuter::findMostRecentTransformOutputFile() const
289 {
290  QString retval;
291  for (int i=0; ; ++i)
292  {
293  QString filename = QString(mLastOutdir + "/TransformParameters.%1.txt").arg(i);
294  if (!QFileInfo(filename).exists())
295  break;
296  retval = filename;
297  }
298  return retval;
299 }
300 
301 Transform3D ElastixExecuter::getFileTransform_ddMd(DataPtr volume)
302 {
303  QString patFolder = mServices->patient()->getActivePatientFolder();
304  CustomMetaImagePtr reader = CustomMetaImage::create(patFolder+"/"+volume->getFilename());
305  Transform3D ddMd = reader->readTransform();
306  return ddMd;
307 }
308 
310 {
311  Transform3D mmMff = this->getAffineResult_mmMff(ok);
312  Transform3D ffMf = this->getFileTransform_ddMd(mFixed);
313  Transform3D mmMm = this->getFileTransform_ddMd(mMoving);
314 
315  return mmMm.inv() * mmMff * ffMf;
316 }
317 
318 Transform3D ElastixExecuter::getAffineResult_mmMff(bool* ok)
319 {
320  QString filename = this->findMostRecentTransformOutputFile();
321  Transform3D mMf = Transform3D::Identity();
322 
323  if (filename.isEmpty())
324  {
325  if (ok)
326  *ok = false;
327 
328  TransformFile file(mLastOutdir+"/moving_M_fixed_registered.cal");
329  mMf = file.read(ok);
330 
331  if (ok && !*ok)
332  {
333  reportWarning("Failed to read registration results.");
334  }
335 
336  return mMf;
337  }
338 
339  while (true)
340  {
341  ElastixParameterFile file(filename);
342 
343  bool useDirectionCosines = file.readParameterBool("UseDirectionCosines");
344  if (useDirectionCosines)
345  {
346  reportWarning("Elastix UseDirectionCosines is not supported. Result is probably wrong.");
347  }
348 
349  QString transformType = file.readParameterString("Transform");
350  if (transformType=="EulerTransform")
351  {
352  if (ok)
353  *ok = true;
354  Transform3D mQf = file.readEulerTransform();
355  // concatenate transforms:
356  mMf = mQf * mMf;
357  }
358  else if (transformType=="AffineTransform")
359  {
360  if (ok)
361  *ok = true;
362  Transform3D mQf = file.readAffineTransform();
363  // concatenate transforms:
364  mMf = mQf * mMf;
365  }
366  else
367  {
368  // accept invalid transforms, but emit warning.
369 // if (ok)
370 // *ok = false;
371  reportWarning(QString("TransformType [%1] is not supported by CustusX. Registration result from %2 ignored.").arg(transformType).arg(filename));
372  }
373 
374  filename = file.readParameterString("InitialTransformParametersFileName");
375  if (filename.isEmpty() || filename=="NoInitialTransform")
376  break;
377  }
378 
379  return mMf;
380 }
381 
382 
384 {
385  if (ok)
386  *ok = false;
387 
388  QString retval;
389  int i=0;
390  for (; ; ++i)
391  {
392  //TODO only mhd supported
393  QString filename = QString(mLastOutdir + "/result.%1.mhd").arg(i);
394  if (!QFileInfo(filename).exists())
395  break;
396  retval = filename;
397  }
398 
399  if (retval.isEmpty())
400  return retval;
401 
402  QString paramFilename = QString(mLastOutdir + "/TransformParameters.%1.txt").arg(i-1);
403  ElastixParameterFile file(paramFilename);
404  QString transform = file.readParameterString("Transform");
405 
406  if ((transform=="BSplineTransform") || (transform=="SplineKernelTransform"))
407  {
408  report(QString("Reading result file %1 created with transform %2").arg(retval).arg(transform));
409  if (ok)
410  *ok = true;
411  return retval;
412  }
413  else
414  {
415  return "";
416  }
417 
418 }
419 
420 
421 // --------------------------------------------------------
422 // --------------------------------------------------------
423 // --------------------------------------------------------
424 
425 
426 
427 
429 {
430  QString transformType = this->readParameterString("Transform");
431  if (transformType!="EulerTransform")
432  reportError("Assert failure: attempting to read EulerTransform");
433 
434  int numberOfParameters = this->readParameterInt("NumberOfParameters");
435  if (numberOfParameters!=6)
436  {
437  reportWarning(QString("Expected 6 Euler parameters, got %1").arg(numberOfParameters));
438  return Transform3D::Identity();
439  }
440  std::vector<double> tp = this->readParameterDoubleVector("TransformParameters");
441  std::vector<double> cor = this->readParameterDoubleVector("CenterOfRotationPoint");
442 
444  Vector3D(tp[0], tp[1], tp[2]),
445  Vector3D(tp[3], tp[4], tp[5]),
446  Vector3D(cor[0], cor[1], cor[2]));
447 
448  return E.toMatrix();
449 }
450 
452 {
453  QString transformType = this->readParameterString("Transform");
454  if (transformType!="AffineTransform")
455  reportError("Assert failure: attempting to read AffineTransform");
456 
457  int numberOfParameters = this->readParameterInt("NumberOfParameters");
458  if (numberOfParameters!=12)
459  {
460  reportWarning(QString("Expected 12 Euler parameters, got %1").arg(numberOfParameters));
461  return Transform3D::Identity();
462  }
463  std::vector<double> tp = this->readParameterDoubleVector("TransformParameters");
464 // std::vector<double> cor = this->readParameterDoubleVector("CenterOfRotationPoint");
465 
466  Transform3D M = Transform3D::Identity();
467 
468  for (int r=0; r<3; ++r)
469  for (int c=0; c<3; ++c)
470  M(r,c) = tp[3*r+c];
471  for (int r=0; r<3; ++r)
472  M(r,3) = tp[9+r];
473 
474  return M;
475 }
476 
477 ElastixParameterFile::ElastixParameterFile(QString filename) : mFile(filename)
478 {
479  if (!mFile.open(QIODevice::ReadOnly | QIODevice::Text))
480  {
481  reportWarning(QString("Can't open ElastiX result file %1").arg(filename));
482  }
483  mText = QString(mFile.readAll());
484 // std::cout << "Loaded text from " << filename << ":\n" << mText << std::endl;
485 }
486 
487 QString ElastixParameterFile::readParameterRawValue(QString key)
488 {
489  QStringList lines = mText.split('\n');
490 
491  QString regexpStr = QString(""
492  "\\s*" // match zero or more whitespace
493  "\\(" // match beginning (
494  "\\s*" // match zero or more whitespace
495  "%1" // key
496  "\\s+" // match one or more whitespace
497  "(" // start grab value
498  "[^\\)]*" // value - anything up to the closing )
499  ")" // end grab value
500  "\\)" // match ending )
501  "[^\\n]*" // rest of line - ignore
502  ).arg(key);
503  QRegExp rx(regexpStr);
504 // std::cout << regexpStr << std::endl;
505 
506  lines.indexOf(rx);
507 
508 // if (lines.indexOf(rx)>=0)
509 // {
510 // std::cout << "FOUND0 " << rx.cap(0) << std::endl;
511 // std::cout << "FOUND1 " << rx.cap(1) << std::endl;
512 // }
513 
514  return rx.cap(1);
515 }
516 
518 {
519  QString retval = this->readParameterRawValue(key);
520  if (retval.startsWith("\"") && retval.endsWith("\""))
521  {
522  retval = retval.replace(0, 1, "");
523  retval = retval.replace(retval.size()-1, 1, "");
524  }
525 
526  return retval;
527 }
528 
530 {
531  QString text = this->readParameterString(key);
532  return (text=="true") ? true : false;
533 }
534 
536 {
537  QString retval = this->readParameterRawValue(key);
538  return retval.toInt();
539 }
540 
541 std::vector<double> ElastixParameterFile::readParameterDoubleVector(QString key)
542 {
543  QString retval = this->readParameterRawValue(key);
544  return convertQString2DoubleVector(retval);
545 }
546 
547 } /* namespace cx */
QString qstring_cast(const T &val)
ElastixParameterFile(QString filename)
bool setInput(QString application, DataPtr fixed, DataPtr moving, QString outdir, QStringList parameterfiles)
DoubleBoundingBox3D transform(const Transform3D &m, const DoubleBoundingBox3D &bb)
void finished()
should be emitted when at the end of postProcessingSlot
QString readParameterString(QString key)
void reportError(QString msg)
Definition: cxLogger.cpp:92
File format for storing a 4x4 matrix.The read/write methods emit error messages if you dont use the o...
Base class for algorithms that wants to time their execution.
boost::shared_ptr< class RegServices > RegServicesPtr
Definition: cxRegServices.h:41
virtual bool isFinished() const
static ElastixEulerTransform create(Vector3D angles_xyz, Vector3D translation, Vector3D centerOfRotation)
int readParameterInt(QString key)
ElastixExecuter(RegServicesPtr services, QObject *parent=NULL)
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
QString timestampSecondsFormatNice()
Definition: cxTime.cpp:47
Transform3D toMatrix() const
QString getNonlinearResultVolume(bool *ok=0)
boost::shared_ptr< class Data > DataPtr
boost::shared_ptr< class CustomMetaImage > CustomMetaImagePtr
void reportWarning(QString msg)
Definition: cxLogger.cpp:91
bool readParameterBool(QString key)
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
void report(QString msg)
Definition: cxLogger.cpp:90
virtual bool isRunning() const
std::vector< double > readParameterDoubleVector(QString key)
void aboutToStart()
emitted at start of execute. Use to perform preprocessing
void setDisplayProcessMessages(bool on)
std::vector< double > convertQString2DoubleVector(const QString &input, bool *ok)
Transform3D getAffineResult_mMf(bool *ok=0)
static CustomMetaImagePtr create(QString filename)
void started(int maxSteps)
emitted at start of run.