From eced0414e93420f93e6c128579ed2e8fc2afa09e Mon Sep 17 00:00:00 2001
From: Matt McCormick <matt.mccormick@kitware.com>
Date: Tue, 4 Jun 2024 16:20:28 -0400
Subject: [PATCH] COMP: Transform API updates

---
 include/itkGridTransform.h                 |  27 +-
 include/itkIRCommon.h                      |  34 +-
 include/itkLegendrePolynomialTransform.h   | 609 +++++++-------
 include/itkLegendrePolynomialTransform.hxx | 906 +++++++++++----------
 include/itkMeshTransform.h                 |  28 +-
 include/itkRBFTransform.h                  |  16 +-
 include/itkRadialDistortionTransform.h     |  14 +-
 include/itkRadialDistortionTransform.hxx   |  23 +-
 src/itkIRCommon.cxx                        |  18 +-
 src/itkRBFTransform.cxx                    |  32 +-
 10 files changed, 857 insertions(+), 850 deletions(-)

diff --git a/include/itkGridTransform.h b/include/itkGridTransform.h
index 05ce386..96421c8 100644
--- a/include/itkGridTransform.h
+++ b/include/itkGridTransform.h
@@ -65,7 +65,7 @@ class GridTransform : public Transform<double, 2, 2>
   typedef Transform<double, 2, 2> Superclass;
 
   // Base inverse transform type:
-  typedef Superclass::InverseTransformType   InverseTransformType;
+  typedef Superclass                         InverseTransformType;
   typedef SmartPointer<InverseTransformType> InverseTransformPointer;
 
   // RTTI:
@@ -233,8 +233,8 @@ class GridTransform : public Transform<double, 2, 2>
   }
 
   // virtual:
-  unsigned int
-  GetNumberOfParameters() const
+  Superclass::NumberOfParametersType
+  GetNumberOfParameters() const override
   {
     return 2 * transform_.grid_.mesh_.size();
   }
@@ -255,7 +255,6 @@ class GridTransform : public Transform<double, 2, 2>
     transform_ = transform;
     GetParameters();
     GetFixedParameters();
-    this->m_Jacobian.SetSize(2, GetNumberOfParameters());
     this->m_FixedParameters[0] = is_inverse ? 1.0 : 0.0;
   }
 
@@ -266,9 +265,8 @@ class GridTransform : public Transform<double, 2, 2>
     return this->m_FixedParameters[0] != 0.0;
   }
 
-  // virtual:
-  const JacobianType &
-  GetJacobian(const InputPointType & point) const
+  void
+  ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const override
   {
     // FIXME: 20061227 -- this function was written and not tested:
 
@@ -281,26 +279,23 @@ class GridTransform : public Transform<double, 2, 2>
 
     unsigned int idx[3];
     double       jac[12];
-    this->m_Jacobian.SetSize(2, GetNumberOfParameters());
-    this->m_Jacobian.Fill(0.0);
+    jacobian.SetSize(2, GetNumberOfParameters());
+    jacobian.Fill(0.0);
     if (transform_.jacobian(point, idx, jac))
     {
       for (unsigned int i = 0; i < 3; i++)
       {
         unsigned int addr = idx[i] * 2;
-        this->m_Jacobian(0, addr) = Su * jac[i * 2];
-        this->m_Jacobian(0, addr + 1) = Su * jac[i * 2 + 1];
-        this->m_Jacobian(1, addr) = Sv * jac[i * 2 + 6];
-        this->m_Jacobian(1, addr + 1) = Sv * jac[i * 2 + 7];
+        jacobian(0, addr) = Su * jac[i * 2];
+        jacobian(0, addr + 1) = Su * jac[i * 2 + 1];
+        jacobian(1, addr) = Sv * jac[i * 2 + 6];
+        jacobian(1, addr + 1) = Sv * jac[i * 2 + 7];
       }
     }
-
-    return this->m_Jacobian;
   }
 
 protected:
   GridTransform()
