CornerDetectorGradient.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 __CornerDetectorGradient_hh__
#define __CornerDetectorGradient_hh__

#include "CornerDetectorBase.hh"
#include <Filter/StructureTensor.hh>
#include <Filter/Gauss.hh>
#include <Base/Math/Matrix2x2.hh>


namespace BIAS {
  /**@cond HIDDEN_SYMBOLS
     @class Feat
     @ingroup g_feature
     @brief  internal class for feature passing only
     @author koeser 09/2004 
  */ 
  class BIASFeatureDetector_EXPORT Feat {
  public:
    int x, y;
    double val;
   
    bool operator<(const Feat& r) const
    { return val>r.val; }
   
    Feat& operator=(const Feat& f)
    { x=f.x; y=f.y; val=f.val; return *this; }
     
  }; // end of class Feat
/** \endcond */
 
  inline std::ostream& operator<<(std::ostream & os, const Feat& f)
  { os << "("<<f.x<<", "<<f.y<<") = "<<f.val; return os; }



  /** @class CornerDetectorGradient
      @ingroup g_feature
      @brief base class for all gradient based corner detectors

      provides  functionality to compute the structure tensor
      | gxx gxy |
      | gyx gyy |
      
      The function detect low pass filters the image, computes gradients and
      the structure tensor and then calls _ComputeCornerness on the structure
      tensor image. This function must be implemented in derived classes, it 
      basically converts the 2x2 structure tensor to a scalar value which
      describes how much cornerlike the region around that pixel is, which 
      is saved in the cornerness image. Values above a threshold are pushed
      into an internal list.
      Afterwards the EnforceMinimumDistance function is called to sort features
      by quality and add them to the final list. New points are accepted only 
      if the have a minimum distance to all previously extracted points.
      In the detect function you may pass a non-empty vector of points,
      which are regarded as "existing corners" and none of the new points may
      lie withing the critical distance around the old ones. This gives the
      possibility to detect new features for a tracker, which has tracked
      some points (where no new features are desired) and lost some features.
      
      @author woelk/koeser 09/2004 */
  template <class StorageType, class CalculationType>
  class BIASFeatureDetector_EXPORT CornerDetectorGradient : 
    public CornerDetectorBase<StorageType>
  {
  public:
    CornerDetectorGradient();
    virtual ~CornerDetectorGradient();

    /** @brief detect corners in a grey image

        see documentation of BIAS::CornerDetectorGradient class for more info
        @param image grey value input image
        @param p all points in p _not_ lying at infinity (p[i][2]!=0) are 
        assumed to be corners of very good quality, in the region around these
        points no new features are detected, leave empty if you just want
        to detect some features in the image
        In any case p is cleared and only new detected points are added before
        returning. Example: assume MaxNum=100, you pass a vector p with 40
        points of which 12 have the third coordinate w=0. You will get back a 
        vector of maximum size 72, it may be smaller if Distance and 
        MinCornerness have too large values
        @param quality after returning, for each point in p there is a 
        quality in quality, no input requirements  */
    virtual int Detect(const Image<StorageType>& image, 
                       std::vector<HomgPoint2D>& p, 
                       std::vector<QUAL>& quality,
                       std::vector<Matrix2x2<double> > *cov = NULL);
  
     /** @brief detect corners given the gradient images

        see documentation of BIAS::CornerDetectorGradient class for more info
        @param gradx image containing gradient x component
        @param grady image containing gradient y component
        @param p all points in p _not_ lying at infinity (p[i][2]!=0) are 
        assumed to be corners of very good quality, in the region around these
        points no new features are detected, leave empty if you just want
        to detect some features in the image.
        In any case p is cleared and only new detected points are added before
        returning. Example: assume MaxNum=100, you pass a vector p with 40
        points of which 12 have the third coordinate w=0. You will get back a 
        vector of maximum size 72, it may be smaller if Distance and 
        MinCornerness have too large values 
        @param quality after returning, for each point in p there is a 
        quality in quality, no input requirements  */
    virtual int Detect(const Image<CalculationType>& gradx, 
                       const Image<CalculationType>& grady, 
                       std::vector<HomgPoint2D>& p, 
                       std::vector<QUAL>& quality,
                       std::vector<Matrix2x2<double> > *cov = NULL);


    /** @brief detect corners avoiding double calculation of cornerness image
      
        see documentation of BIAS::CornerDetectorGradient class for more info
        @param cornerness image with already computed cornerness
        @param p all points in p _not_ lying at infinity (p[i][2]!=0) are 
        assumed to be corners of very good quality, in the region around these
        points no new features are detected, leave empty if you just want
        to detect some features in the image 
        In any case p is cleared and only new detected points are added before
        returning. Example: assume MaxNum=100, you pass a vector p with 40
        points of which 12 have the third coordinate w=0. You will get back a 
        vector of maximum size 72, it may be smaller if Distance and 
        MinCornerness have too large values
        @param quality after returning, for each point in p there is a 
        quality in quality, no input requirements  */
    virtual int DetectFromCornerness(const Image<CalculationType>& cornerness, 
                                     std::vector<HomgPoint2D>& p, 
                                     std::vector<QUAL>& quality);


    /** @brief detect corners avoiding double calculation of cornerness image
        @author Dennis Herzog, June 2005 */
    virtual int DetectFromCornerness(const Image<CalculationType> &cornerness, 
                     const std::vector<Vector2<int> > &region,
                     std::vector<HomgPoint2D> &p, 
                     std::vector<QUAL> &quality);

    /** @brief calculates the cornerness using the gradients gx and gy
        @author woelk 09/2004 */
    virtual int Cornerness(const Image<CalculationType>& gradx, 
                   const Image<CalculationType>& grady, 
                   Image<CalculationType>& cornerness);


    /** @brief rets the min distance of features 
        (simple parameter access function - does nothing else)
        @author Dennis Herzog, 2005-06-22 */    
    inline int GetMinDistance() const
    { return _MinDistance; }
    
    /** @brief sets the min distance of features 
        (simple parameter access function - does nothing else)
        @author Dennis Herzog, 2005-06-22 */    
    inline void SetMinDistance(int mindist)
    { _MinDistance = mindist; }

    /** @brief rets minimum cornernes of detectable features 
        (simple parameter access function - does nothing else)
        @author Dennis Herzog, 2005-07-04 */    
    inline double GetMinCornerness()
    { return _MinCornerness; }

    /** @brief sets minimum cornernes of detectable features 
        @author woelk 02/2006 */    
    virtual inline void SetMinCornerness(double min_cornerness)
    { _MinCornerness = min_cornerness; }

    /** @author woelk 02/2006 */
    inline void SetSubpixelAccuracy(bool subpixel_accuracy)
    { _DoRefinement = subpixel_accuracy; }

    virtual inline void SetStructureTensorHalfWinSize(int halfwinsize)
    { _st.SetHalfWinSize(halfwinsize); }

    virtual inline const Image<CalculationType>& GetCornernessImage() const
    { return _Cornerness; }
    
  protected:    
   
    // exports for dll (cmenk)
    //BIAS_EXPORT_STL_VECTOR( BIASFeatureDetector_EXPORT, class Feat )

    /** calculates the cornerness using the src image
        @author woelk 09/2004 */
    int _CalcCornerness(const Image<StorageType>& src, 
                        Image<CalculationType>& cornerness);

    /** calculates the cornerness using the gradients gx and gy
        @author woelk 09/2004 */
    int _CalcCornerness(const Image<CalculationType>& gradx, 
                        const Image<CalculationType>& grady, 
                        Image<CalculationType>& cornerness);
   
    /** @brief the main function which should be overwritten in derived 
        classes, no implementation in this class

        computes the cornerness from the structure tensor image 
        Uses structure tensor sgxx, _sgxy, _sgyy as input */
    virtual int _ComputeCornerness(Image<CalculationType>& Cornerness);

    /** @brief extracts the local maxima checking a 3x3 region
        into _FeatList */
    void _ExtractLocalMaxima(const Image<CalculationType>& Cornerness);

    /** @brief selects features sorted by quality from _FeatList such that they
        have a min distance as given in the param object

        assumes feature map is zeroed
        @author woelk 09/2004 */
    int _EnforceMinimumDistance(std::vector<HomgPoint2D>& p, 
                                std::vector<QUAL>& quality);

    /** @brief fit quadratic surface around (row;col) to find the maximum 
        cornerness subpixel position and save it in x;y
        @param col column in image where refinement should take place
        @param row row in image where refinement should take place
        @param x output, refined coordinate
        @param y output, refined coordinate
        @param SubPixelCornerness cornerness value at subpixel maximum
        @author koeser 01/2004 */
    bool _RefineCornerPosition(int col, int row, 
                               double &x, double &y, QUAL &SubPixelCornerness);

    /** @brief refines corner position and quality for all corners in list,
        internal Cornerness_ image must be valid 
        @param pos input/output of old/new position
        @param quality holds new cornernesses afterwards
        @author koeser 09/2004 */
    inline void _RefineCorners(std::vector<HomgPoint2D>& pos,
                               std::vector<QUAL>& quality) {
      BIASASSERT(pos.size()==quality.size());
      std::vector<QUAL>::iterator itq = quality.begin();
      for (std::vector<HomgPoint2D>::iterator itp=pos.begin(); 
           itp!=pos.end(); itp++, itq++) 
    _RefineCornerPosition((int)rint((*itp)[0]), (int)rint((*itp)[1]),
                              (*itp)[0], (*itp)[1], *itq);
    }

    /// @brief frees _Cornerness, _sgxx, _sgxy and _sgyy
    void _DeleteInternalMem();

    /// @brief allocates _Cornerness , _sgxx, _sgxy and _sgyy
    void _AllocInternalMem(const int width, const int height);

    /// @brief zeros the feature map 
    inline void _ZeroFeatureMap()
    { memset((void*)_FeatureMap[0], 0, 
             sizeof(bool)*_Cornerness.GetWidth()*_Cornerness.GetHeight()); 
    _NumberFeaturesSoFar = 0;};

    /// @brief fills feature map with every point from p which isnt at infinity
    void _FillFeatureMap(std::vector<HomgPoint2D>& p);

    /** @brief returns the inverse structure tensor, multiplied by noise^2
        @author woelk 03/2006 */
    void GetCov_(const std::vector<HomgPoint2D>& p, 
                 std::vector<Matrix2x2<double> >& cov,
                 const double noise = 5) const;

    /// structure tensor computation object
    StructureTensor<CalculationType, CalculationType> _st;

    /// internal structure tensor images
    Image<CalculationType> _sgxx, _sgxy, _sgyy;

    /// resulting image which tells for each pixel how cornerlike it is
    Image<CalculationType> _Cornerness;

    /// array used to enforce a min distance / neighbor suppression
    bool **_FeatureMap;

    /// temporary list of detected corners
    std::vector<class Feat> _FeatList;

    /// Gauss for pre-smoothing
    Gauss<StorageType, StorageType> _gauss;

    /// initialized with the old features provided by user, then counted up
    /// and compared against _MaxNumFeatures from Base class
    int _NumberFeaturesSoFar;

    // parameter 
    int _MinDistance; ///< minimal distance between points
    double _MinCornerness; ///< threshold to accept a feature
    bool _DoRefinement; ///< tells whether subpixel refinement should be done
  }; // class

} // namespace

#ifndef __APPLE__
// DLL JW
namespace std {
  template class BIASFeatureDetector_EXPORT std::vector<class BIAS::Feat>;
}
#endif

#endif // __CornerDetectorGradient_hh__