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