simplepolymapping.h

#ifndef SIMPLEPOLYMAPPING__H
#define SIMPLEPOLYMAPPING__H

#include <memory> // for auto_ptr

#include "mapping.h"
#include "gtransfo.h"
#include "ccdimage.h"


class SimpleGtransfoMapping : public Mapping
{
 protected:
  bool toFit; 
  unsigned index;
  /* inheritance may also work. Perhaps with some trouble because
     some routines in Mapping and Gtransfo have the same name */
  CountedRef<Gtransfo> transfo; 
  /* to avoid allocation at every call of PosDerivatives. 
     use a pointer for constness */
  std::auto_ptr<GtransfoLin> lin; 

#ifdef STORAGE
  static_assert(std::is_base_of<Gtransfo, Tr>::value,
        "SimpleGtransfoMapping requires a class deriving from Gtransfo");
#endif


 public :

 SimpleGtransfoMapping(const Gtransfo &T, bool ToFit=true) : toFit(ToFit)
    {
      transfo = T.Clone(); 
      lin.reset(new GtransfoLin); // reserve space for derivative computation
    }

  // interface Mapping functions:

  unsigned Npar() const 
  {
    if (toFit) return transfo->Npar();
    else return 0;
  }

  void GetMappingIndices(std::vector<unsigned> &Indices) const
  {
    if (Indices.size() < Npar()) Indices.resize(Npar());
    for (unsigned k=0; k<Npar(); ++k) Indices[k] = index+k;
  }

  void TransformPosAndErrors(const FatPoint &Where,
                 FatPoint &OutPos) const
  {
    transfo->TransformPosAndErrors(Where,OutPos);
  }

  void  PosDerivative(Point &Where, Eigen::Matrix2d &Der, const double & Eps) const
  {
    transfo->Derivative(Where, *lin, Eps);
    Der(0,0) = lin->Coeff(1,0,0);
    // 
    /* This does not work : it was proved by rotating the frame
       see the compilation switch ROTATE_T2 in constrainedpolymodel.cc
    Der(1,0) = lin->Coeff(1,0,1);
    Der(0,1) = lin->Coeff(0,1,0);
    */
    Der(1,0) = lin->Coeff(0,1,0);
    Der(0,1) = lin->Coeff(1,0,1);
    Der(1,1) = lin->Coeff(0,1,1);
  }

  void OffsetParams(const double *Delta)
  {
    transfo->OffsetParams(Delta);
  }

  unsigned Index() const { return index;}

  void  SetIndex(unsigned I) {index=I;}
  
  virtual void ComputeTransformAndDerivatives(const FatPoint &Where,
                          FatPoint &OutPos,
                          Eigen::MatrixX2d &H) const
  {
    TransformPosAndErrors(Where,OutPos);
    transfo->ParamDerivatives(Where, &H(0,0), &H(0,1));      
  }

  const Gtransfo&  Transfo() const {return *transfo;}

};


class SimplePolyMapping : public SimpleGtransfoMapping
{
  // to better condition the 2nd derivative matrix
  Point center;
  // could consider scaling as well.

 public:

  // ! contructor
 SimplePolyMapping(const Point& Center, const GtransfoPoly& P):
  SimpleGtransfoMapping(P), center(Center)
  {
    /* We work in centered coordinates: the routines below center
       coordinates so that the provided initialization has to be
       shifted the other way around: */
    GtransfoLinShift shift(center.x, center.y);
    transfo = new GtransfoPoly(P*shift);
    // check of matrix indexing (once for all)
    MatrixX2d H(3,2);
    assert((&H(1,0) - &H(0,0)) == 1);
  }

  /* We should put the computation of error propagation and
     parameter derivatives into the same Gtransfo routine because
     it could be significantly faster */
  virtual void ComputeTransformAndDerivatives(const FatPoint &Where,
                          FatPoint &OutPos,
                          Eigen::MatrixX2d &H) const
    {
      TransformPosAndErrors(Where,OutPos);
      transfo->ParamDerivatives(Where-center, &H(0,0), &H(0,1));      
    }

  void TransformPosAndErrors(const FatPoint &Where,
                 FatPoint &OutPos) const
  {
    FatPoint where = Where;
    where.x -= center.x;
    where.y -= center.y;
    transfo->TransformPosAndErrors(where,OutPos);
  }




};


#ifdef STORAGE

class SimpleIdentityMapping : public SimpleGtransfoMapping<GtransfoIdentity>
{
 public: 

  virtual void ComputeTransformAndDerivatives(const FatPoint &Where,
                          FatPoint &OutPos,
                          Eigen::MatrixX2d &H) const
  {
    OutPos = Where;
  }

};
#endif


#endif