HarrisDetector.hpp

#ifndef _FDETECT_HARRIS_DETECTOR_HPP_
#define _FDETECT_HARRIS_DETECTOR_HPP_

#include <list>

#include <image/Image.hpp>
#include "image/roi.hpp"

#include "fdetect/Detector.hpp"
#include "fdetect/Descriptor.hpp"
#include "fdetect/InterestFeature.hpp"
#include "fdetect/HarrisDescriptor.hpp"

namespace jafar {
  namespace fdetect {

    class HarrisDetector : public Detector {
    public:
      enum Cornerness {HS_HARRIS, HS_EIGEN} ;
    public:
      HarrisDetector(float scale_factor = 1.0, float derivationSigma = 1.0, float convolutionSigma = 2.0, jafar::fdetect::HarrisDetector::Cornerness cornerness = HS_EIGEN, float threshold = 0.0, bool fixnumber = true, int desirednum = 1000);
    public:
      virtual jafar::fdetect::DetectionResult detectIn(jafar::image::Image const& image);
      virtual bool detectBestPntInRoi(jafar::image::Image const& image,
                                      InterestFeature * featPtr, const jafar::image::ROI * roiPtr = 0 );
    public:
      //         inline double scaleFactor() { return m_scale_factor; }
      float scaleFactor() const;
      float derivationSigma() const;
      float convolutionSigma() const;
      float threshold() const;
      Cornerness cornerness() const;
      int desiredNumber() const;
      int derivationSize() const;
      int convolutionSize() const;
      bool isFixNumber() const;

    private:
      void gaussianDerivativesSA(jafar::image::Image const& image);
      void gaussianDerivatives(jafar::image::Image const& image);
      void gaussianDerivatives(jafar::image::Image const& image, const jafar::image::ROI & roi);
      void gaussianConvolutionSA(jafar::image::Image const& image);
      void gaussianConvolution(jafar::image::Image const& image);
      void gaussianConvolution(jafar::image::Image const& image, const jafar::image::ROI & roi);
      void computeHarrisCurvature (jafar::image::Image const& image);
      void computeEigenCurvature (jafar::image::Image const& image);
      void computeEigenCurvatureWithMax (jafar::image::Image const& image, int pixMax[2], float * scoreMax);
      void cornerCharacteristics(DetectionResult& points, jafar::image::Image const& image);
      //         void cornerCharacteristicsSA(vInterestFeatures& points, jafar::image::Image const& image);
      template<bool autoThreshold, bool useScale>
      DetectionResult extractCorners (jafar::image::Image const& image);
    private:
      float m_scale_factor, m_derivationSigma, m_convolutionSigma, m_threshold;
      Cornerness m_cornerness;
      int m_desirednum, m_derivationSize, m_convolutionSize;
      bool m_fixnumber;
      struct HData {
        float im_x, im_y, im_xx, im_xy, im_yy, im_conv_xx, im_conv_xy, im_conv_yy, im_high_curv, im_low_curv;
      };
      HData* m_data;
    };