-    : Superclass(2, 0)
   {
     this->m_FixedParameters.SetSize(7);
 
diff --git a/include/itkIRCommon.h b/include/itkIRCommon.h
index e6de385..c083731 100644
--- a/include/itkIRCommon.h
+++ b/include/itkIRCommon.h
@@ -113,11 +113,6 @@ extern const char * g_Version;
 //
 typedef itk::Size<2>::SizeValueType image_size_value_t;
 
-//----------------------------------------------------------------
-// array2d
-//
-#define array2d(T) std::vector<std::vector<T>>
-
 //----------------------------------------------------------------
 // mask_t
 //
@@ -1415,7 +1410,7 @@ template <class T>
 std::vector<typename itk::ImageRegionConstIteratorWithIndex<T>>
 DivideROIAmongThreads(typename T::ConstPointer fi,
                       typename T::RegionType   roi,
-                      int                      num_Threads = boost::thread::hardware_concurrency())
+                      int                      num_Threads = std::thread::hardware_concurrency())
 {
   typedef typename itk::ImageRegionConstIteratorWithIndex<T> itex_t;
   typedef typename T::IndexType                              index_t;
@@ -1487,7 +1482,7 @@ template <class T>
 std::vector<typename itk::ImageRegionIteratorWithIndex<T>>
 DivideROIAmongThreads(typename T::Pointer    fi,
                       typename T::RegionType roi,
-                      int                    num_Threads = boost::thread::hardware_concurrency())
+                      int                    num_Threads = std::thread::hardware_concurrency())
 {
   typedef typename itk::ImageRegionIteratorWithIndex<T> itex_t;
   typedef typename T::IndexType                         index_t;
@@ -2003,7 +1998,7 @@ my_metric_mt(double & area,
              const mask_t *                 mi_mask,
 
              const TInterpolator * mi_interpolator,
-             int                   num_Threads = boost::thread::hardware_concurrency())
+             int                   num_Threads = std::thread::hardware_concurrency())
 {
   // return my_metric(area, fi, mi, fi_to_mi, fi_mask, mi_mask, mi_interpolator);
 
@@ -3699,7 +3694,7 @@ make_mosaic(const typename IMG::SpacingType &                       mosaic_sp,
             typename IMG::PixelType background = 255.0,
 
             bool      dont_allocate = false,
-            const int num_threads = boost::thread::hardware_concurrency())
+            const int num_threads = std::thread::hardware_concurrency())
 {
   // NOTE: for backwards compatibility we do not assemble the mosaic mask:
   const bool      assemble_mosaic_mask = false;
@@ -4626,7 +4621,7 @@ make_mask(const mask_t::SpacingType &                             mosaic_sp,
           const std::vector<bool> &                               omit,
           const std::vector<typename transform_t::ConstPointer> & transform,
           const std::vector<mask_t::ConstPointer> &               image_mask,
-          unsigned int                                            num_threads = boost::thread::hardware_concurrency())
+          unsigned int                                            num_threads = std::thread::hardware_concurrency())
 {
   return make_mask_mt<transform_t>(num_threads, mosaic_sp, num_images, omit, transform, image_mask);
 }
