cxElastixExecuter.h

Go to the documentation of this file.

/*=========================================================================
This file is part of CustusX, an Image Guided Therapy Application.
                 
Copyright (c) SINTEF Department of Medical Technology.
All rights reserved.
                 
CustusX is released under a BSD 3-Clause license.
                 
See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
=========================================================================*/


#ifndef CXELASTIXEXECUTER_H_
#define CXELASTIXEXECUTER_H_

#include <QString>
#include "cxForwardDeclarations.h"
#include "cxTransform3D.h"
#include <QObject>
#include <QProcess>
#include <QFile>
#include "cxTimedAlgorithm.h"
#include "cxLogger.h"
#include "cxRegServices.h"
#include "org_custusx_registration_method_commandline_Export.h"

namespace cx
{
class org_custusx_registration_method_commandline_EXPORT ElastixExecuter : public TimedBaseAlgorithm
{
        Q_OBJECT
public:
        ElastixExecuter(RegServicesPtr services, QObject* parent=NULL);
        virtual ~ElastixExecuter();

        void setDisplayProcessMessages(bool on);

        bool setInput(QString application,
                 DataPtr fixed,
                 DataPtr moving,
                 QString outdir,
                 QStringList parameterfiles);
        virtual void execute();
        virtual bool isFinished() const;
    virtual bool isRunning() const;

        Transform3D getAffineResult_mMf(bool* ok = 0) ;
        QString getNonlinearResultVolume(bool* ok = 0);

private slots:
        void processStateChanged(QProcess::ProcessState newState);
        void processError(QProcess::ProcessError error);
        void processFinished(int, QProcess::ExitStatus);
        void processReadyRead();

private:
        QString writeInitTransformToElastixfile(DataPtr fixed, DataPtr moving, QString outdir);
        QString writeInitTransformToCalfile(DataPtr fixed, DataPtr moving, QString outdir);
        QString findMostRecentTransformOutputFile() const;
        Transform3D getAffineResult_mmMff(bool* ok = 0) ;
        Transform3D getFileTransform_ddMd(DataPtr volume);

        QString mLastOutdir;
        QProcess* mProcess;
        DataPtr mFixed;
        DataPtr mMoving;
        RegServicesPtr mServices;
};

class ElastixParameterFile
{
public:
        ElastixParameterFile(QString filename);

        bool isValid() const;
        QString readParameterString(QString key) ;
        bool readParameterBool(QString key) ;
        int readParameterInt(QString key) ;
        std::vector<double> readParameterDoubleVector(QString key) ;
//      void writeParameter(QString key, QStringList value);
        Transform3D readEulerTransform() ;
        Transform3D readAffineTransform();

private:
        QString readParameterRawValue(QString key);
        QFile mFile;
        QString mText;
};

class ElastixEulerTransform
{
public:
        Vector3D mAngles_xyz;
        Vector3D mTranslation;
        Vector3D mCenterOfRotation;

        static ElastixEulerTransform create(Vector3D angles_xyz, Vector3D translation, Vector3D centerOfRotation)
        {
                ElastixEulerTransform retval;
                retval.mAngles_xyz = angles_xyz;
                retval.mTranslation = translation;
                retval.mCenterOfRotation = centerOfRotation;
                return retval;
        }
        static ElastixEulerTransform create(Transform3D M, Vector3D centerOfRotation)
        {
                ElastixEulerTransform retval;
                retval.mAngles_xyz = M.matrix().block<3, 3> (0, 0).eulerAngles(0, 1, 2);
                retval.mCenterOfRotation = centerOfRotation;

                Transform3D R = retval.getRotationMatrix();
                Transform3D C = createTransformTranslate(retval.mCenterOfRotation);
                // solve M = T*Tc*R*Tc.inv() with respect to T:
                Transform3D T = M*C*R.inv()*C.inv();

                retval.mTranslation = T.matrix().block<3, 1> (0, 3);
//              retval.mCenterOfRotation = Vector3D(0,0,0);
                return retval;
        }
        Transform3D toMatrix() const
        {
                Transform3D T = createTransformTranslate(mTranslation);
                Transform3D C = createTransformTranslate(mCenterOfRotation);
                Transform3D R = this->getRotationMatrix();
                Transform3D Q = T*C*R*C.inv();
                return Q;
        }
        static void test()
        {
                std::cout << "==========TEST==============" << std::endl;
                Transform3D M = createTransformRotateX(M_PI/3)*createTransformRotateY(M_PI/4)*createTransformTranslate(Vector3D(0,2,3));

                ElastixEulerTransform E = create(M, Vector3D(30,40,50));
                Transform3D Q = E.toMatrix();

                std::cout << "M\n" << M << std::endl;
                std::cout << "Q\n" << Q << std::endl;
                Transform3D diff = Q*M.inv();
                std::cout << "Q*M.inv\n" << diff << std::endl;
                if (!similar(Transform3D::Identity(), diff))
                        reportError("assertion failure in ElastixEulerTransform");
        }
private:
        ElastixEulerTransform() {}

        Transform3D getRotationMatrix() const
        {
                Eigen::Matrix3d m;
                m =     Eigen::AngleAxisd(mAngles_xyz[0], Eigen::Vector3d::UnitX())
                  * Eigen::AngleAxisd(mAngles_xyz[1], Eigen::Vector3d::UnitY())
                  * Eigen::AngleAxisd(mAngles_xyz[2], Eigen::Vector3d::UnitZ());

                Transform3D R = Transform3D::Identity();
                R.matrix().block<3, 3> (0, 0) = m;
                return R;
        }
};


} /* namespace cx */
#endif /* CXELASTIXEXECUTER_H_ */