ImageAlignment.hh

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003, 2004    (see file CONTACTS for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#ifndef __ImageAlignment__hh__
#define __ImageAlignment__hh__

#include <Base/Common/BIASpragmaStart.hh>

#include <Base/Debug/Debug.hh>
#include <Base/Geometry/HomgPoint2D.hh>
#include <Base/Math/Matrix2x2.hh>
#include <Image/PyramidImage.hh>
#include <Image/TextureTransform.hh>
#include <vector>
#include <Geometry/LocalAffineFrame.hh>


#define D_IMAGEALIGNMENT_PROGRESS   1
#define D_IMAGEALIGNMENT_PARAMETER  2
#define D_IMAGEALIGNMENT_PERPIXEL   4
#define D_IMAGEALIGNMENT_IMAGES     8  

namespace BIAS {


  /** @class ImageAlignment
   *  @ingroup g_tracker
   *  @brief  Inverse Compositional Image Alignment ("Registration")
   *
   * This class implements generic inverse compositional image alignment as 
   * proposed by Baker and Matthews in "Lucas-Kanade 20 years on". The model 
   * for which you estimate the parameters must be of type TextureTransform 
   * (e.g. affine, homography, displacement, camera pose, calibration, ...).
   *
   * The key idea is that the required inversion of the hessian and the 
   * incorporation of the gradients are not recomputed per iteration in this
   * formulation as opposed to the standard KLT.
   *
   *  @author koeser 01/2008 */
  class BIASMatcher2D_EXPORT ImageAlignment : public Debug
  {
  public :
    ImageAlignment();
    ~ImageAlignment() { if (T_!=NULL) delete T_;};

    /** @brief aligns image1 with image2 with respect to parameters params,
     *         scale space tracking is used recursively until maximum 
     *         resolution reached. Requires correct Cov as input.
     */
    int AutoAlign(const PyramidImage<float>& I1, const PyramidImage<float>& I2, 
                  Vector<double>& params, Matrix<double>& Cov);

    /** @brief aligns image1 with image2 with respect to parameters params,
     *         traditional pyramid based alignment is used, mainly for 
     *         comparison reasons. Can be used also if Cov is totally unknown.
     */
    int StrictPyramidAlign(const PyramidImage<float>& I1, 
                           const PyramidImage<float>& I2, 
                           Vector<double>& params, Matrix<double>& Cov);

    /** @brief aligns image1 with image2 with respect to parameters params
     *
     *  if no pyramidlevel is supplied the best pyramid level for each
     *  pixel is chosen automatically using Cov (scale space tracking) */
    int Align(const PyramidImage<float>& I1, const PyramidImage<float>& I2, 
              Vector<double>& params, Matrix<double>& Cov, int PyramidLevel=-1);

    /** @brief pass object type which holds texture transform */
    void SetTextureTransform(const TextureTransform& T) {
      if (T_!=NULL) delete T_;
      T_ = T.Clone();
    }

    /** @brief number of iterations, before the system stops unsuccessfully */
    inline void SetMaxIterations(const int maxiter)        
    { MaxIterations_ = maxiter; }

    /** @brief if a step smaller (norm L2) than this is taken, convergence
        is declared and the system finishes */
    inline void SetUpdateThreshold(const double maxerr)         
    { MaxUpdate_ = maxerr; }

    /** enable brightness variance and offset invariant computation */
    void SetAffineBrightnessInvariance(bool bi) { 
      AffineBrightnessInvariance_ = bi;
    }

    inline int GetMaxIterations() { return MaxIterations_; }

    inline float GetMaxUpdate() { return MaxUpdate_; }
   
    /** @brief choose any rectangular tracking region and resolution
     *
     *  The localaffineframe determines the shape, scale and orientation of the
     *  rectangular region used for tracking. The numberofpixelsis a 
     *  windowsize, e.g. if you pass 9, a 9x9 grid will be used for tracking
     *  If you pass -1, each pixel is used for tracking */
    void SetAlignmentPixelsRectangular(const LocalAffineFrame& LAF,
                                       int numberOfPixels = -1);

    /** @brief return relative affine transform between i1 and i2 if
     *         AffineBrightnessInvariance was active 
     *
     *  i2 = scale * i1 + shift
     */
    void GetLastBrightnessChange(double& scale, double& shift) {
      scale = dev2_ / dev1_;
      shift = mean2_ - scale * mean1_;
    }

    /** @brief enable/disable and set step width for line search: <=0 disables,
     *         0.1 or 0.5 may be good values. 1.0 would mean original step */
    void SetLineSearch(const double& lineSearch) {
      LineSearch_ = lineSearch;
      BIASASSERT(LineSearch_<1.0);
    }
    /** @brief adds a dampening d to hessian's diag
        @param dampening: >=0: d, <0 -d*largest eigvalue */
    void SetDampening(const double& dampening) {
      Dampening_ = dampening;
    }
    /** @brief set expected standard deviation for grey value noise 
        
        This is important if maximum a posteriori-estimation is desired
        and you provide a correctly scaled prior cov for the warp parameters.
        Pixels which have smaller gradient than this are rejected.
     */
    void SetIntensitySigma(const double& sigma) {
      intensitySigma_ = sigma;
    }

  protected:

    /** @brief check for given cov and params, whether control points move more
     *         than pixelAccuracy */
    int CheckConvergence_(const Vector<double>& params, 
                          const Matrix<double>& Cov, 
                          const double& pixelAccuracy = 0.3);

    ///////////////////////////////////////////////////////////////////////////
    //  params
    ///////////////////////////////////////////////////////////////////////////
    
    /// iteration stops after
    int MaxIterations_;

    /// iteration stops if parameter update is less than
    float MaxUpdate_;     

    /// adds a dampening to hessian's diag: >=0: d, <0 d*largest eigvalue
    double Dampening_;  

    /// compute brightness "invariant"
    bool AffineBrightnessInvariance_;

    /// if error is not reduced in one step, linesearch times this step is tried
    /// set to <=0 to disable line search
    double LineSearch_;
    
    /// how much noise is expected on intensity values
    double intensitySigma_;
    ///////////////////////////////////////////////////////////////////////////

    /// the warp model
    TextureTransform *T_;

    /// mean and std deviation of image patch
    double mean1_, mean2_, dev1_, dev2_;

    /// gradients and grey values from image 1
    std::vector<double> TemplateGreyValues_;
    std::vector<Vector2<double> > TemplateGradients_;
    std::vector<double> TemplateScales_;
    std::vector<HomgPoint2D> TemplatePositions_;
    std::vector<HomgPoint2D> ControlPositions_;

    std::vector<HomgPoint2D> ImagePositions_;
    std::vector<Vector<double> > SteepestDescentImages_;
    std::vector<bool> active_;

    /// internal counter for numbering debug images
    int iterationDebug_;
    
  }; // class

} // namespace

#endif