Merge pull request #933 from Yorwba:harris-laplace-detector

* Added files from http://code.opencv.org/issues/1053

The code structure had to be changed because of the
features2d/xfeatures2d split and to avoid changing
multiple modules (gaussian_pyramid doesn't seem ready
for inclusion in the (x)imgproc module).
Except for changes related to this, the patches
remain unchanged, including the outdated
copyright notice at the top of the files.

* Hide Harris-Laplace implementation behind pimpl-creating wrapper.

* Add test for Harris-Laplace feature detector.

* Handle empty images in Harris-Laplace detector.

* Replace HarrisAffineFeatureDetector with more flexible AffineFeature2D

* collapse unnecessary HarrisLaplace class into HarrisLaplaceFeatureDetector

* fold DoG pyramid code into harris_laplace_detector.cpp

* tuck auxiliary functions for Harris-Laplace into anonymous namespace

* use mask to filter keypoints in Harris-Laplace detector

* remove unused value of differentiation scale from calcAffineAdaptation

* make descriptor copy in AffineFeature2D independent of the data type

* Document interface for Harris-Laplace and AffineFeature2D

* Make type conversions explicit and decorate float literals with f in HarrisLaplace and AffineFeature2D

* Replace usage of non-standard log2 in Harris-Laplace feature detector

* Fix virtual overload errors in AffineFeature2D

* Add basic tests for AffineFeature2D and fix what they uncover

* Make type conversions in Harris-Laplace feature detector explicit

* Change license header for Harris-Laplace detector and AffineFeature2D

* Use Matx for small matrices in AffineFeature2D

* Remove redundant attributes of Elliptic_KeyPoint

* Convert Matx22f to Matx22d for inverting in AffineFeature2D

modules/xfeatures2d/doc/xfeatures2d.bib

@@ -111,3 +111,14 @@
   publisher = {{ACM}},
   year      = {2010}
 }
+
+@article{Mikolajczyk2004,
+  title={Scale \& affine invariant interest point detectors},
+  author={Mikolajczyk, Krystian and Schmid, Cordelia},
+  journal={International journal of computer vision},
+  volume={60},
+  number={1},
+  pages={63--86},
+  year={2004},
+  publisher={Springer}
+}

modules/xfeatures2d/include/opencv2/xfeatures2d.hpp

@@ -834,6 +834,97 @@ public:
 
 };
 
+/**
+* @brief Elliptic region around an interest point.
+*/
+class CV_EXPORTS Elliptic_KeyPoint : public KeyPoint
+{
+public:
+    Size_<float> axes; //!< the lengths of the major and minor ellipse axes
+    float si;  //!< the integration scale at which the parameters were estimated
+    Matx23f transf; //!< the transformation between image space and local patch space
+    Elliptic_KeyPoint();
+    Elliptic_KeyPoint(Point2f pt, float angle, Size axes, float size, float si);
+    virtual ~Elliptic_KeyPoint();
+};
+
+/**
+ * @brief Class implementing the Harris-Laplace feature detector as described in @cite Mikolajczyk2004.
+ */
+class CV_EXPORTS_W HarrisLaplaceFeatureDetector : public Feature2D
+{
+public:
+    /**
+     * @brief Creates a new implementation instance.
+     *
+     * @param numOctaves the number of octaves in the scale-space pyramid
+     * @param corn_thresh the threshold for the Harris cornerness measure
+     * @param DOG_thresh the threshold for the Difference-of-Gaussians scale selection
+     * @param maxCorners the maximum number of corners to consider
+     * @param num_layers the number of intermediate scales per octave
+     */
+    CV_WRAP static Ptr<HarrisLaplaceFeatureDetector> create(
+            int numOctaves=6,
+            float corn_thresh=0.01f,
+            float DOG_thresh=0.01f,
+            int maxCorners=5000,
+            int num_layers=4);
+};
+
+/**
+ * @brief Class implementing affine adaptation for key points.
+ *
+ * A @ref FeatureDetector and a @ref DescriptorExtractor are wrapped to augment the
+ * detected points with their affine invariant elliptic region and to compute
+ * the feature descriptors on the regions after warping them into circles.
+ *
+ * The interface is equivalent to @ref Feature2D, adding operations for
+ * @ref Elliptic_KeyPoint "Elliptic_KeyPoints" instead of @ref KeyPoint "KeyPoints".
+ */
+class CV_EXPORTS AffineFeature2D : public Feature2D
+{
+public:
+    /**
+     * @brief Creates an instance wrapping the given keypoint detector and
+     * descriptor extractor.
+     */
+    static Ptr<AffineFeature2D> create(
+        Ptr<FeatureDetector> keypoint_detector,
+        Ptr<DescriptorExtractor> descriptor_extractor);
+
+    /**
+     * @brief Creates an instance where keypoint detector and descriptor
+     * extractor are identical.
+     */
+    static Ptr<AffineFeature2D> create(
+        Ptr<FeatureDetector> keypoint_detector)
+    {
+        return create(keypoint_detector, keypoint_detector);
+    }
+
+    using Feature2D::detect; // overload, don't hide
+    /**
+     * @brief Detects keypoints in the image using the wrapped detector and
+     * performs affine adaptation to augment them with their elliptic regions.
+     */
+    virtual void detect(
+        InputArray image,
+        CV_OUT std::vector<Elliptic_KeyPoint>& keypoints,
+        InputArray mask=noArray() ) = 0;
+
+    using Feature2D::detectAndCompute; // overload, don't hide
+    /**
+     * @brief Detects keypoints and computes descriptors for their surrounding
+     * regions, after warping them into circles.
+     */
+    virtual void detectAndCompute(
+        InputArray image,
+        InputArray mask,
+        CV_OUT std::vector<Elliptic_KeyPoint>& keypoints,
+        OutputArray descriptors,
+        bool useProvidedKeypoints=false ) = 0;
+};
+
 //! @}
 
 }

modules/xfeatures2d/src/affine_feature2d.cpp

