BackwardMapping.hh

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

Copyright (C) 2003-2009    (see file CONTACT 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 __BackwardMapping_hh__
#define __BackwardMapping_hh__

#include <Base/Image/Image.hh>
#include <Image/PyramidImage.hh>
#include <Base/Geometry/HomgPoint2D.hh>
#include <Base/Math/Matrix2x2.hh>
#include <limits>
#include <Image/BackwardMappingLUTStructs.hh>

/// if jacobian gets smaller than this, clip to this value to avoid singularity
/// when computing smoothing kernel
#define JACOBI_THRESH 1e-100
// (1.0/log(2.0));
#define ONE_OVER_LOG2 1.442695

namespace BIAS {


  /** @class BackwardMapping
      @ingroup g_filter
      @brief  Abstract base class to map an image (texture) src to an image
      sink with an arbitrary continuous function (e.g. rotation, homography,
      undistortion, rectification, ...). Provides
      interpolation, border checking and anti-aliasing.

      If you want to learn more about texture handling, a good source is
      Paul Heckbert: "Fundamentals of Texture Mapping and Image Warping",  
      Master's thesis, UCB/CSD 89/516, CS Division, U.C. Berkeley, June 1989
      available at http://www.cs.cmu.edu/~ph/texfund/texfund.pdf

      This class implements a backward mapping and runs over all(*) pixels
      of some output image. For each pixel it calls GetSourceCoordinates
      (which you have to implement in your derived class) to compute
      the coordinates in the source image. The local sampling frequency
      can be computed to avoid aliasing. Either you make an analytic
      implementation or simply use the good (but maybe slow) numerical
      approximation from BackwardMapping.

      For instance if you scale an image by factor 2, you would implement
      GetSourceCoordinates as source_coord = 0.5 * sink_coord and the local
      sampling frequency would be 0.5. Sampling frequencies below 1 wont
      introduce aliasing, while those larger than 1 can.
      Therefore local low pass filtering is simulated by using a pyramid
      of the input image.

      The local sampling frequency of your mapping function can be approximated
      analyzing the jacobian of the function as a R^2->R^2 mapping, see
      HomographyMapping for an example or simply use the numerical solution in
      this class if you are not interested in speed.
      However, using an image pyramid, isotropic scaling is assumed, if you
      resample by 0.3 in x-direction and 5.0 in y direction, both dirs are
      smoothed with the same sigma and you will lose image information. If
      you want to remove such effects, use anisotropic filtering. This
      is the (by far) slowest but best technique. If you use highly nonlinear
      mapping functions, the quality of your result can improve if you
      compute a larger image and reduce size by subsequent subsamping in
      destination space. This technique is called super-sampling.

      To make a pyramid of correct size, Backwardmapping::UpdatePyramidSize
      maps the corners of your sink image and computes the local scale
      at these corners assuming maximum scale at one corner. If your
      derived mapping function has a maximum scale somewhere else (e.g.
      cylindircal rectification), you must override UpdatePyramidSize.

      If possible,bicubic interpolation is done in the image, at borders
      (high pyramid levels) only bilinear. For efficiency (**) or comparison,
      you can disable the usage of pyramids, speeding up the computation
      (bicubic, bilinear, nearest neighbor).

      (*) "all" means: getboundingbox forward-maps src-corners to sink with
      inverse mapping (has to be implmented in derived class !) and these
      corners are clipped against the ROI of sink.

      (**)The first goal of this source code is that it is generic and correct.
      Therefore, in the most inner loop the number of color channels is
      evaluated and weighting factors are computed several times. If you want
      to be fast, add specialized functions for grey, and exploit the
      structure of your problem
      (see e.g. AffineMapping which has a constant Jacobian, ..). You might as
      well consider using the GPU if your mapping is "fixed" (radial
      undistortion with always the same parameters) or use the lookup table
      interfaces if you can afford precomputing your coefficients (e.g. for a
      constant mapping applied hundereds of times)


      OPENMP - open multiprocessing support:
      In the most generic trilinear and anisotropic loop, openmp is enabled
      if it is enabled at all in BIAS (should be by default). So far it is
      not implmented for grey(!) and for bilinear/bicubic mappings (but this is
      only enumerating the right variables in 3 lines of code and some tests).


      @author koeser, 2005

  */
  template <class InputStorageType, class OutputStorageType>
  class BIASImage_EXPORT BackwardMapping {
  public:

    /**
     * @brief Constructor.
     */
    inline
    BackwardMapping()
    {
#ifdef BIAS_DEBUG
      rangeChecked_ = true;
#endif
      aliasing_ = false;
      incomplete_ = false;
      autoPyramidSize_ = true;
      pyramidSize_ = 4;
      interpolationType_ = MapTrilinear;
      offsetX_ = 0;
      offsetY_ = 0;
      pLUT_Nea_ = NULL;
      pLUT_Bil_ = NULL;
      pLUT_Tri_ = NULL;
      superSampling_ = 1.0;
      pConcatenation_ = NULL;
      pyramid_.SetRescaleFactor(2.0);
    }

    /**
     * @brief Destructor.
     */
    inline
    virtual ~BackwardMapping()
    {
      if (pLUT_Nea_)       delete pLUT_Nea_;
      if (pLUT_Bil_)       delete pLUT_Bil_;
      if (pLUT_Tri_)       delete pLUT_Tri_;
      if (pConcatenation_) delete pConcatenation_;
    }

    /** @name User interfaces
        @{
    */

    /** @brief backward mapping with various interpolations

    Sink must be initialized with correct size and background color for invalid
    pixels. This way you can call it several times to create a panorama from
    several source images. For each pixel in sink, GetSourceCoordinates is
    called and the value in source interpolated.
    @param newSink, if true a new sink image is created that is large enough
    to hold the old sink and the new mapped src img.
    @param SuperSampling if > 1.0, intermediate image with isotropically
    enlarged width and height (by factor SuperSampling) is created and
    downsampled for the final image (most useful for nonlinear mappings)
    @return 0 on success, >0 means notification(such as 1 = not all pixels
    could be computed), -1 means failure
    */
    int Map(const Image<InputStorageType>& src,
            Image<OutputStorageType>& sink, InterpolationMethod=MapTrilinear,
            bool newSink=false, double SuperSampling=1.0);

    /** @brief use this for subsequent calls to avoid pyramid recalculation
        when using tri-filtering */
    int Map(Image<OutputStorageType>& sink);

    /** @brief precomputes lookup coordinates and computes displacement map
        int TEXTURE coordinates, according to the size of src (or width,height)
    */
    int GetDisplacementMap(Image<float>& dismap,
                           int width,
                           int height);

    int GetDisplacementMap(Image<float>& dismap,
                           Image<InputStorageType>& src);

    /** @brief precomputes lookup coordinates for accessing src

    If you want to apply the same mapping function over and over again, e.g.
    radial undistortion for a series of images, call this function once and
    MapWithLookupTable every time. sink must be of correct size !
    @return 0 on success, >0 means notification(such as 1 = not all pixels
    could be computed),*/
    int PrepareLookupTableMapping(const Image<InputStorageType>& src,
                                  Image<OutputStorageType>& sink,
                                  InterpolationMethod method,
                                  bool newSink=false);

    /** @brief applies precomputed coordinates in src, fast for repeated
     usages of same mapping function */
    int MapWithLookupTable(const Image<InputStorageType>& src,
                           Image<OutputStorageType>& sink,
                           InterpolationMethod method);

    /** Returns the LUT generated by PrepareLookupTableMapping()
      * as an file.
      **/
    int GetLookupTable(const std::string& filename, InterpolationMethod method);
    int GetLookupTable(float* lutImageRow, InterpolationMethod method);
    /** returns the size of the current LUT.
    * \returns -1 if the current InterpolationMethod has no LUT!
    **/
    int GetLookupTableSize(unsigned int& width, unsigned int& channels, InterpolationMethod method);

    /** Loads the LUT from a file, replaces the call to
     * PrepareLookupTableMapping().
     * \returns the interpolation type in the parameter named <i>method</i>
     * stored in the passed file.
     **/
    int SetLookupTable(const std::string& filename, InterpolationMethod& method);
    int SetLookupTable(const float* lutImageRow, unsigned int width, InterpolationMethod& method);
    /**
        @}
    */

    /** @brief sets pyramid size to pyramid and updates if neccessary */
    void SetPyramidSize(const int newsize);


    /** @brief calculates the bounding box in sink image,
        where to do the backward mapping.

        The resulting coordinates can be negative !
        The br position is the first pixel outside the region(as for a ROI).
        Implement this in a derived class for speeding up calculation.
    */
    virtual int GetBoundingBox(unsigned int /*srcwidth*/,
                               unsigned int /*srcheight*/,
                               unsigned int sinkwidth,
                               unsigned int sinkheight,
                               int &tlx, int &tly,
                               int &brx, int &bry)
    {
      tlx=0; tly=0;
      brx=sinkwidth; bry=sinkheight;
      return 0;
    }


    /** @brief uses corners of sink roi (and few other sample points)
     *      to estimate maximum local scaling (and thus required pyramid level)
    */
    virtual void UpdatePyramidSize(const ROI& ROIsink,
                                   int srcwidth, int srcheight);

    /** wrapper function to allow the lookup implementation to be
     * shared with other algorithms.
     * \attention Do not reimplement this method in any child class
     * implement into GetSourceCoordinates_().
     */
    inline int GetSourceCoordinates(const HomgPoint2D& sink,
                                    HomgPoint2D& source) const {
      if (pConcatenation_==NULL) {
        if (superSampling_==1.0) {
          return GetSourceCoordinates_(sink, source);
        }
        HomgPoint2D sink2(sink);
        sink2[0] = sink2[0] / superSampling_ - 0.5;
        sink2[1] = sink2[1] / superSampling_ - 0.5;
        return GetSourceCoordinates_(sink2, source);
      }
      // have a concatenation
      if (superSampling_==1.0) {
        int result = this->GetSourceCoordinates_(sink, source);
        HomgPoint2D sink2(source);
        if (result<0) return result;
        return pConcatenation_->GetSourceCoordinates(sink2, source);
      } else {
        HomgPoint2D sink2(sink);
        sink2[0] = sink2[0] / superSampling_ - 0.5;
        sink2[1] = sink2[1] / superSampling_ - 0.5;
        int result = this->GetSourceCoordinates_(sink2, source);
        if (result<0) return result;
        sink2 = source;
        return pConcatenation_->GetSourceCoordinates(sink2, source);
      }
    }


    /** @brief get jacobian (including concatenated mappings) */
    int GetJacobian(const HomgPoint2D& sink,
                    Matrix2x2<double>& Jacobian) const {
      if (pConcatenation_==NULL) return  GetJacobian_(sink, Jacobian);
      // evalute concatenated mapping ...
      Matrix2x2<double> tmpJacobian(MatrixIdentity);
      HomgPoint2D source2;
      if (GetSourceCoordinates_(sink, source2)<0) return -1;
      if (pConcatenation_->GetJacobian(source2, tmpJacobian) !=0) return -2;
      if (GetJacobian_(sink, Jacobian)!=0) return -3;
      Jacobian = tmpJacobian * Jacobian;
      return 0;
    }

    /** @brief avoid intermediate image and concatenate two backward mappings

    Imagine that the source image is transformed with pcon into the
    intermediate image, which is then transformed by this into the sink */
    void SetConcatenation(BIAS::BackwardMapping<InputStorageType,
                          OutputStorageType>* pCon) {
      pConcatenation_ = pCon;
    }

    /** @brief generates a siemens star like test image with lots of different
     *         frequencies to test backward mapping
     *
     *  @param testimage initialized image of user defined size which will
     *         be filled
     *  @param highFrequencyCross generates an additional high frequency line
     *  @param dark black image value, e.g. close to 0
     *  @param dark white image value, e.g. close to 255
     *  @param Hinv homography to transform pattern (if desired)
     *  @author koeser 10/2007
     */
    static void GenerateTestImage(Image<InputStorageType>& testimage,
                                  bool highFrequencyCross = true,
                                  InputStorageType dark = 5,
                                  InputStorageType bright = 127,
                                  const Matrix3x3<double>& Hinv =
                                  Matrix3x3<double>(MatrixIdentity));

  protected:

    /** @name geometrical coordinate mapping functions
        @{
    */

    /** @brief actual mapping function which takes a point from sink
        and computes a point in source and a local sampling frequency
        @return 0="ok", +1="out of image", -1="cannot be computed"(e.g. complex, point behind camera, ...)
        The distinction between "out of image" and "cannot be computed" ist useful as some
        methods sample the image before actual mapping to find out about range of local deformations,
        expected bounding box and so on and invalid points cannot be used while those just outside the
        image are still useful.
    */
    virtual int GetSourceCoordinates_(const HomgPoint2D& sink,
                                      HomgPoint2D& source) const;

    /** @brief numeric approximation of jacobian, reimplement in derived class
        to make an analytic jacobian */
    virtual int GetJacobian_(const HomgPoint2D& sink,
                             Matrix2x2<double>& Jacobian) const;

    /** @brief computes pyramid layer for trilinear filtering, e.g.
        localscale 0.5 means best resolution and second best are averaged */
    inline void ComputeLocalPyramidLayer_(const HomgPoint2D& p_sink,
                                          double& localscale) {
      // compute local sampling frequency for antialiasing
      // this is only an approximation, more exact would be SVD:
      // sqrt(s[0]) from J^T * J = U * s * V^T
      Matrix2x2<double> Jacobian;
      GetJacobian(p_sink, Jacobian);
      double
        dx = sqrt(Jacobian[0][0]*Jacobian[0][0]+
                  Jacobian[1][0]*Jacobian[1][0]),
        dy = sqrt(Jacobian[0][1]*Jacobian[0][1]+
                  Jacobian[1][1]*Jacobian[1][1]);
      if (dx<JACOBI_THRESH) dx = JACOBI_THRESH;
      if (dy<JACOBI_THRESH) dy = JACOBI_THRESH;
      // use larger step to select pyramid level
      localscale = (dx > dy)? log(dx)*ONE_OVER_LOG2 : log(dy)*ONE_OVER_LOG2;
      // clip scale to zero since we have no better resolution image than 0
      if (localscale<0.0) localscale = 0.0;
    }

    /**
        @}
    */

    /** @name pixel-value-mapping worker functions
        @{
    */

    /** @brief map image avoiding aliasing

    This is almost correct regarding the sampling theorem, because it uses
    trilinear or anisotropic filtering over a pyramid. The local
    sampling frequency
    is estimated by analyzing the jacobian of the homography as a R^2->R^2
    mapping. This is very slow. If possible, bicubic interpolation is done
    in the image, at borders (high pyramid levels) only bilinear.<br>
    In auto mode pyramid size is checked
    @author koeser 04/2005   */
    int MapTri_(Image<OutputStorageType>& sink);

    /** @brief fast trilinear grey without special features

    assumes one channel grey image, no alpha blending, no concatenation,
    no 2d offset (from auto image size) and no super sampling
    */
    int MapTrilinearGreySimple_(Image<OutputStorageType>& sink);


    /** @brief bilinear, bicubc and nearest neighbor interpolation     */
    int MapBi_(const Image<InputStorageType>& source,
               Image<OutputStorageType>& sink,
               InterpolationMethod interpolationQuality);




    /** @brief worker function, which handles bilinear/bicubic interpolation,
        borders and weights */
    inline int GetImageValue_(const double& xsource,
                              const double& ysource,
                              const Image<InputStorageType>& im,
                              unsigned int channel,
                              double& result);


    /** @brief worker function, which handles bilinear/bicubic interpolation,
        pyramid levels, borders and weights */
    int GetTrilinearImageValue_(const double& xsource,
                                const double& ysource,
                                const double& scale,
                                unsigned int channel,
                                double& result);


    /** @brief use jacobian of mapping function to do anisotropic
        antialiasing */
    int GetAnisotropicImageValue_(const double& xsource,
                                  const double& ysource,
                                  Matrix2x2<double>& Jacobian,
                                  unsigned int channel,
                                  double& thereturnvalue);

    /**
        @}
    */

    /** @brief clip interpolated value to outputstoragetype
        and merge old and new value according to alpha

        This function must be inplemented for diferent template types rather
        than for derived classes, e.g. if bicubic interpolation estimates the
        grey value 257.8 for an unsinged char image, this is rounded and
        clipped here to 255.0 since this is the maximum allowed uchar.
    */
    inline void
    ClipValue_(const double& value, OutputStorageType& newvalue);

    /** @brief clip computed region in dest image to roi
     */
    void ClipBoundingBoxToROICorners_(const BIAS::Image<OutputStorageType>&
                                      sink,
                                      int& tlx,  int& tly,
                                      int& brx,  int& bry) const;

    /** @brief compute maximum number of pyramid levels needed */
    virtual void BuildPyramid_(const Image<InputStorageType>& source,
                               bool forceInit = true,
                               int numlayers = -1);


    /// result in offsetX_,offsetY_
    virtual int CalcCoordOffset_(const Image<OutputStorageType>& sink,
                                 const Image<InputStorageType>& source);

    virtual int ChangeImgSize_(Image<OutputStorageType>& /*sink*/,
                               const Image<InputStorageType>& /*source*/);



    /// this contains a pyramid of the last source image for anti-aliasing
    PyramidImage<InputStorageType> pyramid_;
    /// the roi of the last source image
    ROI ROI_;

    /// pyramid size set by user or computed automatically ?
    bool autoPyramidSize_;
    unsigned pyramidSize_;

    /// which interpolation method are we using
    InterpolationMethod interpolationType_;

    /// if larger than 1, super resolution is active
    double superSampling_;

#ifdef BIAS_DEBUG
    /// The interpolated value is computed in double;
    /// if clipped value is not specialized for your output storage type
    /// (e.g. ushort), this boolean indicates that the generic conversion
    /// (e.g. double->ushort) could have introduced artefacts (65536.0 -> 1).
    bool rangeChecked_;
#endif

    /// set to true if possibility of aliasing is detected during mapping
    bool aliasing_;

    /// set to true if not every pixel could be computed
    bool incomplete_;

    /// lookup-tables for precomputed positions in source, same size as
    /// dst.PixelCount
    BWM_LUT_Entry_Nearest     *pLUT_Nea_;
    BWM_LUT_Entry_Bilinear    *pLUT_Bil_;
    BWM_LUT_Entry_Trilinear   *pLUT_Tri_;

    int                        pLUT_Nea_size;
    int                        pLUT_Bil_size;
    int                        pLUT_Tri_size;

    /// NOT IMPLEMENTED (YET ?)
    BWM_LUT_Entry_Anisotropic *pLUT_Ani_;
    BWM_LUT_Entry_Bicubic     *pLUT_Bic_;

    /// needed for newDist
    int offsetX_,offsetY_;

    /// source image width
    double width_;

    /// source image height
    double height_;

    /// has same size as overlap region used for blending
    BIAS::Image<float> alphaImg_;

    BackwardMapping<InputStorageType, OutputStorageType> *pConcatenation_;

  };

  // --------- inline code which is called for every pixel -----------------

  template <class InputStorageType, class OutputStorageType>
  inline int BackwardMapping<InputStorageType, OutputStorageType>::
  GetImageValue_(const double& x,
                 const double& y,
                 const Image<InputStorageType>& im,
                 unsigned int channel,
                 double& result)
  {
/*
    unsigned int UpperLeftX, UpperLeftY, LowerRightX, LowerRightY;
    im.GetROICorners(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY);

    // compute value only in image (no pyramid), first try bicubic
    if (x < (double)LowerRightX-2.0 && y < (double)LowerRightY-2.0 &&
        x >= double(UpperLeftX)+1.0 && y >= double(UpperLeftY) + 1.0 ) {

      result = im.BicubicInterpolation(x, y, (unsigned short int)channel);

    } else if (x < (double)LowerRightX-1.0 &&
               y < (double)LowerRightY-1.0 &&
               x > double(UpperLeftX) && y > double(UpperLeftY) ) {
      // we are near the border, bilinear
      result = im.BilinearInterpolation(x, y, (unsigned short int)channel);
    } else {
      // we are beyond or AT the border, we could do nearest neighbor here or
      // something else, but for now, skip this pixel
      return -1;
    }
*/
    // compute value only in image (no pyramid, using ROI)
    if (im.GetROI()->CheckBicubicInterpolation(x, y))
    {
      // first try bicubic interpolation
      result = im.BicubicInterpolation(x, y, (unsigned short int)channel);
    }
    else if (im.GetROI()->CheckBilinearInterpolation(x, y))
    {
      // we are near the border, use bilinear interpolation
      result = im.BilinearInterpolation(x, y, (unsigned short int)channel);
    }
    else
    {
      // we are beyond or AT the border, we could do nearest neighbor here or
      // something else, but for now skip this pixel
      return -1;
    }
    return 0;
  }

  // generic implementation tries to use limits
  template <class InputStorageType, class OutputStorageType>
  inline void BackwardMapping<InputStorageType, OutputStorageType>::
  ClipValue_(const double& value, OutputStorageType& newvalue)
  {
#ifdef BIAS_DEBUG
    rangeChecked_ = false;
#endif

    if (value>=double(std::numeric_limits<OutputStorageType>::max()))
      newvalue = std::numeric_limits<OutputStorageType>::max();
    else if (value<=double(std::numeric_limits<OutputStorageType>::min()))
      newvalue = std::numeric_limits<OutputStorageType>::min();
    else newvalue = OutputStorageType(value);
  }

  template <>
  inline void BackwardMapping<unsigned char, unsigned char>::
  ClipValue_(const double& value, unsigned char& newvalue)
  {
    if (value<0.0) newvalue = 0;
    else if (value>255.0) newvalue = 255;
    else newvalue =  ((unsigned char)rint(value));
  }

  template <>
  inline void BackwardMapping<float, float>::
  ClipValue_(const double& value, float& newvalue)
  {
    newvalue = float(value);
  }

  template <>
  inline void BackwardMapping<float, unsigned char>::
  ClipValue_(const double& value, unsigned char& newvalue)
  {
    if (value<0.0) newvalue = 0;
    else if (value>255.0) newvalue = 255;
    else newvalue =  ((unsigned char)rint(value));
  }

  template <>
  inline void BackwardMapping<unsigned char, float>::
  ClipValue_(const double& value, float& newvalue)
  {
    newvalue = float(value);
  }


#ifdef WIN32

  template <>
  void BackwardMapping<float, float>::ClipValue_(const double& value,
                                                 float& newvalue);

  template <>
  void BackwardMapping<unsigned char,
                       unsigned char>::ClipValue_(const double& value,
                                                  unsigned char& newvalue);


  template <>
  void BackwardMapping<unsigned char,
                       float>::ClipValue_(const double& value,
                                          float& newvalue);

  template <>
  void BackwardMapping<float,
                       unsigned char>::ClipValue_(const double& value,
                                                  unsigned char& newvalue);

#endif // WIN32

} // namespace


#endif