/vol/vipdata/irtk/packages/transformation/include/irtkFreeFormTransformation3D.h

/*=========================================================================

  Library   : Image Registration Toolkit (IRTK)
  Module    : $Id: irtkFreeFormTransformation3D.h 2 2008-12-23 12:40:14Z dr $
  Copyright : Imperial College, Department of Computing
              Visual Information Processing (VIP), 2008 onwards
  Date      : $Date: 2008-12-23 12:40:14 +0000 (Tue, 23 Dec 2008) $
  Version   : $Revision: 2 $
  Changes   : $Author: dr $

=========================================================================*/

#ifndef _IRTKFREEFORMTRANSFORMATION3D_H

#define _IRTKFREEFORMTRANSFORMATION3D_H

#include <irtkGeometry.h>

#define FFDLOOKUPTABLESIZE 1000

class irtkFreeFormTransformation3D : public irtkFreeFormTransformation
{

protected:

  int _x;

  int _y;

  int _z;

  double _dx;

  double _dy;

  double _dz;

  double _xaxis[3];

  double _yaxis[3];

  double _zaxis[3];

  irtkPoint _origin;

  irtkMatrix _matL2W;

  irtkMatrix _matW2L;

  double ***_xdata;

  double ***_ydata;

  double ***_zdata;

  static double ***Allocate  (double ***, int, int, int);

  static double ***Deallocate(double ***, int, int, int);

  virtual void UpdateMatrix();

public:

  virtual int GetX() const;

  virtual int GetY() const;

  virtual int GetZ() const;

  virtual int NumberOfDOFs() const;

  virtual void GetSpacing(double &, double &, double &) const;

  virtual void  PutOrientation(double *, double *, double *);

  virtual void  GetOrientation(double *, double *, double *) const;

  virtual void   Put(int, double);

  virtual void   Put(int, int, int, double, double, double);

  virtual double Get(int) const;

  virtual void   Get(int, int, int, double &, double &, double &) const;

  virtual void   PutStatus(int, int, int, _Status);

  virtual void   PutStatus(int, int, int, _Status, _Status, _Status);

  virtual void   PutStatus(int, _Status);

  virtual void   GetStatus(int, int, int, _Status &, _Status &, _Status &);

  virtual _Status GetStatus(int);

  virtual double Bending(double x, double y, double z) = 0;

  virtual void WorldToLattice(double &, double &, double &) const;

  virtual void WorldToLattice(irtkPoint &) const;

  virtual void LatticeToWorld(double &, double &, double &) const;

  virtual void LatticeToWorld(irtkPoint &) const;

  virtual void IndexToLattice(int index, int& i, int& j, int& k) const;

  virtual int  LatticeToIndex(int i, int j, int k) const;

  virtual void  ControlPointLocation(int, double &, double &, double &) const;

  virtual irtkPoint ControlPointLocation(int) const;

  virtual void PutBoundingBox(irtkPoint, irtkPoint);

  virtual void PutBoundingBox(double, double, double,
                              double, double, double);

  virtual void BoundingBox(irtkPoint &, irtkPoint &) const;

  virtual void BoundingBox(double &, double &, double &,
                           double &, double &, double &) const;

  virtual void BoundingBox(int, irtkPoint &, irtkPoint &, double = 1) const = 0;

  virtual void BoundingBox(int, double &, double &, double &,
                           double &, double &, double &, double = 1) const = 0;

  virtual void BoundingBox(irtkGreyImage *, int, int &, int &, int &,
                           int &, int &, int &, double = 1) const = 0;

  virtual double Approximate(double *, double *, double *,
                             double *, double *, double *, int) = 0;

  virtual void Interpolate(double* dxs, double* dys, double* dzs) = 0;

  virtual double Inverse(double &, double &, double &, double, double = 0.01);

  virtual Bool IsIdentity();

  virtual const char *NameOfClass();

};

inline int irtkFreeFormTransformation3D::GetX() const
{
  return _x;
}

inline int irtkFreeFormTransformation3D::GetY() const
{
  return _y;
}

inline int irtkFreeFormTransformation3D::GetZ() const
{
  return _z;
}