    template<bool autoThreshold, bool useScale>
    DetectionResult HarrisDetector::extractCorners (jafar::image::Image const& image)
    {
      struct Corner {
        Corner(float* ihc, float pc) : im_high_corner(ihc), pix_center(pc) { }
        float* im_high_corner;
        float pix_center;
      };

      int nblig = image.height();
      int nbcol = image.width();

      // Variables used when autoThreshold = true
      unsigned int d_number; // nimber of points
      std::list<InterestFeature*> lcorners; // corners list
      DetectionResult vcorners;

      /* Cornerness Type */
      //       switch (params->cornernessType) {
      //         case HS_HARRIS:
      //           pix_center_org = data->im_high_curv;
      //           pix_center_org1 = data->im_low_curv;
      //           break;
      //         case HS_EIGEN:
      //           pix_center_org = data->im_low_curv;
      //           pix_center_org1 = data->im_high_curv;
      //           break;
      //         default:
      //           JFR_PRECOND(false,"ERR_HS_BAD_PARAM");
      //       }

      if(autoThreshold)
        {
          if( m_desirednum == -1) {
            d_number = (int) ((nbcol * nbcol)/400);
          } else {
            d_number = m_desirednum;
          }
          if (d_number > 5000) d_number = 5000;
        }

      int w_size;         /* local window size */
      int half_w_size;    /* local window half size */

      /* Initialisation of some variable use when useScale == true */
      if(useScale)
        {
          w_size = m_derivationSize;
          half_w_size = (int)(0.5*(w_size - 1));
        } else {
        w_size = 5;
        half_w_size = 2;
      }

      /* Computation of reduced first address */
      int redu = ( (m_convolutionSize - 1) / 2 + (m_derivationSize  -1) );

      /********************************************************************/
      /*  Corner Detection Routin */
      /********************************************************************/

      for (int i = redu; i < nblig - redu; i++) {

        HData* p_data_center = m_data + (i * nbcol) + redu;

        HData* pix_win = p_data_center - half_w_size*nbcol - half_w_size;

        HData* p_data_top1 = p_data_center - nbcol;
        HData* p_data_top1_right1 = p_data_top1+1;
        HData* p_data_top1_left1 = p_data_top1-1;

        HData* p_data_down1 = p_data_center + nbcol;
        HData* p_data_down1_right1 = p_data_down1+1;
        HData* p_data_down1_left1 = p_data_down1-1;

        HData* p_data_right1 = p_data_center + 1;
        HData* p_data_left1 = p_data_center - 1;


        for (int j = redu; j < nbcol - redu; j++) {
          if ( (std::min(p_data_center->im_low_curv,p_data_center->im_high_curv) > m_threshold) )
            {
              bool OK_flag = true;
              for (int ii = 0; ii < w_size && OK_flag; ii++) {
                HData* pix_tmp = pix_win + ii * nbcol;
                for (int jj = 0; jj < w_size; jj++) {
                  if (std::min(pix_tmp->im_low_curv,pix_tmp->im_high_curv) >
                      std::min(p_data_center->im_low_curv,p_data_center->im_high_curv)) { OK_flag = false; break; };
                  pix_tmp++;
                }
              }
              if( OK_flag && ( !autoThreshold || lcorners.size() < d_number || ((HarrisDescriptor*)lcorners.back()->descriptor())->low() < std::min(p_data_center->im_high_curv,p_data_center->im_low_curv) ) )
                {
                  float b = 0.125*((p_data_top1_right1->im_high_curv) - (p_data_top1_left1->im_high_curv) +
                                   2.0*(p_data_right1->im_high_curv) - 2.0*(p_data_left1->im_high_curv) +
                                   (p_data_down1_right1->im_high_curv) - (p_data_down1_left1->im_high_curv) );
                  /* [-1 0 1]
                     [-2 0 2]
                     [-1 0 1] */

                  float c = 0.125*((p_data_down1_left1->im_high_curv) - (p_data_top1_left1->im_high_curv) +
                                   2.0*(p_data_down1->im_high_curv) - 2.0*(p_data_top1->im_high_curv) +
                                   (p_data_down1_right1->im_high_curv) - (p_data_top1_right1->im_high_curv));
                  /* [-1 -2 -1]
                     [ 0  0  0]
                     [ 1  2  1] */

                  float d = 0.25*((p_data_top1_left1->im_high_curv) + 2.0*(p_data_left1->im_high_curv) + (p_data_down1_left1->im_high_curv) -
                                  2.0*(p_data_top1->im_high_curv) - 4.0*(p_data_center->im_high_curv) - 2.0*(p_data_down1->im_high_curv) +
                                  (p_data_top1_right1->im_high_curv) + 2.0*(p_data_right1->im_high_curv) + (p_data_down1_right1->im_high_curv));
                  /* [1 -2 1]
                     [2 -4 2]
                     [1 -2 1] */

                  float e = 0.25*(p_data_top1_left1->im_high_curv - p_data_down1_left1->im_high_curv -
                                  p_data_top1_right1->im_high_curv + p_data_down1_right1->im_high_curv);
                  /* [ 1 0 -1]
                     [ 0 0  0]
                     [-1 0  1] */

                  float f = 0.25*((p_data_top1_left1->im_high_curv) - 2.0*(p_data_left1->im_high_curv) + (p_data_down1_left1->im_high_curv)+
                                  2.0*(p_data_top1->im_high_curv) - 4.0*(p_data_center->im_high_curv) + 2.0*(p_data_down1->im_high_curv)+
                                  (p_data_top1_right1->im_high_curv) - 2.0*(p_data_right1->im_high_curv) + (p_data_down1_right1->im_high_curv));
                  /* [ 1  2  1]
                     [-2 -4 -2]
                     [ 1  2  1] */

                  /* Quadratic fit function is
                     F(x,y)=a + b x + c y + d/2 x^2 + e xy + f/2 y^2
                     so, we find subpixel value setting the partial dervative
                     of F to Zero ;
                     [d e][offset_y]+[b]=[0];
                     [e f][offset_x] [c] [0] */

                  float det = d*f - e*e;

                  if (det != 0) {
                    float off_x = (b*e - c*d)/det;
                    float off_y = (c*e - b*f)/det;
                    if ((fabs(off_x) < 1.0)&&(fabs(off_y) < 1.0))
                      {
                        InterestFeature* ip = useScale ? new InterestFeature(j+off_y, i+off_x) : new InterestFeature(j+off_y/2, i+off_x/2);
                        double lHigh, lLow;
                        if( p_data_center->im_high_curv > p_data_center->im_low_curv)
                          {
                            lHigh = p_data_center->im_high_curv ;
                            lLow = p_data_center->im_low_curv;
                          } else {
                          lHigh = p_data_center->im_low_curv;
                          lLow = p_data_center->im_high_curv ;
                        }
                        ip->setDescriptor( new HarrisDescriptor( ip, lHigh, lLow, image ) );
                        ip->setQuality(lLow);
                        if(autoThreshold)
                          {
                            if(lcorners.empty())
                              {
                                lcorners.push_back( ip );
                              } else {
                              if(lcorners.size() >= d_number && lLow > ((HarrisDescriptor*)lcorners.back()->descriptor())->low())
                                { // remove last element, the totalNumber stay equal to d_number
                                  delete lcorners.back();
                                  lcorners.pop_back();
                                }
                              if(lcorners.size() < d_number )
                                {
                                  std::list<InterestFeature*>::iterator it;
                                  if(lLow < ((HarrisDescriptor*)lcorners.back()->descriptor())->low())
                                    {
                                      lcorners.push_back(ip );
                                    } else {
                                    // insert the new corner at his right place
                                    for(it = lcorners.begin(); it != lcorners.end(); it++)
                                      {
                                        if( lLow >= ((HarrisDescriptor*)((*it)->descriptor()))->low())
                                          {
                                            lcorners.insert(it, ip);
                                            break;
                                          }
                                      }
                                  }
                                } else {
                                delete ip;
                              }
                            }
                          } else {
                          vcorners.addFeature(ip);
                        }
                      }
                  }
                }
            }
          p_data_center++;
          pix_win++;
          p_data_top1++;
          p_data_top1_right1++;
          p_data_top1_left1++;

          p_data_down1++;
          p_data_down1_right1++;
          p_data_down1_left1++;

          p_data_right1++;
          p_data_left1++;
        }

      }
      if(autoThreshold) {
        for(std::list<InterestFeature*>::iterator it = lcorners.begin(); it != lcorners.end(); ++it)
          {
            InterestFeature* ip = *it;
            vcorners.addFeature(ip);
          }
      }
      return vcorners;
    }
  }

