[Insight-developers] interrupting filter execution
Stephen R. Aylward
Stephen.Aylward at Kitware.com
Mon Jul 24 22:47:34 EDT 2006
I agree with Jim - base classes shouldn't provide functions that imply
functionality that cannot be provided by the majority of its subclasses.
Conceptually, the subclasses would have to hide the GetInverse
function by making it protected - seems more confusing than the current
state.
An intermediate class can be provided with backward compatibility. I
don't see invasiveness as an issue.
Stephen
Paul Koshevoy wrote:
> Miller, James V (GE, Research) wrote:
>
>>Paul,
>>
>>I'm worried that not all transforms have an inverse (or the inverse is
>>too complicated to evaluate).
>
>
> This is precisely why I need this functionality. I already have such a transform
> -- a bi-variate Legendre polynomial is not analytically invertible for higher
> order polynomials. What I do now is return a pointer to a generic inverse class
> (derived from itk::Transform), templated by the forward transform class. This is
> a very flexible solution.
>
>
>>I think the way to do this is place another class under Transform that
>>is something like InvertibleTransform. This class would have the
>>virtual method to return the inverse of the transform. MatrixOffsetTransformBase
>>and IdentityTransform would inherit from InvertibleTransform.
>
>
> That sounds way more invasive than what I've been doing so far, what's the benefit?
>
>
>>This may not be the best way to do it because a MatrixOffsetTransformBase
>>not have to be invertible (could be a projection). Perhaps we can satisfy
>>this by having a method IsInvertible() on the InvertibleTransform which
>>verifies whether the current transform is actually invertible.
>
>
> How would the IsInvertible() method be better than testing GetInverse() for NULL?
>
> I've attached 4 files to demonstrate how I use the virtual GetInverse() right now.
>
> Paul.
>
>
> ------------------------------------------------------------------------
>
> #ifndef __itkInverseTransform_h
> #define __itkInverseTransform_h
>
> // system includes:
> #include <assert.h>
>
> // ITK includes:
> #include <itkTransform.h>
> #include <itkMacro.h>
>
>
> //----------------------------------------------------------------
> // itk::InverseTransform
> //
> namespace itk
> {
> //----------------------------------------------------------------
> // InverseTransform
> //
> template <class ForwardTransform>
> class InverseTransform :
> public Transform< typename ForwardTransform::ScalarType,
> ForwardTransform::OutputSpaceDimension,
> ForwardTransform::InputSpaceDimension >
> {
> public:
> /** Standard class typedefs. */
> typedef InverseTransform Self;
>
> typedef Transform< typename ForwardTransform::ScalarType,
> ForwardTransform::OutputSpaceDimension,
> ForwardTransform::InputSpaceDimension > Superclass;
>
> typedef SmartPointer< Self > Pointer;
> typedef SmartPointer< const Self > ConstPointer;
>
> /** Base inverse transform type. */
> typedef typename Superclass::InverseTransformType InverseTransformType;
> typedef SmartPointer< InverseTransformType > InverseTransformPointer;
>
> /** New method for creating an object using a factory. */
> itkNewMacro(Self);
>
> /** Run-time type information (and related methods). */
> itkTypeMacro( InverseTransform, Transform );
>
> /** Standard scalar type for this class. */
> typedef typename Superclass::ScalarType ScalarType;
>
> /** Type of the input parameters. */
> typedef typename Superclass::ParametersType ParametersType;
>
> /** Type of the Jacobian matrix. */
> typedef typename Superclass::JacobianType JacobianType;
>
> /** Standard coordinate point type for this class. */
> typedef typename Superclass::InputPointType InputPointType;
> typedef typename Superclass::OutputPointType OutputPointType;
>
> /** Dimension of the domain space. */
> itkStaticConstMacro(InputSpaceDimension,
> unsigned int,
> ForwardTransform::OutputSpaceDimension);
> itkStaticConstMacro(OutputSpaceDimension,
> unsigned int,
> ForwardTransform::InputSpaceDimension);
>
> /** Set the forward transform pointer. */
> void SetForwardTransform(const ForwardTransform * forward)
> { forward_ = forward; }
>
> /** Method to transform a point. */
> virtual OutputPointType TransformPoint(const InputPointType & y) const
> {
> assert(forward_ != NULL);
> return forward_->BackTransformPoint(y);
> }
>
> virtual const JacobianType & GetJacobian(const InputPointType &) const
> {
> itkExceptionMacro(<< "GetJacobian is not implemented "
> "for InverseTransform");
> return this->m_Jacobian;
> };
>
> virtual unsigned int GetNumberOfParameters() const
> { return 0; }
>
> virtual InverseTransformPointer GetInverse() const
> { return const_cast<ForwardTransform *>(forward_); }
>
> protected:
> InverseTransform(): Superclass(0, 0) {}
>
> private:
> // disable default copy constructor and assignment operator:
> InverseTransform(const Self & other);
> const Self & operator = (const Self & t);
>
> // the transform whose inverse we are trying to evaluate:
> const ForwardTransform * forward_;
> };
>
> } // namespace itk
>
> #endif // __itkInverseTransform_h
>
>
> ------------------------------------------------------------------------
>
> #ifndef __itkLegendrePolynomialTransform_h
> #define __itkLegendrePolynomialTransform_h
>
> // system includes:
> #include <iostream>
> #include <assert.h>
>
> // ITK includes:
> #include "itkTransform.h"
> #include "itkExceptionObject.h"
> #include "itkInverseTransform.h"
> #include "itkMacro.h"
>
> // VXL includes:
> #include <vnl/algo/vnl_svd.h>
>
>
> //----------------------------------------------------------------
> // itk::LegendrePolynomialTransform
> //
> // Let
> // A = (u - uc) / Xmax
> // B = (v - vc) / Ymax
> //
> // where uc, vc correspond to the center of the image expressed in
> // the coordinate system of the mosaic.
> //
> // The transform is defined as
> // x(u, v) = Xmax * Sa
> // y(u, v) = Ymax * Sb
> //
> // where
> // Sa = sum(i in [0, N], sum(j in [0, i], a_jk * Pj(A) * Qk(B)));
> // Sb = sum(i in [0, N], sum(j in [0, i], b_jk * Pj(A) * Qk(B)));
> //
> // where k = i - j and (Pj, Qk) are Legendre polynomials
> // of degree (j, k) respectively.
> //
> namespace itk
> {
> template < class TScalar = double, unsigned int N = 2 >
> class LegendrePolynomialTransform :
> public Transform<TScalar, 2, 2>
> {
> public:
> // standard typedefs:
> typedef LegendrePolynomialTransform Self;
> typedef SmartPointer<Self> Pointer;
> typedef SmartPointer<const Self> ConstPointer;
>
> typedef Transform<TScalar, 2, 2> Superclass;
>
> // Base inverse transform type:
> typedef typename Superclass::InverseTransformType InverseTransformType;
> typedef SmartPointer< InverseTransformType > InverseTransformPointer;
>
> // static constant for the degree of the polynomial:
> itkStaticConstMacro(Degree, unsigned int, N);
>
> // static constant for the number of a_jk (or b_jk) coefficients:
> itkStaticConstMacro(CoefficientsPerDimension,
> unsigned int,
> ((N + 1) * (N + 2)) / 2);
>
> // static constant for the length of the parameter vector:
> itkStaticConstMacro(ParameterVectorLength,
> unsigned int,
> (N + 1) * (N + 2));
>
> // RTTI:
> itkTypeMacro(LegendrePolynomialTransform, Transform);
>
> // macro for instantiation through the object factory:
> itkNewMacro(Self);
>
> /** Standard scalar type for this class. */
> typedef typename Superclass::ScalarType ScalarType;
>
> /** Dimension of the domain space. */
> itkStaticConstMacro(InputSpaceDimension, unsigned int, 2);
> itkStaticConstMacro(OutputSpaceDimension, unsigned int, 2);
>
> // shortcuts:
> typedef typename Superclass::ParametersType ParametersType;
> typedef typename Superclass::JacobianType JacobianType;
>
> typedef typename Superclass::InputPointType InputPointType;
> typedef typename Superclass::OutputPointType OutputPointType;
>
> // virtual:
> OutputPointType TransformPoint(const InputPointType & x) const;
>
> // Inverse transformations:
> // If y = Transform(x), then x = BackTransform(y);
> // if no mapping from y to x exists, then an exception is thrown.
> InputPointType BackTransformPoint(const OutputPointType & y) const;
>
> // virtual:
> void SetFixedParameters(const ParametersType & params)
> { this->m_FixedParameters = params; }
>
> // virtual:
> const ParametersType & GetFixedParameters() const
> { return this->m_FixedParameters; }
>
> // virtual:
> void SetParameters(const ParametersType & params)
> { this->m_Parameters = params; }
>
> // virtual:
> const ParametersType & GetParameters() const
> { return this->m_Parameters; }
>
> // virtual:
> unsigned int GetNumberOfParameters() const
> { return ParameterVectorLength; }
>
> // virtual:
> const JacobianType & GetJacobian(const InputPointType & point) const;
>
> // virtual: return an inverse of this transform.
> InverseTransformPointer GetInverse() const
> {
> typedef InverseTransform<Self> InvTransformType;
> typename InvTransformType::Pointer inv = InvTransformType::New();
> inv->SetForwardTransform(this);
> return inv.GetPointer();
> }
>
> // setup the fixed transform parameters:
> void setup(// image bounding box expressed in the physical space:
> const double x_min,
> const double x_max,
> const double y_min,
> const double y_max,
>
> // normalization parameters:
> const double Xmax = 0.0,
> const double Ymax = 0.0)
> {
> double & uc_ = this->m_FixedParameters[0];
> double & vc_ = this->m_FixedParameters[1];
> double & xmax_ = this->m_FixedParameters[2];
> double & ymax_ = this->m_FixedParameters[3];
> double & a00_ = this->m_Parameters[index_a(0, 0)];
> double & b00_ = this->m_Parameters[index_b(0, 0)];
>
> // center of the image:
> double xc = (x_min + x_max) / 2.0;
> double yc = (y_min + y_max) / 2.0;
> uc_ = xc;
> vc_ = yc;
>
> // setup the normalization parameters:
> if (Xmax != 0.0 && Ymax != 0.0)
> {
> xmax_ = Xmax;
> ymax_ = Ymax;
> }
> else
> {
> const double w = x_max - x_min;
> const double h = y_max - y_min;
>
> // -1 : 1
> xmax_ = w / 2.0;
> ymax_ = h / 2.0;
> }
>
> // setup a00, b00 (local translation parameters):
> a00_ = xc / xmax_;
> b00_ = yc / ymax_;
> }
>
> // setup the translation parameters:
> void setup_translation(// translation is expressed in the physical space:
> const double tx_Xmax = 0.0,
> const double ty_Ymax = 0.0)
> {
> // incorporate translation into the (uc, vc) fixed parameters:
> double & uc_ = this->m_FixedParameters[0];
> double & vc_ = this->m_FixedParameters[1];
>
> // FIXME: the signs might be wrong here (20051101):
> uc_ -= tx_Xmax;
> vc_ -= ty_Ymax;
>
> // FIXME: untested:
> /*
> double & a00_ = this->m_Parameters[index_a(0, 0)];
> double & b00_ = this->m_Parameters[index_b(0, 0)];
> a00_ -= tx_Xmax / GetXmax();
> b00_ -= ty_Ymax / GetYmax();
> */
> }
>
> // helper required for numeric inverse transform calculation;
> // evaluate F = T(x), J = dT/dx (another Jacobian):
> void eval(const std::vector<ScalarType> & x,
> std::vector<ScalarType> & F,
> std::vector<std::vector<ScalarType> > & J) const;
>
> // setup a linear system to solve for the parameters of this
> // transform such that it maps points uv to xy:
> void setup_linear_system(const unsigned int start_with_degree,
> const unsigned int degrees_covered,
> const std::vector<InputPointType> & uv,
> const std::vector<OutputPointType> & xy,
> vnl_matrix<double> & M,
> vnl_vector<double> & bx,
> vnl_vector<double> & by) const;
>
> // find the polynomial coefficients such that this transform
> // would map uv to xy:
> void solve_for_parameters(const unsigned int start_with_degree,
> const unsigned int degrees_covered,
> const std::vector<InputPointType> & uv,
> const std::vector<OutputPointType> & xy,
> ParametersType & params) const;
>
> inline void solve_for_parameters(const unsigned int start_with_degree,
> const unsigned int degrees_covered,
> const std::vector<InputPointType> & uv,
> const std::vector<OutputPointType> & xy)
> {
> ParametersType params = GetParameters();
> solve_for_parameters(start_with_degree, degrees_covered, uv, xy, params);
> SetParameters(params);
> }
>
> // calculate the number of coefficients of a given
> // degree range (per dimension):
> inline static unsigned int
> count_coefficients(const unsigned int start_with_degree,
> const unsigned int degrees_covered)
> {
> return
> index_a(0, start_with_degree + degrees_covered) -
> index_a(0, start_with_degree);
> }
>
> // accessors to the warp origin expressed in the mosaic coordinate system:
> inline const double & GetUc() const
> { return this->m_FixedParameters[0]; }
>
> inline const double & GetVc() const
> { return this->m_FixedParameters[1]; }
>
> // accessors to the normalization parameters Xmax, Ymax:
> inline const double & GetXmax() const
> { return this->m_FixedParameters[2]; }
>
> inline const double & GetYmax() const
> { return this->m_FixedParameters[3]; }
>
> // generate a mask of shared parameters:
> static void setup_shared_params_mask(bool shared, std::vector<bool> & mask)
> {
> mask.assign(ParameterVectorLength, shared);
> mask[index_a(0, 0)] = false;
> mask[index_b(0, 0)] = false;
> }
>
> // Convert the j, k indecies associated with a(j, k) coefficient
> // into an index that can be used with the parameters array:
> inline static unsigned int index_a(const unsigned int & j,
> const unsigned int & k)
> { return j + ((j + k) * (j + k + 1)) / 2; }
>
> inline static unsigned int index_b(const unsigned int & j,
> const unsigned int & k)
> { return CoefficientsPerDimension + index_a(j, k); }
>
> // virtual: Generate a platform independent name:
> std::string GetTransformTypeAsString() const
> {
> std::string base = Superclass::GetTransformTypeAsString();
> OStringStream name;
> name << base << '_' << N;
> return name.str();
> }
>
> protected:
> LegendrePolynomialTransform();
>
> // virtual:
> void PrintSelf(std::ostream & s, Indent indent) const;
>
> private:
> // disable default copy constructor and assignment operator:
> LegendrePolynomialTransform(const Self & other);
> const Self & operator = (const Self & t);
>
> // polynomial coefficient accessors:
> inline const double & a(const unsigned int & j,
> const unsigned int & k) const
> { return this->m_Parameters[index_a(j, k)]; }
>
> inline const double & b(const unsigned int & j,
> const unsigned int & k) const
> { return this->m_Parameters[index_b(j, k)]; }
>
> }; // class LegendrePolynomialTransform
>
> } // namespace itk
>
>
> #ifndef ITK_MANUAL_INSTANTIATION
> #include "itkLegendrePolynomialTransform.txx"
> #endif
>
> #endif // __itkLegendrePolynomialTransform_h
>
>
> ------------------------------------------------------------------------
>
> #ifndef _itkLegendrePolynomialTransform_txx
> #define _itkLegendrePolynomialTransform_txx
>
> // local includes:
> #include "itkLegendrePolynomialTransform.h"
> #include "itkNumericInverse.h"
>
> // system includes:
> #include <vector>
> #include <assert.h>
>
> // FIXME:
> extern int BREAK(unsigned int i);
>
>
> namespace itk
> {
> namespace Legendre
> {
> // Legendre polynomials:
> static double P0(const double & /* x */)
> { return 1.0; }
>
> static double P1(const double & x)
> { return x; }
>
> static double P2(const double & x)
> {
> const double xx = x * x;
> return (3.0 * xx - 1.0) / 2.0;
> }
>
> static double P3(const double & x)
> {
> const double xx = x * x;
> return ((5.0 * xx - 3.0) * x) / 2.0;
> }
>
> static double P4(const double & x)
> {
> const double xx = x * x;
> return ((35.0 * xx - 30.0) * xx + 3.0) / 8.0;
> }
>
> static double P5(const double & x)
> {
> const double xx = x * x;
> return (((63.0 * xx - 70.0) * xx + 15.0) * x) / 8.0;
> }
>
> static double P6(const double & x)
> {
> const double xx = x * x;
> return (((231.0 * xx - 315.0) * xx + 105.0) * xx - 5.0) / 16.0;
> }
>
> //----------------------------------------------------------------
> // polynomial_t
> //
> typedef double(*polynomial_t)(const double & x);
>
> //----------------------------------------------------------------
> // A partial table of the Legendre polynomials
> //
> static polynomial_t P[] = { P0, P1, P2, P3, P4, P5, P6 };
>
> // first derivatives of the Legendre polynomials:
> static double dP0(const double & /* x */)
> { return 0.0; }
>
> static double dP1(const double & /* x */)
> { return 1.0; }
>
> static double dP2(const double & x)
> { return 3.0 * x; }
>
> static double dP3(const double & x)
> {
> const double xx = x * x;
> return (15.0 * xx - 3.0) / 2.0;
> }
>
> static double dP4(const double & x)
> {
> const double xx = x * x;
> return ((35.0 * xx - 15.0) * x) / 2.0;
> }
>
> static double dP5(const double & x)
> {
> const double xx = x * x;
> return ((315.0 * xx - 210.0) * xx + 15.0) / 8.0;
> }
>
> static double dP6(const double & x)
> {
> const double xx = x * x;
> return (((693.0 * xx - 630.0) * xx + 105.0) * x) / 8.0;
> }
>
> //----------------------------------------------------------------
> // A partial table of the derivatives of Legendre polynomials
> //
> static polynomial_t dP[] = { dP0, dP1, dP2, dP3, dP4, dP5, dP6 };
> }
>
> //----------------------------------------------------------------
> // LegendrePolynomialTransform
> //
> template <class TScalarType, unsigned int N>
> LegendrePolynomialTransform<TScalarType, N>::
> LegendrePolynomialTransform():
> Superclass(2, ParameterVectorLength)
> {
> // by default, the parameters are initialized for an identity transform:
> // initialize the parameters for an identity transform:
> for (unsigned int i = 0; i < ParameterVectorLength; i++)
> {
> this->m_Parameters[i] = 0.0;
> }
>
> this->m_Parameters[index_a(1, 0)] = 1.0;
> this->m_Parameters[index_b(0, 1)] = 1.0;
>
> // allocate some space for the fixed parameters:
> this->m_FixedParameters.SetSize(4);
> setup(-1, 1, -1, 1);
>
> this->m_Parameters[index_a(0, 0)] = GetUc() / GetXmax();
> this->m_Parameters[index_b(0, 0)] = GetVc() / GetYmax();
>
> // zero-out the Jacobian:
> for (unsigned int i = 0; i < ParameterVectorLength; i++)
> {
> this->m_Jacobian(0, i) = 0.0;
> this->m_Jacobian(1, i) = 0.0;
> }
> }
>
> //----------------------------------------------------------------
> // TransformPoint
> //
> template <class TScalarType, unsigned int N>
> typename LegendrePolynomialTransform<TScalarType, N>::OutputPointType
> LegendrePolynomialTransform<TScalarType, N>::
> TransformPoint(const InputPointType & x) const
> {
> const double & uc = this->GetUc();
> const double & vc = this->GetVc();
> const double & Xmax = this->GetXmax();
> const double & Ymax = this->GetYmax();
>
> const ScalarType & u = x[0];
> const ScalarType & v = x[1];
>
> const double A = (u - uc) / Xmax;
> const double B = (v - vc) / Ymax;
>
> double P[N + 1];
> double Q[N + 1];
> for (unsigned int i = 0; i <= N; i++)
> {
> P[i] = Legendre::P[i](A);
> Q[i] = Legendre::P[i](B);
> }
>
> double Sa = 0.0;
> double Sb = 0.0;
>
> for (unsigned int i = 0; i <= N; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> unsigned int k = i - j;
> double PjQk = P[j] * Q[k];
> Sa += a(j, k) * PjQk;
> Sb += b(j, k) * PjQk;
> }
> }
>
> OutputPointType y;
> y[0] = Xmax * Sa;
> y[1] = Ymax * Sb;
>
> return y;
> }
>
> //----------------------------------------------------------------
> // BackTransformPoint
> //
> template <class TScalarType, unsigned int N>
> typename LegendrePolynomialTransform<TScalarType, N>::InputPointType
> LegendrePolynomialTransform<TScalarType, N>::
> BackTransformPoint(const OutputPointType & y) const
> {
> NumericInverse<LegendrePolynomialTransform<ScalarType, N> > inverse(*this);
>
> std::vector<ScalarType> vy(2);
> std::vector<ScalarType> vx(2);
> vy[0] = y[0];
> vy[1] = y[1];
>
> // initialize x: first guess - x is close to y:
> vx = vy;
> const double & uc = this->GetUc();
> const double & vc = this->GetVc();
>
> vx[0] += uc;
> vx[1] += vc;
> bool ok = inverse.transform(vy, vx, true);
> if (!ok)
> {
> // FIXME:
> vx[0] = vy[0] + uc;
> vx[1] = vy[1] + vc;
> ok = inverse.transform(vy, vx, true);
> }
> assert(ok);
>
> OutputPointType x;
> x[0] = vx[0];
> x[1] = vx[1];
>
> return x;
> }
>
> //----------------------------------------------------------------
> // GetJacobian
> //
> template <class TScalarType, unsigned int N>
> const typename LegendrePolynomialTransform<TScalarType, N>::JacobianType &
> LegendrePolynomialTransform<TScalarType, N>::
> GetJacobian(const InputPointType & x) const
> {
> const double & uc = this->GetUc();
> const double & vc = this->GetVc();
> const double & Xmax = this->GetXmax();
> const double & Ymax = this->GetYmax();
>
> const ScalarType & u = x[0];
> const ScalarType & v = x[1];
>
> const double A = (u - uc) / Xmax;
> const double B = (v - vc) / Ymax;
>
> double P[N + 1];
> double Q[N + 1];
> for (unsigned int i = 0; i <= N; i++)
> {
> P[i] = Legendre::P[i](A);
> Q[i] = Legendre::P[i](B);
> }
>
> // derivatives with respect to a_jk, b_jk:
> for (unsigned int i = 0; i <= N; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> unsigned int k = i - j;
>
> double PjQk = P[j] * Q[k];
> this->m_Jacobian(0, index_a(j, k)) = Xmax * PjQk;
> this->m_Jacobian(1, index_b(j, k)) = Ymax * PjQk;
> }
> }
>
> // done:
> return this->m_Jacobian;
> }
>
> //----------------------------------------------------------------
> // eval
> //
> template <class TScalarType, unsigned int N>
> void
> LegendrePolynomialTransform<TScalarType, N>::
> eval(const std::vector<ScalarType> & x,
> std::vector<ScalarType> & F,
> std::vector<std::vector<ScalarType> > & J) const
> {
> const double & uc = this->GetUc();
> const double & vc = this->GetVc();
> const double & Xmax = this->GetXmax();
> const double & Ymax = this->GetYmax();
>
> const ScalarType & u = x[0];
> const ScalarType & v = x[1];
>
> const double A = (u - uc) / Xmax;
> const double B = (v - vc) / Ymax;
>
> double P[N + 1];
> double Q[N + 1];
> for (unsigned int i = 0; i <= N; i++)
> {
> P[i] = Legendre::P[i](A);
> Q[i] = Legendre::P[i](B);
> }
>
> double Sa = 0.0;
> double Sb = 0.0;
>
> // derivatives with respect to a_jk, b_jk:
> for (unsigned int i = 0; i <= N; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> unsigned int k = i - j;
> double PjQk = P[j] * Q[k];
> Sa += a(j, k) * PjQk;
> Sb += b(j, k) * PjQk;
> }
> }
>
> F[0] = Xmax * Sa;
> F[1] = Ymax * Sb;
>
> // derivatives with respect to u:
> double dSa_du = 0.0;
> double dSb_du = 0.0;
>
> for (unsigned int i = 1; i <= N; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> unsigned int k = i - j;
> double dPjQk = Legendre::dP[j](A) * Q[k];
> dSa_du += a(j, k) * dPjQk;
> dSb_du += b(j, k) * dPjQk;
> }
> }
>
> // derivatives with respect to v:
> double dSa_dv = 0.0;
> double dSb_dv = 0.0;
>
> for (unsigned int i = 1; i <= N; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> unsigned int k = i - j;
> double dQkPj = Legendre::dP[k](B) * P[j];
> dSa_dv += a(j, k) * dQkPj;
> dSb_dv += b(j, k) * dQkPj;
> }
> }
>
> // dx/du:
> J[0][0] = dSa_du;
>
> // dx/dv:
> J[0][1] = Xmax / Ymax * dSa_dv;
>
> // dy/du:
> J[1][0] = Ymax / Xmax * dSb_du;
>
> // dy/dv:
> J[1][1] = dSb_dv;
> }
>
> //----------------------------------------------------------------
> // setup_linear_system
> //
> template <class TScalarType, unsigned int N>
> void
> LegendrePolynomialTransform<TScalarType, N>::
> setup_linear_system(const unsigned int start_with_degree,
> const unsigned int degrees_covered,
> const std::vector<InputPointType> & uv,
> const std::vector<OutputPointType> & xy,
> vnl_matrix<double> & M,
> vnl_vector<double> & bx,
> vnl_vector<double> & by) const
> {
> if (degrees_covered == 0) return;
>
> const double & uc = this->GetUc();
> const double & vc = this->GetVc();
> const double & Xmax = this->GetXmax();
> const double & Ymax = this->GetYmax();
>
> static double P[N + 1];
> static double Q[N + 1];
>
> const unsigned int & num_points = uv.size();
> assert(num_points == xy.size());
>
> const unsigned int offset =
> index_a(0, start_with_degree);
>
> const unsigned int extent =
> index_a(0, start_with_degree + degrees_covered) - offset;
>
> M.set_size(num_points, extent);
> bx.set_size(num_points);
> by.set_size(num_points);
>
> for (unsigned int row = 0; row < num_points; row++)
> {
> const ScalarType & u = uv[row][0];
> const ScalarType & v = uv[row][1];
> const ScalarType & x = xy[row][0];
> const ScalarType & y = xy[row][1];
>
> const double A = (u - uc) / Xmax;
> const double B = (v - vc) / Ymax;
>
> bx[row] = x / Xmax;
> by[row] = y / Ymax;
>
> for (unsigned int i = 0; i < start_with_degree + degrees_covered; i++)
> {
> P[i] = Legendre::P[i](A);
> Q[i] = Legendre::P[i](B);
> }
>
> for (unsigned int i = 0; i < start_with_degree; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> const unsigned int k = i - j;
> double PjQk = P[j] * Q[k];
> bx[row] -= a(j, k) * PjQk;
> by[row] -= b(j, k) * PjQk;
> }
> }
>
> for (unsigned int i = start_with_degree;
> i < start_with_degree + degrees_covered;
> i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> const unsigned int k = i - j;
> const unsigned int col = index_a(j, k) - offset;
>
> M[row][col] = P[j] * Q[k];
> }
> }
> }
> }
>
> //----------------------------------------------------------------
> // solve_for_parameters
> //
> template <class TScalarType, unsigned int N>
> void
> LegendrePolynomialTransform<TScalarType, N>::
> solve_for_parameters(const unsigned int start_with_degree,
> const unsigned int degrees_covered,
> const std::vector<InputPointType> & uv,
> const std::vector<OutputPointType> & xy,
> ParametersType & params) const
> {
> if (degrees_covered == 0) return;
>
> vnl_matrix<double> M;
> vnl_vector<double> bx;
> vnl_vector<double> by;
>
> setup_linear_system(start_with_degree, degrees_covered, uv, xy, M, bx, by);
> vnl_svd<double> svd(M);
>
> vnl_vector<double> xa = svd.solve(bx);
> vnl_vector<double> xb = svd.solve(by);
>
> unsigned int offset = index_a(0, start_with_degree);
> for (unsigned int i = start_with_degree;
> i < start_with_degree + degrees_covered;
> i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> const unsigned int k = i - j;
> const unsigned int a_jk = index_a(j, k);
> const unsigned int b_jk = index_b(j, k);
>
> params[a_jk] = xa[a_jk - offset];
> params[b_jk] = xb[a_jk - offset];
> }
> }
> }
>
> //----------------------------------------------------------------
> // PrintSelf
> //
> template <class TScalarType, unsigned int N>
> void
> LegendrePolynomialTransform<TScalarType, N>::
> PrintSelf(std::ostream & os, Indent indent) const
> {
> Superclass::PrintSelf(os,indent);
>
> for (unsigned int i = 0; i <= N; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> unsigned int k = i - j;
> os << indent << "a(" << j << ", " << k << ") = " << a(j, k)
> << std::endl;
> }
> }
> for (unsigned int i = 0; i <= N; i++)
> {
> for (unsigned int j = 0; j <= i; j++)
> {
> unsigned int k = i - j;
> os << indent << "b(" << j << ", " << k << ") = " << b(j, k)
> << std::endl;
> }
> }
> os << indent << "uc = " << GetUc() << std::endl
> << indent << "vc = " << GetVc() << std::endl
> << indent << "Xmax = " << GetXmax() << std::endl
> << indent << "Ymax = " << GetYmax() << std::endl;
> #if 0
> for (unsigned int i = 0; i < ParameterVectorLength; i++)
> {
> os << indent << "params[" << i << "] = " << this->m_Parameters[i] << ';'
> << std::endl;
> }
> #endif
> }
>
> } // namespace itk
>
>
> #endif // _itkLegendrePolynomialTransform_txx
>
>
> ------------------------------------------------------------------------
>
> #ifndef __itkNumericInverse_h
> #define __itkNumericInverse_h
>
> // system includes:
> #include <iostream>
> #include <cassert>
> #include <vector>
>
> // ITK includes:
> #include <itkTransform.h>
> #include <itkMacro.h>
>
> namespace help
> {
> //----------------------------------------------------------------
> // nonlinear_system_evaluator_t
> //
> template <class ScalarType>
> class nonlinear_system_evaluator_t
> {
> public:
> virtual ~nonlinear_system_evaluator_t() {}
>
> // evaluate the nonlinear system of equations
> // and its Jacobian (dF/dx) at a given point:
> virtual void
> operator() (const std::vector<ScalarType> & x,
> std::vector<ScalarType> & F,
> std::vector<std::vector<ScalarType> > & J) const = 0;
> };
>
> //----------------------------------------------------------------
> // ludcmp
> //
> template <class ScalarType>
> bool
> ludcmp(std::vector<std::vector<ScalarType> > & a,
> std::vector<unsigned int> & indx,
> ScalarType & d)
> {
> // shortcut:
> const unsigned int n = a.size();
>
> const ScalarType TINY = ScalarType(1e-20);
>
> unsigned int imax = 0;
> ScalarType big, dum, temp;
>
> std::vector<ScalarType> vv(n);
> d = ScalarType(1);
>
> for (unsigned int i = 0; i < n; i++)
> {
> big = ScalarType(0);
> for (unsigned int j = 0; j < n; j++)
> {
> if ((temp = fabs(a[i][j])) > big) big = temp;
> }
>
> if (big == ScalarType(0)) return false;
>
> vv[i] = ScalarType(1) / big;
> }
>
> for (unsigned int j = 0; j < n; j++)
> {
> for (unsigned int i = 0; i < j; i++)
> {
> ScalarType sum = a[i][j];
> for (unsigned int k = 0; k < i; k++)
> {
> sum -= a[i][k] * a[k][j];
> }
> a[i][j] = sum;
> }
>
> big = ScalarType(0);
> for (unsigned int i = j; i < n; i++)
> {
> ScalarType sum = a[i][j];
> for (unsigned int k = 0; k < j; k++)
> {
> sum -= a[i][k] * a[k][j];
> }
> a[i][j] = sum;
>
> if ((dum = vv[i] * fabs(sum)) >= big)
> {
> big = dum;
> imax = i;
> }
> }
>
> if (j != imax)
> {
> for (unsigned int k = 0; k < n; k++)
> {
> dum = a[imax][k];
> a[imax][k] = a[j][k];
> a[j][k] = dum;
> }
>
> d = -d;
> vv[imax] = vv[j];
> }
>
> indx[j] = imax;
> if (a[j][j] == ScalarType(0))
> {
> a[j][j] = TINY;
> }
>
> if (j != n)
> {
> dum = ScalarType(1) / a[j][j];
> for (unsigned int i = j + 1; i < n; i++)
> {
> a[i][j] *= dum;
> }
> }
> }
>
> return true;
> }
>
> //----------------------------------------------------------------
> // lubksb
> //
> template <class ScalarType>
> void
> lubksb(const std::vector<std::vector<ScalarType> > & a,
> const std::vector<unsigned int> & indx,
> std::vector<ScalarType> & b)
> {
> // shortcut:
> const unsigned int n = a.size();
>
> unsigned int ii = 0;
> for (unsigned int i = 0; i < n; i++)
> {
> const unsigned int & ip = indx[i];
> ScalarType sum = b[ip];
> b[ip] = b[i];
>
> if (ii)
> {
> for (unsigned int j = ii; j <= i - 1; j++)
> {
> sum -= a[i][j] * b[j];
> }
> }
> else if (sum)
> {
> ii = i;
> }
>
> b[i] = sum;
> }
>
> for (int i = n - 1; i >= 0; i--)
> {
> ScalarType sum = b[i];
> for (unsigned int j = i + 1; j < n; j++)
> {
> sum -= a[i][j] * b[j];
> }
>
> b[i] = sum / a[i][i];
> }
> }
>
> //----------------------------------------------------------------
> // mnewt
> //
> template <class ScalarType>
> bool
> mnewt(const nonlinear_system_evaluator_t<ScalarType> & usrfun,
> std::vector<ScalarType> & x,
> const unsigned int & ntrial,
> const ScalarType & tolx,
> const ScalarType & tolf)
> {
> // shortcut:
> const unsigned int n = x.size();
>
> std::vector<unsigned int> indx(n);
> std::vector<ScalarType> p(n);
> std::vector<ScalarType> fvec(n);
>
> std::vector<std::vector<ScalarType> > fjac(n);
> for (unsigned int i = 0; i < n; i++) fjac[i].resize(n);
>
> for (unsigned int k = 0; k < ntrial; k++)
> {
> usrfun(x, fvec, fjac);
>
> ScalarType errf = ScalarType(0);
> for (unsigned int i = 0; i < n; i++)
> {
> errf += fabs(fvec[i]);
> }
> if (errf <= tolf) break;
>
> for (unsigned int i = 0; i < n; i++)
> {
> p[i] = -fvec[i];
> }
>
> ScalarType d;
> if (!ludcmp(fjac, indx, d)) return false;
> lubksb(fjac, indx, p);
>
> ScalarType errx = ScalarType(0);
> for (unsigned int i = 0; i < n; i++)
> {
> errx += fabs(p[i]);
> x[i] += p[i];
> }
>
> if (errx <= tolx) break;
> }
>
> return true;
> }
>
> } // namespace help
>
> namespace itk
> {
> //----------------------------------------------------------------
> // NumericInverse
> //
> template <class TTransform>
> class NumericInverse :
> public help::nonlinear_system_evaluator_t<typename TTransform::ScalarType>
> {
> public:
> typedef TTransform TransformType;
> typedef typename TTransform::ScalarType ScalarType;
>
> NumericInverse(const TransformType & transform):
> transform_(transform)
> {}
>
> // virtual:
> void operator() (const std::vector<ScalarType> & x,
> std::vector<ScalarType> & F,
> std::vector<std::vector<ScalarType> > & J) const
> {
> transform_.eval(x, F, J);
>
> const unsigned int n = x.size();
> for (unsigned int i = 0; i < n; i++)
> {
> F[i] -= y_[i];
> }
> }
>
> // If y = Transform(x), then x = BackTransform(y).
> // given y, find x:
> bool transform(const std::vector<ScalarType> & y,
> std::vector<ScalarType> & x,
> bool x_is_initialized = false) const
> {
> y_ = y;
> if (!x_is_initialized) x = y;
> return mnewt(*this, x, 50, 1e-17, 1e-17);
> }
>
> private:
> // the transform whose inverse we are trying to evaluate:
> const TransformType & transform_;
>
> // the point for which we are tryying to find the inverse mapping:
> mutable std::vector<ScalarType> y_;
> };
>
> } // namespace itk
>
> #endif // __itkNumericInverse_h
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers
--
=============================================================
Stephen R. Aylward, Ph.D.
Chief Medical Scientist
Kitware, Inc.
http://www.kitware.com
More information about the Insight-developers
mailing list