diff --git a/modules/xfeatures2d/doc/xfeatures2d.bib b/modules/xfeatures2d/doc/xfeatures2d.bib index 93d5b804c..9222663a5 100644 --- a/modules/xfeatures2d/doc/xfeatures2d.bib +++ b/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} +} diff --git a/modules/xfeatures2d/include/opencv2/xfeatures2d.hpp b/modules/xfeatures2d/include/opencv2/xfeatures2d.hpp index 1185de033..93c16feff 100644 --- a/modules/xfeatures2d/include/opencv2/xfeatures2d.hpp +++ b/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_ 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 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 create( + Ptr keypoint_detector, + Ptr descriptor_extractor); + + /** + * @brief Creates an instance where keypoint detector and descriptor + * extractor are identical. + */ + static Ptr create( + Ptr 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& 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& keypoints, + OutputArray descriptors, + bool useProvidedKeypoints=false ) = 0; +}; + //! @} } diff --git a/modules/xfeatures2d/src/affine_feature2d.cpp b/modules/xfeatures2d/src/affine_feature2d.cpp new file mode 100644 index 000000000..60f43bf57 --- /dev/null +++ b/modules/xfeatures2d/src/affine_feature2d.cpp @@ -0,0 +1,681 @@ +// 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& keypoints, std::vector& affRegions); +void calcAffineCovariantDescriptors( const Ptr& dextractor, const Mat& img, std::vector& 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 (y, x); + M(0, 1) = M(1, 0) = dxy.at (y, x); + M(1, 1) = dy2.at (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_ 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 (cyPr - 1 + j, cxPr - 1 + t); + float dy2 = Lym2smooth.at (cyPr - 1 + j, cxPr - 1 + t); + float dxy = Lxmysmooth.at (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 (0, 0)) * 3 * si; + ax2 = 1.f / std::abs(uVal.at (1, 0)) * 3 * si; + phi = float(atan(uV.at (1, 0) / uV.at (0, 0)) * (180) / CV_PI); + keypoint.axes = Size_ (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 (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 (0, 0) = sqrt(eigVal.at (0, 0)); + float eval2 = eigVal.at (1, 0) = sqrt(eigVal.at (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 (0, 0); + float uval2 = uVal.at (1, 0); + + if (std::abs(uval1) < std::abs(uval2)) + { + uVal.at (0, 0) = 1; + uVal.at (1, 0) = uval2 / uval1; + } else + { + uVal.at (1, 0) = 1; + uVal.at (0, 0) = uval1 / uval2; + } + + Mat D = Mat::diag(uVal); + //U normalized + U = Mat(uVec * D * uVinv); + + return max(std::abs(uVal.at (0, 0)), std::abs(uVal.at (1, 0))) / min( + std::abs(uVal.at (0, 0)), std::abs(uVal.at (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 (0, 0)); + double eval2 = std::abs(eval.at (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 & keypoints, + std::vector & affRegions) +{ + for (size_t i = 0; i < keypoints.size(); ++i) + { + KeyPoint kp = keypoints[i]; + Elliptic_KeyPoint ex(kp.pt, 0, Size_ (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& dextractor, const Mat& img, + std::vector& 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::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_ 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 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 keypoint_detector, + Ptr 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& keypoints, InputArray mask); + void detectAndCompute(InputArray image, InputArray mask, std::vector& keypoints, OutputArray descriptors, bool useProvidedKeypoints); + void detectAndCompute(InputArray image, InputArray mask, std::vector& keypoints, OutputArray descriptors, bool useProvidedKeypoints); + int descriptorSize() const; + int descriptorType() const; + int defaultNorm() const; +private: + Ptr m_keypoint_detector; + Ptr m_descriptor_extractor; +}; + +Ptr AffineFeature2D::create( + Ptr keypoint_detector, + Ptr descriptor_extractor) +{ + return makePtr(keypoint_detector, descriptor_extractor); +} + +void AffineFeature2D_Impl::detect( + InputArray image, + std::vector& keypoints, + InputArray mask) +{ + std::vector 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& keypoints, + OutputArray descriptors, + bool useProvidedKeypoints) +{ + if(!useProvidedKeypoints) + { + std::vector 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& 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_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(); +} + +} +} diff --git a/modules/xfeatures2d/src/ellipticKeyPoint.cpp b/modules/xfeatures2d/src/ellipticKeyPoint.cpp new file mode 100644 index 000000000..d71b00266 --- /dev/null +++ b/modules/xfeatures2d/src/ellipticKeyPoint.cpp @@ -0,0 +1,22 @@ +// 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() { +} + +} +} diff --git a/modules/xfeatures2d/src/harris_lapace_detector.cpp b/modules/xfeatures2d/src/harris_lapace_detector.cpp new file mode 100644 index 000000000..a76dfd37a --- /dev/null +++ b/modules/xfeatures2d/src/harris_lapace_detector.cpp @@ -0,0 +1,629 @@ +// 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 layers; + Octave(); + Octave(std::vector layers); + virtual ~Octave(); + std::vector getLayers(); + Mat getLayerAt(int i); + }; + + class DOGOctave + { + public: + std::vector layers; + + DOGOctave(); + DOGOctave(std::vector layers); + virtual ~DOGOctave(); + std::vector getLayers(); + Mat getLayerAt(int i); + }; + +private: + std::vector octaves; + std::vector 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 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 _layers) : layers(_layers) {} + +/** + * Return layers of the Octave + */ +std::vector 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 _layers) : layers(_layers) {} + +Pyramid::DOGOctave::~DOGOctave() +{ +} + +std::vector 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& keypoints, InputArray mask=noArray() ); + + int numOctaves; + float corn_thresh; + float DOG_thresh; + int maxCorners; + int num_layers; +}; + +Ptr HarrisLaplaceFeatureDetector::create( + int numOctaves, + float corn_thresh, + float DOG_thresh, + int maxCorners, + int num_layers) +{ + return makePtr(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& 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_ 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 (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 (row, col); + float dy2f = Lym2smooth.at (row, col); + float dxyf = Lxmysmooth.at (row, col); + float det = dx2f * dy2f - dxyf * dxyf; + float tr = dx2f + dy2f; + float cornerness = det - (0.04f * tr * tr); + cornern_mat.at (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 (y, x); + if (val != 0 && val == corn_dilate.at (y, x)) + { + + float curVal = curDOG.at (y, x); + float prevVal = prevDOG.at (y, x); + float succVal = succDOG.at (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(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); +} + + +} +} diff --git a/modules/xfeatures2d/test/test_features2d.cpp b/modules/xfeatures2d/test/test_features2d.cpp index d42218cc4..d0b470d9e 100644 --- a/modules/xfeatures2d/test/test_features2d.cpp +++ b/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 */