@@ -4640,7 +4635,7 @@ template <class transform_t>
 mask_t::Pointer
 make_mask(const std::vector<typename transform_t::ConstPointer> & transform,
           const std::vector<typename mask_t::ConstPointer> &      image_mask,
-          unsigned int                                            num_threads = boost::thread::hardware_concurrency())
+          unsigned int                                            num_threads = std::thread::hardware_concurrency())
 {
   const unsigned int num_images = transform.size();
   return make_mask<transform_t>(
@@ -5159,8 +5154,8 @@ save_rgb(const char *             fn_save,
 inline static translate_transform_t::Pointer
 inverse(const translate_transform_t * transform)
 {
-  translate_transform_t::Pointer inv = NULL;
-  if (transform != NULL)
+  translate_transform_t::Pointer inv;
+  if (transform != nullptr)
   {
     inv = translate_transform_t::New();
     inv->SetOffset(neg(transform->GetOffset()));
@@ -6300,8 +6295,9 @@ approx_transform(const pnt2d_t &          tile_min,
                  const unsigned int       version = 1,
                  const bool               refine = false)
 {
-  base_transform_t::Pointer inverse = forward->GetInverse();
-  assert(inverse.GetPointer() != NULL);
+  base_transform_t::Pointer inverse;
+  forward->GetInverse(inverse);
+  assert(inverse.GetPointer() != nullptr);
 
   // calculate the shift:
   pnt2d_t center =
@@ -6434,8 +6430,9 @@ solve_for_transform(const pnt2d_t &                          tile_min,
                     const unsigned int                       version = 1,
                     const bool                               refine = false)
 {
-  base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile_0->GetInverse();
-  assert(tile_to_mosaic.GetPointer() != NULL);
+  base_transform_t::Pointer tile_to_mosaic;
+  mosaic_to_tile_0->GetInverse(tile_to_mosaic);
+  assert(tile_to_mosaic.GetPointer() != nullptr);
 
   // calculate the shift:
   pnt2d_t center =
@@ -6467,7 +6464,8 @@ solve_for_transform(const T *                                tile,
                     const unsigned int                       version = 1,
                     const bool                               refine = false)
 {
-  base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile_0->GetInverse();
+  base_transform_t::Pointer tile_to_mosaic;
+  mosaic_to_tile_0->GetInverse(tile_to_mosaic);
   assert(tile_to_mosaic.GetPointer() != NULL);
 
   typename T::SizeType    sz = tile->GetLargestPossibleRegion().GetSize();
diff --git a/include/itkLegendrePolynomialTransform.h b/include/itkLegendrePolynomialTransform.h
index 5dc98e0..65eb639 100644
--- a/include/itkLegendrePolynomialTransform.h
+++ b/include/itkLegendrePolynomialTransform.h
@@ -48,361 +48,388 @@
 
 //----------------------------------------------------------------
 // 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>
+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 Superclass                         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)
   {
-  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, // mosaic
-			      const std::vector<OutputPointType> & xy,// tile
-			      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)
+    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:
+  IdentifierType
+  GetNumberOfParameters() const override
+  {
+    return ParameterVectorLength;
+  }
+
+  void
+  ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const override;
+
+  // 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)
     {
-      mask.assign(ParameterVectorLength, shared);
-      mask[index_a(0, 0)] = false;
-      mask[index_b(0, 0)] = false;
+      xmax_ = Xmax;
+      ymax_ = Ymax;
     }
-    
-    // Convert the j, k indeces 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
+    else
     {
-      std::string base = Superclass::GetTransformTypeAsString();
-      std::ostringstream name;
-      name << base << '_' << N;
-      return name.str();
+      const double w = x_max - x_min;
+      const double h = y_max - y_min;
+
+      // -1 : 1
+      xmax_ = w / 2.0;
+      ymax_ = h / 2.0;
     }
-    
-  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
-  
+
+    // 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, // mosaic
+                       const std::vector<OutputPointType> & xy, // tile
+                       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 indeces 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();
+    std::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.hxx"
+#  include "itkLegendrePolynomialTransform.hxx"
 #endif
 
 
 //----------------------------------------------------------------
 // setup_transform
-// 
+//
 template <class transform_t>
 typename transform_t::Pointer
-setup_transform(const itk::Point<double, 2> & bbox_min,
-		const itk::Point<double, 2> & bbox_max)
+setup_transform(const itk::Point<double, 2> & bbox_min, const itk::Point<double, 2> & bbox_max)
 {
   double w = bbox_max[0] - bbox_min[0];
   double h = bbox_max[1] - bbox_min[1];
   double Umax = w / 2.0;
   double Vmax = h / 2.0;
-  
+
   typename transform_t::Pointer t = transform_t::New();
-  t->setup(bbox_min[0],
-	   bbox_max[0],
-	   bbox_min[1],
-	   bbox_max[1],
-	   Umax,
-	   Vmax);
-  
+  t->setup(bbox_min[0], bbox_max[0], bbox_min[1], bbox_max[1], Umax, Vmax);
+
   return t;
 }
 
 
 //----------------------------------------------------------------
 // setup_transform
-// 
+//
 template <class transform_t, class image_t>
 typename transform_t::Pointer
 setup_transform(const image_t * image)
 {
   typedef typename image_t::IndexType index_t;
-  
+
   index_t i00;
   i00[0] = 0;
   i00[1] = 0;
-  
+
   typename image_t::PointType origin;
   image->TransformIndexToPhysicalPoint(i00, origin);
-  
+
   index_t i11;
   i11[0] = 1;
   i11[1] = 1;
-  
+
   typename image_t::PointType spacing;
   image->TransformIndexToPhysicalPoint(i11, spacing);
   spacing[0] -= origin[0];
   spacing[1] -= origin[1];
-  
-  typename image_t::SizeType sz =
-    image->GetLargestPossibleRegion().GetSize();
-  
+
+  typename image_t::SizeType sz = image->GetLargestPossibleRegion().GetSize();
+
   typename image_t::PointType bbox_min = origin;
   typename image_t::PointType bbox_max;
   bbox_max[0] = bbox_min[0] + spacing[0] * double(sz[0]);
   bbox_max[1] = bbox_min[1] + spacing[1] * double(sz[1]);
-  
+
   return setup_transform<transform_t>(bbox_min, bbox_max);
 }
 
diff --git a/include/itkLegendrePolynomialTransform.hxx b/include/itkLegendrePolynomialTransform.hxx
index 8fe8946..00df6fd 100644
--- a/include/itkLegendrePolynomialTransform.hxx
+++ b/include/itkLegendrePolynomialTransform.hxx
@@ -43,504 +43,508 @@
 
 namespace itk
 {
-  namespace Legendre
+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 };
+} // namespace Legendre
+
+//----------------------------------------------------------------
+// LegendrePolynomialTransform
+//
+template <class TScalarType, unsigned int N>
+LegendrePolynomialTransform<TScalarType, N>::LegendrePolynomialTransform()
+{
+  // 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++)
   {
-    // 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 };
+    this->m_Parameters[i] = 0.0;
   }
-  
-  //----------------------------------------------------------------
-  // LegendrePolynomialTransform
-  // 
-  template <class TScalarType, unsigned int N>
-  LegendrePolynomialTransform<TScalarType, N>::
-  LegendrePolynomialTransform():
-    Superclass(2, ParameterVectorLength)
+
+  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();
+}
+
+//----------------------------------------------------------------
+// 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++)
   {
-    // 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;
-    }
+    P[i] = Legendre::P[i](A);
+    Q[i] = Legendre::P[i](B);
   }
-  
-  //----------------------------------------------------------------
-  // TransformPoint
-  // 
-  template <class TScalarType, unsigned int N>
-  typename LegendrePolynomialTransform<TScalarType, N>::OutputPointType
-  LegendrePolynomialTransform<TScalarType, N>::
-  TransformPoint(const InputPointType & x) const
+
+  double Sa = 0.0;
+  double Sb = 0.0;
+
+  for (unsigned int i = 0; i <= N; i++)
   {
-    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++)
+    for (unsigned int j = 0; j <= i; j++)
     {
-      P[i] = Legendre::P[i](A);
-      Q[i] = Legendre::P[i](B);
+      unsigned int k = i - j;
+      double       PjQk = P[j] * Q[k];
+      Sa += a(j, k) * PjQk;
+      Sb += b(j, k) * PjQk;
     }
-    
-    double Sa = 0.0;
-    double Sb = 0.0;
-    
-    for (unsigned int i = 0; i <= N; i++)
+  }
+
+  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)
+  {
+    itk::ExceptionObject e(__FILE__, __LINE__);
+    e.SetDescription("could not perform a numeric inverse transformation "
+                     "for a Legendre polynomial transform");
+    throw e;
+  }
+
+  OutputPointType x;
+  x[0] = vx[0];
+  x[1] = vx[1];
+
+  return x;
+}
+
+template <class TScalarType, unsigned int N>
+void
+LegendrePolynomialTransform<TScalarType, N>::ComputeJacobianWithRespectToParameters(const InputPointType & x,
+                                                                                    JacobianType & jacobian) 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);
+  }
+
+  // zero-out the Jacobian:
+  for (unsigned int i = 0; i < ParameterVectorLength; i++)
+  {
+    jacobian(0, i) = 0.0;
+    jacobian(1, i) = 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++)
     {
-      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;
-      }
+      unsigned int k = i - j;
+
+      double PjQk = P[j] * Q[k];
+      jacobian(0, index_a(j, k)) = Xmax * PjQk;
+      jacobian(1, index_b(j, k)) = Ymax * 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
+}
+
+//----------------------------------------------------------------
+// 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++)
   {
-    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)
+    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++)
     {
-      itk::ExceptionObject e(__FILE__, __LINE__);
-      e.SetDescription("could not perform a numeric inverse transformation "
-		       "for a Legendre polynomial transform");
-      throw e;
+      unsigned int k = i - j;
+      double       PjQk = P[j] * Q[k];
+      Sa += a(j, k) * PjQk;
+      Sb += b(j, k) * PjQk;
     }