  // namespace fdetect_v2 {

  //   /**
  //    * This class is the implements a detector of Harris points as described in : "A Combined
  //    * Corner and Edge Detector" par C. Harris et M. Stephens en 1988.
  //    * \ingroup fdetect
  //    */
  //   class HarrisDetector : public fdetect_v2::Detector<fdetect_v2::HarrisDescriptor> 
  //  {
  //  public:
  //    enum Cornerness {HS_HARRIS, HS_EIGEN} ;
  //  public:
  //    /**
  //     * @param scale_factor the scale used to search the interest point, it usually set to 1.0,
  //     *        but it might prove usefull to have it set to a different value if you know that
  //     *        there is a scale difference between the two images you are trying to match.
  //     * @param derivationSigma the sigma used by the derivation kernel
  //     * @param convolutionSigma the sigma used by the convolution kernel
  //     * @param cornerness the type of cornerness (either using eigen values or the "harris" cornerness)
  //     * @param threshold use to select point that are correctly different from their neighbourhood
  //     * @param fixnumber set it to true if you want a given number of point which have the maximum
  //     *        cornerness (note that the threshold is still used to avoid selecting points which won't
  //     *        be good enought)
  //     * @param desirednum the number of points (if set to -1 and fixnumber is true, an automatic number
  //     *        of points is used)
  //     */
  //    HarrisDetector(float scale_factor = 1.0, float derivationSigma = 1.0, float convolutionSigma = 2.0, jafar::fdetect_v2::HarrisDetector::Cornerness cornerness = HS_EIGEN, float threshold = 0.0, bool fixnumber = true, int desirednum = 1000);
  //  public:
  //    virtual fdetect_v2::DetectionResult<fdetect_v2::HarrisDescriptor> detectIn(image::Image const& image);
  //    virtual bool detectBestPntInRoi(jafar::image::Image const& image,
  //                                    InterestFeature<fdetect_v2::HarrisDescriptor> * featPtr, 
  //                                    const jafar::image::ROI * roiPtr = 0 );
  //  public:
  //    //         inline double scaleFactor() { return m_scale_factor; }
  //     float scaleFactor() const;
  //     float derivationSigma() const;
  //     float convolutionSigma() const;
  //     float threshold() const;
  //     Cornerness cornerness() const;
  //     int desiredNumber() const;
  //     int derivationSize() const;
  //     int convolutionSize() const;
  //     bool isFixNumber() const;

