ffd.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */
 
 
#ifndef AIMS_REGISTRATION_FFD_H
#define AIMS_REGISTRATION_FFD_H
 
#include <cartodata/volume/volume.h>                  // Volume
#include <aims/io/writer.h>                           // aims::Reader / Writer
#include <aims/math/bspline.h>                        // aims::TabulBSpline
#include <aims/resampling/nearestneighborresampler.h>
#include <soma-io/transformation/transformation.h>    // soma::Transformation
 
class Graph;
 
namespace aims {
 
//==========================================================================
  //   FFD TRANSFORMATION
  //==========================================================================
 
  class FfdTransformation : public soma::Transformation3d
  {
  public:
    FfdTransformation( int dimX = 0, int dimY = 1, int dimZ = 1,
                       float sizeX = 1., float sizeY = 1., float sizeZ = 1. );
    template <typename T>
    FfdTransformation( int dimX, int dimY, int dimZ,
                       const carto::rc_ptr<carto::Volume<T> > & test_volume );
    template <typename T>
    FfdTransformation( int dimX, int dimY, int dimZ,
                       const carto::Volume<T> & test_volume );
 
    FfdTransformation( const FfdTransformation & other );
    FfdTransformation( const carto::rc_ptr<carto::Volume<Point3df> > & other );
    FfdTransformation & operator=( const FfdTransformation & other );
 
    operator const carto::VolumeRef<Point3df> & () const
    { return _ctrlPointDelta; }
    operator carto::VolumeRef<Point3df> & ()
    { return _ctrlPointDelta; }
    operator const carto::rc_ptr<carto::Volume<Point3df> > & () const
    { return _ctrlPointDelta; }
    operator carto::rc_ptr<carto::Volume<Point3df> > & ()
    { return _ctrlPointDelta; }
    carto::rc_ptr<carto::Volume<Point3df> > & ctrlDeformations()
    { return _ctrlPointDelta; }
 
    bool isIdentity() const CARTO_OVERRIDE { return false; }
 
    //--- Control Knots ------------------------------------------------------
    Point3df     getCtrlKnot( int nx, int ny, int nz ) const;
    void         updateCtrlKnot( int nx, int ny, int nz, const Point3df & newCtrlKnot );
    void         updateAllCtrlKnot(
      const carto::rc_ptr<carto::Volume<Point3df> > & newCtrlKnotGrid );
    void         updateAllCtrlKnotFromDeformation(
      const carto::rc_ptr<carto::Volume<Point3df> > & newDeformationGrid );
 
    //--- Modify -------------------------------------------------------------
    void         increaseResolution( const Point3d & addKnots );
    // void         inverseTransform();
    // void         estimateLocalDisplacement( const Point3df & voxelSize);
 
    //--- Deformation --------------------------------------------------------
    Point3dd     deformation( const Point3dd& p_mm ) const;
    Point3dd     ffdCoord( const Point3dd& p_mm ) const
    { return mmToSplineVox(p_mm); }
 
    //--- Parameters ---------------------------------------------------------
    int          dimX() const { return _dimx; }
    int          dimY() const { return _dimy; }
    int          dimZ() const { return _dimz; }
    int getSizeX() const { return _dimx; }
    int getSizeY() const { return _dimy; }
    int getSizeZ() const { return _dimz; }
    std::vector<int> getSize() const
    { std::vector<int> s( 3 ); s[0] = _dimx; s[1] = _dimy; s[2] = _dimz;
      return s; }
    float        sizeX() const { return _vsx; }
    float        sizeY() const { return _vsy; }
    float        sizeZ() const { return _vsz; }
    std::vector<float> getVoxelSize() const
    { std::vector<float> v( 3 ); v[0] = _vsx; v[1] = _vsy; v[2] = _vsz;
      return v; }
    bool isFlat( int i ) const
    {
      switch(i)
      {
        case 0:  return _flatx;
        case 1:  return _flaty;
        case 2:  return _flatz;
        default: return false;
      }
    }
    bool isXFlat() const { return _flatx; }
    bool isYFlat() const { return _flaty; }
    bool isZFlat() const { return _flatz; }
    virtual void updateDimensions();
 