+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html
+
+/*
+ * Functions to perform affine adaptation of keypoint and to calculate descriptors of elliptic regions
+ */
+
+#include "precomp.hpp"
+#include "opencv2/imgproc.hpp"
+#include "opencv2/imgproc/types_c.h"
+
+namespace {
+
+using namespace cv;
+using namespace cv::xfeatures2d;
+
+/*
+* Functions to perform affine adaptation of circular keypoint
+*/
+void calcAffineCovariantRegions(const Mat& image, const std::vector<KeyPoint>& keypoints, std::vector<Elliptic_KeyPoint>& affRegions);
+void calcAffineCovariantDescriptors( const Ptr<DescriptorExtractor>& dextractor, const Mat& img, std::vector<Elliptic_KeyPoint>& affRegions, Mat& descriptors );
+
+void calcSecondMomentMatrix(const Mat & dx2, const Mat & dxy, const Mat & dy2, Point p, Matx22f& M);
+bool calcAffineAdaptation(const Mat & image, Elliptic_KeyPoint& keypoint);
+float selIntegrationScale(const Mat & image, float si, Point c);
+float selDifferentiationScale(const Mat & image, Mat & Lxm2smooth, Mat & Lxmysmooth, Mat & Lym2smooth, float si, Point c);
+float calcSecondMomentSqrt(const Mat & dx2, const Mat & dxy, const Mat & dy2, Point p, Matx22f& Mk);
+float normMaxEval(Matx22f & U, Mat& uVal, Mat& uVect);
+
+/*
+ * Calculates second moments matrix in point p
+ */
+void calcSecondMomentMatrix(const Mat & dx2, const Mat & dxy, const Mat & dy2, Point p, Matx22f & M)
+{
+    int x = p.x;
+    int y = p.y;
+
+    M(0, 0) = dx2.at<float> (y, x);
+    M(0, 1) = M(1, 0) = dxy.at<float> (y, x);
+    M(1, 1) = dy2.at<float> (y, x);
+}
+
+/*
+ * Performs affine adaptation
+ */
+bool calcAffineAdaptation(const Mat & fimage, Elliptic_KeyPoint & keypoint)
+{
+    Matx23f transf; /*Transformation matrix*/
+    Matx21f   size; /*Image size after transformation*/
+    Matx21f      c; /*Transformed point*/
+    Matx21f      p; /*Image point*/
+
+    Matx22f U(1.f, 0.f, 0.f, 1.f); /*Normalization matrix*/
+
+    Mat warpedImg, Lxm2smooth, Lym2smooth, Lxmysmooth, img_roi;
+    Matx22f Mk;
+    float Qinv = 1, q, si = keypoint.si;
+    bool divergence = false, convergence = false;
+    int i = 0;
+
+    //Coordinates in image
+    int py = (int) keypoint.pt.y;
+    int px = (int) keypoint.pt.x;
+
+    //Roi coordinates
+    int roix, roiy;
+
+    //Coordinates in U-trasformation
+    int cx = px;
+    int cy = py;
+    int cxPr = cx;
+    int cyPr = cy;
+
+    float radius = keypoint.size / 2 * 1.4f;
+    float half_width, half_height;
+
+    Rect roi;
+    float ax1, ax2;
+    float phi = 0;
+    ax1 = ax2 = keypoint.size / 2;
+    Mat drawImg;
+
+    //Affine adaptation
+    while (i <= 10 && !divergence && !convergence)
+    {
+        //Transformation matrix
+        transf = Matx23f(
+            U(0,0), U(0,1), 0.f,
+            U(1,0), U(1,1), 0.f
+        );
+        keypoint.transf = transf;
+
+        Size_<float> boundingBox;
+
+        float ac_b2 = float(determinant(U));
+        boundingBox.width  = ceil(U(1, 1)/ac_b2 * 3 * si*1.4f );
+        boundingBox.height = ceil(U(0, 0)/ac_b2 * 3 * si*1.4f );
+
+        //Create window around interest point
+        half_width = std::min((float) std::min(fimage.cols - px-1, px), boundingBox.width);
+        half_height = std::min((float) std::min(fimage.rows - py-1, py), boundingBox.height);
+        roix = max(px - (int) boundingBox.width, 0);
+        roiy = max(py - (int) boundingBox.height, 0);
+        roi = Rect(roix, roiy, px - roix + int(half_width)+1, py - roiy + int(half_height)+1);
+
+        //create ROI
+        img_roi = fimage(roi);
+
+
+        //Point within the ROI
+        p(0, 0) = float(px - roix);
+        p(1, 0) = float(py - roiy);
+
+        if (half_width <= 0 || half_height <= 0)
+            return divergence;
+
+        //Find coordinates of square's angles to find size of warped ellipse's bounding box
+        float u00 = U(0, 0);
+        float u01 = U(0, 1);
+        float u10 = U(1, 0);
+        float u11 = U(1, 1);
+
+        float minx = u01 * img_roi.rows < 0 ? u01 * img_roi.rows : 0;
+        float miny = u10 * img_roi.cols < 0 ? u10 * img_roi.cols : 0;
+        float maxx = (u00 * img_roi.cols > u00 * img_roi.cols + u01 * img_roi.rows ? u00
+                * img_roi.cols : u00 * img_roi.cols + u01 * img_roi.rows) - minx;
+        float maxy = (u11 * img_roi.rows > u10 * img_roi.cols + u11 * img_roi.rows ? u11
+                * img_roi.rows : u10 * img_roi.cols + u11 * img_roi.rows) - miny;
+
+        //Shift
+        transf(0, 2) = -minx;
+        transf(1, 2) = -miny;
+
+        /*float min_width = minx >= 0 ? u00 * img_roi.cols - u01 * img_roi.rows : u00 * img_roi.cols
+                + u01 * img_roi.rows;
+        float min_height = miny >= 0 ? u11 * img_roi.rows - u10 * img_roi.cols : u10 * img_roi.cols
+                + u11 * img_roi.rows;*/
+
+        if (maxx >=  2*radius+1 && maxy >=  2*radius+1)
+        {
+            //Size of normalized window must be 2*radius
+            //Transformation
+            Mat warpedImgRoi;
+            warpAffine(img_roi, warpedImgRoi, transf, Size(int(maxx), int(maxy)),INTER_AREA, BORDER_REPLICATE);
+
+            //Point in U-Normalized coordinates
+            c = U * p;
+            cx = int(c(0, 0) - minx);
+            cy = int(c(1, 0) - miny);
+
+            if (warpedImgRoi.rows > 2 * radius+1 && warpedImgRoi.cols > 2 * radius+1)
+            {
+                //Cut around normalized patch
+                roix = std::max(cx - int(ceil(radius)), 0);
+                roiy = std::max(cy - int(ceil(radius)), 0);
+                roi = Rect(roix, roiy,
+                        cx - roix + std::min(int(ceil(radius)), warpedImgRoi.cols - cx-1)+1,
+                        cy - roiy + std::min(int(ceil(radius)), warpedImgRoi.rows - cy-1)+1);
+                warpedImg = warpedImgRoi(roi);
+
+                //Coordinates in cutted ROI
+                cx = cx - roix;
+                cy = cy - roiy;
+            } else
+                warpedImgRoi.copyTo(warpedImg);
+
+            //Integration Scale selection
+            si = selIntegrationScale(warpedImg, si, Point(cx, cy));
+            //Differentation scale selection
+            selDifferentiationScale(warpedImg, Lxm2smooth, Lxmysmooth, Lym2smooth, si,
+                    Point(cx, cy));
+
+            //Spatial Localization
+            cxPr = cx; //Previous iteration point in normalized window
+            cyPr = cy;
+
+            float cornMax = 0;
+            for (int j = 0; j < 3; j++)
+            {
+                for (int t = 0; t < 3; t++)
+                {
+                    float dx2 = Lxm2smooth.at<float> (cyPr - 1 + j, cxPr - 1 + t);
+                    float dy2 = Lym2smooth.at<float> (cyPr - 1 + j, cxPr - 1 + t);
+                    float dxy = Lxmysmooth.at<float> (cyPr - 1 + j, cxPr - 1 + t);
+                    float det = dx2 * dy2 - dxy * dxy;
+                    float tr = dx2 + dy2;
+                    float cornerness = det - (0.04f * tr*tr);
+                    if (cornerness > cornMax)
+                    {
+                        cornMax = cornerness;
+                        cx = cxPr - 1 + t;
+                        cy = cyPr - 1 + j;
+                    }
+                }
+            }
+
+            //Transform point in image coordinates
+            p(0, 0) = float(px);
+            p(1, 0) = float(py);
+            //Displacement vector
+            c(0, 0) = float(cx - cxPr);
+            c(1, 0) = float(cy - cyPr);
+            //New interest point location in image
+            p = p + Matx22f(Matx22d(U).inv()) * c;
+            px = int(p(0, 0));
+            py = int(p(1, 0));
+
+            q = calcSecondMomentSqrt(Lxm2smooth, Lxmysmooth, Lym2smooth, Point(cx, cy), Mk);
+
+            float ratio = 1 - q;
+
+            //if ratio == 1 means q == 0 and one axes equals to 0
+            if (!isnan(ratio) && ratio != 1)
+            {
+                //Update U matrix
+                U = U * Mk;
+
+                Mat uVal, uV;
+                eigen(U, uVal, uV);
+
+                Qinv = normMaxEval(U, uVal, uV);
+
+                //Keypoint doesn't converge
+                if (Qinv >= 6)
+                    divergence = true;
+
+                //Keypoint converges
+                else if (ratio <= 0.05f)
+                {
+                    convergence = true;
+
+                    //Set transformation matrix
+                    transf = Matx23f(
+                        U(0,0), U(0,1), 0.f,
+                        U(1,0), U(1,1), 0.f
+                    );
+                    keypoint.transf = transf;
+
+                    ax1 = 1.f / std::abs(uVal.at<float> (0, 0)) * 3 * si;
+                    ax2 = 1.f / std::abs(uVal.at<float> (1, 0)) * 3 * si;
+                    phi = float(atan(uV.at<float> (1, 0) / uV.at<float> (0, 0)) * (180) / CV_PI);
+                    keypoint.axes = Size_<float> (ax1, ax2);
+                    keypoint.angle = phi;
+                    keypoint.pt = Point2f( (float) px, (float) py);
+                    keypoint.si = si;
+                    keypoint.size = 2 * 3 * si;
+
+                } else
+                    radius = 3 * si * 1.4f;
+
+            } else divergence = true;
+
+        } else divergence = true;
+
+        ++i;
+    }
+
+    return convergence;
+}
+
+/*
+ * Selects the integration scale that maximize LoG in point c
+ */
+float selIntegrationScale(const Mat & image, float si, Point c)
+{
+    Mat Lap, L;
+    int cx = c.x;
+    int cy = c.y;
+    float maxLap = 0;
+    float maxsx = si;
+    int gsize;
+    float sigma, sigma_prev = 0;
+
+    image.copyTo(L);
+    /* Search best integration scale between previous and successive layer
+     */
+    for (float u = 0.7f; u <= 1.41f; u += 0.1f)
+    {
+        float sik = u * si;
+        sigma = sqrt(powf(sik, 2) - powf(sigma_prev, 2));
+
+        gsize = int(ceil(sigma * 3)) * 2 + 1;
+
+        GaussianBlur(L, L, Size(gsize, gsize), sigma);
+        sigma_prev = sik;
+
+        Laplacian(L, Lap, CV_32F, 3);
+
+        float lapVal = sik * sik * std::abs(Lap.at<float> (cy, cx));
+
+        if (u == 0.7f)
+            maxLap = lapVal;
+
+        if (lapVal >= maxLap)
+        {
+            maxLap = lapVal;
+            maxsx = sik;
+        }
+    }
+    return maxsx;
+}
+
+/*
+ * Calculates second moments matrix square root
+ */
+float calcSecondMomentSqrt(const Mat & dx2, const Mat & dxy, const Mat & dy2, Point p, Matx22f & Mk)
+{
+    Mat V, eigVal, Vinv, D;
+    Matx22f M;
+
+    calcSecondMomentMatrix(dx2, dxy, dy2, p, M);
+
+    /* *
+     * M = V * D * V.inv()
+     * V has eigenvectors as columns
+     * D is a diagonal Matrix with eigenvalues as elements
+     * V.inv() is the inverse of V
+     * */
+
+    eigen(M, eigVal, V);
+    V = V.t();
+    Vinv = V.inv();
+
+    float eval1 = eigVal.at<float> (0, 0) = sqrt(eigVal.at<float> (0, 0));
+    float eval2 = eigVal.at<float> (1, 0) = sqrt(eigVal.at<float> (1, 0));
+
+    D = Mat::diag(eigVal);
+
+    //square root of M
+    Mk = Mat(V * D * Vinv);
+    //return q isotropic measure
+    return min(eval1, eval2) / max(eval1, eval2);
+}
+
+float normMaxEval(Matx22f & U, Mat & uVal, Mat & uVec)
+{
+    /* *
+     * Decomposition:
+     * U = V * D * V.inv()
+     * V has eigenvectors as columns
+     * D is a diagonal Matrix with eigenvalues as elements
+     * V.inv() is the inverse of V
+     * */
+    uVec = uVec.t();
+    Mat uVinv = uVec.inv();
+
+    //Normalize min eigenvalue to 1 to expand patch in the direction of min eigenvalue of U.inv()
+    float uval1 = uVal.at<float> (0, 0);
+    float uval2 = uVal.at<float> (1, 0);
+
+    if (std::abs(uval1) < std::abs(uval2))
+    {
+        uVal.at<float> (0, 0) = 1;
+        uVal.at<float> (1, 0) = uval2 / uval1;
+    } else
+    {
+        uVal.at<float> (1, 0) = 1;
+        uVal.at<float> (0, 0) = uval1 / uval2;
+    }
+
+    Mat D = Mat::diag(uVal);
+    //U normalized
+    U = Mat(uVec * D * uVinv);
+
+    return max(std::abs(uVal.at<float> (0, 0)), std::abs(uVal.at<float> (1, 0))) / min(
+            std::abs(uVal.at<float> (0, 0)), std::abs(uVal.at<float> (1, 0))); //define the direction of warping
+}
+
+/*
+ * Selects diffrentiation scale
+ */
+float selDifferentiationScale(const Mat & img, Mat & Lxm2smooth, Mat & Lxmysmooth,
+        Mat & Lym2smooth, float si, Point c)
+{
+    float s = 0.5f;
+    float sdk = s * si;
+    float sigma_prev = 0, sigma;
+
+    Mat L, dx2, dxy, dy2;
+
+    double qMax = 0;
+
+    //Gaussian kernel size
+    int gsize;
+    Size ksize;
+
+    img.copyTo(L);
+
+    while (s <= 0.751f)
+    {
+        Matx22f M;
+        float sd = s * si;
+
+        //Smooth previous smoothed image L
+        sigma = sqrt(powf(sd, 2) - powf(sigma_prev, 2));
+
+        gsize = int(ceil(sigma * 3)) * 2 + 1;
+
+        GaussianBlur(L, L, Size(gsize, gsize), sigma);
+
+        sigma_prev = sd;
+
+        //X and Y derivatives
+        Mat Lx, Ly;
+        Sobel(L, Lx, L.depth(), 1, 0, 1);
+        Lx = Lx * sd;
+        Sobel(L, Ly, L.depth(), 0, 1, 1);
+        Ly = Ly * sd;
+
+        //Size of gaussian kernel
+        gsize = int(ceil(si * 3)) * 2 + 1;
+        ksize = Size(gsize, gsize);
+
+        Mat Lxm2 = Lx.mul(Lx);
+        GaussianBlur(Lxm2, dx2, ksize, si);
+
+        Mat Lym2 = Ly.mul(Ly);
+        GaussianBlur(Lym2, dy2, ksize, si);
+
+        Mat Lxmy = Lx.mul(Ly);
+        GaussianBlur(Lxmy, dxy, ksize, si);
+
+        calcSecondMomentMatrix(dx2, dxy, dy2, Point(c.x, c.y), M);
+
+        //calc eigenvalues
+        Mat eval;
+        eigen(M, eval);
+        double eval1 = std::abs(eval.at<float> (0, 0));
+        double eval2 = std::abs(eval.at<float> (1, 0));
+        double q = min(eval1, eval2) / max(eval1, eval2);
+
+        if (q >= qMax)
+        {
+            qMax = q;
+            sdk = sd;
+            dx2.copyTo(Lxm2smooth);
+            dxy.copyTo(Lxmysmooth);
+            dy2.copyTo(Lym2smooth);
+
+        }
+        s += 0.05f;
+    }
+
+    return sdk;
+}
+
+void calcAffineCovariantRegions(const Mat & image, const std::vector<KeyPoint> & keypoints,
+        std::vector<Elliptic_KeyPoint> & affRegions)
+{
+    for (size_t i = 0; i < keypoints.size(); ++i)
+    {
+        KeyPoint kp = keypoints[i];
+        Elliptic_KeyPoint ex(kp.pt, 0, Size_<float> (kp.size / 2, kp.size / 2), kp.size,
+                kp.size / 6);
+
+        if (calcAffineAdaptation(image, ex))
+            affRegions.push_back(ex);
+    }
+    //Erase similar keypoint
+    float maxDiff = 4;
+    Mat colorimg;
+    for (size_t i = 0; i < affRegions.size(); i++)
+    {
+        Elliptic_KeyPoint kp1 = affRegions[i];
+        for (size_t j = i+1; j < affRegions.size(); j++){
+
+            Elliptic_KeyPoint kp2 = affRegions[j];
+
+            if(norm(kp1.pt-kp2.pt)<=maxDiff){
+                float phi1, phi2;
+                Size axes1, axes2;
+                float si1, si2;
+                phi1 = kp1.angle;
+                phi2 = kp2.angle;
+                axes1 = kp1.axes;
+                axes2 = kp2.axes;
+                si1 = kp1.si;
+                si2 = kp2.si;
+                if(std::abs(phi1-phi2)<15 && std::max(si1,si2)/std::min(si1,si2)<1.4f && axes1.width-axes2.width<5 && axes1.height-axes2.height<5){
+                    affRegions.erase(affRegions.begin()+j);
+                    j--;
+                }
+            }
+        }
+    }
+}
+
+void calcAffineCovariantDescriptors(const Ptr<DescriptorExtractor>& dextractor, const Mat& img,
+        std::vector<Elliptic_KeyPoint>& affRegions, Mat& descriptors)
+{
+
+    assert(!affRegions.empty());
+    int descriptorSize = dextractor->descriptorSize();
+    int descriptorType = dextractor->descriptorType();
+    descriptors.create(Size(descriptorSize, int(affRegions.size())), descriptorType);
+    descriptors.setTo(0);
+
+    int i = 0;
+
+    for (std::vector<Elliptic_KeyPoint>::iterator it = affRegions.begin(); it < affRegions.end(); ++it)
+    {
+        Point p = it->pt;
+
+        Matx21f size;
+        size(0, 0) = size(1, 0) = it->size;
+
+        //U matrix
+        Matx23f transf = it->transf;
+        Matx22f U(
+            transf(0,0), transf(0,1),
+            transf(1,0), transf(1,1)
+        );
+
+        float radius = it->size / 2;
+        float si = it->si;
+
+        Size_<float> boundingBox;
+
+        float ac_b2 = float(determinant(U));
+        boundingBox.width  = ceil(U(1, 1)/ac_b2 * 3 * si );
+        boundingBox.height = ceil(U(0, 0)/ac_b2 * 3 * si );
+
+        //Create window around interest point
+        float half_width = std::min((float) std::min(img.cols - p.x-1, p.x), boundingBox.width);
+        float half_height = std::min((float) std::min(img.rows - p.y-1, p.y), boundingBox.height);
+        int roix = max(p.x - (int) boundingBox.width, 0);
+        int roiy = max(p.y - (int) boundingBox.height, 0);
+        Rect roi = Rect(roix, roiy, p.x - roix + int(half_width)+1, p.y - roiy + int(half_height)+1);
+
+        Mat img_roi = img(roi);
+
+        size(0, 0) = float(img_roi.cols);
+        size(1, 0) = float(img_roi.rows);
+
+        size = U * size;
+
+        Mat transfImgRoi, transfImg;
+        warpAffine(img_roi, transfImgRoi, transf, Size(int(ceil(size(0, 0))), int(ceil(size(1, 0)))),
+                INTER_AREA, BORDER_DEFAULT);
+
+        Matx21f c; //Transformed point
+        Matx21f pt; //Image point
+        //Point within the Roi
+        pt(0, 0) = float(p.x - roix);
+        pt(1, 0) = float(p.y - roiy);
+
+        //Point in U-Normalized coordinates
+        c = U * pt;
+        float cx = c(0, 0);
+        float cy = c(1, 0);
+
+        //Cut around point to have patch of 2*keypoint->size
+
+        roix = std::max(int(ceil(cx - radius)), 0);
+        roiy = std::max(int(ceil(cy - radius)), 0);
+
+        roi = Rect(roix, roiy, int(ceil(std::min(cx - roix + radius, size(0, 0)))),
+                int(ceil(std::min(cy - roiy + radius, size(1, 0)))));
+        transfImg = transfImgRoi(roi);
+
+        cx = c(0, 0) - roix;
+        cy = c(1, 0) - roiy;
+
+        Mat tmpDesc;
+        KeyPoint kp(Point(int(cx), int(cy)), it->size);
+
+        std::vector<KeyPoint> k(1, kp);
+
+        transfImg.convertTo(transfImg, CV_8U);
+        dextractor->compute(transfImg, k, tmpDesc);
+
+        tmpDesc.row(0).copyTo(descriptors.row(i));
+
+        i++;
+
+    }
+
+}
+
+} // anonymous namespace
+
+namespace cv
+{
+namespace xfeatures2d
+{
+class AffineFeature2D_Impl : public AffineFeature2D
+{
+public:
+    AffineFeature2D_Impl(
+        Ptr<FeatureDetector> keypoint_detector,
+        Ptr<DescriptorExtractor> descriptor_extractor
+    ) : m_keypoint_detector(keypoint_detector)
+      , m_descriptor_extractor(descriptor_extractor) {}
+protected:
+    using Feature2D::detect; // overload, don't hide
+    void detect(InputArray image, std::vector<Elliptic_KeyPoint>& keypoints, InputArray mask);
+    void detectAndCompute(InputArray image, InputArray mask, std::vector<Elliptic_KeyPoint>& keypoints, OutputArray descriptors, bool useProvidedKeypoints);
+    void detectAndCompute(InputArray image, InputArray mask, std::vector<KeyPoint>& keypoints, OutputArray descriptors, bool useProvidedKeypoints);
+    int descriptorSize() const;
+    int descriptorType() const;
+    int defaultNorm() const;
+private:
+    Ptr<FeatureDetector> m_keypoint_detector;
+    Ptr<DescriptorExtractor> m_descriptor_extractor;
+};
+
+Ptr<AffineFeature2D> AffineFeature2D::create(
+    Ptr<FeatureDetector> keypoint_detector,
+    Ptr<DescriptorExtractor> descriptor_extractor)
+{
+    return makePtr<AffineFeature2D_Impl>(keypoint_detector, descriptor_extractor);
+}
+
+void AffineFeature2D_Impl::detect(
+    InputArray image,
+    std::vector<Elliptic_KeyPoint>& keypoints,
+    InputArray mask)
+{
+    std::vector<KeyPoint> non_elliptic_keypoints;
+    m_keypoint_detector->detect(image, non_elliptic_keypoints, mask);
+    Mat fimage;
+    image.getMat().convertTo(fimage, CV_32F, 1.f/255);
+    calcAffineCovariantRegions(fimage, non_elliptic_keypoints, keypoints);
+}
+
+void AffineFeature2D_Impl::detectAndCompute(
+        InputArray image,
+        InputArray mask,
+        std::vector<Elliptic_KeyPoint>& keypoints,
+        OutputArray descriptors,
+        bool useProvidedKeypoints)
+{
+    if(!useProvidedKeypoints)
+    {
+        std::vector<KeyPoint> non_elliptic_keypoints;
+        m_keypoint_detector->detect(image, non_elliptic_keypoints, mask);
+        Mat fimage;
+        image.getMat().convertTo(fimage, CV_32F, 1.f/255);
+        calcAffineCovariantRegions(fimage, non_elliptic_keypoints, keypoints);
+    }
+    if(descriptors.needed())calcAffineCovariantDescriptors(m_descriptor_extractor, image.getMat(), keypoints, descriptors.getMatRef());
+}
+
+void AffineFeature2D_Impl::detectAndCompute(
+        InputArray image,
+        InputArray mask,
+        std::vector<KeyPoint>& keypoints,
+        OutputArray descriptors,
+        bool useProvidedKeypoints)
+{
+    if(!useProvidedKeypoints)
+    {
+        m_keypoint_detector->detect(image, keypoints, mask);
+    }
+    if(descriptors.needed()) {
+        Mat fimage;
+        image.getMat().convertTo(fimage, CV_32F, 1.f/255);
+        std::vector<Elliptic_KeyPoint> elliptic_keypoints;
+        calcAffineCovariantRegions(fimage, keypoints, elliptic_keypoints);
+        calcAffineCovariantDescriptors(m_descriptor_extractor, image.getMat(), elliptic_keypoints, descriptors.getMatRef());
+    }
+}
+
+int AffineFeature2D_Impl::descriptorSize() const
+{
+    return m_descriptor_extractor->descriptorSize();
+}
+
+int AffineFeature2D_Impl::descriptorType() const
+{
+    return m_descriptor_extractor->descriptorType();
+}
+
+int AffineFeature2D_Impl::defaultNorm() const
+{
+    return m_descriptor_extractor->defaultNorm();
+}
+
+}
+}