  //  private:
  //    void gaussianDerivativesSA(jafar::image::Image const& image);
  //    void gaussianDerivatives(jafar::image::Image const& image);
  //    void gaussianDerivatives(jafar::image::Image const& image, const jafar::image::ROI & roi);
  //    void gaussianConvolutionSA(jafar::image::Image const& image);
  //    void gaussianConvolution(jafar::image::Image const& image);
  //    void gaussianConvolution(jafar::image::Image const& image, const jafar::image::ROI & roi);
  //    void computeHarrisCurvature (jafar::image::Image const& image);
  //    void computeEigenCurvature (jafar::image::Image const& image);
  //    void computeEigenCurvatureWithMax (jafar::image::Image const& image, int pixMax[2], float * scoreMax);
  //    void cornerCharacteristics(DetectionResult<HarrisDescriptor>& points, jafar::image::Image const& image);
  //    //         void cornerCharacteristicsSA(vInterestFeatures& points, jafar::image::Image const& image);
  //    template<bool autoThreshold, bool useScale>
  //    DetectionResult<fdetect_v2::HarrisDescriptor> extractCorners (jafar::image::Image const& image);
  //  private:
  //    float m_scale_factor, m_derivationSigma, m_convolutionSigma, m_threshold;
  //    Cornerness m_cornerness;
  //    int m_desirednum, m_derivationSize, m_convolutionSize;
  //    bool m_fixnumber;
  //    struct HData {
  //      float im_x, im_y, im_xx, im_xy, im_yy, im_conv_xx, im_conv_xy, im_conv_yy, im_high_curv, im_low_curv;
  //    };
  //    HData* m_data;
  //   };

  //   template<bool autoThreshold, bool useScale>
  //   DetectionResult<fdetect_v2::HarrisDescriptor> HarrisDetector::extractCorners (jafar::image::Image const& image)
  //   {
  //    using namespace jafar::fdetect_v2;
  //     struct Corner 
  //    {
  //       Corner(float* ihc, float pc) : im_high_corner(ihc), pix_center(pc) { }
  //       float* im_high_corner;
  //       float pix_center;
  //     };

  //     int nblig = image.height();
  //     int nbcol = image.width();

  //    // Variables used when autoThreshold = true
  //     unsigned int d_number; // nimber of points
  //     std::list< InterestFeature<HarrisDescriptor>* > lcorners; // corners list
  //     DetectionResult<HarrisDescriptor> vcorners;

  //     /* Cornerness Type */
  //    //       switch (params->cornernessType) {
  //    //         case HS_HARRIS:
  //    //           pix_center_org = data->im_high_curv;
  //    //           pix_center_org1 = data->im_low_curv;
  //    //           break;
  //    //         case HS_EIGEN:
  //    //           pix_center_org = data->im_low_curv;
  //    //           pix_center_org1 = data->im_high_curv;
  //    //           break;
  //    //         default:
  //    //           JFR_PRECOND(false,"ERR_HS_BAD_PARAM");
  //    //       }

