fmri_inference.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 fmri_inference.h Tools for performing ICA, including rewriting images as
 * matrices. This file should really be fmri inference helpers or something
 * along that line, because it also include GLM helpers.
 *
 *****************************************************************************/
#ifndef ICA_HELPERS_H
#define ICA_HELPERS_H

#include <Eigen/Dense>
#include <Eigen/IterativeSolvers>
#include <string>
#include <vector>

#include "npltypes.h"
#include "mrimage.h"
#include "utility.h"
#include "macros.h"

namespace npl {

class MatrixReorg;

VectorXd onDiskSVD(const MatrixReorg& A,
        int rank, size_t poweriters, double varthresh, double cvarthresh,
        MatrixXd* U=NULL, MatrixXd* V=NULL);

void gicaCreateMatrices(size_t tcat, size_t scat, vector<std::string> masks,
        vector<std::string> inputs, std::string prefix, double maxmem, bool normts,
        bool verbose);

void gicaReduceFull(std::string inpref, std::string outpref, double varthresh,
        double cvarthresh, size_t maxrank, bool verbose);

void gicaReduceProb(std::string inpref, std::string outpref, double varthresh,
        double cvarthresh, size_t rank, size_t poweriters, bool verbose);

void gicaTemporalICA(std::string reorgpref, std::string reducepref,
        std::string outpref, bool verbose);

void gicaSpatialICA(std::string reorgpref, std::string reducepref,
        std::string outpref, bool verbose);

VectorXd covSVD(const MatrixReorg& A, double varthresh, double cvarthresh,
        size_t maxrank, MatrixXd* U, MatrixXd* V);

class MatMap
{
public:
    MatMap() : mat(NULL, 0, 0)
    {
    };

    MatMap(std::string filename, size_t rows, size_t cols) : mat(NULL, 0, 0)
    {
        create(filename, rows, cols);
    };

    MatMap(std::string filename, bool writeable=false) : mat(NULL, 0, 0)
    {
        open(filename, writeable);
    };

    ~MatMap() { close(); };

    void open(std::string filename, bool writeable = false);

    void create(std::string filename, size_t newrows, size_t newcols);

    void close()
    {
        datamap.close();
    };

    bool isopen()
    {
        return datamap.isopen();
    };

    const size_t& rows() const { return m_rows; };
    const size_t& cols() const { return m_cols; };

    Eigen::Map<MatrixXd> mat;
private:
    MemMap datamap;
    size_t m_rows;
    size_t m_cols;
};

class MatrixReorg
{
public:

    MatrixReorg(std::string prefix="", size_t maxd=(1<<30), bool verbose=true);

    int createMats(size_t timeblocks, size_t spaceblocks,
            const std::vector<std::string>& masknames,
            const std::vector<std::string>& filenames, bool normts = true);

    int checkMats();

//  inline int nwide() const { return m_outrows.size(); };
    inline int ntall() const { return m_outcols.size(); };

    inline const vector<int>& tallMatCols() const
    {
        return m_outcols;
    };

    inline int tallMatRows() const
    {
        return m_totalrows;
    };

    inline std::string tallMatName(size_t ii) const
    {
        return m_prefix+"_tall_"+std::to_string(ii);
    };

    inline std::string inColMaskName(size_t ii) const
    {
        return m_prefix+"_mask_"+std::to_string(ii)+".nii.gz";
    };

    void preMult(Eigen::Ref<MatrixXd> out, const Eigen::Ref<const MatrixXd> in,
            bool transpose = false) const;

    void postMult(Eigen::Ref<MatrixXd> out, const Eigen::Ref<const MatrixXd> in,
            bool transpose = false) const;

    inline int rows() const { return m_totalrows; };
    inline int cols() const { return m_totalcols; };

    std::string mask_name(size_t ii) const
    {
        return m_prefix+"_mask_"+std::to_string(ii)+".nii.gz";
    };

    std::string info_name() const
    {
        return m_prefix+".info";
    };
    std::string tall_name(size_t ii) const
    {
        return m_prefix+"_tall_"+std::to_string(ii);
    };

    private:
    int m_totalrows;
    int m_totalcols;
    std::string m_prefix;
    bool m_verbose;
//  vector<int> m_outrows;
    vector<int> m_outcols;
    size_t m_maxdoubles;
};

void fmriGLM(ptr<const MRImage> fmri, const MatrixXd& X,
        ptr<MRImage> bimg, ptr<MRImage> Timg, ptr<MRImage> pimg);


void fmriBandPass(ptr<MRImage> inimg, double cuton, double cutoff);

ptr<MRImage> regressOut(ptr<const MRImage> inimg, const MatrixXd& X);

MatrixXd extractLabelAVG(ptr<const MRImage> fmri,
        ptr<const MRImage> labelmap);

MatrixXd extractLabelPCA(ptr<const MRImage> fmri,
        ptr<const MRImage> labelmap, size_t outsz);

MatrixXd extractLabelICA(ptr<const MRImage> fmri,
        ptr<const MRImage> labelmap, size_t outsz);

} // NPL

#endif //ICA_HELPERS_H