Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /branches/vis15/src/lib/common/image.cxx
ViewVC logotype

View of /branches/vis15/src/lib/common/image.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3884 - (download) (as text) (annotate)
Fri May 20 13:35:25 2016 UTC (3 years, 1 month ago) by jhr
File size: 10216 byte(s)
working on merge
/*! \file image.cxx
 *
 * \author John Reppy
 */

/*
 * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
 *
 * COPYRIGHT (c) 2016 The University of Chicago
 * All rights reserved.
 */

#include <cstdlib>
#include "diderot/image.hxx"

namespace diderot {

/******************** 1D images ********************/

static void load_transform_1d (const Nrrd *nin, double &s, double &t)
{
  // compute the offset to the first space axis
    int base = nin->dim - nin->spaceDim;

  // Image axis Scaling and Rotation
    s = nin->axis[base+0].spaceDirection[0];

  // Image location
    t = nin->spaceOrigin[0];
}

template <typename REAL>
__details::image1d<REAL>::~image1d ()
{
    free (this->_data);
}

template <typename REAL>
bool image1d<REAL>::load (struct WorldBase *wrld, std::string const &name)
{
    Nrrd *nin = load_nrrd_file (wrld, name);
    if (nin == nullptr)
        return true;
    else {
        bool sts = this->load (wrld, nin);
        nrrdNix (nin);
        return sts;
    }
}

template <typename REAL>
bool image1d<REAL>::load (struct WorldBase *wrld, const Nrrd *nin)
{
    this->_img = std::make_shared<__details::image1d<REAL>>();

  // compute the offset to the first space axis
    int base = nin->dim - nin->spaceDim;

// FIXME: we should also be checking the shape of the voxels and the type of the samples.
  // check that the image format is what we expect
    if (nin->spaceDim != 1) {
        error (wrld, "nrrd has unexpected dimension %d, expected 1\n", nin->spaceDim);
        return true;
    }
    else if (nin->axis[base].kind != nrrdKindSpace) {
        error (wrld, "nrrd has unexpected axis structure: %s\n",
            airEnumStr(nrrdKind, nin->axis[base].kind));
        return true;
    }

    this->_img->_dim     = 1;
    this->_img->_size[0] = nin->axis[base+0].size;
    this->_img->_data    = nin->data;
    this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);

//printf("SetImage \"%s\": space dim = 1, axes = <%d>\n",
//name, img->size[base+0]);

  // from the Nrrd file, we load the scaling and translation
    double s, t;

    load_transform_1d (nin, s, t);

    this->_img->_s = static_cast<REAL>(1.0 / s);
    this->_img->_t = static_cast<REAL>(-t / s);

    return false;
}

/******************** 2D images ********************/

static void load_transform_2d (const Nrrd *nin, mat2x2<double> &m, real2<double> &t)
{
#define M(i,j)  m[(i)*2+(j)]

  // compute the offset to the first space axis
    int base = nin->dim - nin->spaceDim;

  // Image axis Scaling and Rotation
    M(0,0) = nin->axis[base+0].spaceDirection[0];
    M(1,0) = nin->axis[base+0].spaceDirection[1];
    M(0,1) = nin->axis[base+1].spaceDirection[0];
    M(1,1) = nin->axis[base+1].spaceDirection[1];

  // Image location
    t[0] = nin->spaceOrigin[0];
    t[1] = nin->spaceOrigin[1];

#undef M
}

inline double det2 (double a, double b, double c, double d)
{
    return (a*d - b*c);
}

/*! \brief compute the inverse of \arg m, storing the result in \arg i.
 *  \param m the matrix to invert
 *  \param i the inverted matrix
 */
inline void invert (mat2x2<double> &i, mat2x2<double> &m)
{
        double scale = 1.0 / det2(m[0], m[1], m[2], m[3]);

        i[0] =  scale * m[3];
        i[1] = -scale * m[1];
        i[2] = -scale * m[2];
        i[3] =  scale * m[0];
}

