registration.h

Go to the documentation of this file.

/******************************************************************************
 * Copyright 2014 Micah C Chambers (micahc.vt@gmail.com)
 *
 * NPL is free software: you can redistribute it and/or modify it under the
 * terms of the BSD 2-Clause License available in LICENSE or at
 * http://opensource.org/licenses/BSD-2-Clause
 *
 * @file registration.h Tools for registering 2D or 3D images.
 *
 *****************************************************************************/

#ifndef REGISTRATION_H
#define REGISTRATION_H

#include "mrimage.h"
#include "accessors.h"
#include "iterators.h"
#include <Eigen/Dense>

#include <memory>

namespace npl {

class RigidInfoComp
{
public:

    RigidInfoComp(bool mindiff);

    int valueGrad(const Eigen::VectorXd& params, double& val,
            Eigen::VectorXd& grad);

    int grad(const Eigen::VectorXd& params, Eigen::VectorXd& grad);

    int value(const Eigen::VectorXd& params, double& val);

    bool m_compdiff;

    Metric m_metric;

    void setBins(size_t nbins, size_t krad);

    void setFixed(ptr<const MRImage> fixed);

    ptr<const MRImage> getFixed() { return m_fixed; };

    void setMoving(ptr<const MRImage> moving);

    ptr<const MRImage> getMoving() { return m_moving; };

    size_t nparam() { return 6; };

private:

    ptr<const MRImage> m_fixed;
    ptr<const MRImage> m_moving;
    ptr<MRImage> m_dmoving;

    int m_bins;

    int m_krad;

    // for interpolating moving image, and iterating fixed
    LinInterp3DView<double> m_move_get;
    LinInterp3DView<double> m_dmove_get;
    NDConstIter<double> m_fit;

    double m_center[3];

    NDArrayStore<1, double> m_pdfmove;
    NDArrayStore<1, double> m_pdffix;

    NDArrayStore<2, double> m_pdfjoint;
    NDArrayStore<3, double> m_dpdfjoint;
    NDArrayStore<2, double> m_dpdfmove;

    vector<double> m_gradHjoint;
    vector<double> m_gradHmove;

    double m_Hfix;
    double m_Hmove;
    double m_Hjoint;

    double m_rangemove[2];
    double m_rangefix[2];
    double m_wmove;
    double m_wfix;
};

class RigidCorrComp
{
public:

    RigidCorrComp(bool mindiff);

    int valueGrad(const Eigen::VectorXd& params, double& val,
            Eigen::VectorXd& grad);

    int grad(const Eigen::VectorXd& params, Eigen::VectorXd& grad);

    int value(const Eigen::VectorXd& params, double& val);

    void setFixed(ptr<const MRImage> fixed);

    ptr<const MRImage> getFixed() { return m_fixed; };

    void setMoving(ptr<const MRImage> moving);

    ptr<const MRImage> getMoving() { return m_moving; };

    size_t nparam() { return 6; };

    bool m_compdiff;
private:

    ptr<const MRImage> m_fixed;
    ptr<const MRImage> m_moving;
    ptr<MRImage> m_dmoving;

    double m_center[3];

};

class ProbDistCorrInfoComp
{
public:

    ProbDistCorrInfoComp(bool mindiff);

    int valueGrad(const Eigen::VectorXd& params, double& val,
            Eigen::VectorXd& grad);

    int grad(const Eigen::VectorXd& params, Eigen::VectorXd& grad);

    int value(const Eigen::VectorXd& params, double& val);

    bool m_compdiff;

    double m_tps_reg;

    double m_jac_reg;

    Metric m_metric;

    size_t nparam() { return m_deform ? m_deform->elements() : 0; };

    void initialize(ptr<const MRImage> fixed, ptr<const MRImage> moving,
            size_t m_krad, size_t nbins, double space, int dir);

    ptr<const MRImage> getFixed() { return m_fixed; };

    ptr<const MRImage> getMoving() { return m_moving; };

    ptr<MRImage> getDeform() { return m_deform; };

private:
    int metric(double& val, VectorXd& grad);

    int metric(double& val);

    /* Variables:
     *
     * m_bins, m_krad
     *
     */

    int m_dir;

    ptr<const MRImage> m_fixed;
    ptr<const MRImage> m_moving;
    ptr<MRImage> m_dmoving;
    ptr<MRImage> m_deform;

    ptr<MRImage> m_fixed_cache;

    void updateCaches();

    int m_bins;

    int m_krad;

    double m_knotspace;

    // for interpolating moving image, and iterating fixed
    VectorXd gradbuff;

    NDArrayStore<1, double> m_pdffix;

    NDArrayStore<1, double> m_pdfmove;

    NDArrayStore<2, double> m_pdfjoint;

    MRImageStore<5, double> m_dpdfjoint;

    MRImageStore<4, double> m_dpdfmove;

    MRImageStore<3, double> m_gradHjoint;

    MRImageStore<3, double> m_gradHmove;

    double m_Hfix;