-    
-    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
+
+  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++)
   {
-    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++)
+    for (unsigned int j = 0; j <= i; j++)
     {
-      P[i] = Legendre::P[i](A);
-      Q[i] = Legendre::P[i](B);
+      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 a_jk, b_jk:
-    for (unsigned int i = 0; i <= N; i++)
+  }
+
+  // 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++)
     {
-      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;
-      }
+      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;
     }
-    
-    // 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
+
+  // 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 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 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;
-    
-    double P[N + 1];
-    double Q[N + 1];
-    for (unsigned int i = 0; i <= N; i++)
+
+    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);
     }
-    
-    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 i = 0; i < start_with_degree; 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;
+        const unsigned int k = i - j;
+        double             PjQk = P[j] * Q[k];
+        bx[row] -= a(j, k) * PjQk;
+        by[row] -= b(j, k) * PjQk;
       }
     }
-    
-    // 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 i = start_with_degree; i < start_with_degree + degrees_covered; 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;
+        const unsigned int k = i - j;
+        const unsigned int col = index_a(j, k) - offset;
+
+        M[row][col] = P[j] * Q[k];
       }
     }
-    
-    // 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
+}
+
+//----------------------------------------------------------------
+// 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);
+
+  // use SVD to find the inverse of M:
+  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++)
   {
-    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++)
+    for (unsigned int j = 0; j <= i; j++)
     {
-      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];
-	}
-      }
+      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];
     }
   }