template <typename REAL>
__details::image2d<REAL>::~image2d ()
{
    free (this->_data);
}

template <typename REAL>
bool image2d<REAL>::load (struct WorldBase *wrld, std::string const &name)
{
    Nrrd *nin = load_nrrd_file (wrld, name);
    if (nin == nullptr)
        return true;
    else {
        bool sts = this->load (wrld, nin);
        nrrdNix (nin);
        return sts;
    }
}

template <typename REAL>
bool image2d<REAL>::load (struct WorldBase *wrld, const Nrrd *nin)
{
    this->_img = std::make_shared<__details::image2d<REAL>>();

  // compute the offset to the first space axis
    int base = nin->dim - nin->spaceDim;

  // check that the image format is what we expect
    if (nin->spaceDim != 2) {
        error (wrld, "nrrd has unexpected dimension %d, expected 2\n", nin->spaceDim);
        return true;
    }
    else if ((nin->axis[base].kind != nrrdKindSpace)
    ||       (nin->axis[base+1].kind != nrrdKindSpace)) {
        error (wrld, "nrrd has unexpected axis structure: %s %s\n",
            airEnumStr(nrrdKind, nin->axis[base].kind),
            airEnumStr(nrrdKind, nin->axis[base+1].kind));
        return true;
    }

    this->_img->_dim	 = 2;
    this->_img->_size[0] = nin->axis[base+0].size;
    this->_img->_size[1] = nin->axis[base+1].size;
    this->_img->_data	 = nin->data;
    this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);

  // from the Nrrd file, we load the affine image-to-world transform matrix
    mat2x2<double> m;              	// rotation and scaling
    real2<double> t;                	// translation part of the transform
    load_transform_2d (nin, m, t);

  // compute inverse of m, which is the transform from the world basis to the image basis
    mat2x2<double> mInv;
    invert (mInv, m);

  // copy inverted matrix into the image data structure
    for (int i = 0;  i < 2;  i++) {
        for (int j = 0;  j < 2;  j++) {
	    REAL r = static_cast<REAL>(mInv[2*i+j]);
            this->_img->_w2i[2*i+j] = r;
	    this->_img->_w2iT[2*j+i] = r;  	// transpose
        }
    }

  // transform the translation vector: inv([M t]) = [inv(M) -inv(M)t]
    for (int i = 0;  i < 2;  i++) {
	double sum = 0;
	for (int j = 0;  j < 2;  j++) {
	    sum += mInv[2*i+j] * t[j];
	}
	this->_img->_tVec[i] = static_cast<REAL>(-sum);
    }

    return false;
}

/******************** 3D images ********************/

static void load_transform_3d (const Nrrd *nin, mat3x3<double> &m, real3<double> &t)
{
#define M(i,j)  m[3*(i)+(j)]

  // compute the offset to the first space axis
    int base = nin->dim - nin->spaceDim;

  // Image axis Scaling and Rotation
    M(0,0) = nin->axis[base+0].spaceDirection[0];
    M(1,0) = nin->axis[base+0].spaceDirection[1];
    M(2,0) = nin->axis[base+0].spaceDirection[2];
    M(0,1) = nin->axis[base+1].spaceDirection[0];
    M(1,1) = nin->axis[base+1].spaceDirection[1];
    M(2,1) = nin->axis[base+1].spaceDirection[2];
    M(0,2) = nin->axis[base+2].spaceDirection[0];
    M(1,2) = nin->axis[base+2].spaceDirection[1];
    M(2,2) = nin->axis[base+2].spaceDirection[2];

  // Image location
    t[0] = nin->spaceOrigin[0];
    t[1] = nin->spaceOrigin[1];
    t[2] = nin->spaceOrigin[2];

#undef M
}

/*! \brief compute the inverse of \arg m, storing the result in \arg i.
 *  \param m the matrix to invert
 *  \param i the inverted matrix
 */