    double m_Hmove;

    double m_Hjoint;

    double m_rangefix[2];

    double m_wfix;
};


class DistCorrInfoComp
{
public:

    DistCorrInfoComp(bool mindiff);

    int valueGrad(const Eigen::VectorXd& params, double& val,
            Eigen::VectorXd& grad);

    int grad(const Eigen::VectorXd& params, Eigen::VectorXd& grad);

    int value(const Eigen::VectorXd& params, double& val);

    int m_dir;

    bool m_compdiff;

    double m_tps_reg;

    double m_jac_reg;

    Metric m_metric;

    size_t nparam() { return m_deform ? m_deform->elements() : 0; };

    void setBins(size_t nbins, size_t krad);

    void initializeKnots(double space);

    void setFixed(ptr<const MRImage> fixed);

    ptr<const MRImage> getFixed() { return m_fixed; };

    void setMoving(ptr<const MRImage> moving, int dir = -1);

    ptr<const MRImage> getMoving() { return m_moving; };

    ptr<MRImage> getDeform() { return m_deform; };

private:
    int metric(double& val, VectorXd& grad);

    int metric(double& val);

    /* Variables:
     *
     * m_bins, m_krad
     *
     */

    ptr<const MRImage> m_fixed;
    ptr<const MRImage> m_moving;
    ptr<MRImage> m_dmoving;
    ptr<MRImage> m_deform;

    ptr<MRImage> m_move_cache;
    ptr<MRImage> m_dmove_cache;
    ptr<MRImage> m_corr_cache;

    void updateCaches();

    int m_bins;

    int m_krad;

    double m_knotspace;

    // for interpolating moving image, and iterating fixed
    VectorXd gradbuff;

    NDArrayStore<1, double> m_pdfmove;

    NDArrayStore<1, double> m_pdffix;

    NDArrayStore<2, double> m_pdfjoint;

    MRImageStore<5, double> m_dpdfjoint;

    MRImageStore<4, double> m_dpdfmove;

    MRImageStore<3, double> m_gradHjoint;

    MRImageStore<3, double> m_gradHmove;

    double m_Hfix;

    double m_Hmove;

    double m_Hjoint;

    double m_rangemove[2];

    double m_rangefix[2];

    double m_wmove;

    double m_wfix;
};


struct Rigid3DTrans
{
    Vector3d rotation; //Rx, Ry, Rz
    Vector3d shift;
    Vector3d center;

    bool ras_coord;

    Rigid3DTrans() {
        ras_coord = true;
        rotation.setZero();
        shift.setZero();
        center.setZero();
    };

    void invert();

    Matrix3d rotMatrix();

    void setRotation(const Matrix3d& rot);

    void toRASCoords(ptr<const MRImage> in);

    void toIndexCoords(ptr<const MRImage> in, bool forcegridcenter);

};

ptr<MRImage> motionCorrect(ptr<const MRImage> input, size_t ref);

Rigid3DTrans corReg3D(ptr<const MRImage> fixed, ptr<const MRImage> moving,
        const std::vector<double>& sigmas);

Rigid3DTrans informationReg3D(ptr<const MRImage> fixed,
        ptr<const MRImage> moving, const std::vector<double>& sigmas,
        size_t nbins = 128, size_t binradius = 4, std::string metric = "MI",
        double stopx = 0.001);

ptr<MRImage> infoDistCor(ptr<const MRImage> fixed, ptr<const MRImage> moving,
        bool otsu, int dir, double bspace, double jac, double tps,
        const std::vector<double>& sigmas,
        size_t nbins = 128, size_t binradius = 4, string metric = "MI",
        size_t hist = 8, double stopx = 1e-5, double beta = 0.5);

int cor3DDerivTest(double step, double tol, ptr<const MRImage> in1,
        ptr<const MRImage> in2);

int information3DDerivTest(double step, double tol,
        ptr<const MRImage> in1, ptr<const MRImage> in2);

int distcorDerivTest(double step, double tol,
        shared_ptr<const MRImage> in1, shared_ptr<const MRImage> in2,
        double regj = 0, double regt = 0);

inline
std::ostream& operator<< (std::ostream& stream, const Rigid3DTrans& rigid)
{
    stream << "Rigid3DTrans ";
    if(rigid.ras_coord)
        stream << "(In RAS)\n";
    else
        stream << "(In Index)\n";

    stream << "Rotation: ";
    for(size_t ii=0; ii<2; ii++)
        stream << rigid.rotation[ii] << ", ";
    stream << rigid.rotation[2] << "\n";

    stream << "Center: ";
    for(size_t ii=0; ii<2; ii++)
        stream << rigid.center[ii] << ", ";
    stream << rigid.center[2] << "\n";

    stream << "Shift : ";
    for(size_t ii=0; ii<2; ii++)
        stream << rigid.shift[ii] << ", ";
    stream << rigid.shift[2] << "\n";

    return stream;
}

} // npl
#endif  //REGISTRATION_H