-  
-  //----------------------------------------------------------------
-  // 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
+}
+
+//----------------------------------------------------------------
+// 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++)
   {
-    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);
-    
-    // use SVD to find the inverse of M:
-    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++)
     {
-      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];
-      }
+      unsigned int k = i - j;
+      os << indent << "a(" << j << ", " << k << ") = " << a(j, k) << std::endl;
     }
   }
-  
-  //----------------------------------------------------------------
-  // PrintSelf
-  // 
-  template <class TScalarType, unsigned int N>
-  void
-  LegendrePolynomialTransform<TScalarType, N>::
-  PrintSelf(std::ostream & os, Indent indent) const
+  for (unsigned int i = 0; i <= N; i++)
   {
-    Superclass::PrintSelf(os,indent);
-    
-    for (unsigned int i = 0; i <= N; i++)
+    for (unsigned int j = 0; j <= i; j++)
     {
-      for (unsigned int j = 0; j <= i; j++)
-      {
-	unsigned int k = i - j;
-	os << indent << "a(" << j << ", " << k << ") = " << a(j, k)
-	   << std::endl;
-      }
+      unsigned int k = i - j;
+      os << indent << "b(" << j << ", " << k << ") = " << b(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;
+  }
+  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++)
     {
@@ -548,8 +552,8 @@ namespace itk
 	 << std::endl;
     }
 #endif
-  }
-  
+}
+
 } // namespace itk
 
 
diff --git a/include/itkMeshTransform.h b/include/itkMeshTransform.h
index 26feee2..94cb905 100644
--- a/include/itkMeshTransform.h
+++ b/include/itkMeshTransform.h
@@ -65,7 +65,7 @@ class MeshTransform : public Transform<double, 2, 2>
   typedef Transform<double, 2, 2> Superclass;
 
   // Base inverse transform type:
-  typedef Superclass::InverseTransformType   InverseTransformType;
+  typedef Superclass                         InverseTransformType;
   typedef SmartPointer<InverseTransformType> InverseTransformPointer;
 
   // RTTI:
@@ -248,9 +248,8 @@ class MeshTransform : public Transform<double, 2, 2>
     return this->m_Parameters;
   }
 
-  // virtual:
-  unsigned int
-  GetNumberOfParameters() const
+  IdentifierType
+  GetNumberOfParameters() const override
   {
     return 4 * transform_.grid_.mesh_.size();
   }
@@ -271,7 +270,6 @@ class MeshTransform : public Transform<double, 2, 2>
     transform_ = transform;
     GetParameters();
     GetFixedParameters();
-    this->m_Jacobian.SetSize(2, GetNumberOfParameters());
     this->m_FixedParameters[0] = is_inverse ? 1.0 : 0.0;
   }
 
@@ -282,9 +280,8 @@ class MeshTransform : public Transform<double, 2, 2>
     return this->m_FixedParameters[0] != 0.0;
   }
 