inline int irtkFreeFormTransformation3D::NumberOfDOFs() const
{
  return 3*_x*_y*_z;
}

inline double irtkFreeFormTransformation3D::Get(int index) const
{
  int i, j, k;

  if (index >= 3*_x*_y*_z) {
    cerr << "irtkFreeFormTransformation3D::Get: No such dof" << endl;
    exit(1);
  }
  if (index < _x*_y*_z) {
    i = index/(_y*_z);
    j = index%(_y*_z)/_z;
    k = index%(_y*_z)%_z;
    return _xdata[k][j][i];
  } else {
    if (index < 2*_x*_y*_z) {
      index -= _x*_y*_z;
      i = index/(_y*_z);
      j = index%(_y*_z)/_z;
      k = index%(_y*_z)%_z;
      return _ydata[k][j][i];
    } else {
      index -= 2*_x*_y*_z;
      i = index/(_y*_z);
      j = index%(_y*_z)/_z;
      k = index%(_y*_z)%_z;
      return _zdata[k][j][i];
    }
  }
}

inline void irtkFreeFormTransformation3D::Get(int i, int j, int k, double &x, double &y, double &z) const
{
  if ((i < 0) || (i >= _x) || (j < 0) || (j >= _y) || (k < 0) || (k >= _z)) {
    cerr << "irtkFreeFormTransformation3D::Get: No such dof" << endl;
    exit(1);
  }
  x = _xdata[k][j][i];
  y = _ydata[k][j][i];
  z = _zdata[k][j][i];
}

inline void irtkFreeFormTransformation3D::Put(int index, double x)
{
  irtkPoint p;
  int i, j, k;

  if (index >= 3*_x*_y*_z) {
    cerr << "irtkFreeFormTransformation3D::Put: No such dof" << endl;
    exit(1);
  }
  if (index < _x*_y*_z) {
    i = index/(_y*_z);
    j = index%(_y*_z)/_z;
    k = index%(_y*_z)%_z;
    _xdata[k][j][i] = x;
  } else {
    if (index < 2*_x*_y*_z) {
      index -= _x*_y*_z;
      i = index/(_y*_z);
      j = index%(_y*_z)/_z;
      k = index%(_y*_z)%_z;
      _ydata[k][j][i] = x;
    } else {
      index -= 2*_x*_y*_z;
      i = index/(_y*_z);
      j = index%(_y*_z)/_z;
      k = index%(_y*_z)%_z;
      _zdata[k][j][i] = x;
    }
  }
}

inline void irtkFreeFormTransformation3D::Put(int i, int j, int k, double x, double y, double z)
{
  if ((i < 0) || (i >= _x) || (j < 0) || (j >= _y) || (k < 0) || (k >= _z)) {
    cerr << "irtkFreeFormTransformation3D::Put: No such dof" << endl;
    exit(1);
  }
  _xdata[k][j][i] = x;
  _ydata[k][j][i] = y;
  _zdata[k][j][i] = z;
}

inline void irtkFreeFormTransformation3D::GetSpacing(double &dx, double &dy, double &dz) const
{
  dx = _dx;
  dy = _dy;
  dz = _dz;
}

inline void irtkFreeFormTransformation3D::PutOrientation(double *xaxis, double *yaxis, double *zaxis)
{
  _xaxis[0] = xaxis[0];
  _xaxis[1] = xaxis[1];
  _xaxis[2] = xaxis[2];
  _yaxis[0] = yaxis[0];
  _yaxis[1] = yaxis[1];
  _yaxis[2] = yaxis[2];
  _zaxis[0] = zaxis[0];
  _zaxis[1] = zaxis[1];
  _zaxis[2] = zaxis[2];

  this->UpdateMatrix();
}

inline void irtkFreeFormTransformation3D::GetOrientation(double *xaxis, double *yaxis, double *zaxis) const
{
  xaxis[0] = _xaxis[0];
  xaxis[1] = _xaxis[1];
  xaxis[2] = _xaxis[2];
  yaxis[0] = _yaxis[0];
  yaxis[1] = _yaxis[1];
  yaxis[2] = _yaxis[2];
  zaxis[0] = _zaxis[0];
  zaxis[1] = _zaxis[1];
  zaxis[2] = _zaxis[2];
}