modules/xfeatures2d/src/ellipticKeyPoint.cpp

+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html
+
+#include "precomp.hpp"
+
+namespace cv {
+namespace xfeatures2d {
+
+Elliptic_KeyPoint::Elliptic_KeyPoint(Point2f _pt, float _angle, Size _axes, float _size, float _si) :
+    KeyPoint(_pt,_size,_angle), axes(_axes), si(_si) {
+}
+
+Elliptic_KeyPoint::Elliptic_KeyPoint(){
+
+}
+
+Elliptic_KeyPoint::~Elliptic_KeyPoint() {
+}
+
+}
+}

modules/xfeatures2d/src/harris_lapace_detector.cpp

+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html
+
+#include "precomp.hpp"
+
+namespace {
+
+using namespace cv;
+
+bool sort_func(KeyPoint kp1, KeyPoint kp2)
+{
+    return (kp1.response > kp2.response);
+}
+
+class Pyramid
+{
+
+protected:
+    class Octave
+    {
+    public:
+        std::vector<Mat> layers;
+        Octave();
+        Octave(std::vector<Mat> layers);
+        virtual ~Octave();
+        std::vector<Mat> getLayers();
+        Mat getLayerAt(int i);
+    };
+
+    class DOGOctave
+    {
+    public:
+        std::vector<Mat> layers;
+
+        DOGOctave();
+        DOGOctave(std::vector<Mat> layers);
+        virtual ~DOGOctave();
+        std::vector<Mat> getLayers();
+        Mat getLayerAt(int i);
+    };
+
+private:
+    std::vector<Octave> octaves;
+    std::vector<DOGOctave> DOG_octaves;
+    void build(const Mat& img, bool DOG);
+public:
+    class Params
+    {
+    public:
+        int octavesN;
+        int layersN;
+        float sigma0;
+        int omin;
+        float step;
+        Params();
+        Params(int octavesN, int layersN, float sigma0, int omin);
+        void clear();
+    };
+    Params params;
+
+    Pyramid();
+    Pyramid(const Mat& img, int octavesN, int layersN = 2, float sigma0 = 1, int omin = 0,
+            bool DOG = false);
+    Mat getLayer(int octave, int layer);
+    Mat getDOGLayer(int octave, int layer);
+    float getSigma(int octave, int layer);
+    float getSigma(int layer);
+
+    virtual ~Pyramid();
+    Params getParams();
+    void clear();
+    bool empty();
+};
+
+Pyramid::Pyramid()
+{
+}
+
+/**
+ * Pyramid class constructor
+ * octavesN_: number of octaves
+ * layersN_: number of layers before subsampling layer
+ * sigma0_: starting sigma (depends on detector's type, i.e. SIFT sigma0 = 1.6, Harris sigma0 = 1)
+ * omin_: if omin<0 an octave is added before first octave. In this octave the image size is doubled
+ * _DOG: if true, a DOG pyramid is build
+ */
+Pyramid::Pyramid(const Mat & img, int octavesN_, int layersN_, float sigma0_, int omin_, bool _DOG) :
+    params(octavesN_, layersN_, sigma0_, omin_)
+{
+    build(img, _DOG);
+}
+
+/**
+ * Build gaussian pyramid with layersN_ + 3 layers and 2^(1/layersN_) step between layers
+ * each octave is downsampled of a factor of 2
+ */
+
+void Pyramid::build(const Mat& img, bool DOG)
+{
+
+    Size ksize(0, 0);
+    int gsize;
+
+    Size imgSize = img.size();
+    int minSize = MIN(imgSize.width, imgSize.height);
+    int octavesN = MIN(params.octavesN, int(floor(log((double) minSize)/log(2))));
+    float sigma0 = params.sigma0;
+    float sigma = sigma0;
+    int layersN = params.layersN + 3;
+    int omin = params.omin;
+    float k = params.step;
+
+    /*layer to downsample*/
+    int down_lay = int(1 / log(k));
+
+    int octave, layer;
+    double sigmaN = 0.5;
+
+    std::vector<Mat> layers, DOG_layers;
+    /* standard deviation of current layer*/
+    float sigma_curr = sigma;
+    /* standard deviation of previous layer*/
+    float sigma_prev = sigma;
+
+    if (omin < 0)
+    {
+        omin = -1;
+        Mat tmp_img;
+        Mat blurred_img;
+        gsize = int(ceil(sigmaN * 3)) * 2 + 1;
+        GaussianBlur(img, blurred_img, Size(gsize,gsize), sigmaN);
+        resize(blurred_img, tmp_img, ksize, 2, 2, INTER_AREA);
+        layers.push_back(tmp_img);
+
+        for (layer = 1; layer < layersN; layer++)
+        {
+            sigma_curr = getSigma(layer);
+            sigma = sqrt(powf(sigma_curr, 2) - powf(sigma_prev, 2));
+            Mat prev_lay = layers[layer - 1], curr_lay, DOG_lay;
+            /* smoothing is applied on previous layer so sigma_curr^2 = sigma^2 + sigma_prev^2 */
+            gsize = int(ceil(sigma * 3)) * 2 + 1;
+            GaussianBlur(prev_lay, curr_lay, Size(gsize,gsize), sigma);
+            layers.push_back(curr_lay);
+            if (DOG)
+            {
+                absdiff(curr_lay, prev_lay, DOG_lay);
+                DOG_layers.push_back(DOG_lay);
+            }
+            sigma_prev = sigma_curr;
+
+        }
+        Octave tmp_oct(layers);
+        octaves.push_back(tmp_oct);
+        layers.clear();
+
+        if (DOG)
+        {
+            DOGOctave tmp_DOG_Oct(DOG_layers);
+            DOG_octaves.push_back(tmp_DOG_Oct);
+            DOG_layers.clear();
+        }
+
+    }
+
+    /* Presmoothing on first layer */
+    float sb = float(sigmaN) / powf(2.0f, (float) omin);
+    sigma = sigma0;
+    if (sigma0 > sb)
+        sigma = sqrt(sigma0 * sigma0 - sb * sb);
+
+    /*1° step on image*/
+    Mat tmpImg;
+    gsize = int(ceil(sigma * 3)) * 2 + 1;
+    GaussianBlur(img, tmpImg, Size(gsize,gsize), sigma);
+    layers.push_back(tmpImg);
+
+    /*for every octave build layers*/
+    sigma_prev = sigma;
+
+    for (octave = 0; octave < octavesN; octave++)
+    {
+        for (layer = 1; layer < layersN; layer++)
+        {
+            sigma_curr = getSigma(layer);
+            sigma = sqrt(powf(sigma_curr, 2) - powf(sigma_prev, 2));
+
+            Mat prev_lay = layers[layer - 1], curr_lay, DOG_lay;
+            gsize = int(ceil(sigma * 3)) * 2 + 1;
+            GaussianBlur(prev_lay, curr_lay, Size(gsize,gsize), sigma);
+            layers.push_back(curr_lay);
+
+            if (DOG)
+            {
+                absdiff(curr_lay, prev_lay, DOG_lay);
+                DOG_layers.push_back(DOG_lay);
+            }
+            sigma_prev = sigma_curr;
+        }
+
+        Mat resized_lay;
+        resize(layers[down_lay], resized_lay, ksize, 1.0f / 2, 1.0f / 2, INTER_AREA);
+
+        Octave tmp_oct(layers);
+        octaves.push_back(tmp_oct);
+        if (DOG)
+        {
+            DOGOctave tmp_DOG_Oct(DOG_layers);
+            DOG_octaves.push_back(tmp_DOG_Oct);
+            DOG_layers.clear();
+        }
+        sigma_curr = sigma_prev = sigma0;
+        layers.clear();
+        layers.push_back(resized_lay);
+
+    }
+
+}
+
+/**
+ * Return layer at indicated octave and layer numbers
+ */
+Mat Pyramid::getLayer(int octave, int layer)
+{
+    return octaves[octave].getLayerAt(layer);
+}
+
+/**
+ * Return DOG layer at indicated octave and layer numbers
+ */
+Mat Pyramid::getDOGLayer(int octave, int layer)
+{
+    CV_Assert(!DOG_octaves.empty());
+    return DOG_octaves[octave].getLayerAt(layer);
+}
+
+/**
+ * Return sigma value of indicated octave and layer
+ */
+float Pyramid::getSigma(int octave, int layer)
+{
+
+    return powf(2.0f, float(octave)) * powf(params.step, float(layer)) * params.sigma0;
+}
+
+/**
+ * Return sigma value of indicated layer
+ * sigma value of layer is the same at each octave
+ * i.e. sigma of first layer at each octave is sigma0
+ */
+float Pyramid::getSigma(int layer)
+{
+
+    return powf(params.step, float(layer)) * params.sigma0;
+}
+
+/**
+ * Destructor of Pyramid class
+ */
+Pyramid::~Pyramid()
+{
+    clear();
+}
+
+/**
+ * Clear octaves and params
+ */
+void Pyramid::clear()
+{
+    octaves.clear();
+    params.clear();
+}
+
+/**
+ * Empty Pyramid
+ * @return
+ */
+bool Pyramid::empty()
+{
+    return octaves.empty();
+}
+
+Pyramid::Params::Params()
+{
+}
+
+/**
+ * Params for Pyramid class
+ *
+ */
+Pyramid::Params::Params(int octavesN_, int layersN_, float sigma0_, int omin_) :
+    octavesN(octavesN_), layersN(layersN_), sigma0(sigma0_), omin(omin_)
+{
+    CV_Assert(layersN > 0 && octavesN_>0);
+    step = powf(2, 1.0f / layersN);
+}
+
+/**
+ * Returns Pyramid's params
+ */
+Pyramid::Params Pyramid::getParams()
+{
+    return params;
+}
+
+/**
+ * Set to zero all params
+ */
+void Pyramid::Params::clear()
+{
+    octavesN = 0;
+    layersN = 0;
+    sigma0 = 0;
+    omin = 0;
+    step = 0;
+}
+
+/**
+ * Create an Octave with layers
+ */
+Pyramid::Octave::Octave(std::vector<Mat> _layers) : layers(_layers) {}
+
+/**
+ * Return layers of the Octave
+ */
+std::vector<Mat> Pyramid::Octave::getLayers()
+{
+    return layers;
+}
+
+Pyramid::Octave::Octave()
+{
+}
+
+/**
+ * Return the Octave's layer at index i
+ */
+Mat Pyramid::Octave::getLayerAt(int i)
+{
+    CV_Assert(i < (int) layers.size());
+    return layers[i];
+}
+
+Pyramid::Octave::~Octave()
+{
+}
+
+Pyramid::DOGOctave::DOGOctave()
+{
+}
+
+Pyramid::DOGOctave::DOGOctave(std::vector<Mat> _layers) : layers(_layers) {}
+
+Pyramid::DOGOctave::~DOGOctave()
+{
+}
+
+std::vector<Mat> Pyramid::DOGOctave::getLayers()
+{
+    return layers;
+}
+
+Mat Pyramid::DOGOctave::getLayerAt(int i)
+{
+    CV_Assert(i < (int) layers.size());
+    return layers[i];
+}
+
+} // anonymous namespace
+
+namespace cv
+{
+namespace xfeatures2d
+{
+
+/*
+ *  HarrisLaplaceFeatureDetector_Impl
+ */
+class HarrisLaplaceFeatureDetector_Impl : public HarrisLaplaceFeatureDetector
+{
+public:
+    HarrisLaplaceFeatureDetector_Impl(
+        int numOctaves=6,
+        float corn_thresh=0.01,
+        float DOG_thresh=0.01,
+        int maxCorners=5000,
+        int num_layers=4
+    );
+    virtual void read( const FileNode& fn );
+    virtual void write( FileStorage& fs ) const;
+
+protected:
+    void detect( InputArray image, std::vector<KeyPoint>& keypoints, InputArray mask=noArray() );
+
+    int numOctaves;
+    float corn_thresh;
+    float DOG_thresh;
+    int maxCorners;
+    int num_layers;
+};
+
+Ptr<HarrisLaplaceFeatureDetector> HarrisLaplaceFeatureDetector::create(
+    int numOctaves,
+    float corn_thresh,
+    float DOG_thresh,
+    int maxCorners,
+    int num_layers)
+{
+    return makePtr<HarrisLaplaceFeatureDetector_Impl>(numOctaves, corn_thresh, DOG_thresh, maxCorners, num_layers);
+}
+
+HarrisLaplaceFeatureDetector_Impl::HarrisLaplaceFeatureDetector_Impl(
+    int _numOctaves,
+    float _corn_thresh,
+    float _DOG_thresh,
+    int _maxCorners,
+    int _num_layers
+) :
+    numOctaves(_numOctaves),
+    corn_thresh(_corn_thresh),
+    DOG_thresh(_DOG_thresh),
+    maxCorners(_maxCorners),
+    num_layers(_num_layers)
+{
+    CV_Assert(num_layers == 2 || num_layers==4);
+}
+
+void HarrisLaplaceFeatureDetector_Impl::read (const FileNode& fn)
+{
+    numOctaves = fn["numOctaves"];
+    corn_thresh = fn["corn_thresh"];
+    DOG_thresh = fn["DOG_thresh"];
+    maxCorners = fn["maxCorners"];
+    num_layers = fn["num_layers"];
+}
+
+void HarrisLaplaceFeatureDetector_Impl::write (FileStorage& fs) const
+{
+    fs << "numOctaves" << numOctaves;
+    fs << "corn_thresh" << corn_thresh;
+    fs << "DOG_thresh" << DOG_thresh;
+    fs << "maxCorners" << maxCorners;
+    fs << "num_layers" << num_layers;
+}
+
+/*
+ * Detect method
+ * The method detect Harris corners on scale space as described in
+ * "K. Mikolajczyk and C. Schmid.
+ * Scale & affine invariant interest point detectors.
+ * International Journal of Computer Vision, 2004"
+ */
+void HarrisLaplaceFeatureDetector_Impl::detect(InputArray img, std::vector<KeyPoint>& keypoints, InputArray msk )
+{
+    Mat image = img.getMat();
+    if( image.empty() )
+    {
+        keypoints.clear();
+        return;
+    }
+    Mat mask = msk.getMat();
+    if( !mask.empty() )
+    {
+        CV_Assert(mask.type() == CV_8UC1);
+        CV_Assert(mask.size == image.size);
+    }
+    Mat_<float> dx2, dy2, dxy;
+    Mat Lx, Ly;
+    float si, sd;
+    int gsize;
+    Mat fimage;
+    image.convertTo(fimage, CV_32F, 1.f/255);
+    /*Build gaussian pyramid*/
+    Pyramid pyr(fimage, numOctaves, num_layers, 1, -1, true);
+    keypoints = std::vector<KeyPoint> (0);
+
+    /*Find Harris corners on each layer*/
+    for (int octave = 0; octave <= numOctaves; octave++)
+    {
+        for (int layer = 1; layer <= num_layers; layer++)
+        {
+            if (octave == 0)
+                layer = num_layers;
+
+            Mat Lxm2smooth, Lxmysmooth, Lym2smooth;
+
+            si = powf(2.f, layer / (float) num_layers);
+            sd = si * 0.7f;
+
+            Mat curr_layer;
+            if (num_layers == 4)
+            {
+                if (layer == 1)
+                {
+                    Mat tmp = pyr.getLayer(octave - 1, num_layers - 1);
+                    resize(tmp, curr_layer, Size(0, 0), 0.5, 0.5, INTER_AREA);
+
+                } else
+                    curr_layer = pyr.getLayer(octave, layer - 2);
+            } else /*if num_layer==2*/
+            {
+
+                curr_layer = pyr.getLayer(octave, layer - 1);
+            }
+
+            /*Calculates second moment matrix*/
+
+            /*Derivatives*/
+            Sobel(curr_layer, Lx, CV_32F, 1, 0, 1);
+            Sobel(curr_layer, Ly, CV_32F, 0, 1, 1);
+
+            /*Normalization*/
+            Lx = Lx * sd;
+            Ly = Ly * sd;
+
+            Mat Lxm2 = Lx.mul(Lx);
+            Mat Lym2 = Ly.mul(Ly);
+            Mat Lxmy = Lx.mul(Ly);
+
+            gsize = int(ceil(si * 3)) * 2 + 1;
+
+            /*Convolution*/
+            GaussianBlur(Lxm2, Lxm2smooth, Size(gsize, gsize), si, si, BORDER_REPLICATE);
+            GaussianBlur(Lym2, Lym2smooth, Size(gsize, gsize), si, si, BORDER_REPLICATE);
+            GaussianBlur(Lxmy, Lxmysmooth, Size(gsize, gsize), si, si, BORDER_REPLICATE);
+
+            Mat cornern_mat(curr_layer.size(), CV_32F);
+
+            /*Calculates cornerness in each pixel of the image*/
+            for (int row = 0; row < curr_layer.rows; row++)
+            {
+                for (int col = 0; col < curr_layer.cols; col++)
+                {
+                    float dx2f = Lxm2smooth.at<float> (row, col);
+                    float dy2f = Lym2smooth.at<float> (row, col);
+                    float dxyf = Lxmysmooth.at<float> (row, col);
+                    float det = dx2f * dy2f - dxyf * dxyf;
+                    float tr = dx2f + dy2f;
+                    float cornerness = det - (0.04f * tr * tr);
+                    cornern_mat.at<float> (row, col) = cornerness;
+                }
+            }
+
+            double maxVal = 0;
+            Mat corn_dilate;
+
+            /*Find max cornerness value and rejects all corners that are lower than a threshold*/
+            minMaxLoc(cornern_mat, 0, &maxVal, 0, 0);
+            threshold(cornern_mat, cornern_mat, maxVal * corn_thresh, 0, THRESH_TOZERO);
+            dilate(cornern_mat, corn_dilate, Mat());
+
+            Size imgsize = curr_layer.size();
+
+            /*Verify for each of the initial points whether the DoG attains a maximum at the scale of the point*/
+            Mat prevDOG, curDOG, succDOG;
+            prevDOG = pyr.getDOGLayer(octave, layer - 1);
+            curDOG = pyr.getDOGLayer(octave, layer);
+            succDOG = pyr.getDOGLayer(octave, layer + 1);
+
+            for (int y = 1; y < imgsize.height - 1; y++)
+            {
+                for (int x = 1; x < imgsize.width - 1; x++)
+                {
+                    float val = cornern_mat.at<float> (y, x);
+                    if (val != 0 && val == corn_dilate.at<float> (y, x))
+                    {
+
+                        float curVal = curDOG.at<float> (y, x);
+                        float prevVal =  prevDOG.at<float> (y, x);
+                        float succVal = succDOG.at<float> (y, x);
+
+                        KeyPoint kp(
+                                Point2f(x * powf(2.0f, (float) octave - 1) + powf(2.0f, (float) octave - 1) / 2,
+                                        y * powf(2.0f, (float) octave - 1) + powf(2.0f, (float) octave - 1) / 2),
+                                3 * powf(2.0f, (float) octave - 1) * si * 2, 0, val, octave);
+
+                        if(!mask.empty() && mask.at<unsigned char>(int(kp.pt.y), int(kp.pt.x)) == 0)
+                        {
+                            // ignore keypoints where mask is zero
+                            continue;
+                        }
+
+                        /*Check whether keypoint size is inside the image*/
+                        float start_kp_x = kp.pt.x - kp.size / 2;
+                        float start_kp_y = kp.pt.y - kp.size / 2;
+                        float end_kp_x = start_kp_x + kp.size;
+                        float end_kp_y = start_kp_y + kp.size;
+
+                        if (curVal > prevVal && curVal > succVal && curVal >= DOG_thresh
+                                && start_kp_x > 0 && start_kp_y > 0 && end_kp_x < image.cols
+                                && end_kp_y < image.rows)
+                            keypoints.push_back(kp);
+
+                    }
+                }
+            }
+
+        }
+
+    }
+
+    /*Sort keypoints in decreasing cornerness order*/
+    sort(keypoints.begin(), keypoints.end(), sort_func);
+    for (size_t i = 1; i < keypoints.size(); i++)
+    {
+        float max_diff = powf(2, keypoints[i].octave + 1.f / 2);
+
+        if (keypoints[i].response == keypoints[i - 1].response && norm(
+                keypoints[i].pt - keypoints[i - 1].pt) <= max_diff)
+        {
+
+            float x = (keypoints[i].pt.x + keypoints[i - 1].pt.x) / 2;
+            float y = (keypoints[i].pt.y + keypoints[i - 1].pt.y) / 2;
+
+            keypoints[i].pt = Point2f(x, y);
+            --i;
+            keypoints.erase(keypoints.begin() + i);
+
+        }
+    }
+
+    /*Select strongest keypoints*/
+    if (maxCorners > 0 && maxCorners < (int) keypoints.size())
+        keypoints.resize(maxCorners);
+}
+
+
+}
+}

modules/xfeatures2d/test/test_features2d.cpp

@@ -997,6 +997,24 @@ TEST( Features2d_Detector_STAR, regression )
     test.safe_run();
 }
 
+TEST( Features2d_Detector_Harris_Laplace, regression )
+{
+    CV_FeatureDetectorTest test( "detector-harris-laplace", HarrisLaplaceFeatureDetector::create() );
+    test.safe_run();
+}
+
+TEST( Features2d_Detector_Harris_Laplace_Affine_Keypoint_Invariance, regression )
+{
+    CV_FeatureDetectorTest test( "detector-harris-laplace", AffineFeature2D::create(HarrisLaplaceFeatureDetector::create()));
+    test.safe_run();
+}
+
+TEST( Features2d_Detector_Harris_Laplace_Affine, regression )
+{
+    CV_FeatureDetectorTest test( "detector-harris-laplace-affine", AffineFeature2D::create(HarrisLaplaceFeatureDetector::create()));
+    test.safe_run();
+}
+
 /*
  * Descriptors
  */