-  // virtual:
-  const JacobianType &
-  GetJacobian(const InputPointType & point) const
+  void
+  ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const override
   {
     // FIXME: 20061227 -- this function was written and not tested:
 
@@ -297,26 +294,23 @@ class MeshTransform : public Transform<double, 2, 2>
 
     unsigned int idx[3];
     double       jac[12];
-    this->m_Jacobian.SetSize(2, GetNumberOfParameters());
-    this->m_Jacobian.Fill(0.0);
+    jacobian.SetSize(2, GetNumberOfParameters());
+    jacobian.Fill(0.0);
     if (transform_.jacobian(point, idx, jac))
     {
       for (unsigned int i = 0; i < 3; i++)
       {
         unsigned int addr = idx[i] * 2;
-        this->m_Jacobian(0, addr) = Su * jac[i * 2];
-        this->m_Jacobian(0, addr + 1) = Su * jac[i * 2 + 1];
-        this->m_Jacobian(1, addr) = Sv * jac[i * 2 + 6];
-        this->m_Jacobian(1, addr + 1) = Sv * jac[i * 2 + 7];
+        jacobian(0, addr) = Su * jac[i * 2];
+        jacobian(0, addr + 1) = Su * jac[i * 2 + 1];
+        jacobian(1, addr) = Sv * jac[i * 2 + 6];
+        jacobian(1, addr + 1) = Sv * jac[i * 2 + 7];
       }
     }
-
-    return this->m_Jacobian;
   }
 
 protected:
   MeshTransform()
-    : Superclass(2, 0)
   {
     this->m_FixedParameters.SetSize(8);
 
diff --git a/include/itkRBFTransform.h b/include/itkRBFTransform.h
index 4b60199..36c0c43 100644
--- a/include/itkRBFTransform.h
+++ b/include/itkRBFTransform.h
@@ -80,7 +80,7 @@ class RBFTransform : public Transform<double, 2, 2>
   typedef Transform<double, 2, 2> Superclass;
 
   // Base inverse transform type:
-  typedef Superclass::InverseTransformType   InverseTransformType;
+  typedef Superclass                         InverseTransformType;
   typedef SmartPointer<InverseTransformType> InverseTransformPointer;
 
   // RTTI:
@@ -141,9 +141,8 @@ class RBFTransform : public Transform<double, 2, 2>
     return this->m_Parameters;
   }
 
-  // virtual:
-  unsigned int
-  GetNumberOfParameters() const
+  IdentifierType
+  GetNumberOfParameters() const override
   {
     return this->m_Parameters.size();
   }
@@ -154,8 +153,8 @@ class RBFTransform : public Transform<double, 2, 2>
 
   // setup the transform parameters:
   void
-  setup( // image bounding box expressed in the image space,
-         // defines transform normalization parameters:
+  setup(                              // image bounding box expressed in the image space,
+                                      // defines transform normalization parameters:
     const OutputPointType & tile_min, // tile space
     const OutputPointType & tile_max, // tile space
 
@@ -164,9 +163,8 @@ class RBFTransform : public Transform<double, 2, 2>
     const InputPointType *  uv,      // mosaic space
     const OutputPointType * xy);     // tile space
 
-  // virtual:
-  const JacobianType &
-  GetJacobian(const InputPointType & point) const;
+  void
+  ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const override;
 
 #if 0
     // helper required for numeric inverse transform calculation;
diff --git a/include/itkRadialDistortionTransform.h b/include/itkRadialDistortionTransform.h
index dde7fbe..4145f7c 100644
--- a/include/itkRadialDistortionTransform.h
+++ b/include/itkRadialDistortionTransform.h
@@ -77,8 +77,8 @@ class RadialDistortionTransform : public Transform<TScalar, 2, 2>
   typedef Transform<TScalar, 2, 2> Superclass;
 
   // Base inverse transform type:
-  typedef typename Superclass::InverseTransformType InverseTransformType;
-  typedef SmartPointer<InverseTransformType>        InverseTransformPointer;
+  typedef Superclass                         InverseTransformType;
+  typedef SmartPointer<InverseTransformType> InverseTransformPointer;
 
   // static constant for the number of polynomial coefficients:
   itkStaticConstMacro(Nk, unsigned int, N);
@@ -137,16 +137,14 @@ class RadialDistortionTransform : public Transform<TScalar, 2, 2>
     return this->m_Parameters;
   }
 
-  // virtual: mumber of parameters that define this transform:
-  unsigned int
-  GetNumberOfParameters() const
+  IdentifierType
+  GetNumberOfParameters() const override
   {
     return N + 2;
   }
 
-  // virtual:
-  const JacobianType &
-  GetJacobian(const InputPointType & point) const;
+  void
+  ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const override;
 
   // virtual: return an inverse of this transform.
   InverseTransformPointer