inline void invert (mat3x3<double> &i, mat3x3<double> &m)
{
	double det =
		( m[0]*m[4]*m[8]
		+ m[3]*m[7]*m[2]
		+ m[6]*m[1]*m[5]
		- m[6]*m[4]*m[2]
		- m[3]*m[1]*m[8]
		- m[0]*m[7]*m[5]);
        double scale = 1.0 / det;

        i[0] = scale * det2 (m[4], m[5], m[7], m[8]);
        i[1] = scale * det2 (m[2], m[1], m[8], m[7]);
        i[2] = scale * det2 (m[1], m[2], m[4], m[5]);
        i[3] = scale * det2 (m[5], m[3], m[8], m[6]);
        i[4] = scale * det2 (m[0], m[2], m[6], m[8]);
        i[5] = scale * det2 (m[2], m[0], m[5], m[3]);
        i[6] = scale * det2 (m[3], m[4], m[6], m[7]);
        i[7] = scale * det2 (m[1], m[0], m[7], m[6]);
        i[8] = scale * det2 (m[0], m[1], m[3], m[4]);
}

template <typename REAL>
__details::image3d<REAL>::~image3d ()
{
    free (this->_data);
}

template <typename REAL>
bool image3d<REAL>::load (struct WorldBase *wrld, std::string const &name)
{
    Nrrd *nin = load_nrrd_file (wrld, name);
    if (nin == nullptr)
        return true;
    else {
        bool sts = this->load (wrld, nin);
        nrrdNix (nin);
        return sts;
    }
}

template <typename REAL>
bool image3d<REAL>::load (struct WorldBase *wrld, const Nrrd *nin)
{
    this->_img = std::make_shared<__details::image3d<REAL>>();

  // compute the offset to the first space axis
    int base = nin->dim - nin->spaceDim;

  // check that the image format is what we expect
    if (nin->spaceDim != 3) {
        error (wrld, "nrrd has unexpected dimension %d, expected 3\n", nin->spaceDim);
        return true;
    }
    else if ((nin->axis[base].kind != nrrdKindSpace)
    ||       (nin->axis[base+1].kind != nrrdKindSpace)
    ||       (nin->axis[base+2].kind != nrrdKindSpace)) {
        error (wrld, "nrrd has unexpected axis structure: %s %s %s\n",
            airEnumStr(nrrdKind, nin->axis[base].kind),
            airEnumStr(nrrdKind, nin->axis[base+1].kind),
            airEnumStr(nrrdKind, nin->axis[base+2].kind));
        return true;
    }

    this->_img->_dim     = 3;
    this->_img->_size[0] = nin->axis[base+0].size;
    this->_img->_size[1] = nin->axis[base+1].size;
    this->_img->_size[2] = nin->axis[base+2].size;
    this->_img->_data    = nin->data;
    this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);

  // from the Nrrd file, we load the affine image-to-world transform matrix
    mat3x3<double> m;              	// rotation and scaling
    real3<double> t;                   	// translation part of the transform
    load_transform_3d (nin, m, t);

  // compute inverse of m, which is the transform from the world basis to the image basis
    mat3x3<double> mInv;
    invert (mInv, m);

  // copy results into the image data structure
    for (int i = 0;  i < 3;  i++) {
        for (int j = 0;  j < 3;  j++) {
	    REAL r = static_cast<REAL>(mInv[2*i+j]);
            this->_img->_w2i[3*i+j]  = r;
            this->_img->_w2iT[3*j+i] = r;	// transpose
        }
    }

  // transform the translation vector: inv([M t]) = [inv(M) -inv(M)t]
    for (int i = 0;  i < 3;  i++) {
	double sum = 0;
	for (int j = 0;  j < 3;  j++) {
	    sum += mInv[3*i+j] * t[j];
	}
	this->_img->_tVec[i] = static_cast<REAL>(-sum);
    }

    return false;
}

} // namespace Diderot

void foo (std::string name)
{
    diderot::image1d<float> img1;
    diderot::image2d<float> img2;
    diderot::image3d<float> img3;

    img1.load(nullptr, name);
    img2.load(nullptr, name);
    img3.load(nullptr, name);

    float x = img2.translate()[0];

}

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0