    //--- Debug --------------------------------------------------------------
    void  printControlPointsGrid() const;
    void  writeDebugCtrlKnots( const std::string & filename ) const;
    void  writeDebugDeformations( const std::string & filename,
                                  int dimX, int dimY, int dimZ,
                                  float sizeX, float sizeY, float sizeZ ) const;
 
    //--- Output -------------------------------------------------------------
    void   write( const std::string & filename ) const;
 
  protected:
    //--- Protected methods --------------------------------------------------
    Point3dd     splineVoxToMm( const Point3dd& p ) const;
    Point3dd     mmToSplineVox( const Point3dd& p ) const;
    void         updateGridResolution(
      const carto::rc_ptr<carto::Volume<Point3df> > & newGrid );
 
    virtual Point3dd _deformation( const Point3dd& p_mm ) const = 0;
 
    carto::VolumeRef<Point3df>  _ctrlPointDelta;
    int _dimx, _dimy, _dimz;
    float _vsx, _vsy, _vsz;
    bool _flatx, _flaty, _flatz;
  };
 
  template <typename T>
  inline
  FfdTransformation::FfdTransformation(
    int dimX, int dimY, int dimZ,
    const carto::rc_ptr<carto::Volume<T> > & test_volume ):
    _ctrlPointDelta( dimX, dimY, dimZ ),
    _dimx( dimX ), _dimy( dimY ), _dimz( dimZ ),
    _flatx( dimX == 1 ), _flaty( dimY == 1 ), _flatz( dimZ == 1 )
  {
    _ctrlPointDelta = Point3df(0., 0., 0.);
    _ctrlPointDelta.setVoxelSize(
      _flatx ? test_volume->getVoxelSize()[0] : double(test_volume->getSizeX() - 1) / double(dimX - 1) * test_volume->getVoxelSize()[0],
      _flaty ? test_volume->getVoxelSize()[1] : double(test_volume->getSizeY() - 1) / double(dimY - 1) * test_volume->getVoxelSize()[1],
      _flatz ? test_volume->getVoxelSize()[2] : double(test_volume->getSizeZ() - 1) / double(dimZ - 1) * test_volume->getSizeZ()
    );
    updateDimensions();
  }
 
  template <typename T>
  inline
  FfdTransformation::FfdTransformation( int dimX, int dimY, int dimZ,
                                        const carto::Volume<T> & test_volume ):
    _ctrlPointDelta( dimX, dimY, dimZ ),
    _dimx( dimX ), _dimy( dimY ), _dimz( dimZ ),
    _flatx( dimX == 1 ), _flaty( dimY == 1 ), _flatz( dimZ == 1 )
  {
    _ctrlPointDelta = Point3df(0., 0., 0.);
    std::vector<float> vs = test_volume.getVoxelSize();
    _ctrlPointDelta.setVoxelSize(
      _flatx ? vs[0] : double(test_volume.getSizeX() - 1) / double(dimX - 1) * vs[0],
      _flaty ? vs[1] : double(test_volume.getSizeY() - 1) / double(dimY - 1) * vs[1],
      _flatz ? vs[2] : double(test_volume.getSizeZ() - 1) / double(dimZ - 1) * vs[2]
    );
    updateDimensions();
  }
 
  inline Point3dd
  FfdTransformation::splineVoxToMm( const Point3dd & p ) const
  {
    return Point3dd( p[0] * sizeX(),
                     p[1] * sizeY(),
                     p[2] * sizeZ() );
  }
 
  inline Point3dd
  FfdTransformation::mmToSplineVox( const Point3dd & p ) const
  {
    return Point3dd( p[0] / sizeX(),
                     p[1] / sizeY(),
                     p[2] / sizeZ() );
  }
 
  //==========================================================================
  //   FFD TRANSFORMATION
  //==========================================================================
 
  class SplineFfd: public FfdTransformation
  {
   public:
    //--- Constructor --------------------------------------------------------
    SplineFfd( int dimX = 0, int dimY = 1, int dimZ = 1,
                  float sizeX = 1., float sizeY = 1., float sizeZ = 1. );
    template <typename T>
    SplineFfd( int dimX, int dimY, int dimZ,
                  const carto::rc_ptr<carto::Volume<T> > & test_volume );
    template <typename T>
    SplineFfd( int dimX, int dimY, int dimZ,
                  const carto::Volume<T> & test_volume );
 