diff --git a/include/itkRadialDistortionTransform.hxx b/include/itkRadialDistortionTransform.hxx
index 5784f7c..67bf53f 100644
--- a/include/itkRadialDistortionTransform.hxx
+++ b/include/itkRadialDistortionTransform.hxx
@@ -44,7 +44,6 @@ namespace itk
 //
 template <class TScalar, unsigned int N>
 RadialDistortionTransform<TScalar, N>::RadialDistortionTransform()
-  : Superclass(2, N + 2) // there are k0, k1,.. kN-1, tx, ty parameters:
 {
   // by default, the parameters are initialized for an identity transform:
   this->m_Parameters[0] = 1.0;
@@ -139,12 +138,10 @@ RadialDistortionTransform<TScalar, N>::BackTransformPoint(const OutputPointType
   return x;
 }
 
-//----------------------------------------------------------------
-// GetJacobian
-//
 template <class TScalar, unsigned int N>
-const typename RadialDistortionTransform<TScalar, N>::JacobianType &
-RadialDistortionTransform<TScalar, N>::GetJacobian(const InputPointType & x) const
+void
+RadialDistortionTransform<TScalar, N>::ComputeJacobianWithRespectToParameters(const InputPointType & x,
+                                                                              JacobianType &         jacobian) const
 {
   const double & ac = this->m_FixedParameters[0];
   const double & bc = this->m_FixedParameters[1];
@@ -168,8 +165,8 @@ RadialDistortionTransform<TScalar, N>::GetJacobian(const InputPointType & x) con
   double RRmax2n = double(1);
   for (unsigned int n = 0; n < N; n++)
   {
-    this->m_Jacobian(0, n) = A * RRmax2n;
-    this->m_Jacobian(1, n) = B * RRmax2n;
+    jacobian(0, n) = A * RRmax2n;
+    jacobian(1, n) = B * RRmax2n;
     RRmax2n *= RRmax2;
   }
 
@@ -192,13 +189,11 @@ RadialDistortionTransform<TScalar, N>::GetJacobian(const InputPointType & x) con
 
   double S = this->m_Parameters[0] + RRmax2 * P;
 
-  this->m_Jacobian(0, N) = Rmax * S + (double(2) * A2 / Rmax) * Q;     // dx/dta
-  this->m_Jacobian(1, N + 1) = Rmax * S + (double(2) * B2 / Rmax) * Q; // dy/dtb
-
-  this->m_Jacobian(0, N + 1) = ((double(2) * A * B) / Rmax) * Q; // dx/dtb
-  this->m_Jacobian(1, N) = this->m_Jacobian(0, N + 1);           // dy/dta
+  jacobian(0, N) = Rmax * S + (double(2) * A2 / Rmax) * Q;     // dx/dta
+  jacobian(1, N + 1) = Rmax * S + (double(2) * B2 / Rmax) * Q; // dy/dtb
 
-  return this->m_Jacobian;
+  jacobian(0, N + 1) = ((double(2) * A * B) / Rmax) * Q; // dx/dtb
+  jacobian(1, N) = jacobian(0, N + 1);                   // dy/dta
 }
 
 //----------------------------------------------------------------
diff --git a/src/itkIRCommon.cxx b/src/itkIRCommon.cxx
index 7d9c335..c6d8bb7 100644
--- a/src/itkIRCommon.cxx
+++ b/src/itkIRCommon.cxx
@@ -579,8 +579,9 @@ calc_tile_mosaic_bbox(const base_transform_t * mosaic_to_tile,
     return true;
   }
 
-  base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse();
-  if (tile_to_mosaic.GetPointer() == NULL)
+  base_transform_t::Pointer tile_to_mosaic;
+  mosaic_to_tile->GetInverse(tile_to_mosaic);
+  if (tile_to_mosaic.GetPointer() == nullptr)
   {
     return false;
   }
@@ -681,7 +682,8 @@ find_inverse(const pnt2d_t &          tile_min, // tile space
              const double             min_error_sqrd,
              const unsigned int       pick_up_pace_steps)
 {
-  base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse();
+  base_transform_t::Pointer tile_to_mosaic;
+  mosaic_to_tile->GetInverse(tile_to_mosaic);
   return find_inverse(tile_min,
                       tile_max,
                       mosaic_to_tile,
@@ -980,8 +982,9 @@ generate_landmarks_v1(const pnt2d_t &          tile_min,
 {
   // #define DEBUG_LANDMARKS
 
-  base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse();
-  if (tile_to_mosaic.GetPointer() == NULL)
+  base_transform_t::Pointer tile_to_mosaic;
+  mosaic_to_tile->GetInverse(tile_to_mosaic);
+  if (tile_to_mosaic.GetPointer() == nullptr)
   {
     return false;
   }
@@ -1185,8 +1188,9 @@ generate_landmarks_v2(const pnt2d_t &          tile_min,
 {
   // #define DEBUG_LANDMARKS
 
-  base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse();
-  if (tile_to_mosaic.GetPointer() == NULL)
+  base_transform_t::Pointer tile_to_mosaic;
+  mosaic_to_tile->GetInverse(tile_to_mosaic);
+  if (tile_to_mosaic.GetPointer() == nullptr)
   {
     return false;
   }
diff --git a/src/itkRBFTransform.cxx b/src/itkRBFTransform.cxx
index 33e3c98..3cbc5a9 100644
--- a/src/itkRBFTransform.cxx
+++ b/src/itkRBFTransform.cxx
@@ -50,12 +50,9 @@ namespace itk
 // RBFTransform::RBFTransform
 //
 RBFTransform::RBFTransform()
-  : Superclass(2, 0)
 {
   this->m_FixedParameters.SetSize(4);
   this->m_Parameters.SetSize(6);
-  this->m_Jacobian.SetSize(2, 6);
-  this->m_Jacobian.Fill(0);
 
   double & Xmax = this->m_FixedParameters[0];
   double & Ymax = this->m_FixedParameters[1];
@@ -160,8 +157,8 @@ RBFTransform::GetInverse() const
 // RBFTransform::setup
 //
 void
-RBFTransform::setup( // image bounding box expressed in the image space,
-                     // defines transform normalization parameters:
+RBFTransform::setup(                // image bounding box expressed in the image space,
+                                    // defines transform normalization parameters:
   const OutputPointType & tile_min, // tile space
   const OutputPointType & tile_max, // tile space
 
@@ -172,8 +169,6 @@ RBFTransform::setup( // image bounding box expressed in the image space,
 {
   this->m_FixedParameters.SetSize(4 + num_pts * 2);
   this->m_Parameters.SetSize(6 + num_pts * 2);
-  this->m_Jacobian.SetSize(2, 6 + num_pts * 2);
-  this->m_Jacobian.Fill(0);
 
   // setup the normalization parameters:
   double & Xmax = this->m_FixedParameters[0];
@@ -277,8 +272,9 @@ RBFTransform::setup( // image bounding box expressed in the image space,
 //----------------------------------------------------------------
 // RBFTransform::GetJacobian
 //
-const RBFTransform::JacobianType &
-RBFTransform::GetJacobian(const InputPointType & point) const
+void
+RBFTransform::ComputeJacobianWithRespectToParameters(const InputPointType & point,
+                                                     JacobianType &         jacobian) const override
 {
   // shortcuts:
   const double & Xmax = GetXmax();
@@ -292,13 +288,13 @@ RBFTransform::GetJacobian(const InputPointType & point) const
   const double A = (u - uc) / Xmax;
   const double B = (v - vc) / Ymax;
 
-  this->m_Jacobian(0, index_a(0)) = Xmax;
-  this->m_Jacobian(0, index_a(1)) = Xmax * A;
-  this->m_Jacobian(0, index_a(2)) = Xmax * B;
+  jacobian(0, index_a(0)) = Xmax;
+  jacobian(0, index_a(1)) = Xmax * A;
+  jacobian(0, index_a(2)) = Xmax * B;
 
-  this->m_Jacobian(1, index_b(0)) = Ymax;
-  this->m_Jacobian(1, index_b(1)) = Ymax * A;
-  this->m_Jacobian(1, index_b(2)) = Ymax * B;
+  jacobian(1, index_b(0)) = Ymax;
+  jacobian(1, index_b(1)) = Ymax * A;
+  jacobian(1, index_b(2)) = Ymax * B;
 
   const unsigned int num_pts = num_points();
   const double *     uv_vec = &(this->m_FixedParameters[index_uv(0)]);
@@ -313,11 +309,9 @@ RBFTransform::GetJacobian(const InputPointType & point) const
     double lR = R == 0 ? 0 : log(R);
     double Q = R * lR;
 
-    this->m_Jacobian(0, index_f(i)) = Xmax * Q;
-    this->m_Jacobian(1, index_g(i)) = Ymax * Q;
+    jacobian(0, index_f(i)) = Xmax * Q;
+    jacobian(1, index_g(i)) = Ymax * Q;
   }
-
-  return this->m_Jacobian;
 }
 
 #if 0