  //     if(autoThreshold)
  //    {
  //      if( m_desirednum == -1) 
  //      {
  //        d_number = (int) ((nbcol * nbcol)/400);
  //      } 
  //      else 
  //      {
  //        d_number = m_desirednum;
  //      }
  //      if (d_number > 5000) 
  //        d_number = 5000;
  //    }

  //     int w_size;         /* local window size */
  //     int half_w_size;    /* local window half size */

  //     /* Initialisation of some variable use when useScale == true */
  //     if(useScale)
  //    {
  //      w_size = m_derivationSize;
  //      half_w_size = (int)(0.5*(w_size - 1));
  //    } 
  //    else 
  //    {
  //       w_size = 5;
  //       half_w_size = 2;
  //     }
      
  //     /* Computation of reduced first address */
  //     int redu = ( (m_convolutionSize - 1) / 2 + (m_derivationSize  -1) );

  //     /********************************************************************/
  //     /*  Corner Detection Routin */
  //     /********************************************************************/

  //     for (int i = redu; i < nblig - redu; i++) 
  //    {
  //       HData* p_data_center = m_data + (i * nbcol) + redu;

  //       HData* pix_win = p_data_center - half_w_size*nbcol - half_w_size;

  //       HData* p_data_top1 = p_data_center - nbcol;
  //       HData* p_data_top1_right1 = p_data_top1+1;
  //       HData* p_data_top1_left1 = p_data_top1-1;

  //       HData* p_data_down1 = p_data_center + nbcol;
  //       HData* p_data_down1_right1 = p_data_down1+1;
  //       HData* p_data_down1_left1 = p_data_down1-1;

  //       HData* p_data_right1 = p_data_center + 1;
  //       HData* p_data_left1 = p_data_center - 1;


  //       for (int j = redu; j < nbcol - redu; j++) 
  //      {
  //         if ( (std::min(p_data_center->im_low_curv,p_data_center->im_high_curv) > m_threshold) )
  //        {
  //          bool OK_flag = true;
  //          for (int ii = 0; ii < w_size && OK_flag; ii++) 
  //          {
  //            HData* pix_tmp = pix_win + ii * nbcol;
  //            for (int jj = 0; jj < w_size; jj++) 
  //            {
  //              if (std::min(pix_tmp->im_low_curv,pix_tmp->im_high_curv) >
  //                  std::min(p_data_center->im_low_curv,p_data_center->im_high_curv)) 
  //              { 
  //                OK_flag = false; 
  //                break; 
  //              };
  //              pix_tmp++;
  //            }
  //          }
  //          if( OK_flag && ( !autoThreshold || 
  //                           lcorners.size() < d_number || 
  //                           (lcorners.back()->descriptor()).low() < std::min(p_data_center->im_high_curv,p_data_center->im_low_curv) ) )
  //          {
  //            float b = 0.125*((p_data_top1_right1->im_high_curv) - (p_data_top1_left1->im_high_curv) +
  //                             2.0*(p_data_right1->im_high_curv) - 2.0*(p_data_left1->im_high_curv) +
  //                             (p_data_down1_right1->im_high_curv) - (p_data_down1_left1->im_high_curv) );
  //            /* [-1 0 1]
  //               [-2 0 2]
  //               [-1 0 1] */
              
  //            float c = 0.125*((p_data_down1_left1->im_high_curv) - (p_data_top1_left1->im_high_curv) +
  //                             2.0*(p_data_down1->im_high_curv) - 2.0*(p_data_top1->im_high_curv) +
  //                             (p_data_down1_right1->im_high_curv) - (p_data_top1_right1->im_high_curv));
  //            /* [-1 -2 -1]
  //               [ 0  0  0]
  //               [ 1  2  1] */
              
  //            float d = 0.25*((p_data_top1_left1->im_high_curv) + 2.0*(p_data_left1->im_high_curv) + (p_data_down1_left1->im_high_curv) -
  //                            2.0*(p_data_top1->im_high_curv) - 4.0*(p_data_center->im_high_curv) - 2.0*(p_data_down1->im_high_curv) +
  //                            (p_data_top1_right1->im_high_curv) + 2.0*(p_data_right1->im_high_curv) + (p_data_down1_right1->im_high_curv));
  //            /* [1 -2 1]
  //               [2 -4 2]
  //               [1 -2 1] */
              