    SplineFfd( const SplineFfd & other );
    SplineFfd( const carto::rc_ptr<carto::Volume<Point3df> > & other );
    SplineFfd & operator=( const SplineFfd & other );
 
    virtual void updateDimensions();
    //--- Deformation --------------------------------------------------------
    Point3dd     deformation( const Point3dd& p_mm ) const;
    double spline3( double x ) const { return _spline(x); }
 
  private:
    virtual Point3dd transformDouble( double x, double y, double z ) const;
    virtual Point3dd _deformation( const Point3dd& p_mm ) const
    { return deformation_private(p_mm); };
    Point3dd deformation_private( const Point3dd& p_mm ) const;
 
    const aims::TabulBSpline _spline;
    std::vector<int> _mirrorcoefvecx;
    std::vector<int> _mirrorcoefvecy;
    std::vector<int> _mirrorcoefvecz;
    int *_mirrorcoefx;
    int *_mirrorcoefy;
    int *_mirrorcoefz;
  };
 
  template <typename T>
  inline
  SplineFfd::SplineFfd( int dimX, int dimY, int dimZ,
                        const carto::rc_ptr<carto::Volume<T> > & test_volume ):
    FfdTransformation( dimX, dimY, dimZ, test_volume ),
    _spline(3, 0)
  {
  }
 
  template <typename T>
  inline
  SplineFfd::SplineFfd( int dimX, int dimY, int dimZ,
                        const carto::Volume<T> & test_volume ):
    FfdTransformation( dimX, dimY, dimZ, test_volume ),
    _spline(3, 0)
  {
  }
 
 
  //==========================================================================
  //   FFD TRILINEAR RESAMPLED TRANSFORMATION
  //==========================================================================
 
  class TrilinearFfd: public FfdTransformation
  {
   public:
    //--- Constructor --------------------------------------------------------
    TrilinearFfd( int dimX = 0, int dimY = 1, int dimZ = 1,
                  float sizeX = 1., float sizeY = 1., float sizeZ = 1. );
    template <typename T>
    TrilinearFfd( int dimX, int dimY, int dimZ,
                  const carto::rc_ptr<carto::Volume<T> > & test_volume );
 
    TrilinearFfd( const TrilinearFfd & other );
    TrilinearFfd( const carto::rc_ptr<carto::Volume<Point3df> > & other );
    TrilinearFfd & operator=( const TrilinearFfd & other );
 
    //--- Deformation --------------------------------------------------------
    Point3dd deformation( const Point3dd& p_mm ) const;
  private:
    virtual Point3dd transformDouble( double x, double y, double z ) const;
    virtual Point3dd _deformation( const Point3dd& p_mm ) const
    { return deformation_private(p_mm); };
    Point3dd deformation_private( const Point3dd& p_mm ) const;
  };
 
  template <typename T>
  inline
  TrilinearFfd::TrilinearFfd(
    int dimX, int dimY, int dimZ,
    const carto::rc_ptr<carto::Volume<T> > & test_volume ):
    FfdTransformation( dimX, dimY, dimZ, test_volume )
  {
  }
 
//============================================================================
//   FFD READER/WRITER
//============================================================================
 
  template<>
  class Writer<aims::FfdTransformation>
    : public Writer<carto::Volume<Point3df> >
  {
    typedef Writer<carto::Volume<Point3df> > base;
  public:
    inline Writer(): base() {}
    inline Writer( const std::string& filename,
                   carto::Object options = carto::none() ):
      base( filename, options ) {}
    virtual ~Writer() {}
    virtual bool write( const aims::FfdTransformation & obj,
                        bool ascii = false,
                        const std::string* format = 0 );
  };
 
} // namespace aims
 
 
// ---
 
namespace soma
{
 
  template <> inline std::string
  DataTypeCode<aims::FfdTransformation>::objectType()
  {
    return "FFDTransformation";
  }
 
  template <> inline std::string
  DataTypeCode<aims::FfdTransformation>::dataType()
  {
    return "VOID";
  }
 
  template <> inline std::string
  DataTypeCode<aims::FfdTransformation>::name()
  {
    return "FFDTransformation";
  }
 
}
 
 
#endif