inline void irtkFreeFormTransformation3D::PutBoundingBox(double x1, double y1, double z1, double x2, double y2, double z2)
{
  // Initialize control point domain
  _origin._x = (x2 + x1) / 2.0;
  _origin._y = (y2 + y1) / 2.0;
  _origin._z = (z2 + z1) / 2.0;

  double a = x1 * _xaxis[0] + y1 * _xaxis[1] + z1 * _xaxis[2];
  double b = x1 * _yaxis[0] + y1 * _yaxis[1] + z1 * _yaxis[2];
  double c = x1 * _zaxis[0] + y1 * _zaxis[1] + z1 * _zaxis[2];
  x1 = a;
  y1 = b;
  z1 = c;
  a = x2 * _xaxis[0] + y2 * _xaxis[1] + z2 * _xaxis[2];
  b = x2 * _yaxis[0] + y2 * _yaxis[1] + z2 * _yaxis[2];
  c = x2 * _zaxis[0] + y2 * _zaxis[1] + z2 * _zaxis[2];
  x2 = a;
  y2 = b;
  z2 = c;

  // Initialize control point spacing
  _dx = (x2 - x1) / (_x - 1);
  _dy = (y2 - y1) / (_y - 1);
  if (z2 > z1) {
    _dz = (z2 - z1) / (_z - 1);
  } else {
    _dz = 1;
  }

  this->UpdateMatrix();
}

inline void irtkFreeFormTransformation3D::PutBoundingBox(irtkPoint p1, irtkPoint p2)
{
  this->PutBoundingBox(p1._x, p1._y, p1._z, p2._x, p2._y, p2._z);
}

inline void irtkFreeFormTransformation3D::WorldToLattice(double &x, double &y, double &z) const
{
  double a, b, c;

  // Pre-multiply point with transformation matrix
  a = _matW2L(0, 0)*x+_matW2L(0, 1)*y+_matW2L(0, 2)*z+_matW2L(0, 3);
  b = _matW2L(1, 0)*x+_matW2L(1, 1)*y+_matW2L(1, 2)*z+_matW2L(1, 3);
  c = _matW2L(2, 0)*x+_matW2L(2, 1)*y+_matW2L(2, 2)*z+_matW2L(2, 3);

  // Copy result back
  x = a;
  y = b;
  z = c;
}

inline void irtkFreeFormTransformation3D::WorldToLattice(irtkPoint &p) const
{
  this->WorldToLattice(p._x, p._y, p._z);
}

inline void irtkFreeFormTransformation3D::LatticeToWorld(double &x, double &y, double &z) const
{
  double a, b, c;

  // Pre-multiply point with transformation matrix
  a = _matL2W(0, 0)*x+_matL2W(0, 1)*y+_matL2W(0, 2)*z+_matL2W(0, 3);
  b = _matL2W(1, 0)*x+_matL2W(1, 1)*y+_matL2W(1, 2)*z+_matL2W(1, 3);
  c = _matL2W(2, 0)*x+_matL2W(2, 1)*y+_matL2W(2, 2)*z+_matL2W(2, 3);

  // Copy result back
  x = a;
  y = b;
  z = c;
}

inline void irtkFreeFormTransformation3D::LatticeToWorld(irtkPoint &p) const
{
  this->LatticeToWorld(p._x, p._y, p._z);
}

inline void irtkFreeFormTransformation3D::IndexToLattice(int index, int& i, int& j, int& k) const
{

  if (index >= _x*_y*_z) {
    index -= _x*_y*_z;
    if (index >= _x*_y*_z) {
      index -= _x*_y*_z;
    }
  }
  i = index/(_y*_z);
  j = index%(_y*_z)/_z;
  k = index%(_y*_z)%_z;
}

inline int irtkFreeFormTransformation3D::LatticeToIndex(int i, int j, int k) const
{
  return i * _y * _z + j * _z + k;
}

inline const char *irtkFreeFormTransformation3D::NameOfClass()
{
  return "irtkFreeFormTransformation3D";
}

#endif

Generated on Mon Mar 23 16:55:53 2009 for IRTK by Doxygen 1.5.4