  //            float e = 0.25*(p_data_top1_left1->im_high_curv - p_data_down1_left1->im_high_curv -
  //                            p_data_top1_right1->im_high_curv + p_data_down1_right1->im_high_curv);
  //            /* [ 1 0 -1]
  //               [ 0 0  0]
  //               [-1 0  1] */
              
  //            float f = 0.25*((p_data_top1_left1->im_high_curv) - 2.0*(p_data_left1->im_high_curv) + (p_data_down1_left1->im_high_curv)+
  //                            2.0*(p_data_top1->im_high_curv) - 4.0*(p_data_center->im_high_curv) + 2.0*(p_data_down1->im_high_curv)+
  //                            (p_data_top1_right1->im_high_curv) - 2.0*(p_data_right1->im_high_curv) + (p_data_down1_right1->im_high_curv));
  //            /* [ 1  2  1]
  //               [-2 -4 -2]
  //               [ 1  2  1] */
              
  //            /* Quadratic fit function is
  //               F(x,y)=a + b x + c y + d/2 x^2 + e xy + f/2 y^2
  //               so, we find subpixel value setting the partial dervative
  //               of F to Zero ;
  //               [d e][offset_y]+[b]=[0];
  //               [e f][offset_x] [c] [0] */
              
  //            float det = d*f - e*e;
              
  //            if (det != 0) 
  //            {
  //              float off_x = (b*e - c*d)/det;
  //              float off_y = (c*e - b*f)/det;
  //              if ((fabs(off_x) < 1.0)&&(fabs(off_y) < 1.0))
  //              {
  //                InterestFeature<HarrisDescriptor>* ip = useScale ? new InterestFeature<HarrisDescriptor>(j+off_y, i+off_x) : new InterestFeature<HarrisDescriptor>(j+off_y/2, i+off_x/2);
  //                double lHigh, lLow;
  //                if( p_data_center->im_high_curv > p_data_center->im_low_curv)
  //                {
  //                  lHigh = p_data_center->im_high_curv ;
  //                  lLow = p_data_center->im_low_curv;
  //                } 
  //                else 
  //                {
  //                  lHigh = p_data_center->im_low_curv;
  //                  lLow = p_data_center->im_high_curv ;
  //                }
  //                HarrisDescriptor hd( lHigh, lLow, ip->u(), ip->v(), image );
  //                ip->setDescriptor( hd );
  //                ip->setQuality(lLow);
  //                if(autoThreshold)
  //                {
  //                  if(lcorners.empty())
  //                  {
  //                    lcorners.push_back( ip );
  //                  } 
  //                  else 
  //                  {
  //                    if(lcorners.size() >= d_number && lLow > (lcorners.back()->descriptor()).low())
  //                    { // remove last element, the totalNumber stay equal to d_number
  //                      delete lcorners.back();
  //                      lcorners.pop_back();
  //                    }
  //                    if(lcorners.size() < d_number )
  //                    {
  //                      std::list< InterestFeature<HarrisDescriptor>* >::iterator it;
  //                      if(lLow < (lcorners.back()->descriptor()).low())
  //                      {
  //                        lcorners.push_back(ip );
  //                      } 
  //                      else 
  //                      {
  //                        // insert the new corner at his right place
  //                        for(it = lcorners.begin(); it != lcorners.end(); it++)
  //                        {
  //                          if( lLow >= ((*it)->descriptor()).low() )
  //                          {
  //                            lcorners.insert(it, ip);
  //                            break;
  //                          }
  //                        }
  //                      }
  //                    } 
  //                    else 
  //                    {
  //                      delete ip;
  //                    }
  //                  }
  //                } 
  //                else 
  //                {
  //                  vcorners.push_back(ip);
  //                }
  //              }
  //            }
  //          }
  //        }
  //         p_data_center++;
  //         pix_win++;
  //         p_data_top1++;
  //         p_data_top1_right1++;
  //         p_data_top1_left1++;

  //         p_data_down1++;
  //         p_data_down1_right1++;
  //         p_data_down1_left1++;

  //         p_data_right1++;
  //         p_data_left1++;
  //       }

  //     }
  //     if(autoThreshold) 
  //    {
  //       for(std::list< InterestFeature<HarrisDescriptor>* >::const_iterator it = lcorners.begin(); 
  //          it != lcorners.end(); 
  //          ++it)
  //      {
  //        vcorners.push_back(*it);
  //      }
  //     }
  //     return vcorners;
  //   }
  // }
}

#endif