Merge pull request #26405 from kaingwade:rename_features2d

Rename features2d #26405

This PR renames the module _features2d_ to _features_ as one of the Big OpenCV Cleanup #25007. 
Related PR: opencv/opencv_contrib: [#3820](https://github.com/opencv/opencv_contrib/pull/3820) opencv/ci-gha-workflow: [#192](https://github.com/opencv/ci-gha-workflow/pull/192)

modules/features/src/affine_feature.cpp

@@ -0,0 +1,360 @@
+// 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.
+//
+// This file is based on code issued with the following license.
+/*********************************************************************
+* Software License Agreement (BSD License)
+*
+*  Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+*  Copyright (C) 2008-2013, Willow Garage Inc., all rights reserved.
+*  Copyright (C) 2013, Evgeny Toropov, all rights reserved.
+*  Third party copyrights are property of their respective owners.
+*
+*  Redistribution and use in source and binary forms, with or without
+*  modification, are permitted provided that the following conditions
+*  are met:
+*
+*   * Redistributions of source code must retain the above copyright
+*     notice, this list of conditions and the following disclaimer.
+*   * Redistributions in binary form must reproduce the above
+*     copyright notice, this list of conditions and the following
+*     disclaimer in the documentation and/or other materials provided
+*     with the distribution.
+*   * The name of the copyright holders may not be used to endorse
+*     or promote products derived from this software without specific
+*     prior written permission.
+*
+*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+*  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+*  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+*  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+*  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+*  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+*  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+*  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+*  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+*  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+*  POSSIBILITY OF SUCH DAMAGE.
+*********************************************************************/
+
+/*
+ Guoshen Yu, Jean-Michel Morel, ASIFT: An Algorithm for Fully Affine
+ Invariant Comparison,  Image Processing On Line, 1 (2011), pp. 11-38.
+ https://doi.org/10.5201/ipol.2011.my-asift
+ */
+
+#include "precomp.hpp"
+#include <iostream>
+namespace cv {
+
+class AffineFeature_Impl CV_FINAL : public AffineFeature
+{
+public:
+    explicit AffineFeature_Impl(const Ptr<Feature2D>& backend,
+            int maxTilt, int minTilt, float tiltStep, float rotateStepBase);
+
+    int descriptorSize() const CV_OVERRIDE
+    {
+        return backend_->descriptorSize();
+    }
+
+    int descriptorType() const CV_OVERRIDE
+    {
+        return backend_->descriptorType();
+    }
+
+    int defaultNorm() const CV_OVERRIDE
+    {
+        return backend_->defaultNorm();
+    }
+
+    void detectAndCompute(InputArray image, InputArray mask, std::vector<KeyPoint>& keypoints,
+            OutputArray descriptors, bool useProvidedKeypoints=false) CV_OVERRIDE;
+
+    void setViewParams(const std::vector<float>& tilts, const std::vector<float>& rolls) CV_OVERRIDE;
+    void getViewParams(std::vector<float>& tilts, std::vector<float>& rolls) const CV_OVERRIDE;
+
+protected:
+    void splitKeypointsByView(const std::vector<KeyPoint>& keypoints_,
+            std::vector< std::vector<KeyPoint> >& keypointsByView) const;
+
+    const Ptr<Feature2D> backend_;
+    int maxTilt_;
+    int minTilt_;
+    float tiltStep_;
+    float rotateStepBase_;
+
+    // Tilt factors.
+    std::vector<float> tilts_;
+    // Roll factors.
+    std::vector<float> rolls_;
+
+private:
+    AffineFeature_Impl(const AffineFeature_Impl &); // copy disabled
+    AffineFeature_Impl& operator=(const AffineFeature_Impl &); // assign disabled
+};
+
+AffineFeature_Impl::AffineFeature_Impl(const Ptr<FeatureDetector>& backend,
+        int maxTilt, int minTilt, float tiltStep, float rotateStepBase)
+    : backend_(backend), maxTilt_(maxTilt), minTilt_(minTilt), tiltStep_(tiltStep), rotateStepBase_(rotateStepBase)
+{
+    int i = minTilt_;
+    if( i == 0 )
+    {
+        tilts_.push_back(1);
+        rolls_.push_back(0);
+        i++;
+    }
+    float tilt = 1;
+    for( ; i <= maxTilt_; i++ )
+    {
+        tilt *= tiltStep_;
+        float rotateStep = rotateStepBase_ / tilt;
+        int rollN = cvFloor(180.0f / rotateStep);
+        if( rollN * rotateStep == 180.0f )
+            rollN--;
+        for( int j = 0; j <= rollN; j++ )
+        {
+            tilts_.push_back(tilt);
+            rolls_.push_back(rotateStep * j);
+        }
+    }
+}
+
+void AffineFeature_Impl::setViewParams(const std::vector<float>& tilts,
+        const std::vector<float>& rolls)
+{
+    CV_Assert(tilts.size() == rolls.size());
+    tilts_ = tilts;
+    rolls_ = rolls;
+}
+
+void AffineFeature_Impl::getViewParams(std::vector<float>& tilts,
+        std::vector<float>& rolls) const
+{
+    tilts = tilts_;
+    rolls = rolls_;
+}
+
+void AffineFeature_Impl::splitKeypointsByView(const std::vector<KeyPoint>& keypoints_,
+        std::vector< std::vector<KeyPoint> >& keypointsByView) const
+{
+    for( size_t i = 0; i < keypoints_.size(); i++ )
+    {
+        const KeyPoint& kp = keypoints_[i];
+        CV_Assert( kp.class_id >= 0 && kp.class_id < (int)tilts_.size() );
+        keypointsByView[kp.class_id].push_back(kp);
+    }
+}
+
+class skewedDetectAndCompute : public ParallelLoopBody
+{
+public:
+    skewedDetectAndCompute(
+        const std::vector<float>& _tilts,
+        const std::vector<float>& _rolls,
+        std::vector< std::vector<KeyPoint> >& _keypointsCollection,
+        std::vector<Mat>& _descriptorCollection,
+        const Mat& _image,
+        const Mat& _mask,
+        const bool _do_keypoints,
+        const bool _do_descriptors,
+        const Ptr<Feature2D>& _backend)
+        : tilts(_tilts),
+          rolls(_rolls),
+          keypointsCollection(_keypointsCollection),
+          descriptorCollection(_descriptorCollection),
+          image(_image),
+          mask(_mask),
+          do_keypoints(_do_keypoints),
+          do_descriptors(_do_descriptors),
+          backend(_backend) {}
+
+    void operator()( const cv::Range& range ) const CV_OVERRIDE
+    {
+        CV_TRACE_FUNCTION();
+
+        const int begin = range.start;
+        const int end = range.end;
+
+        for( int a = begin; a < end; a++ )
+        {
+            Mat warpedImage, warpedMask;
+            Matx23f pose, invPose;
+            affineSkew(tilts[a], rolls[a], warpedImage, warpedMask, pose);
+            invertAffineTransform(pose, invPose);
+
+            std::vector<KeyPoint> wKeypoints;
+            Mat wDescriptors;
+            if( !do_keypoints )
+            {
+                const std::vector<KeyPoint>& keypointsInView = keypointsCollection[a];
+                if( keypointsInView.size() == 0 ) // when there are no keypoints in this affine view
+                    continue;
+
+                std::vector<Point2f> pts_, pts;
+                KeyPoint::convert(keypointsInView, pts_);
+                transform(pts_, pts, pose);
+                wKeypoints.resize(keypointsInView.size());
+                for( size_t wi = 0; wi < wKeypoints.size(); wi++ )
+                {
+                    wKeypoints[wi] = keypointsInView[wi];
+                    wKeypoints[wi].pt = pts[wi];
+                }
+            }
+            backend->detectAndCompute(warpedImage, warpedMask, wKeypoints, wDescriptors, !do_keypoints);
+            if( do_keypoints )
+            {
+                // KeyPointsFilter::runByPixelsMask( wKeypoints, warpedMask );
+                if( wKeypoints.size() == 0 )
+                {
+                    keypointsCollection[a].clear();
+                    continue;
+                }
+                std::vector<Point2f> pts_, pts;
+                KeyPoint::convert(wKeypoints, pts_);
+                transform(pts_, pts, invPose);
+
+                keypointsCollection[a].resize(wKeypoints.size());
+                for( size_t wi = 0; wi < wKeypoints.size(); wi++ )
+                {
+                    keypointsCollection[a][wi] = wKeypoints[wi];
+                    keypointsCollection[a][wi].pt = pts[wi];
+                    keypointsCollection[a][wi].class_id = a;
+                }
+            }
+            if( do_descriptors )
+                wDescriptors.copyTo(descriptorCollection[a]);
+        }
+    }
+private:
+    void affineSkew(float tilt, float phi,
+            Mat& warpedImage, Mat& warpedMask, Matx23f& pose) const
+    {
+        int h = image.size().height;
+        int w = image.size().width;
+        Mat rotImage;
+
+        Mat mask0;
+        if( mask.empty() )
+            mask0 = Mat(h, w, CV_8UC1, 255);
+        else
+            mask0 = mask;
+        pose = Matx23f(1,0,0,
+                       0,1,0);
+
+        if( phi == 0 )
+            image.copyTo(rotImage);
+        else
+        {
+            phi = phi * (float)CV_PI / 180;
+            float s = std::sin(phi);
+            float c = std::cos(phi);
+            Matx22f A(c, -s, s, c);
+            Matx<float, 4, 2> corners(0, 0, (float)w, 0, (float)w,(float)h, 0, (float)h);
+            Mat tf(corners * A.t());
+            Mat tcorners;
+            tf.convertTo(tcorners, CV_32S);
+            Rect rect = boundingRect(tcorners);
+            h = rect.height; w = rect.width;
+            pose = Matx23f(c, -s, -(float)rect.x,
+                        s,  c, -(float)rect.y);
+            warpAffine(image, rotImage, pose, Size(w, h), INTER_LINEAR, BORDER_REPLICATE, Scalar(), cv::ALGO_HINT_ACCURATE);
+        }
+        if( tilt == 1 )
+            warpedImage = rotImage;
+        else
+        {
+            float s = 0.8f * sqrt(tilt * tilt - 1);
+            GaussianBlur(rotImage, rotImage, Size(0, 0), s, 0.01);
+            resize(rotImage, warpedImage, Size(0, 0), 1.0/tilt, 1.0, INTER_NEAREST);
+            pose(0, 0) /= tilt;
+            pose(0, 1) /= tilt;
+            pose(0, 2) /= tilt;
+        }
+        if( phi != 0 || tilt != 1 )
+            warpAffine(mask0, warpedMask, pose, warpedImage.size(), INTER_NEAREST, BORDER_CONSTANT, Scalar(), cv::ALGO_HINT_ACCURATE);
+        else
+            warpedMask = mask0;
+    }
+
+
+    const std::vector<float>& tilts;
+    const std::vector<float>& rolls;
+    std::vector< std::vector<KeyPoint> >& keypointsCollection;
+    std::vector<Mat>& descriptorCollection;
+    const Mat& image;
+    const Mat& mask;
+    const bool do_keypoints;
+    const bool do_descriptors;
+    const Ptr<Feature2D>& backend;
+};
+
+void AffineFeature_Impl::detectAndCompute(InputArray _image, InputArray _mask,
+        std::vector<KeyPoint>& keypoints,
+        OutputArray _descriptors,
+        bool useProvidedKeypoints)
+{
+    CV_TRACE_FUNCTION();
+
+    bool do_keypoints = !useProvidedKeypoints;
+    bool do_descriptors = _descriptors.needed();
+    Mat image = _image.getMat(), mask = _mask.getMat();
+    Mat descriptors;
+
+    if( (!do_keypoints && !do_descriptors) || _image.empty() )
+        return;
+
+    std::vector< std::vector<KeyPoint> > keypointsCollection(tilts_.size());
+    std::vector< Mat > descriptorCollection(tilts_.size());
+
+    if( do_keypoints )
+        keypoints.clear();
+    else
+        splitKeypointsByView(keypoints, keypointsCollection);
+
+    parallel_for_(Range(0, (int)tilts_.size()), skewedDetectAndCompute(tilts_, rolls_, keypointsCollection, descriptorCollection,
+        image, mask, do_keypoints, do_descriptors, backend_));
+
+    if( do_keypoints )
+        for( size_t i = 0; i < keypointsCollection.size(); i++ )
+        {
+            const std::vector<KeyPoint>& keys = keypointsCollection[i];
+            keypoints.insert(keypoints.end(), keys.begin(), keys.end());
+        }
+
+    if( do_descriptors )
+    {
+        _descriptors.create((int)keypoints.size(), backend_->descriptorSize(), backend_->descriptorType());
+        descriptors = _descriptors.getMat();
+        int iter = 0;
+        for( size_t i = 0; i < descriptorCollection.size(); i++ )
+        {
+            const Mat& descs = descriptorCollection[i];
+            if( descs.empty() )
+                continue;
+            Mat roi(descriptors, Rect(0, iter, descriptors.cols, descs.rows));
+            descs.copyTo(roi);
+            iter += descs.rows;
+        }
+    }
+}
+
+
+Ptr<AffineFeature> AffineFeature::create(const Ptr<Feature2D>& backend,
+                                         int maxTilt, int minTilt, float tiltStep, float rotateStepBase)
+{
+    CV_Assert(minTilt < maxTilt);
+    CV_Assert(tiltStep > 0);
+    CV_Assert(rotateStepBase > 0);
+    return makePtr<AffineFeature_Impl>(backend, maxTilt, minTilt, tiltStep, rotateStepBase);
+}
+
+String AffineFeature::getDefaultName() const
+{
+    return (Feature2D::getDefaultName() + ".AffineFeature");
+}
+
+} // namespace

modules/features/src/annoy.cpp

@@ -0,0 +1,264 @@
+// 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"
+#include "../3rdparty/annoy/annoylib.h"
+#include <opencv2/core/utils/logger.hpp>
+
+
+namespace cv
+{
+
+struct Random
+{
+    static const uint64 default_seed = 0xffffffff;
+    #if __cplusplus < 201103L
+    typedef uint64 seed_type;
+    #endif
+
+    RNG rng;
+
+    Random(uint64 seed = default_seed)
+    {
+        rng.state = seed;
+    }
+
+    inline int flip()
+    {
+        // Draw random 0 or 1
+        return rng.next() & 1;
+    }
+
+    inline size_t index(size_t n)
+    {
+        // Draw random integer between 0 and n-1 where n is at most the number of data points you have
+        return rng(unsigned(n));
+    }
+
+    inline void set_seed(uint64 seed)
+    {
+        rng.state = seed;
+    }
+};
+
+template <typename DataType, typename DistanceType>
+class ANNIndexImpl : public ANNIndex
+{
+public:
+    ANNIndexImpl(int dimension) : dim(dimension)
+    {
+        index = makePtr<::cvannoy::AnnoyIndex<int, DataType, DistanceType, Random, ::cvannoy::AnnoyIndexSingleThreadedBuildPolicy>>(dimension);
+    }
+
+    void addItems(InputArray _dataset) CV_OVERRIDE
+    {
+        CV_Assert(!_dataset.empty());
+
+        Mat features = _dataset.getMat();
+        CV_Assert(features.cols == dim);
+        CV_Assert(features.type() == cv::DataType<DataType>::type);
+
+        int num = features.rows;
+        char* msg = nullptr;
+        if (!index->add_item(0, features.ptr<DataType>(0), &msg))
+        {
+            if (msg)
+            {
+                String errorMsg = msg;
+                free(msg);
+                CV_Error(Error::StsError, errorMsg);
+            }
+            else
+            {
+                CV_Error(Error::StsError, "Fail to add an item.");
+            }
+        }
+        for (int i = 1; i < num; ++i)
+            index->add_item(i, features.ptr<DataType>(i));
+    }
+
+    void build(int trees) CV_OVERRIDE
+    {
+        if (index->get_n_items() <= 0)
+            CV_Error(Error::StsError, "No items added. Please add items before building the index.");
+
+        if (trees <= 0)
+            trees = -1;
+
+        char* msg = nullptr;
+        if (!index->build(trees, -1, &msg))
+        {
+            if (msg)
+            {
+                String errorMsg = msg;
+                free(msg);
+                CV_Error(Error::StsError, errorMsg);
+            }
+            else
+            {
+                CV_Error(Error::StsError, "Fail to build the index.");
+            }
+        }
+    }
+
+    void knnSearch(InputArray _query, OutputArray _indices, OutputArray _dists, int knn, int search_k) CV_OVERRIDE
+    {
+        CV_Assert(!_query.empty() && _query.isContinuous());
+        Mat query = _query.getMat(), indices, dists;
+        CV_Assert(query.type() == cv::DataType<DataType>::type);
+        CV_Assert(knn > 0 && knn <= index->get_n_items());
+
+        int numQuery = query.rows;
+        if (_indices.needed())
+        {
+            indices = _indices.getMat();
+            if (!indices.isContinuous() || indices.type() != CV_32S ||
+                indices.rows != numQuery || indices.cols != knn)
+            {
+                if (!indices.isContinuous())
+                    _indices.release();
+                _indices.create(numQuery, knn, CV_32S);
+                indices = _indices.getMat();
+            }
+        }
+        else
+            indices.create(numQuery, knn, CV_32S);
+
+        if (_dists.needed())
+        {
+            dists = _dists.getMat();
+            if (!dists.isContinuous() || dists.type() != cv::DataType<DataType>::type ||
+                dists.rows != numQuery || dists.cols != knn)
+            {
+                if (!_dists.isContinuous())
+                    _dists.release();
+                _dists.create(numQuery, knn, cv::DataType<DataType>::type);
+                dists = _dists.getMat();
+            }
+        }
+        else
+            dists.create(numQuery, knn, cv::DataType<DataType>::type);
+
+        auto processBatch = [&](const Range& range)
+        {
+            std::vector<int> nns;
+            std::vector<DataType> distances;
+
+            for (int i = range.start; i < range.end; ++i)
+            {
+                index->get_nns_by_vector(query.ptr<DataType>(i), knn, search_k, &nns, &distances);
+
+                std::copy(nns.begin(), nns.end(), indices.ptr<int>(i));
+                std::copy(distances.begin(), distances.end(), dists.ptr<DataType>(i));
+
+                nns.clear();
+                distances.clear();
+            }
+        };
+
+        parallel_for_(Range(0, numQuery), processBatch);
+    }
+
+    void save(const String &filename, bool prefault) CV_OVERRIDE
+    {
+        char* msg = nullptr;
+        if (!index->save(filename.c_str(), prefault, &msg))
+        {
+            if (msg)
+            {
+                String errorMsg = msg;
+                free(msg);
+                CV_Error(Error::StsError, errorMsg);
+            }
+            else
+            {
+                CV_Error(Error::StsError, "Fail to save the index.");
+            }
+        }
+    }
+
+    void load(const String &filename, bool prefault) CV_OVERRIDE
+    {
+        char* msg = nullptr;
+        if (!index->load(filename.c_str(), prefault, &msg))
+        {
+            if (msg)
+            {
+                String errorMsg = msg;
+                free(msg);
+                CV_Error(Error::StsError, errorMsg);
+            }
+            else
+            {
+                CV_Error(Error::StsError, "Fail to load the index.");
+            }
+        }
+    }
+
+    int getTreeNumber() CV_OVERRIDE
+    {
+        return index->get_n_trees();
+    }
+
+    int getItemNumber() CV_OVERRIDE
+    {
+        return index->get_n_items();
+    }
+
+    bool setOnDiskBuild(const String &filename) CV_OVERRIDE
+    {
+        char* msg = nullptr;
+        if (index->on_disk_build(filename.c_str(), &msg))
+            return true;
+        else
+        {
+            if (msg)
+            {
+                String errorMsg = msg;
+                CV_LOG_ERROR(NULL, errorMsg);
+                free(msg);
+            }
+            else
+            {
+                CV_LOG_ERROR(NULL, "Cannot set build on disk.");
+            }
+            return false;
+        }
+    }
+
+    void setSeed(int seed) CV_OVERRIDE
+    {
+        index->set_seed(static_cast<uint32_t>(seed));
+    }
+
+private:
+    int dim;
+    Ptr<::cvannoy::AnnoyIndex<int, DataType, DistanceType, Random, ::cvannoy::AnnoyIndexSingleThreadedBuildPolicy>> index;
+};
+
+Ptr<ANNIndex> ANNIndex::create(int dim, ANNIndex::Distance distType)
+{
+    switch (distType)
+    {
+        case ANNIndex::DIST_EUCLIDEAN:
+            return makePtr<ANNIndexImpl<float, ::cvannoy::Euclidean>>(dim);
+            break;
+        case ANNIndex::DIST_MANHATTAN:
+            return makePtr<ANNIndexImpl<float, ::cvannoy::Manhattan>>(dim);
+            break;
+        case ANNIndex::DIST_ANGULAR:
+            return makePtr<ANNIndexImpl<float, ::cvannoy::Angular>>(dim);
+            break;
+        case ANNIndex::DIST_HAMMING:
+            return makePtr<ANNIndexImpl<uchar, ::cvannoy::Hamming>>(dim);
+            break;
+        case ANNIndex::DIST_DOTPRODUCT:
+            return makePtr<ANNIndexImpl<float, ::cvannoy::DotProduct>>(dim);
+            break;
+        default:
+            CV_Error(Error::StsBadArg, "Unknown/unsupported distance type");
+    }
+};
+
+}

modules/features/src/blobdetector.cpp

@@ -0,0 +1,490 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+#include <iterator>
+#include <limits>
+
+#include <opencv2/core/utils/logger.hpp>
+
+// Requires CMake flag: DEBUG_opencv_features=ON
+//#define DEBUG_BLOB_DETECTOR
+
+#ifdef DEBUG_BLOB_DETECTOR
+#include "opencv2/highgui.hpp"
+#endif
+
+namespace cv
+{
+
+class CV_EXPORTS_W SimpleBlobDetectorImpl : public SimpleBlobDetector
+{
+public:
+
+  explicit SimpleBlobDetectorImpl(const SimpleBlobDetector::Params &parameters = SimpleBlobDetector::Params());
+
+  virtual void read( const FileNode& fn ) CV_OVERRIDE;
+  virtual void write( FileStorage& fs ) const CV_OVERRIDE;
+
+  void setParams(const SimpleBlobDetector::Params& _params ) CV_OVERRIDE {
+    SimpleBlobDetectorImpl::validateParameters(_params);
+    params = _params;
+  }
+
+  SimpleBlobDetector::Params getParams() const CV_OVERRIDE { return params; }
+
+  static void validateParameters(const SimpleBlobDetector::Params& p)
+  {
+      if (p.thresholdStep <= 0)
+          CV_Error(Error::StsBadArg, "thresholdStep>0");
+
+      if (p.minThreshold > p.maxThreshold || p.minThreshold < 0)
+          CV_Error(Error::StsBadArg, "0<=minThreshold<=maxThreshold");
+
+      if (p.minDistBetweenBlobs <=0 )
+          CV_Error(Error::StsBadArg, "minDistBetweenBlobs>0");
+
+      if (p.minArea > p.maxArea || p.minArea <=0)
+          CV_Error(Error::StsBadArg, "0<minArea<=maxArea");
+
+      if (p.minCircularity > p.maxCircularity || p.minCircularity <= 0)
+          CV_Error(Error::StsBadArg, "0<minCircularity<=maxCircularity");
+
+      if (p.minInertiaRatio > p.maxInertiaRatio || p.minInertiaRatio <= 0)
+          CV_Error(Error::StsBadArg, "0<minInertiaRatio<=maxInertiaRatio");
+
+      if (p.minConvexity > p.maxConvexity || p.minConvexity <= 0)
+          CV_Error(Error::StsBadArg, "0<minConvexity<=maxConvexity");
+  }
+
+protected:
+  struct CV_EXPORTS Center
+  {
+      Point2d location;
+      double radius;
+      double confidence;
+  };
+
+  virtual void detect( InputArray image, std::vector<KeyPoint>& keypoints, InputArray mask=noArray() ) CV_OVERRIDE;
+  virtual void findBlobs(InputArray image, InputArray binaryImage, std::vector<Center> &centers,
+                         std::vector<std::vector<Point> > &contours, std::vector<Moments> &moments) const;
+  virtual const std::vector<std::vector<Point> >& getBlobContours() const CV_OVERRIDE;
+
+  Params params;
+  std::vector<std::vector<Point> > blobContours;
+};
+
+/*
+*  SimpleBlobDetector
+*/
+SimpleBlobDetector::Params::Params()
+{
+    thresholdStep = 10;
+    minThreshold = 50;
+    maxThreshold = 220;
+    minRepeatability = 2;
+    minDistBetweenBlobs = 10;
+
+    filterByColor = true;
+    blobColor = 0;
+
+    filterByArea = true;
+    minArea = 25;
+    maxArea = 5000;
+
+    filterByCircularity = false;
+    minCircularity = 0.8f;
+    maxCircularity = std::numeric_limits<float>::max();
+
+    filterByInertia = true;
+    //minInertiaRatio = 0.6;
+    minInertiaRatio = 0.1f;
+    maxInertiaRatio = std::numeric_limits<float>::max();
+
+    filterByConvexity = true;
+    //minConvexity = 0.8;
+    minConvexity = 0.95f;
+    maxConvexity = std::numeric_limits<float>::max();
+
+    collectContours = false;
+}
+
+void SimpleBlobDetector::Params::read(const cv::FileNode& fn )
+{
+    thresholdStep = fn["thresholdStep"];
+    minThreshold = fn["minThreshold"];
+    maxThreshold = fn["maxThreshold"];
+
+    minRepeatability = (size_t)(int)fn["minRepeatability"];
+    minDistBetweenBlobs = fn["minDistBetweenBlobs"];
+
+    filterByColor = (int)fn["filterByColor"] != 0 ? true : false;
+    blobColor = (uchar)(int)fn["blobColor"];
+
+    filterByArea = (int)fn["filterByArea"] != 0 ? true : false;
+    minArea = fn["minArea"];
+    maxArea = fn["maxArea"];
+
+    filterByCircularity = (int)fn["filterByCircularity"] != 0 ? true : false;
+    minCircularity = fn["minCircularity"];
+    maxCircularity = fn["maxCircularity"];
+
+    filterByInertia = (int)fn["filterByInertia"] != 0 ? true : false;
+    minInertiaRatio = fn["minInertiaRatio"];
+    maxInertiaRatio = fn["maxInertiaRatio"];
+
+    filterByConvexity = (int)fn["filterByConvexity"] != 0 ? true : false;
+    minConvexity = fn["minConvexity"];
+    maxConvexity = fn["maxConvexity"];
+
+    collectContours = (int)fn["collectContours"] != 0 ? true : false;
+}
+
+void SimpleBlobDetector::Params::write(cv::FileStorage& fs) const
+{
+    fs << "thresholdStep" << thresholdStep;
+    fs << "minThreshold" << minThreshold;
+    fs << "maxThreshold" << maxThreshold;
+
+    fs << "minRepeatability" << (int)minRepeatability;
+    fs << "minDistBetweenBlobs" << minDistBetweenBlobs;
+
+    fs << "filterByColor" << (int)filterByColor;
+    fs << "blobColor" << (int)blobColor;
+
+    fs << "filterByArea" << (int)filterByArea;
+    fs << "minArea" << minArea;
+    fs << "maxArea" << maxArea;
+
+    fs << "filterByCircularity" << (int)filterByCircularity;
+    fs << "minCircularity" << minCircularity;
+    fs << "maxCircularity" << maxCircularity;
+
+    fs << "filterByInertia" << (int)filterByInertia;
+    fs << "minInertiaRatio" << minInertiaRatio;
+    fs << "maxInertiaRatio" << maxInertiaRatio;
+
+    fs << "filterByConvexity" << (int)filterByConvexity;
+    fs << "minConvexity" << minConvexity;
+    fs << "maxConvexity" << maxConvexity;
+
+    fs << "collectContours" << (int)collectContours;
+}
+
+SimpleBlobDetectorImpl::SimpleBlobDetectorImpl(const SimpleBlobDetector::Params &parameters) :
+params(parameters)
+{
+}
+
+void SimpleBlobDetectorImpl::read( const cv::FileNode& fn )
+{
+    SimpleBlobDetector::Params rp;
+    rp.read(fn);
+    SimpleBlobDetectorImpl::validateParameters(rp);
+    params = rp;
+}
+
+void SimpleBlobDetectorImpl::write( cv::FileStorage& fs ) const
+{
+    writeFormat(fs);
+    params.write(fs);
+}
+
+void SimpleBlobDetectorImpl::findBlobs(InputArray _image, InputArray _binaryImage, std::vector<Center> &centers,
+                                       std::vector<std::vector<Point> > &contoursOut, std::vector<Moments> &momentss) const
+{
+    CV_INSTRUMENT_REGION();
+
+    Mat image = _image.getMat(), binaryImage = _binaryImage.getMat();
+    CV_UNUSED(image);
+    centers.clear();
+    contoursOut.clear();
+    momentss.clear();
+
+    std::vector < std::vector<Point> > contours;
+    findContours(binaryImage, contours, RETR_LIST, CHAIN_APPROX_NONE);
+
+#ifdef DEBUG_BLOB_DETECTOR
+    Mat keypointsImage;
+    cvtColor(binaryImage, keypointsImage, COLOR_GRAY2RGB);
+
+    Mat contoursImage;
+    cvtColor(binaryImage, contoursImage, COLOR_GRAY2RGB);
+    drawContours( contoursImage, contours, -1, Scalar(0,255,0) );
+    imshow("contours", contoursImage );
+#endif
+
+    for (size_t contourIdx = 0; contourIdx < contours.size(); contourIdx++)
+    {
+        Center center;
+        center.confidence = 1;
+        Moments moms = moments(contours[contourIdx]);
+        if (params.filterByArea)
+        {
+            double area = moms.m00;
+            if (area < params.minArea || area >= params.maxArea)
+                continue;
+        }
+
+        if (params.filterByCircularity)
+        {
+            double area = moms.m00;
+            double perimeter = arcLength(contours[contourIdx], true);
+            double ratio = 4 * CV_PI * area / (perimeter * perimeter);
+            if (ratio < params.minCircularity || ratio >= params.maxCircularity)
+                continue;
+        }
+
+        if (params.filterByInertia)
+        {
+            double denominator = std::sqrt(std::pow(2 * moms.mu11, 2) + std::pow(moms.mu20 - moms.mu02, 2));
+            const double eps = 1e-2;
+            double ratio;
+            if (denominator > eps)
+            {
+                double cosmin = (moms.mu20 - moms.mu02) / denominator;
+                double sinmin = 2 * moms.mu11 / denominator;
+                double cosmax = -cosmin;
+                double sinmax = -sinmin;
+
+                double imin = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmin - moms.mu11 * sinmin;
+                double imax = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmax - moms.mu11 * sinmax;
+                ratio = imin / imax;
+            }
+            else
+            {
+                ratio = 1;
+            }
+
+            if (ratio < params.minInertiaRatio || ratio >= params.maxInertiaRatio)
+                continue;
+
+            center.confidence = ratio * ratio;
+        }
+
+        if (params.filterByConvexity)
+        {
+            std::vector < Point > hull;
+            convexHull(contours[contourIdx], hull);
+            double area = moms.m00;
+            double hullArea = contourArea(hull);
+            if (fabs(hullArea) < DBL_EPSILON)
+                continue;
+            double ratio = area / hullArea;
+            if (ratio < params.minConvexity || ratio >= params.maxConvexity)
+                continue;
+        }
+
+        if(moms.m00 == 0.0)
+            continue;
+        center.location = Point2d(moms.m10 / moms.m00, moms.m01 / moms.m00);
+
+        if (params.filterByColor)
+        {
+            if (binaryImage.at<uchar> (cvRound(center.location.y), cvRound(center.location.x)) != params.blobColor)
+                continue;
+        }
+
+        //compute blob radius
+        {
+            std::vector<double> dists;
+            for (size_t pointIdx = 0; pointIdx < contours[contourIdx].size(); pointIdx++)
+            {
+                Point2d pt = contours[contourIdx][pointIdx];
+                dists.push_back(norm(center.location - pt));
+            }
+            std::sort(dists.begin(), dists.end());
+            center.radius = (dists[(dists.size() - 1) / 2] + dists[dists.size() / 2]) / 2.;
+        }
+
+        centers.push_back(center);
+        if (params.collectContours)
+        {
+            contoursOut.push_back(contours[contourIdx]);
+            momentss.push_back(moms);
+        }
+
+#ifdef DEBUG_BLOB_DETECTOR
+        circle( keypointsImage, center.location, 1, Scalar(0,0,255), 1 );
+#endif
+    }
+#ifdef DEBUG_BLOB_DETECTOR
+    imshow("bk", keypointsImage );
+    waitKey();
+#endif
+}
+
+void SimpleBlobDetectorImpl::detect(InputArray image, std::vector<cv::KeyPoint>& keypoints, InputArray mask)
+{
+    CV_INSTRUMENT_REGION();
+
+    keypoints.clear();
+    blobContours.clear();
+
+    CV_Assert(params.minRepeatability != 0);
+    Mat grayscaleImage;
+    if (image.channels() == 3 || image.channels() == 4)
+        cvtColor(image, grayscaleImage, COLOR_BGR2GRAY);
+    else
+        grayscaleImage = image.getMat();
+
+    if (grayscaleImage.type() != CV_8UC1) {
+        CV_Error(Error::StsUnsupportedFormat, "Blob detector only supports 8-bit images!");
+    }
+
+    CV_CheckGT(params.thresholdStep, 0.0f, "");
+    if (params.minThreshold + params.thresholdStep >= params.maxThreshold)
+    {
+        // https://github.com/opencv/opencv/issues/6667
+        CV_LOG_ONCE_INFO(NULL, "SimpleBlobDetector: params.minDistBetweenBlobs is ignored for case with single threshold");
+        CV_CheckEQ(params.minRepeatability, 1u, "Incompatible parameters for case with single threshold");
+    }
+
+    std::vector < std::vector<Center> > centers;
+    std::vector<Moments> momentss;
+    for (double thresh = params.minThreshold; thresh < params.maxThreshold; thresh += params.thresholdStep)
+    {
+        Mat binarizedImage;
+        threshold(grayscaleImage, binarizedImage, thresh, 255, THRESH_BINARY);
+
+        std::vector < Center > curCenters;
+        std::vector<std::vector<Point> > curContours;
+        std::vector<Moments> curMomentss;
+        findBlobs(grayscaleImage, binarizedImage, curCenters, curContours, curMomentss);
+        std::vector < std::vector<Center> > newCenters;
+        std::vector<std::vector<Point> > newContours;
+        std::vector<Moments> newMomentss;
+        for (size_t i = 0; i < curCenters.size(); i++)
+        {
+            bool isNew = true;
+            for (size_t j = 0; j < centers.size(); j++)
+            {
+                double dist = norm(centers[j][ centers[j].size() / 2 ].location - curCenters[i].location);
+                isNew = dist >= params.minDistBetweenBlobs && dist >= centers[j][ centers[j].size() / 2 ].radius && dist >= curCenters[i].radius;
+                if (!isNew)
+                {
+                    centers[j].push_back(curCenters[i]);
+
+                    size_t k = centers[j].size() - 1;
+                    while( k > 0 && curCenters[i].radius < centers[j][k-1].radius )
+                    {
+                        centers[j][k] = centers[j][k-1];
+                        k--;
+                    }
+
+                    if (params.collectContours)
+                    {
+                        if (curCenters[i].confidence > centers[j][k].confidence
+                            || (curCenters[i].confidence == centers[j][k].confidence && curMomentss[i].m00 > momentss[j].m00))
+                        {
+                            blobContours[j] = curContours[i];
+                            momentss[j] = curMomentss[i];
+                        }
+                    }
+                    centers[j][k] = curCenters[i];
+
+                    break;
+                }
+            }
+            if (isNew)
+            {
+                newCenters.push_back(std::vector<Center> (1, curCenters[i]));
+                if (params.collectContours)
+                {
+                    newContours.push_back(curContours[i]);
+                    newMomentss.push_back(curMomentss[i]);
+                }
+            }
+        }
+        std::copy(newCenters.begin(), newCenters.end(), std::back_inserter(centers));
+        if (params.collectContours)
+        {
+            std::copy(newContours.begin(), newContours.end(), std::back_inserter(blobContours));
+            std::copy(newMomentss.begin(), newMomentss.end(), std::back_inserter(momentss));
+        }
+    }
+
+    for (size_t i = 0; i < centers.size(); i++)
+    {
+        if (centers[i].size() < params.minRepeatability)
+            continue;
+        Point2d sumPoint(0, 0);
+        double normalizer = 0;
+        for (size_t j = 0; j < centers[i].size(); j++)
+        {
+            sumPoint += centers[i][j].confidence * centers[i][j].location;
+            normalizer += centers[i][j].confidence;
+        }
+        sumPoint *= (1. / normalizer);
+        KeyPoint kpt(sumPoint, (float)(centers[i][centers[i].size() / 2].radius) * 2.0f);
+        keypoints.push_back(kpt);
+    }
+
+    if (!mask.empty())
+    {
+        if (params.collectContours)
+        {
+            KeyPointsFilter::runByPixelsMask2VectorPoint(keypoints, blobContours, mask.getMat());
+        }
+        else
+        {
+            KeyPointsFilter::runByPixelsMask(keypoints, mask.getMat());
+        }
+    }
+}
+
+const std::vector<std::vector<Point> >& SimpleBlobDetectorImpl::getBlobContours() const {
+    return blobContours;
+}
+
+Ptr<SimpleBlobDetector> SimpleBlobDetector::create(const SimpleBlobDetector::Params& params)
+{
+    SimpleBlobDetectorImpl::validateParameters(params);
+    return makePtr<SimpleBlobDetectorImpl>(params);
+}
+
+String SimpleBlobDetector::getDefaultName() const
+{
+    return (Feature2D::getDefaultName() + ".SimpleBlobDetector");
+}
+
+}

modules/features/src/draw.cpp

@@ -0,0 +1,279 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+const int draw_shift_bits = 4;
+const int draw_multiplier = 1 << draw_shift_bits;
+
+namespace cv
+{
+
+/*
+ * Functions to draw keypoints and matches.
+ */
+static inline void _drawKeypoint( InputOutputArray img, const KeyPoint& p, const Scalar& color, DrawMatchesFlags flags )
+{
+    CV_Assert( !img.empty() );
+    Point center( cvRound(p.pt.x * draw_multiplier), cvRound(p.pt.y * draw_multiplier) );
+
+    if( !!(flags & DrawMatchesFlags::DRAW_RICH_KEYPOINTS) )
+    {
+        int radius = cvRound(p.size/2 * draw_multiplier); // KeyPoint::size is a diameter
+
+        // draw the circles around keypoints with the keypoints size
+        circle( img, center, radius, color, 1, LINE_AA, draw_shift_bits );
+
+        // draw orientation of the keypoint, if it is applicable
+        if( p.angle != -1 )
+        {
+            float srcAngleRad = p.angle*(float)CV_PI/180.f;
+            Point orient( cvRound(cos(srcAngleRad)*radius ),
+                          cvRound(sin(srcAngleRad)*radius )
+                        );
+            line( img, center, center+orient, color, 1, LINE_AA, draw_shift_bits );
+        }
+#if 0
+        else
+        {
+            // draw center with R=1
+            int radius = 1 * draw_multiplier;
+            circle( img, center, radius, color, 1, LINE_AA, draw_shift_bits );
+        }
+#endif
+    }
+    else
+    {
+        // draw center with R=3
+        int radius = 3 * draw_multiplier;
+        circle( img, center, radius, color, 1, LINE_AA, draw_shift_bits );
+    }
+}
+
+void drawKeypoints( InputArray image, const std::vector<KeyPoint>& keypoints, InputOutputArray outImage,
+                    const Scalar& _color, DrawMatchesFlags flags )
+{
+    CV_INSTRUMENT_REGION();
+
+    if( !(flags & DrawMatchesFlags::DRAW_OVER_OUTIMG) )
+    {
+        if (image.type() == CV_8UC3 || image.type() == CV_8UC4)
+        {
+            image.copyTo(outImage);
+        }
+        else if( image.type() == CV_8UC1 )
+        {
+            cvtColor( image, outImage, COLOR_GRAY2BGR );
+        }
+        else
+        {
+            CV_Error( Error::StsBadArg, "Incorrect type of input image: " + typeToString(image.type()) );
+        }
+    }
+
+    RNG& rng=theRNG();
+    bool isRandColor = _color == Scalar::all(-1);
+
+    CV_Assert( !outImage.empty() );
+    std::vector<KeyPoint>::const_iterator it = keypoints.begin(),
+                                     end = keypoints.end();
+    for( ; it != end; ++it )
+    {
+        Scalar color = isRandColor ? Scalar( rng(256), rng(256), rng(256), 255 ) : _color;
+        _drawKeypoint( outImage, *it, color, flags );
+    }
+}
+
+static void _prepareImage(InputArray src, const Mat& dst)
+{
+    CV_CheckType(src.type(), src.type() == CV_8UC1 || src.type() == CV_8UC3 || src.type() == CV_8UC4, "Unsupported source image");
+    CV_CheckType(dst.type(), dst.type() == CV_8UC3 || dst.type() == CV_8UC4, "Unsupported destination image");
+    const int src_cn = src.channels();
+    const int dst_cn = dst.channels();
+
+    if (src_cn == dst_cn)
+        src.copyTo(dst);
+    else if (src_cn == 1)
+        cvtColor(src, dst, dst_cn == 3 ? COLOR_GRAY2BGR : COLOR_GRAY2BGRA);
+    else if (src_cn == 3 && dst_cn == 4)
+        cvtColor(src, dst, COLOR_BGR2BGRA);
+    else if (src_cn == 4 && dst_cn == 3)
+        cvtColor(src, dst, COLOR_BGRA2BGR);
+    else
+        CV_Error(Error::StsInternal, "");
+}
+
+static void _prepareImgAndDrawKeypoints( InputArray img1, const std::vector<KeyPoint>& keypoints1,
+                                         InputArray img2, const std::vector<KeyPoint>& keypoints2,
+                                         InputOutputArray _outImg, Mat& outImg1, Mat& outImg2,
+                                         const Scalar& singlePointColor, DrawMatchesFlags flags )
+{
+    Mat outImg;
+    Size img1size = img1.size(), img2size = img2.size();
+    Size size( img1size.width + img2size.width, MAX(img1size.height, img2size.height) );
+    if( !!(flags & DrawMatchesFlags::DRAW_OVER_OUTIMG) )
+    {
+        outImg = _outImg.getMat();
+        if( size.width > outImg.cols || size.height > outImg.rows )
+            CV_Error( Error::StsBadSize, "outImg has size less than need to draw img1 and img2 together" );
+        outImg1 = outImg( Rect(0, 0, img1size.width, img1size.height) );
+        outImg2 = outImg( Rect(img1size.width, 0, img2size.width, img2size.height) );
+    }
+    else
+    {
+        const int cn1 = img1.channels(), cn2 = img2.channels();
+        const int out_cn = std::max(3, std::max(cn1, cn2));
+        _outImg.create(size, CV_MAKETYPE(img1.depth(), out_cn));
+        outImg = _outImg.getMat();
+        outImg = Scalar::all(0);
+        outImg1 = outImg( Rect(0, 0, img1size.width, img1size.height) );
+        outImg2 = outImg( Rect(img1size.width, 0, img2size.width, img2size.height) );
+
+        _prepareImage(img1, outImg1);
+        _prepareImage(img2, outImg2);
+    }
+
+    // draw keypoints
+    if( !(flags & DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS) )
+    {
+        Mat _outImg1 = outImg( Rect(0, 0, img1size.width, img1size.height) );
+        drawKeypoints( _outImg1, keypoints1, _outImg1, singlePointColor, flags | DrawMatchesFlags::DRAW_OVER_OUTIMG );
+
+        Mat _outImg2 = outImg( Rect(img1size.width, 0, img2size.width, img2size.height) );
+        drawKeypoints( _outImg2, keypoints2, _outImg2, singlePointColor, flags | DrawMatchesFlags::DRAW_OVER_OUTIMG );
+    }
+}
+
+static inline void _drawMatch( InputOutputArray outImg, InputOutputArray outImg1, InputOutputArray outImg2 ,
+                          const KeyPoint& kp1, const KeyPoint& kp2, const Scalar& matchColor, DrawMatchesFlags flags,
+                          const int matchesThickness )
+{
+    RNG& rng = theRNG();
+    bool isRandMatchColor = matchColor == Scalar::all(-1);
+    Scalar color = isRandMatchColor ? Scalar( rng(256), rng(256), rng(256), 255 ) : matchColor;
+
+    _drawKeypoint( outImg1, kp1, color, flags );
+    _drawKeypoint( outImg2, kp2, color, flags );
+
+    Point2f pt1 = kp1.pt,
+            pt2 = kp2.pt,
+            dpt2 = Point2f( std::min(pt2.x+outImg1.size().width, float(outImg.size().width-1)), pt2.y );
+
+    line( outImg,
+          Point(cvRound(pt1.x*draw_multiplier), cvRound(pt1.y*draw_multiplier)),
+          Point(cvRound(dpt2.x*draw_multiplier), cvRound(dpt2.y*draw_multiplier)),
+          color, matchesThickness, LINE_AA, draw_shift_bits );
+}
+
+void drawMatches( InputArray img1, const std::vector<KeyPoint>& keypoints1,
+                  InputArray img2, const std::vector<KeyPoint>& keypoints2,
+                  const std::vector<DMatch>& matches1to2, InputOutputArray outImg,
+                  const Scalar& matchColor, const Scalar& singlePointColor,
+                  const std::vector<char>& matchesMask, DrawMatchesFlags flags )
+{
+    drawMatches( img1, keypoints1,
+                 img2, keypoints2,
+                 matches1to2, outImg,
+                 1, matchColor,
+                 singlePointColor, matchesMask,
+                 flags);
+}
+
+void drawMatches( InputArray img1, const std::vector<KeyPoint>& keypoints1,
+                  InputArray img2, const std::vector<KeyPoint>& keypoints2,
+                  const std::vector<DMatch>& matches1to2, InputOutputArray outImg,
+                  const int matchesThickness, const Scalar& matchColor,
+                  const Scalar& singlePointColor, const std::vector<char>& matchesMask,
+                  DrawMatchesFlags flags )
+{
+    if( !matchesMask.empty() && matchesMask.size() != matches1to2.size() )
+        CV_Error( Error::StsBadSize, "matchesMask must have the same size as matches1to2" );
+
+    Mat outImg1, outImg2;
+    _prepareImgAndDrawKeypoints( img1, keypoints1, img2, keypoints2,
+                                 outImg, outImg1, outImg2, singlePointColor, flags );
+
+    // draw matches
+    for( size_t m = 0; m < matches1to2.size(); m++ )
+    {
+        if( matchesMask.empty() || matchesMask[m] )
+        {
+            int i1 = matches1to2[m].queryIdx;
+            int i2 = matches1to2[m].trainIdx;
+            CV_Assert(i1 >= 0 && i1 < static_cast<int>(keypoints1.size()));
+            CV_Assert(i2 >= 0 && i2 < static_cast<int>(keypoints2.size()));
+
+            const KeyPoint &kp1 = keypoints1[i1], &kp2 = keypoints2[i2];
+            _drawMatch( outImg, outImg1, outImg2, kp1, kp2, matchColor, flags, matchesThickness );
+        }
+    }
+}
+
+
+void drawMatches( InputArray img1, const std::vector<KeyPoint>& keypoints1,
+                  InputArray img2, const std::vector<KeyPoint>& keypoints2,
+                  const std::vector<std::vector<DMatch> >& matches1to2, InputOutputArray outImg,
+                  const Scalar& matchColor, const Scalar& singlePointColor,
+                  const std::vector<std::vector<char> >& matchesMask, DrawMatchesFlags flags )
+{
+    if( !matchesMask.empty() && matchesMask.size() != matches1to2.size() )
+        CV_Error( Error::StsBadSize, "matchesMask must have the same size as matches1to2" );
+
+    Mat outImg1, outImg2;
+    _prepareImgAndDrawKeypoints( img1, keypoints1, img2, keypoints2,
+                                 outImg, outImg1, outImg2, singlePointColor, flags );
+
+    // draw matches
+    for( size_t i = 0; i < matches1to2.size(); i++ )
+    {
+        for( size_t j = 0; j < matches1to2[i].size(); j++ )
+        {
+            int i1 = matches1to2[i][j].queryIdx;
+            int i2 = matches1to2[i][j].trainIdx;
+            if( matchesMask.empty() || matchesMask[i][j] )
+            {
+                const KeyPoint &kp1 = keypoints1[i1], &kp2 = keypoints2[i2];
+                _drawMatch( outImg, outImg1, outImg2, kp1, kp2, matchColor, flags, 1 );
+            }
+        }
+    }
+}
+}

modules/features/src/dynamic.cpp

@@ -0,0 +1,47 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+ //
+ //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+ //
+ //  By downloading, copying, installing or using the software you agree to this license.
+ //  If you do not agree to this license, do not download, install,
+ //  copy or use the software.
+ //
+ //
+ //                           License Agreement
+ //                For Open Source Computer Vision Library
+ //
+ // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+ // Copyright (C) 2009-2010, Willow Garage Inc., all rights reserved.
+ // Third party copyrights are property of their respective owners.
+ //
+ // Redistribution and use in source and binary forms, with or without modification,
+ // are permitted provided that the following conditions are met:
+ //
+ //   * Redistribution's of source code must retain the above copyright notice,
+ //     this list of conditions and the following disclaimer.
+ //
+ //   * Redistribution's in binary form must reproduce the above copyright notice,
+ //     this list of conditions and the following disclaimer in the documentation
+ //     and/or other materials provided with the distribution.
+ //
+ //   * The name of the copyright holders may not be used to endorse or promote products
+ //     derived from this software without specific prior written permission.
+ //
+ // This software is provided by the copyright holders and contributors "as is" and
+ // any express or implied warranties, including, but not limited to, the implied
+ // warranties of merchantability and fitness for a particular purpose are disclaimed.
+ // In no event shall the Intel Corporation or contributors be liable for any direct,
+ // indirect, incidental, special, exemplary, or consequential damages
+ // (including, but not limited to, procurement of substitute goods or services;
+ // loss of use, data, or profits; or business interruption) however caused
+ // and on any theory of liability, whether in contract, strict liability,
+ // or tort (including negligence or otherwise) arising in any way out of
+ // the use of this software, even if advised of the possibility of such damage.
+ //
+ //M*/
+
+#include "precomp.hpp"
+namespace cv
+{
+
+}

modules/features/src/evaluation.cpp

@@ -0,0 +1,572 @@
+//*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+#include <limits>
+
+using namespace cv;
+
+template<typename _Tp> static int solveQuadratic(_Tp a, _Tp b, _Tp c, _Tp& x1, _Tp& x2)
+{
+    if( a == 0 )
+    {
+        if( b == 0 )
+        {
+            x1 = x2 = 0;
+            return c == 0;
+        }
+        x1 = x2 = -c/b;
+        return 1;
+    }
+
+    _Tp d = b*b - 4*a*c;
+    if( d < 0 )
+    {
+        x1 = x2 = 0;
+        return 0;
+    }
+    if( d > 0 )
+    {
+        d = std::sqrt(d);
+        double s = 1/(2*a);
+        x1 = (-b - d)*s;
+        x2 = (-b + d)*s;
+        if( x1 > x2 )
+            std::swap(x1, x2);
+        return 2;
+    }
+    x1 = x2 = -b/(2*a);
+    return 1;
+}
+
+//for android ndk
+#undef _S
+static inline Point2f applyHomography( const Mat_<double>& H, const Point2f& pt )
+{
+    double z = H(2,0)*pt.x + H(2,1)*pt.y + H(2,2);
+    if( z )
+    {
+        double w = 1./z;
+        return Point2f( (float)((H(0,0)*pt.x + H(0,1)*pt.y + H(0,2))*w), (float)((H(1,0)*pt.x + H(1,1)*pt.y + H(1,2))*w) );
+    }
+    return Point2f( std::numeric_limits<float>::max(), std::numeric_limits<float>::max() );
+}
+
+static inline void linearizeHomographyAt( const Mat_<double>& H, const Point2f& pt, Mat_<double>& A )
+{
+    A.create(2,2);
+    double p1 = H(0,0)*pt.x + H(0,1)*pt.y + H(0,2),
+           p2 = H(1,0)*pt.x + H(1,1)*pt.y + H(1,2),
+           p3 = H(2,0)*pt.x + H(2,1)*pt.y + H(2,2),
+           p3_2 = p3*p3;
+    if( p3 )
+    {
+        A(0,0) = H(0,0)/p3 - p1*H(2,0)/p3_2; // fxdx
+        A(0,1) = H(0,1)/p3 - p1*H(2,1)/p3_2; // fxdy
+
+        A(1,0) = H(1,0)/p3 - p2*H(2,0)/p3_2; // fydx
+        A(1,1) = H(1,1)/p3 - p2*H(2,1)/p3_2; // fydx
+    }
+    else
+        A.setTo(Scalar::all(std::numeric_limits<double>::max()));
+}
+
+class EllipticKeyPoint
+{
+public:
+    EllipticKeyPoint();
+    EllipticKeyPoint( const Point2f& _center, const Scalar& _ellipse );
+
+    static void convert( const std::vector<KeyPoint>& src, std::vector<EllipticKeyPoint>& dst );
+    static void convert( const std::vector<EllipticKeyPoint>& src, std::vector<KeyPoint>& dst );
+
+    static Mat_<double> getSecondMomentsMatrix( const Scalar& _ellipse );
+    Mat_<double> getSecondMomentsMatrix() const;
+
+    void calcProjection( const Mat_<double>& H, EllipticKeyPoint& projection ) const;
+    static void calcProjection( const std::vector<EllipticKeyPoint>& src, const Mat_<double>& H, std::vector<EllipticKeyPoint>& dst );
+
+    Point2f center;
+    Scalar ellipse; // 3 elements a, b, c: ax^2+2bxy+cy^2=1
+    Size_<float> axes; // half length of ellipse axes
+    Size_<float> boundingBox; // half sizes of bounding box which sides are parallel to the coordinate axes
+};
+
+EllipticKeyPoint::EllipticKeyPoint()
+{
+    *this = EllipticKeyPoint(Point2f(0,0), Scalar(1, 0, 1) );
+}
+
+EllipticKeyPoint::EllipticKeyPoint( const Point2f& _center, const Scalar& _ellipse )
+{
+    center = _center;
+    ellipse = _ellipse;
+
+    double a = ellipse[0], b = ellipse[1], c = ellipse[2];
+    double ac_b2 = a*c - b*b;
+    double x1, x2;
+    solveQuadratic(1., -(a+c), ac_b2, x1, x2);
+    axes.width = (float)(1/sqrt(x1));
+    axes.height = (float)(1/sqrt(x2));
+
+    boundingBox.width = (float)sqrt(ellipse[2]/ac_b2);
+    boundingBox.height = (float)sqrt(ellipse[0]/ac_b2);
+}
+
+Mat_<double> EllipticKeyPoint::getSecondMomentsMatrix( const Scalar& _ellipse )
+{
+    Mat_<double> M(2, 2);
+    M(0,0) = _ellipse[0];
+    M(1,0) = M(0,1) = _ellipse[1];
+    M(1,1) = _ellipse[2];
+    return M;
+}
+
+Mat_<double> EllipticKeyPoint::getSecondMomentsMatrix() const
+{
+    return getSecondMomentsMatrix(ellipse);
+}
+
+void EllipticKeyPoint::calcProjection( const Mat_<double>& H, EllipticKeyPoint& projection ) const
+{
+    Point2f dstCenter = applyHomography(H, center);
+
+    Mat_<double> invM; invert(getSecondMomentsMatrix(), invM);
+    Mat_<double> Aff; linearizeHomographyAt(H, center, Aff);
+    Mat_<double> dstM; invert(Aff*invM*Aff.t(), dstM);
+
+    projection = EllipticKeyPoint( dstCenter, Scalar(dstM(0,0), dstM(0,1), dstM(1,1)) );
+}
+
+void EllipticKeyPoint::convert( const std::vector<KeyPoint>& src, std::vector<EllipticKeyPoint>& dst )
+{
+    CV_INSTRUMENT_REGION();
+
+    if( !src.empty() )
+    {
+        dst.resize(src.size());
+        for( size_t i = 0; i < src.size(); i++ )
+        {
+            float rad = src[i].size/2;
+            CV_Assert( rad );
+            float fac = 1.f/(rad*rad);
+            dst[i] = EllipticKeyPoint( src[i].pt, Scalar(fac, 0, fac) );
+        }
+    }
+}
+
+void EllipticKeyPoint::convert( const std::vector<EllipticKeyPoint>& src, std::vector<KeyPoint>& dst )
+{
+    CV_INSTRUMENT_REGION();
+
+    if( !src.empty() )
+    {
+        dst.resize(src.size());
+        for( size_t i = 0; i < src.size(); i++ )
+        {
+            Size_<float> axes = src[i].axes;
+            float rad = sqrt(axes.height*axes.width);
+            dst[i] = KeyPoint(src[i].center, 2*rad );
+        }
+    }
+}
+
+void EllipticKeyPoint::calcProjection( const std::vector<EllipticKeyPoint>& src, const Mat_<double>& H, std::vector<EllipticKeyPoint>& dst )
+{
+    if( !src.empty() )
+    {
+        CV_Assert( !H.empty() && H.cols == 3 && H.rows == 3);
+        dst.resize(src.size());
+        std::vector<EllipticKeyPoint>::const_iterator srcIt = src.begin();
+        std::vector<EllipticKeyPoint>::iterator       dstIt = dst.begin();
+        for( ; srcIt != src.end() && dstIt != dst.end(); ++srcIt, ++dstIt )
+            srcIt->calcProjection(H, *dstIt);
+    }
+}
+
+static void filterEllipticKeyPointsByImageSize( std::vector<EllipticKeyPoint>& keypoints, const Size& imgSize )
+{
+    if( !keypoints.empty() )
+    {
+        std::vector<EllipticKeyPoint> filtered;
+        filtered.reserve(keypoints.size());
+        std::vector<EllipticKeyPoint>::const_iterator it = keypoints.begin();
+        for( int i = 0; it != keypoints.end(); ++it, i++ )
+        {
+            if( it->center.x + it->boundingBox.width < imgSize.width &&
+                it->center.x - it->boundingBox.width > 0 &&
+                it->center.y + it->boundingBox.height < imgSize.height &&
+                it->center.y - it->boundingBox.height > 0 )
+                filtered.push_back(*it);
+        }
+        keypoints.assign(filtered.begin(), filtered.end());
+    }
+}
+
+struct IntersectAreaCounter
+{
+    IntersectAreaCounter( float _dr, int _minx,
+                          int _miny, int _maxy,
+                          const Point2f& _diff,
+                          const Scalar& _ellipse1, const Scalar& _ellipse2 ) :
+               dr(_dr), bua(0), bna(0), minx(_minx), miny(_miny), maxy(_maxy),
+               diff(_diff), ellipse1(_ellipse1), ellipse2(_ellipse2) {}
+    IntersectAreaCounter( const IntersectAreaCounter& counter, Split )
+    {
+        *this = counter;
+        bua = 0;
+        bna = 0;
+    }
+
+    void operator()( const BlockedRange& range )
+    {
+        CV_Assert( miny < maxy );
+        CV_Assert( dr > FLT_EPSILON );
+
+        int temp_bua = bua, temp_bna = bna;
+        for( int i = range.begin(); i != range.end(); i++ )
+        {
+            float rx1 = minx + i*dr;
+            float rx2 = rx1 - diff.x;
+            for( float ry1 = (float)miny; ry1 <= (float)maxy; ry1 += dr )
+            {
+                float ry2 = ry1 - diff.y;
+                //compute the distance from the ellipse center
+                float e1 = (float)(ellipse1[0]*rx1*rx1 + 2*ellipse1[1]*rx1*ry1 + ellipse1[2]*ry1*ry1);
+                float e2 = (float)(ellipse2[0]*rx2*rx2 + 2*ellipse2[1]*rx2*ry2 + ellipse2[2]*ry2*ry2);
+                //compute the area
+                if( e1<1 && e2<1 ) temp_bna++;
+                if( e1<1 || e2<1 ) temp_bua++;
+            }
+        }
+        bua = temp_bua;
+        bna = temp_bna;
+    }
+
+    void join( IntersectAreaCounter& ac )
+    {
+        bua += ac.bua;
+        bna += ac.bna;
+    }
+
+    float dr;
+    int bua, bna;
+
+    int minx;
+    int miny, maxy;
+
+    Point2f diff;
+    Scalar ellipse1, ellipse2;
+
+};
+
+struct SIdx
+{
+    SIdx() : S(-1), i1(-1), i2(-1) {}
+    SIdx(float _S, int _i1, int _i2) : S(_S), i1(_i1), i2(_i2) {}
+    float S;
+    int i1;
+    int i2;
+
+    bool operator<(const SIdx& v) const { return S > v.S; }
+
+    struct UsedFinder
+    {
+        UsedFinder(const SIdx& _used) : used(_used) {}
+        const SIdx& used;
+        bool operator()(const SIdx& v) const { return  (v.i1 == used.i1 || v.i2 == used.i2); }
+        UsedFinder& operator=(const UsedFinder&) = delete;
+        // To avoid -Wdeprecated-copy warning, copy constructor is needed.
+        UsedFinder(const UsedFinder&) = default;
+    };
+};
+
+static void computeOneToOneMatchedOverlaps( const std::vector<EllipticKeyPoint>& keypoints1, const std::vector<EllipticKeyPoint>& keypoints2t,
+                                            bool commonPart, std::vector<SIdx>& overlaps, float minOverlap )
+{
+    CV_Assert( minOverlap >= 0.f );
+    overlaps.clear();
+    if( keypoints1.empty() || keypoints2t.empty() )
+        return;
+
+    overlaps.clear();
+    overlaps.reserve(cvRound(keypoints1.size() * keypoints2t.size() * 0.01));
+
+    for( size_t i1 = 0; i1 < keypoints1.size(); i1++ )
+    {
+        EllipticKeyPoint kp1 = keypoints1[i1];
+        float maxDist = sqrt(kp1.axes.width*kp1.axes.height),
+              fac = 30.f/maxDist;
+        if( !commonPart )
+            fac=3;
+
+        maxDist = maxDist*4;
+        fac = 1.f/(fac*fac);
+
+        EllipticKeyPoint keypoint1a = EllipticKeyPoint( kp1.center, Scalar(fac*kp1.ellipse[0], fac*kp1.ellipse[1], fac*kp1.ellipse[2]) );
+
+        for( size_t i2 = 0; i2 < keypoints2t.size(); i2++ )
+        {
+            EllipticKeyPoint kp2 = keypoints2t[i2];
+            Point2f diff = kp2.center - kp1.center;
+
+            if( norm(diff) < maxDist )
+            {
+                EllipticKeyPoint keypoint2a = EllipticKeyPoint( kp2.center, Scalar(fac*kp2.ellipse[0], fac*kp2.ellipse[1], fac*kp2.ellipse[2]) );
+                //find the largest eigenvalue
+                int maxx =  (int)ceil(( keypoint1a.boundingBox.width > (diff.x+keypoint2a.boundingBox.width)) ?
+                                     keypoint1a.boundingBox.width : (diff.x+keypoint2a.boundingBox.width));
+                int minx = (int)floor((-keypoint1a.boundingBox.width < (diff.x-keypoint2a.boundingBox.width)) ?
+                                    -keypoint1a.boundingBox.width : (diff.x-keypoint2a.boundingBox.width));
+
+                int maxy =  (int)ceil(( keypoint1a.boundingBox.height > (diff.y+keypoint2a.boundingBox.height)) ?
+                                     keypoint1a.boundingBox.height : (diff.y+keypoint2a.boundingBox.height));
+                int miny = (int)floor((-keypoint1a.boundingBox.height < (diff.y-keypoint2a.boundingBox.height)) ?
+                                    -keypoint1a.boundingBox.height : (diff.y-keypoint2a.boundingBox.height));
+                int mina = (maxx-minx) < (maxy-miny) ? (maxx-minx) : (maxy-miny) ;
+
+                //compute the area
+                float dr = (float)mina/50.f;
+                int N = (int)floor((float)(maxx - minx) / dr);
+                IntersectAreaCounter ac( dr, minx, miny, maxy, diff, keypoint1a.ellipse, keypoint2a.ellipse );
+                parallel_reduce( BlockedRange(0, N+1), ac );
+                if( ac.bna > 0 )
+                {
+                    float ov =  (float)ac.bna / (float)ac.bua;
+                    if( ov >= minOverlap )
+                        overlaps.push_back(SIdx(ov, (int)i1, (int)i2));
+                }
+            }
+        }
+    }
+
+    std::sort( overlaps.begin(), overlaps.end() );
+
+    typedef std::vector<SIdx>::iterator It;
+
+    It pos = overlaps.begin();
+    It end = overlaps.end();
+
+    while(pos != end)
+    {
+        It prev = pos++;
+        end = std::remove_if(pos, end, SIdx::UsedFinder(*prev));
+    }
+    overlaps.erase(pos, overlaps.end());
+}
+
+static void calculateRepeatability( const Mat& img1, const Mat& img2, const Mat& H1to2,
+                                    const std::vector<KeyPoint>& _keypoints1, const std::vector<KeyPoint>& _keypoints2,
+                                    float& repeatability, int& correspondencesCount,
+                                    Mat* thresholdedOverlapMask=0  )
+{
+    std::vector<EllipticKeyPoint> keypoints1, keypoints2, keypoints1t, keypoints2t;
+    EllipticKeyPoint::convert( _keypoints1, keypoints1 );
+    EllipticKeyPoint::convert( _keypoints2, keypoints2 );
+
+    // calculate projections of key points
+    EllipticKeyPoint::calcProjection( keypoints1, H1to2, keypoints1t );
+    Mat H2to1; invert(H1to2, H2to1);
+    EllipticKeyPoint::calcProjection( keypoints2, H2to1, keypoints2t );
+
+    float overlapThreshold;
+    bool ifEvaluateDetectors = thresholdedOverlapMask == 0;
+    if( ifEvaluateDetectors )
+    {
+        overlapThreshold = 1.f - 0.4f;
+
+        // remove key points from outside of the common image part
+        Size sz1 = img1.size(), sz2 = img2.size();
+        filterEllipticKeyPointsByImageSize( keypoints1, sz1 );
+        filterEllipticKeyPointsByImageSize( keypoints1t, sz2 );
+        filterEllipticKeyPointsByImageSize( keypoints2, sz2 );
+        filterEllipticKeyPointsByImageSize( keypoints2t, sz1 );
+    }
+    else
+    {
+        overlapThreshold = 1.f - 0.5f;
+
+        thresholdedOverlapMask->create( (int)keypoints1.size(), (int)keypoints2t.size(), CV_8UC1 );
+        thresholdedOverlapMask->setTo( Scalar::all(0) );
+    }
+    size_t size1 = keypoints1.size(), size2 = keypoints2t.size();
+    size_t minCount = MIN( size1, size2 );
+
+    // calculate overlap errors
+    std::vector<SIdx> overlaps;
+    computeOneToOneMatchedOverlaps( keypoints1, keypoints2t, ifEvaluateDetectors, overlaps, overlapThreshold/*min overlap*/ );
+
+    correspondencesCount = -1;
+    repeatability = -1.f;
+    if( overlaps.empty() )
+        return;
+
+    if( ifEvaluateDetectors )
+    {
+        // regions one-to-one matching
+        correspondencesCount = (int)overlaps.size();
+        repeatability = minCount ? (float)correspondencesCount / minCount : -1;
+    }
+    else
+    {
+        for( size_t i = 0; i < overlaps.size(); i++ )
+        {
+            int y = overlaps[i].i1;
+            int x = overlaps[i].i2;
+            thresholdedOverlapMask->at<uchar>(y,x) = 1;
+        }
+    }
+}
+
+void cv::evaluateFeatureDetector( const Mat& img1, const Mat& img2, const Mat& H1to2,
+                              std::vector<KeyPoint>* _keypoints1, std::vector<KeyPoint>* _keypoints2,
+                              float& repeatability, int& correspCount,
+                              const Ptr<FeatureDetector>& _fdetector )
+{
+    CV_INSTRUMENT_REGION();
+
+    Ptr<FeatureDetector> fdetector(_fdetector);
+    std::vector<KeyPoint> *keypoints1, *keypoints2, buf1, buf2;
+    keypoints1 = _keypoints1 != 0 ? _keypoints1 : &buf1;
+    keypoints2 = _keypoints2 != 0 ? _keypoints2 : &buf2;
+
+    if( (keypoints1->empty() || keypoints2->empty()) && !fdetector )
+        CV_Error( Error::StsBadArg, "fdetector must not be empty when keypoints1 or keypoints2 is empty" );
+
+    if( keypoints1->empty() )
+        fdetector->detect( img1, *keypoints1 );
+    if( keypoints2->empty() )
+        fdetector->detect( img2, *keypoints2 );
+
+    calculateRepeatability( img1, img2, H1to2, *keypoints1, *keypoints2, repeatability, correspCount );
+}
+
+struct DMatchForEvaluation : public DMatch
+{
+    uchar isCorrect;
+    DMatchForEvaluation( const DMatch &dm ) : DMatch( dm ), isCorrect(0) {}
+};
+
+static inline float recall( int correctMatchCount, int correspondenceCount )
+{
+    return correspondenceCount ? (float)correctMatchCount / (float)correspondenceCount : -1;
+}
+
+static inline float precision( int correctMatchCount, int falseMatchCount )
+{
+    return correctMatchCount + falseMatchCount ? (float)correctMatchCount / (float)(correctMatchCount + falseMatchCount) : -1;
+}
+
+void cv::computeRecallPrecisionCurve( const std::vector<std::vector<DMatch> >& matches1to2,
+                                      const std::vector<std::vector<uchar> >& correctMatches1to2Mask,
+                                      std::vector<Point2f>& recallPrecisionCurve )
+{
+    CV_INSTRUMENT_REGION();
+
+    CV_Assert( matches1to2.size() == correctMatches1to2Mask.size() );
+
+    std::vector<DMatchForEvaluation> allMatches;
+    int correspondenceCount = 0;
+    for( size_t i = 0; i < matches1to2.size(); i++ )
+    {
+        for( size_t j = 0; j < matches1to2[i].size(); j++ )
+        {
+            DMatchForEvaluation match = matches1to2[i][j];
+            match.isCorrect = correctMatches1to2Mask[i][j] ;
+            allMatches.push_back( match );
+            correspondenceCount += match.isCorrect != 0 ? 1 : 0;
+        }
+    }
+
+    std::sort( allMatches.begin(), allMatches.end() );
+
+    int correctMatchCount = 0, falseMatchCount = 0;
+    recallPrecisionCurve.resize( allMatches.size() );
+    for( size_t i = 0; i < allMatches.size(); i++ )
+    {
+        if( allMatches[i].isCorrect )
+            correctMatchCount++;
+        else
+            falseMatchCount++;
+
+        float r = recall( correctMatchCount, correspondenceCount );
+        float p =  precision( correctMatchCount, falseMatchCount );
+        recallPrecisionCurve[i] = Point2f(1-p, r);
+    }
+}
+
+float cv::getRecall( const std::vector<Point2f>& recallPrecisionCurve, float l_precision )
+{
+    CV_INSTRUMENT_REGION();
+
+    int nearestPointIndex = getNearestPoint( recallPrecisionCurve, l_precision );
+
+    float recall = -1.f;
+
+    if( nearestPointIndex >= 0 )
+        recall = recallPrecisionCurve[nearestPointIndex].y;
+
+    return recall;
+}
+
+int cv::getNearestPoint( const std::vector<Point2f>& recallPrecisionCurve, float l_precision )
+{
+    CV_INSTRUMENT_REGION();
+
+    int nearestPointIndex = -1;
+
+    if( l_precision >= 0 && l_precision <= 1 )
+    {
+        float minDiff = FLT_MAX;
+        for( size_t i = 0; i < recallPrecisionCurve.size(); i++ )
+        {
+            float curDiff = std::fabs(l_precision - recallPrecisionCurve[i].x);
+            if( curDiff <= minDiff )
+            {
+                nearestPointIndex = (int)i;
+                minDiff = curDiff;
+            }
+        }
+    }
+
+    return nearestPointIndex;
+}

modules/features/src/fast.avx2.cpp

@@ -0,0 +1,184 @@
+/* This is FAST corner detector, contributed to OpenCV by the author, Edward Rosten.
+   Below is the original copyright and the references */
+
+/*
+Copyright (c) 2006, 2008 Edward Rosten
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+    *Redistributions of source code must retain the above copyright
+     notice, this list of conditions and the following disclaimer.
+
+    *Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+
+    *Neither the name of the University of Cambridge nor the names of
+     its contributors may be used to endorse or promote products derived
+     from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+The references are:
+ * Machine learning for high-speed corner detection,
+   E. Rosten and T. Drummond, ECCV 2006
+ * Faster and better: A machine learning approach to corner detection
+   E. Rosten, R. Porter and T. Drummond, PAMI, 2009
+*/
+
+#include "precomp.hpp"
+#include "fast.hpp"
+#include "opencv2/core/hal/intrin.hpp"
+
+namespace cv
+{
+namespace opt_AVX2
+{
+
+class FAST_t_patternSize16_AVX2_Impl CV_FINAL: public FAST_t_patternSize16_AVX2
+{
+public:
+    FAST_t_patternSize16_AVX2_Impl(int _cols, int _threshold, bool _nonmax_suppression, const int* _pixel):
+                                   cols(_cols), nonmax_suppression(_nonmax_suppression), pixel(_pixel)
+    {
+        //patternSize = 16
+        t256c = (char)_threshold;
+        threshold = std::min(std::max(_threshold, 0), 255);
+    }
+
+    virtual void process(int &j, const uchar* &ptr, uchar* curr, int* cornerpos, int &ncorners) CV_OVERRIDE
+    {
+        static const __m256i delta256 = _mm256_broadcastsi128_si256(_mm_set1_epi8((char)(-128))), K16_256 = _mm256_broadcastsi128_si256(_mm_set1_epi8((char)8));
+        const __m256i t256 = _mm256_broadcastsi128_si256(_mm_set1_epi8(t256c));
+        for (; j < cols - 32 - 3; j += 32, ptr += 32)
+        {
+            __m256i m0, m1;
+            __m256i v0 = _mm256_loadu_si256((const __m256i*)ptr);
+
+            __m256i v1 = _mm256_xor_si256(_mm256_subs_epu8(v0, t256), delta256);
+            v0 = _mm256_xor_si256(_mm256_adds_epu8(v0, t256), delta256);
+
+            __m256i x0 = _mm256_sub_epi8(_mm256_loadu_si256((const __m256i*)(ptr + pixel[0])), delta256);
+            __m256i x1 = _mm256_sub_epi8(_mm256_loadu_si256((const __m256i*)(ptr + pixel[4])), delta256);
+            __m256i x2 = _mm256_sub_epi8(_mm256_loadu_si256((const __m256i*)(ptr + pixel[8])), delta256);
+            __m256i x3 = _mm256_sub_epi8(_mm256_loadu_si256((const __m256i*)(ptr + pixel[12])), delta256);
+
+            m0 = _mm256_and_si256(_mm256_cmpgt_epi8(x0, v0), _mm256_cmpgt_epi8(x1, v0));
+            m1 = _mm256_and_si256(_mm256_cmpgt_epi8(v1, x0), _mm256_cmpgt_epi8(v1, x1));
+            m0 = _mm256_or_si256(m0, _mm256_and_si256(_mm256_cmpgt_epi8(x1, v0), _mm256_cmpgt_epi8(x2, v0)));
+            m1 = _mm256_or_si256(m1, _mm256_and_si256(_mm256_cmpgt_epi8(v1, x1), _mm256_cmpgt_epi8(v1, x2)));
+            m0 = _mm256_or_si256(m0, _mm256_and_si256(_mm256_cmpgt_epi8(x2, v0), _mm256_cmpgt_epi8(x3, v0)));
+            m1 = _mm256_or_si256(m1, _mm256_and_si256(_mm256_cmpgt_epi8(v1, x2), _mm256_cmpgt_epi8(v1, x3)));
+            m0 = _mm256_or_si256(m0, _mm256_and_si256(_mm256_cmpgt_epi8(x3, v0), _mm256_cmpgt_epi8(x0, v0)));
+            m1 = _mm256_or_si256(m1, _mm256_and_si256(_mm256_cmpgt_epi8(v1, x3), _mm256_cmpgt_epi8(v1, x0)));
+            m0 = _mm256_or_si256(m0, m1);
+
+            unsigned int mask = _mm256_movemask_epi8(m0); //unsigned is important!
+            if (mask == 0){
+                continue;
+            }
+            if ((mask & 0xffff) == 0)
+            {
+                j -= 16;
+                ptr -= 16;
+                continue;
+            }
+
+            __m256i c0 = _mm256_setzero_si256(), c1 = c0, max0 = c0, max1 = c0;
+            for (int k = 0; k < 25; k++)
+            {
+                __m256i x = _mm256_xor_si256(_mm256_loadu_si256((const __m256i*)(ptr + pixel[k])), delta256);
+                m0 = _mm256_cmpgt_epi8(x, v0);
+                m1 = _mm256_cmpgt_epi8(v1, x);
+
+                c0 = _mm256_and_si256(_mm256_sub_epi8(c0, m0), m0);
+                c1 = _mm256_and_si256(_mm256_sub_epi8(c1, m1), m1);
+
+                max0 = _mm256_max_epu8(max0, c0);
+                max1 = _mm256_max_epu8(max1, c1);
+            }
+
+            max0 = _mm256_max_epu8(max0, max1);
+            unsigned int m = _mm256_movemask_epi8(_mm256_cmpgt_epi8(max0, K16_256));
+
+            for (int k = 0; m > 0 && k < 32; k++, m >>= 1)
+                if (m & 1)
+                {
+                    cornerpos[ncorners++] = j + k;
+                    if (nonmax_suppression)
+                    {
+                        short d[25];
+                        for (int q = 0; q < 25; q++)
+                            d[q] = (short)(ptr[k] - ptr[k + pixel[q]]);
+                        v_int16x8 q0 = v_setall_s16(-1000), q1 = v_setall_s16(1000);
+                        for (int q = 0; q < 16; q += 8)
+                        {
+                            v_int16x8 v0_ = v_load(d + q + 1);
+                            v_int16x8 v1_ = v_load(d + q + 2);
+                            v_int16x8 a = v_min(v0_, v1_);
+                            v_int16x8 b = v_max(v0_, v1_);
+                            v0_ = v_load(d + q + 3);
+                            a = v_min(a, v0_);
+                            b = v_max(b, v0_);
+                            v0_ = v_load(d + q + 4);
+                            a = v_min(a, v0_);
+                            b = v_max(b, v0_);
+                            v0_ = v_load(d + q + 5);
+                            a = v_min(a, v0_);
+                            b = v_max(b, v0_);
+                            v0_ = v_load(d + q + 6);
+                            a = v_min(a, v0_);
+                            b = v_max(b, v0_);
+                            v0_ = v_load(d + q + 7);
+                            a = v_min(a, v0_);
+                            b = v_max(b, v0_);
+                            v0_ = v_load(d + q + 8);
+                            a = v_min(a, v0_);
+                            b = v_max(b, v0_);
+                            v0_ = v_load(d + q);
+                            q0 = v_max(q0, v_min(a, v0_));
+                            q1 = v_min(q1, v_max(b, v0_));
+                            v0_ = v_load(d + q + 9);
+                            q0 = v_max(q0, v_min(a, v0_));
+                            q1 = v_min(q1, v_max(b, v0_));
+                        }
+                        q0 = v_max(q0, v_sub(v_setzero_s16(), q1));
+                        curr[j + k] = (uchar)(v_reduce_max(q0) - 1);
+                    }
+                }
+        }
+        _mm256_zeroupper();
+    }
+
+    virtual ~FAST_t_patternSize16_AVX2_Impl() CV_OVERRIDE {}
+
+private:
+    int cols;
+    char t256c;
+    int threshold;
+    bool nonmax_suppression;
+    const int* pixel;
+};
+
+Ptr<FAST_t_patternSize16_AVX2> FAST_t_patternSize16_AVX2::getImpl(int _cols, int _threshold, bool _nonmax_suppression, const int* _pixel)
+{
+    return Ptr<FAST_t_patternSize16_AVX2>(new FAST_t_patternSize16_AVX2_Impl(_cols, _threshold, _nonmax_suppression, _pixel));
+}
+
+}
+}

modules/features/src/fast.cpp

@@ -0,0 +1,557 @@
+/* This is FAST corner detector, contributed to OpenCV by the author, Edward Rosten.
+   Below is the original copyright and the references */
+
+/*
+Copyright (c) 2006, 2008 Edward Rosten
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+    *Redistributions of source code must retain the above copyright
+     notice, this list of conditions and the following disclaimer.
+
+    *Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+
+    *Neither the name of the University of Cambridge nor the names of
+     its contributors may be used to endorse or promote products derived
+     from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+The references are:
+ * Machine learning for high-speed corner detection,
+   E. Rosten and T. Drummond, ECCV 2006
+ * Faster and better: A machine learning approach to corner detection
+   E. Rosten, R. Porter and T. Drummond, PAMI, 2009
+*/
+
+#include "precomp.hpp"
+#include "fast.hpp"
+#include "fast_score.hpp"
+#include "opencl_kernels_features.hpp"
+#include "hal_replacement.hpp"
+#include "opencv2/core/hal/intrin.hpp"
+#include "opencv2/core/utils/buffer_area.private.hpp"
+
+namespace cv
+{
+
+template<int patternSize>
+void FAST_t(InputArray _img, std::vector<KeyPoint>& keypoints, int threshold, bool nonmax_suppression)
+{
+    Mat img = _img.getMat();
+    const int K = patternSize/2, N = patternSize + K + 1;
+    int i, j, k, pixel[25];
+    makeOffsets(pixel, (int)img.step, patternSize);
+
+#if CV_SIMD128
+    const int quarterPatternSize = patternSize/4;
+    v_uint8x16 delta = v_setall_u8(0x80), t = v_setall_u8((char)threshold), K16 = v_setall_u8((char)K);
+#if CV_TRY_AVX2
+    Ptr<opt_AVX2::FAST_t_patternSize16_AVX2> fast_t_impl_avx2;
+    if(CV_CPU_HAS_SUPPORT_AVX2)
+        fast_t_impl_avx2 = opt_AVX2::FAST_t_patternSize16_AVX2::getImpl(img.cols, threshold, nonmax_suppression, pixel);
+#endif
+
+#endif
+
+    keypoints.clear();
+
+    threshold = std::min(std::max(threshold, 0), 255);
+
+    uchar threshold_tab[512];
+    for( i = -255; i <= 255; i++ )
+        threshold_tab[i+255] = (uchar)(i < -threshold ? 1 : i > threshold ? 2 : 0);
+
+    uchar* buf[3] = { 0 };
+    int* cpbuf[3] = { 0 };
+    utils::BufferArea area;
+    for (unsigned idx = 0; idx < 3; ++idx)
+    {
+        area.allocate(buf[idx], img.cols);
+        area.allocate(cpbuf[idx], img.cols + 1);
+    }
+    area.commit();
+
+    for (unsigned idx = 0; idx < 3; ++idx)
+    {
+        memset(buf[idx], 0, img.cols);
+    }
+
+    for(i = 3; i < img.rows-2; i++)
+    {
+        const uchar* ptr = img.ptr<uchar>(i) + 3;
+        uchar* curr = buf[(i - 3)%3];
+        int* cornerpos = cpbuf[(i - 3)%3] + 1; // cornerpos[-1] is used to store a value
+        memset(curr, 0, img.cols);
+        int ncorners = 0;
+
+        if( i < img.rows - 3 )
+        {
+            j = 3;
+#if CV_SIMD128
+            {
+                if( patternSize == 16 )
+                {
+#if CV_TRY_AVX2
+                    if (fast_t_impl_avx2)
+                        fast_t_impl_avx2->process(j, ptr, curr, cornerpos, ncorners);
+#endif
+                    //vz if (j <= (img.cols - 27)) //it doesn't make sense using vectors for less than 8 elements
+                    {
+                        for (; j < img.cols - 16 - 3; j += 16, ptr += 16)
+                        {
+                            v_uint8x16 v = v_load(ptr);
+                            v_int8x16 v0 = v_reinterpret_as_s8(v_xor(v_add(v, t), delta));
+                            v_int8x16 v1 = v_reinterpret_as_s8(v_xor(v_sub(v, t), delta));
+
+                            v_int8x16 x0 = v_reinterpret_as_s8(v_sub_wrap(v_load(ptr + pixel[0]), delta));
+                            v_int8x16 x1 = v_reinterpret_as_s8(v_sub_wrap(v_load(ptr + pixel[quarterPatternSize]), delta));
+                            v_int8x16 x2 = v_reinterpret_as_s8(v_sub_wrap(v_load(ptr + pixel[2*quarterPatternSize]), delta));
+                            v_int8x16 x3 = v_reinterpret_as_s8(v_sub_wrap(v_load(ptr + pixel[3*quarterPatternSize]), delta));
+
+                            v_int8x16 m0, m1;
+                            m0 = v_and(v_lt(v0, x0), v_lt(v0, x1));
+                            m1 = v_and(v_lt(x0, v1), v_lt(x1, v1));
+                            m0 = v_or(m0, v_and(v_lt(v0, x1), v_lt(v0, x2)));
+                            m1 = v_or(m1, v_and(v_lt(x1, v1), v_lt(x2, v1)));
+                            m0 = v_or(m0, v_and(v_lt(v0, x2), v_lt(v0, x3)));
+                            m1 = v_or(m1, v_and(v_lt(x2, v1), v_lt(x3, v1)));
+                            m0 = v_or(m0, v_and(v_lt(v0, x3), v_lt(v0, x0)));
+                            m1 = v_or(m1, v_and(v_lt(x3, v1), v_lt(x0, v1)));
+                            m0 = v_or(m0, m1);
+
+                            if( !v_check_any(m0) )
+                                continue;
+                            if( !v_check_any(v_combine_low(m0, m0)) )
+                            {
+                                j -= 8;
+                                ptr -= 8;
+                                continue;
+                            }
+
+                            v_int8x16 c0 = v_setzero_s8();
+                            v_int8x16 c1 = v_setzero_s8();
+                            v_uint8x16 max0 = v_setzero_u8();
+                            v_uint8x16 max1 = v_setzero_u8();
+                            for( k = 0; k < N; k++ )
+                            {
+                                v_int8x16 x = v_reinterpret_as_s8(v_xor(v_load((ptr + pixel[k])), delta));
+                                m0 = v_lt(v0, x);
+                                m1 = v_lt(x, v1);
+
+                                c0 = v_and(v_sub_wrap(c0, m0), m0);
+                                c1 = v_and(v_sub_wrap(c1, m1), m1);
+
+                                max0 = v_max(max0, v_reinterpret_as_u8(c0));
+                                max1 = v_max(max1, v_reinterpret_as_u8(c1));
+                            }
+
+                            max0 = v_lt(K16, v_max(max0, max1));
+                            unsigned int m = v_signmask(v_reinterpret_as_s8(max0));
+
+                            for( k = 0; m > 0 && k < 16; k++, m >>= 1 )
+                            {
+                                if( m & 1 )
+                                {
+                                    cornerpos[ncorners++] = j+k;
+                                    if(nonmax_suppression)
+                                    {
+                                        short d[25];
+                                        for (int _k = 0; _k < 25; _k++)
+                                            d[_k] = (short)(ptr[k] - ptr[k + pixel[_k]]);
+
+                                        v_int16x8 a0, b0, a1, b1;
+                                        a0 = b0 = a1 = b1 = v_load(d + 8);
+                                        for(int shift = 0; shift < 8; ++shift)
+                                        {
+                                            v_int16x8 v_nms = v_load(d + shift);
+                                            a0 = v_min(a0, v_nms);
+                                            b0 = v_max(b0, v_nms);
+                                            v_nms = v_load(d + 9 + shift);
+                                            a1 = v_min(a1, v_nms);
+                                            b1 = v_max(b1, v_nms);
+                                        }
+                                        curr[j + k] = (uchar)(v_reduce_max(v_max(v_max(a0, a1), v_sub(v_setzero_s16(), v_min(b0, b1)))) - 1);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+#endif
+            for( ; j < img.cols - 3; j++, ptr++ )
+            {
+                int v = ptr[0];
+                const uchar* tab = &threshold_tab[0] - v + 255;
+                int d = tab[ptr[pixel[0]]] | tab[ptr[pixel[8]]];
+
+                if( d == 0 )
+                    continue;
+
+                d &= tab[ptr[pixel[2]]] | tab[ptr[pixel[10]]];
+                d &= tab[ptr[pixel[4]]] | tab[ptr[pixel[12]]];
+                d &= tab[ptr[pixel[6]]] | tab[ptr[pixel[14]]];
+
+                if( d == 0 )
+                    continue;
+
+                d &= tab[ptr[pixel[1]]] | tab[ptr[pixel[9]]];
+                d &= tab[ptr[pixel[3]]] | tab[ptr[pixel[11]]];
+                d &= tab[ptr[pixel[5]]] | tab[ptr[pixel[13]]];
+                d &= tab[ptr[pixel[7]]] | tab[ptr[pixel[15]]];
+
+                if( d & 1 )
+                {
+                    int vt = v - threshold, count = 0;
+
+                    for( k = 0; k < N; k++ )
+                    {
+                        int x = ptr[pixel[k]];
+                        if(x < vt)
+                        {
+                            if( ++count > K )
+                            {
+                                cornerpos[ncorners++] = j;
+                                if(nonmax_suppression)
+                                    curr[j] = (uchar)cornerScore<patternSize>(ptr, pixel, threshold);
+                                break;
+                            }
+                        }
+                        else
+                            count = 0;
+                    }
+                }
+
+                if( d & 2 )
+                {
+                    int vt = v + threshold, count = 0;
+
+                    for( k = 0; k < N; k++ )
+                    {
+                        int x = ptr[pixel[k]];
+                        if(x > vt)
+                        {
+                            if( ++count > K )
+                            {
+                                cornerpos[ncorners++] = j;
+                                if(nonmax_suppression)
+                                    curr[j] = (uchar)cornerScore<patternSize>(ptr, pixel, threshold);
+                                break;
+                            }
+                        }
+                        else
+                            count = 0;
+                    }
+                }
+            }
+        }
+
+        cornerpos[-1] = ncorners;
+
+        if( i == 3 )
+            continue;
+
+        const uchar* prev = buf[(i - 4 + 3)%3];
+        const uchar* pprev = buf[(i - 5 + 3)%3];
+        cornerpos = cpbuf[(i - 4 + 3)%3] + 1; // cornerpos[-1] is used to store a value
+        ncorners = cornerpos[-1];
+
+        for( k = 0; k < ncorners; k++ )
+        {
+            j = cornerpos[k];
+            int score = prev[j];
+            if( !nonmax_suppression ||
+               (score > prev[j+1] && score > prev[j-1] &&
+                score > pprev[j-1] && score > pprev[j] && score > pprev[j+1] &&
+                score > curr[j-1] && score > curr[j] && score > curr[j+1]) )
+            {
+                keypoints.push_back(KeyPoint((float)j, (float)(i-1), 7.f, -1, (float)score));
+            }
+        }
+    }
+}
+
+#ifdef HAVE_OPENCL
+template<typename pt>
+struct cmp_pt
+{
+    bool operator ()(const pt& a, const pt& b) const { return a.y < b.y || (a.y == b.y && a.x < b.x); }
+};
+
+static bool ocl_FAST( InputArray _img, std::vector<KeyPoint>& keypoints,
+                     int threshold, bool nonmax_suppression, int maxKeypoints )
+{
+    UMat img = _img.getUMat();
+    if( img.cols < 7 || img.rows < 7 )
+        return false;
+    size_t globalsize[] = { (size_t)img.cols-6, (size_t)img.rows-6 };
+
+    ocl::Kernel fastKptKernel("FAST_findKeypoints", ocl::features::fast_oclsrc);
+    if (fastKptKernel.empty())
+        return false;
+
+    UMat kp1(1, maxKeypoints*2+1, CV_32S);
+
+    UMat ucounter1(kp1, Rect(0,0,1,1));
+    ucounter1.setTo(Scalar::all(0));
+
+    if( !fastKptKernel.args(ocl::KernelArg::ReadOnly(img),
+                            ocl::KernelArg::PtrReadWrite(kp1),
+                            maxKeypoints, threshold).run(2, globalsize, 0, true))
+        return false;
+
+    Mat mcounter;
+    ucounter1.copyTo(mcounter);
+    int i, counter = mcounter.at<int>(0);
+    counter = std::min(counter, maxKeypoints);
+
+    keypoints.clear();
+
+    if( counter == 0 )
+        return true;
+
+    if( !nonmax_suppression )
+    {
+        Mat m;
+        kp1(Rect(0, 0, counter*2+1, 1)).copyTo(m);
+        const Point* pt = (const Point*)(m.ptr<int>() + 1);
+        for( i = 0; i < counter; i++ )
+            keypoints.push_back(KeyPoint((float)pt[i].x, (float)pt[i].y, 7.f, -1, 1.f));
+    }
+    else
+    {
+        UMat kp2(1, maxKeypoints*3+1, CV_32S);
+        UMat ucounter2 = kp2(Rect(0,0,1,1));
+        ucounter2.setTo(Scalar::all(0));
+
+        ocl::Kernel fastNMSKernel("FAST_nonmaxSupression", ocl::features::fast_oclsrc);
+        if (fastNMSKernel.empty())
+            return false;
+
+        size_t globalsize_nms[] = { (size_t)counter };
+        if( !fastNMSKernel.args(ocl::KernelArg::PtrReadOnly(kp1),
+                                ocl::KernelArg::PtrReadWrite(kp2),
+                                ocl::KernelArg::ReadOnly(img),
+                                counter, counter).run(1, globalsize_nms, 0, true))
+            return false;
+
+        Mat m2;
+        kp2(Rect(0, 0, counter*3+1, 1)).copyTo(m2);
+        Point3i* pt2 = (Point3i*)(m2.ptr<int>() + 1);
+        int newcounter = std::min(m2.at<int>(0), counter);
+
+        std::sort(pt2, pt2 + newcounter, cmp_pt<Point3i>());
+
+        for( i = 0; i < newcounter; i++ )
+            keypoints.push_back(KeyPoint((float)pt2[i].x, (float)pt2[i].y, 7.f, -1, (float)pt2[i].z));
+    }
+
+    return true;
+}
+#endif
+
+
+
+static inline int hal_FAST(cv::Mat& src, std::vector<KeyPoint>& keypoints, int threshold, bool nonmax_suppression, FastFeatureDetector::DetectorType type)
+{
+    if (threshold > 20)
+        return CV_HAL_ERROR_NOT_IMPLEMENTED;
+
+    cv::Mat scores(src.size(), src.type());
+
+    int error = cv_hal_FAST_dense(src.data, src.step, scores.data, scores.step, src.cols, src.rows, type);
+
+    if (error != CV_HAL_ERROR_OK)
+        return error;
+
+    cv::Mat suppressedScores(src.size(), src.type());
+
+    if (nonmax_suppression)
+    {
+        error = cv_hal_FAST_NMS(scores.data, scores.step, suppressedScores.data, suppressedScores.step, scores.cols, scores.rows);
+
+        if (error != CV_HAL_ERROR_OK)
+            return error;
+    }
+    else
+    {
+        suppressedScores = scores;
+    }
+
+    if (!threshold && nonmax_suppression) threshold = 1;
+
+    cv::KeyPoint kpt(0, 0, 7.f, -1, 0);
+
+    unsigned uthreshold = (unsigned) threshold;
+
+    int ofs = 3;
+
+    int stride = (int)suppressedScores.step;
+    const unsigned char* pscore = suppressedScores.data;
+
+    keypoints.clear();
+
+    for (int y = ofs; y + ofs < suppressedScores.rows; ++y)
+    {
+        kpt.pt.y = (float)(y);
+        for (int x = ofs; x + ofs < suppressedScores.cols; ++x)
+        {
+            unsigned score = pscore[y * stride + x];
+            if (score > uthreshold)
+            {
+                kpt.pt.x = (float)(x);
+                kpt.response = (nonmax_suppression != 0) ? (float)((int)score - 1) : 0.f;
+                keypoints.push_back(kpt);
+            }
+        }
+    }
+
+    return CV_HAL_ERROR_OK;
+}
+
+void FAST(InputArray _img, std::vector<KeyPoint>& keypoints, int threshold, bool nonmax_suppression, FastFeatureDetector::DetectorType type)
+{
+    CV_INSTRUMENT_REGION();
+
+    CV_OCL_RUN(_img.isUMat() && type == FastFeatureDetector::TYPE_9_16,
+               ocl_FAST(_img, keypoints, threshold, nonmax_suppression, 10000));
+
+    cv::Mat img = _img.getMat();
+    CALL_HAL(fast_dense, hal_FAST, img, keypoints, threshold, nonmax_suppression, type);
+
+    size_t keypoints_count;
+    CALL_HAL(fast, cv_hal_FAST, img.data, img.step, img.cols, img.rows,
+             (uchar*)(keypoints.data()), &keypoints_count, threshold, nonmax_suppression, type);
+
+    switch(type) {
+    case FastFeatureDetector::TYPE_5_8:
+        FAST_t<8>(_img, keypoints, threshold, nonmax_suppression);
+        break;
+    case FastFeatureDetector::TYPE_7_12:
+        FAST_t<12>(_img, keypoints, threshold, nonmax_suppression);
+        break;
+    case FastFeatureDetector::TYPE_9_16:
+        FAST_t<16>(_img, keypoints, threshold, nonmax_suppression);
+        break;
+    }
+}
+
+
+class FastFeatureDetector_Impl CV_FINAL : public FastFeatureDetector
+{
+public:
+    FastFeatureDetector_Impl( int _threshold, bool _nonmaxSuppression, FastFeatureDetector::DetectorType _type )
+    : threshold(_threshold), nonmaxSuppression(_nonmaxSuppression), type(_type)
+    {}
+
+    void read( const FileNode& fn) CV_OVERRIDE
+    {
+      // if node is empty, keep previous value
+      if (!fn["threshold"].empty())
+        fn["threshold"] >> threshold;
+      if (!fn["nonmaxSuppression"].empty())
+        fn["nonmaxSuppression"] >> nonmaxSuppression;
+      if (!fn["type"].empty())
+        fn["type"] >> type;
+    }
+    void write( FileStorage& fs) const CV_OVERRIDE
+    {
+      if(fs.isOpened())
+      {
+        fs << "name" << getDefaultName();
+        fs << "threshold" << threshold;
+        fs << "nonmaxSuppression" << nonmaxSuppression;
+        fs << "type" << type;
+      }
+    }
+
+    void detect( InputArray _image, std::vector<KeyPoint>& keypoints, InputArray _mask ) CV_OVERRIDE
+    {
+        CV_INSTRUMENT_REGION();
+
+        if(_image.empty())
+        {
+            keypoints.clear();
+            return;
+        }
+
+        Mat mask = _mask.getMat(), grayImage;
+        UMat ugrayImage;
+        _InputArray gray = _image;
+        if( _image.type() != CV_8U )
+        {
+            _OutputArray ogray = _image.isUMat() ? _OutputArray(ugrayImage) : _OutputArray(grayImage);
+            cvtColor( _image, ogray, COLOR_BGR2GRAY );
+            gray = ogray;
+        }
+        FAST( gray, keypoints, threshold, nonmaxSuppression, type );
+        KeyPointsFilter::runByPixelsMask( keypoints, mask );
+    }
+
+    void set(int prop, double value)
+    {
+        if(prop == THRESHOLD)
+            threshold = cvRound(value);
+        else if(prop == NONMAX_SUPPRESSION)
+            nonmaxSuppression = value != 0;
+        else if(prop == FAST_N)
+            type = static_cast<FastFeatureDetector::DetectorType>(cvRound(value));
+        else
+            CV_Error(Error::StsBadArg, "");
+    }
+
+    double get(int prop) const
+    {
+        if(prop == THRESHOLD)
+            return threshold;
+        if(prop == NONMAX_SUPPRESSION)
+            return nonmaxSuppression;
+        if(prop == FAST_N)
+            return static_cast<int>(type);
+        CV_Error(Error::StsBadArg, "");
+        return 0;
+    }
+
+    void setThreshold(int threshold_) CV_OVERRIDE { threshold = threshold_; }
+    int getThreshold() const CV_OVERRIDE { return threshold; }
+
+    void setNonmaxSuppression(bool f) CV_OVERRIDE { nonmaxSuppression = f; }
+    bool getNonmaxSuppression() const CV_OVERRIDE { return nonmaxSuppression; }
+
+    void setType(FastFeatureDetector::DetectorType type_) CV_OVERRIDE{ type = type_; }
+    FastFeatureDetector::DetectorType getType() const CV_OVERRIDE{ return type; }
+
+    int threshold;
+    bool nonmaxSuppression;
+    FastFeatureDetector::DetectorType type;
+};
+
+Ptr<FastFeatureDetector> FastFeatureDetector::create( int threshold, bool nonmaxSuppression, FastFeatureDetector::DetectorType type )
+{
+    return makePtr<FastFeatureDetector_Impl>(threshold, nonmaxSuppression, type);
+}
+
+String FastFeatureDetector::getDefaultName() const
+{
+    return (Feature2D::getDefaultName() + ".FastFeatureDetector");
+}
+
+}

modules/features/src/fast.hpp

@@ -0,0 +1,62 @@
+/* This is FAST corner detector, contributed to OpenCV by the author, Edward Rosten.
+   Below is the original copyright and the references */
+
+/*
+Copyright (c) 2006, 2008 Edward Rosten
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+    *Redistributions of source code must retain the above copyright
+     notice, this list of conditions and the following disclaimer.
+
+    *Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+
+    *Neither the name of the University of Cambridge nor the names of
+     its contributors may be used to endorse or promote products derived
+     from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+The references are:
+ * Machine learning for high-speed corner detection,
+   E. Rosten and T. Drummond, ECCV 2006
+ * Faster and better: A machine learning approach to corner detection
+   E. Rosten, R. Porter and T. Drummond, PAMI, 2009
+*/
+
+#ifndef OPENCV_FEATURES_FAST_HPP
+#define OPENCV_FEATURES_FAST_HPP
+
+namespace cv
+{
+namespace opt_AVX2
+{
+#if CV_TRY_AVX2
+class FAST_t_patternSize16_AVX2
+{
+public:
+    static Ptr<FAST_t_patternSize16_AVX2> getImpl(int _cols, int _threshold, bool _nonmax_suppression, const int* _pixel);
+    virtual void process(int &j, const uchar* &ptr, uchar* curr, int* cornerpos, int &ncorners) = 0;
+    virtual ~FAST_t_patternSize16_AVX2() {}
+};
+#endif
+}
+}
+#endif

modules/features/src/fast_score.cpp

@@ -0,0 +1,366 @@
+/* This is FAST corner detector, contributed to OpenCV by the author, Edward Rosten.
+   Below is the original copyright and the references */
+
+/*
+Copyright (c) 2006, 2008 Edward Rosten
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+    *Redistributions of source code must retain the above copyright
+     notice, this list of conditions and the following disclaimer.
+
+    *Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+
+    *Neither the name of the University of Cambridge nor the names of
+     its contributors may be used to endorse or promote products derived
+     from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+The references are:
+ * Machine learning for high-speed corner detection,
+   E. Rosten and T. Drummond, ECCV 2006
+ * Faster and better: A machine learning approach to corner detection
+   E. Rosten, R. Porter and T. Drummond, PAMI, 2009
+*/
+
+#include "fast_score.hpp"
+#include "opencv2/core/hal/intrin.hpp"
+#define VERIFY_CORNERS 0
+
+namespace cv {
+
+void makeOffsets(int pixel[25], int rowStride, int patternSize)
+{
+    static const int offsets16[][2] =
+    {
+        {0,  3}, { 1,  3}, { 2,  2}, { 3,  1}, { 3, 0}, { 3, -1}, { 2, -2}, { 1, -3},
+        {0, -3}, {-1, -3}, {-2, -2}, {-3, -1}, {-3, 0}, {-3,  1}, {-2,  2}, {-1,  3}
+    };
+
+    static const int offsets12[][2] =
+    {
+        {0,  2}, { 1,  2}, { 2,  1}, { 2, 0}, { 2, -1}, { 1, -2},
+        {0, -2}, {-1, -2}, {-2, -1}, {-2, 0}, {-2,  1}, {-1,  2}
+    };
+
+    static const int offsets8[][2] =
+    {
+        {0,  1}, { 1,  1}, { 1, 0}, { 1, -1},
+        {0, -1}, {-1, -1}, {-1, 0}, {-1,  1}
+    };
+
+    const int (*offsets)[2] = patternSize == 16 ? offsets16 :
+                              patternSize == 12 ? offsets12 :
+                              patternSize == 8  ? offsets8  : 0;
+
+    CV_Assert(pixel && offsets);
+
+    int k = 0;
+    for( ; k < patternSize; k++ )
+        pixel[k] = offsets[k][0] + offsets[k][1] * rowStride;
+    for( ; k < 25; k++ )
+        pixel[k] = pixel[k - patternSize];
+}
+
+#if VERIFY_CORNERS
+static void testCorner(const uchar* ptr, const int pixel[], int K, int N, int threshold) {
+    // check that with the computed "threshold" the pixel is still a corner
+    // and that with the increased-by-1 "threshold" the pixel is not a corner anymore
+    for( int delta = 0; delta <= 1; delta++ )
+    {
+        int v0 = std::min(ptr[0] + threshold + delta, 255);
+        int v1 = std::max(ptr[0] - threshold - delta, 0);
+        int c0 = 0, c1 = 0;
+
+        for( int k = 0; k < N; k++ )
+        {
+            int x = ptr[pixel[k]];
+            if(x > v0)
+            {
+                if( ++c0 > K )
+                    break;
+                c1 = 0;
+            }
+            else if( x < v1 )
+            {
+                if( ++c1 > K )
+                    break;
+                c0 = 0;
+            }
+            else
+            {
+                c0 = c1 = 0;
+            }
+        }
+        CV_Assert( (delta == 0 && std::max(c0, c1) > K) ||
+                   (delta == 1 && std::max(c0, c1) <= K) );
+    }
+}
+#endif
+
+template<>
+int cornerScore<16>(const uchar* ptr, const int pixel[], int threshold)
+{
+    const int K = 8, N = K*3 + 1;
+    int k, v = ptr[0];
+    short d[N];
+    for( k = 0; k < N; k++ )
+        d[k] = (short)(v - ptr[pixel[k]]);
+
+#if CV_SIMD128
+    if (true)
+    {
+        v_int16x8 q0 = v_setall_s16(-1000), q1 = v_setall_s16(1000);
+        for (k = 0; k < 16; k += 8)
+        {
+            v_int16x8 v0 = v_load(d + k + 1);
+            v_int16x8 v1 = v_load(d + k + 2);
+            v_int16x8 a = v_min(v0, v1);
+            v_int16x8 b = v_max(v0, v1);
+            v0 = v_load(d + k + 3);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 4);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 5);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 6);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 7);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 8);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k);
+            q0 = v_max(q0, v_min(a, v0));
+            q1 = v_min(q1, v_max(b, v0));
+            v0 = v_load(d + k + 9);
+            q0 = v_max(q0, v_min(a, v0));
+            q1 = v_min(q1, v_max(b, v0));
+        }
+        q0 = v_max(q0, v_sub(v_setzero_s16(), q1));
+        threshold = v_reduce_max(q0) - 1;
+    }
+    else
+#endif
+    {
+
+        int a0 = threshold;
+        for( k = 0; k < 16; k += 2 )
+        {
+            int a = std::min((int)d[k+1], (int)d[k+2]);
+            a = std::min(a, (int)d[k+3]);
+            if( a <= a0 )
+                continue;
+            a = std::min(a, (int)d[k+4]);
+            a = std::min(a, (int)d[k+5]);
+            a = std::min(a, (int)d[k+6]);
+            a = std::min(a, (int)d[k+7]);
+            a = std::min(a, (int)d[k+8]);
+            a0 = std::max(a0, std::min(a, (int)d[k]));
+            a0 = std::max(a0, std::min(a, (int)d[k+9]));
+        }
+
+        int b0 = -a0;
+        for( k = 0; k < 16; k += 2 )
+        {
+            int b = std::max((int)d[k+1], (int)d[k+2]);
+            b = std::max(b, (int)d[k+3]);
+            b = std::max(b, (int)d[k+4]);
+            b = std::max(b, (int)d[k+5]);
+            if( b >= b0 )
+                continue;
+            b = std::max(b, (int)d[k+6]);
+            b = std::max(b, (int)d[k+7]);
+            b = std::max(b, (int)d[k+8]);
+
+            b0 = std::min(b0, std::max(b, (int)d[k]));
+            b0 = std::min(b0, std::max(b, (int)d[k+9]));
+        }
+
+        threshold = -b0 - 1;
+    }
+
+#if VERIFY_CORNERS
+    testCorner(ptr, pixel, K, N, threshold);
+#endif
+    return threshold;
+}
+
+template<>
+int cornerScore<12>(const uchar* ptr, const int pixel[], int threshold)
+{
+    const int K = 6, N = K*3 + 1;
+    int k, v = ptr[0];
+    short d[N + 4];
+    for( k = 0; k < N; k++ )
+        d[k] = (short)(v - ptr[pixel[k]]);
+#if CV_SIMD128
+    for( k = 0; k < 4; k++ )
+        d[N+k] = d[k];
+#endif
+
+#if CV_SIMD128
+    if (true)
+    {
+        v_int16x8 q0 = v_setall_s16(-1000), q1 = v_setall_s16(1000);
+        for (k = 0; k < 16; k += 8)
+        {
+            v_int16x8 v0 = v_load(d + k + 1);
+            v_int16x8 v1 = v_load(d + k + 2);
+            v_int16x8 a = v_min(v0, v1);
+            v_int16x8 b = v_max(v0, v1);
+            v0 = v_load(d + k + 3);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 4);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 5);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k + 6);
+            a = v_min(a, v0);
+            b = v_max(b, v0);
+            v0 = v_load(d + k);
+            q0 = v_max(q0, v_min(a, v0));
+            q1 = v_min(q1, v_max(b, v0));
+            v0 = v_load(d + k + 7);
+            q0 = v_max(q0, v_min(a, v0));
+            q1 = v_min(q1, v_max(b, v0));
+        }
+        q0 = v_max(q0, v_sub(v_setzero_s16(), q1));
+        threshold = v_reduce_max(q0) - 1;
+    }
+    else
+#endif
+    {
+        int a0 = threshold;
+        for( k = 0; k < 12; k += 2 )
+        {
+            int a = std::min((int)d[k+1], (int)d[k+2]);
+            if( a <= a0 )
+                continue;
+            a = std::min(a, (int)d[k+3]);
+            a = std::min(a, (int)d[k+4]);
+            a = std::min(a, (int)d[k+5]);
+            a = std::min(a, (int)d[k+6]);
+            a0 = std::max(a0, std::min(a, (int)d[k]));
+            a0 = std::max(a0, std::min(a, (int)d[k+7]));
+        }
+
+        int b0 = -a0;
+        for( k = 0; k < 12; k += 2 )
+        {
+            int b = std::max((int)d[k+1], (int)d[k+2]);
+            b = std::max(b, (int)d[k+3]);
+            b = std::max(b, (int)d[k+4]);
+            if( b >= b0 )
+                continue;
+            b = std::max(b, (int)d[k+5]);
+            b = std::max(b, (int)d[k+6]);
+
+            b0 = std::min(b0, std::max(b, (int)d[k]));
+            b0 = std::min(b0, std::max(b, (int)d[k+7]));
+        }
+
+        threshold = -b0-1;
+    }
+#if VERIFY_CORNERS
+    testCorner(ptr, pixel, K, N, threshold);
+#endif
+    return threshold;
+}
+
+template<>
+int cornerScore<8>(const uchar* ptr, const int pixel[], int threshold)
+{
+    const int K = 4, N = K * 3 + 1;
+    int k, v = ptr[0];
+    short d[N];
+    for (k = 0; k < N; k++)
+        d[k] = (short)(v - ptr[pixel[k]]);
+
+#if CV_SIMD128 \
+    && (!defined(CV_SIMD128_CPP) || (!defined(__GNUC__) || __GNUC__ != 5))  // "movdqa" bug on "v_load(d + 1)" line (Ubuntu 16.04 + GCC 5.4)
+    if (true)
+    {
+        v_int16x8 v0 = v_load(d + 1);
+        v_int16x8 v1 = v_load(d + 2);
+        v_int16x8 a = v_min(v0, v1);
+        v_int16x8 b = v_max(v0, v1);
+        v0 = v_load(d + 3);
+        a = v_min(a, v0);
+        b = v_max(b, v0);
+        v0 = v_load(d + 4);
+        a = v_min(a, v0);
+        b = v_max(b, v0);
+        v0 = v_load(d);
+        v_int16x8 q0 = v_min(a, v0);
+        v_int16x8 q1 = v_max(b, v0);
+        v0 = v_load(d + 5);
+        q0 = v_max(q0, v_min(a, v0));
+        q1 = v_min(q1, v_max(b, v0));
+        q0 = v_max(q0, v_sub(v_setzero_s16(), q1));
+        threshold = v_reduce_max(q0) - 1;
+    }
+    else
+#endif
+    {
+        int a0 = threshold;
+        for( k = 0; k < 8; k += 2 )
+        {
+            int a = std::min((int)d[k+1], (int)d[k+2]);
+            if( a <= a0 )
+                continue;
+            a = std::min(a, (int)d[k+3]);
+            a = std::min(a, (int)d[k+4]);
+            a0 = std::max(a0, std::min(a, (int)d[k]));
+            a0 = std::max(a0, std::min(a, (int)d[k+5]));
+        }
+
+        int b0 = -a0;
+        for( k = 0; k < 8; k += 2 )
+        {
+            int b = std::max((int)d[k+1], (int)d[k+2]);
+            b = std::max(b, (int)d[k+3]);
+            if( b >= b0 )
+                continue;
+            b = std::max(b, (int)d[k+4]);
+
+            b0 = std::min(b0, std::max(b, (int)d[k]));
+            b0 = std::min(b0, std::max(b, (int)d[k+5]));
+        }
+
+        threshold = -b0-1;
+    }
+
+#if VERIFY_CORNERS
+    testCorner(ptr, pixel, K, N, threshold);
+#endif
+    return threshold;
+}
+
+} // namespace cv

modules/features/src/fast_score.hpp

@@ -0,0 +1,62 @@
+/* This is FAST corner detector, contributed to OpenCV by the author, Edward Rosten.
+   Below is the original copyright and the references */
+
+/*
+Copyright (c) 2006, 2008 Edward Rosten
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+    *Redistributions of source code must retain the above copyright
+     notice, this list of conditions and the following disclaimer.
+
+    *Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+
+    *Neither the name of the University of Cambridge nor the names of
+     its contributors may be used to endorse or promote products derived
+     from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+The references are:
+ * Machine learning for high-speed corner detection,
+   E. Rosten and T. Drummond, ECCV 2006
+ * Faster and better: A machine learning approach to corner detection
+   E. Rosten, R. Porter and T. Drummond, PAMI, 2009
+*/
+
+#ifndef __OPENCV_FEATURES_2D_FAST_HPP__
+#define __OPENCV_FEATURES_2D_FAST_HPP__
+
+#ifdef __cplusplus
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+void makeOffsets(int pixel[25], int row_stride, int patternSize);
+
+template<int patternSize>
+int cornerScore(const uchar* ptr, const int pixel[], int threshold);
+
+}
+
+#endif
+#endif

modules/features/src/feature2d.cpp

@@ -0,0 +1,224 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                          License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+using std::vector;
+
+Feature2D::~Feature2D() {}
+
+/*
+ * Detect keypoints in an image.
+ * image        The image.
+ * keypoints    The detected keypoints.
+ * mask         Mask specifying where to look for keypoints (optional). Must be a char
+ *              matrix with non-zero values in the region of interest.
+ */
+void Feature2D::detect( InputArray image,
+                        std::vector<KeyPoint>& keypoints,
+                        InputArray mask )
+{
+    CV_INSTRUMENT_REGION();
+
+    if( image.empty() )
+    {
+        keypoints.clear();
+        return;
+    }
+    detectAndCompute(image, mask, keypoints, noArray(), false);
+}
+
+
+void Feature2D::detect( InputArrayOfArrays images,
+                        std::vector<std::vector<KeyPoint> >& keypoints,
+                        InputArrayOfArrays masks )
+{
+    CV_INSTRUMENT_REGION();
+
+    int nimages = (int)images.total();
+
+    if (!masks.empty())
+    {
+        CV_Assert(masks.total() == (size_t)nimages);
+    }
+
+    keypoints.resize(nimages);
+
+    if (images.isMatVector())
+    {
+       for (int i = 0; i < nimages; i++)
+       {
+           detect(images.getMat(i), keypoints[i], masks.empty() ? noArray() : masks.getMat(i));
+       }
+    }
+    else
+    {
+        // assume UMats
+        for (int i = 0; i < nimages; i++)
+        {
+            detect(images.getUMat(i), keypoints[i], masks.empty() ? noArray() : masks.getUMat(i));
+        }
+    }
+
+
+}
+
+/*
+ * Compute the descriptors for a set of keypoints in an image.
+ * image        The image.
+ * keypoints    The input keypoints. Keypoints for which a descriptor cannot be computed are removed.
+ * descriptors  Copmputed descriptors. Row i is the descriptor for keypoint i.
+ */
+void Feature2D::compute( InputArray image,
+                         std::vector<KeyPoint>& keypoints,
+                         OutputArray descriptors )
+{
+    CV_INSTRUMENT_REGION();
+
+    if( image.empty() )
+    {
+        descriptors.release();
+        return;
+    }
+    detectAndCompute(image, noArray(), keypoints, descriptors, true);
+}
+
+void Feature2D::compute( InputArrayOfArrays images,
+                         std::vector<std::vector<KeyPoint> >& keypoints,
+                         OutputArrayOfArrays descriptors )
+{
+    CV_INSTRUMENT_REGION();
+
+    if( !descriptors.needed() )
+        return;
+
+    int nimages = (int)images.total();
+
+    CV_Assert( keypoints.size() == (size_t)nimages );
+    // resize descriptors to appropriate size and compute
+    if (descriptors.isMatVector())
+    {
+        vector<Mat>& vec = *(vector<Mat>*)descriptors.getObj();
+        vec.resize(nimages);
+        for (int i = 0; i < nimages; i++)
+        {
+            compute(images.getMat(i), keypoints[i], vec[i]);
+        }
+    }
+    else if (descriptors.isUMatVector())
+    {
+        vector<UMat>& vec = *(vector<UMat>*)descriptors.getObj();
+        vec.resize(nimages);
+        for (int i = 0; i < nimages; i++)
+        {
+            compute(images.getUMat(i), keypoints[i], vec[i]);
+        }
+    }
+    else
+    {
+        CV_Error(Error::StsBadArg, "descriptors must be vector<Mat> or vector<UMat>");
+    }
+}
+
+
+/* Detects keypoints and computes the descriptors */
+void Feature2D::detectAndCompute( InputArray, InputArray,
+                                  std::vector<KeyPoint>&,
+                                  OutputArray,
+                                  bool )
+{
+    CV_INSTRUMENT_REGION();
+
+    CV_Error(Error::StsNotImplemented, "");
+}
+
+void Feature2D::write( const String& fileName ) const
+{
+    FileStorage fs(fileName, FileStorage::WRITE);
+    write(fs);
+}
+
+void Feature2D::read( const String& fileName )
+{
+    FileStorage fs(fileName, FileStorage::READ);
+    read(fs.root());
+}
+
+void Feature2D::write( FileStorage&) const
+{
+}
+
+void Feature2D::read( const FileNode&)
+{
+}
+
+int Feature2D::descriptorSize() const
+{
+    return 0;
+}
+
+int Feature2D::descriptorType() const
+{
+    return CV_32F;
+}
+
+int Feature2D::defaultNorm() const
+{
+    int tp = descriptorType();
+    return tp == CV_8U ? NORM_HAMMING : NORM_L2;
+}
+
+// Return true if detector object is empty
+bool Feature2D::empty() const
+{
+    return true;
+}
+
+String Feature2D::getDefaultName() const
+{
+    return "Feature2D";
+}
+
+}

modules/features/src/gftt.cpp

@@ -0,0 +1,185 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+class GFTTDetector_Impl CV_FINAL : public GFTTDetector
+{
+public:
+    GFTTDetector_Impl( int _nfeatures, double _qualityLevel,
+                      double _minDistance, int _blockSize, int _gradientSize,
+                      bool _useHarrisDetector, double _k )
+        : nfeatures(_nfeatures), qualityLevel(_qualityLevel), minDistance(_minDistance),
+        blockSize(_blockSize), gradSize(_gradientSize), useHarrisDetector(_useHarrisDetector), k(_k)
+    {
+    }
+
+    void read( const FileNode& fn) CV_OVERRIDE
+    {
+      // if node is empty, keep previous value
+      if (!fn["nfeatures"].empty())
+        fn["nfeatures"] >> nfeatures;
+      if (!fn["qualityLevel"].empty())
+        fn["qualityLevel"] >> qualityLevel;
+      if (!fn["minDistance"].empty())
+        fn["minDistance"] >> minDistance;
+      if (!fn["blockSize"].empty())
+        fn["blockSize"] >> blockSize;
+      if (!fn["gradSize"].empty())
+        fn["gradSize"] >> gradSize;
+      if (!fn["useHarrisDetector"].empty())
+        fn["useHarrisDetector"] >> useHarrisDetector;
+      if (!fn["k"].empty())
+        fn["k"] >> k;
+    }
+    void write( FileStorage& fs) const CV_OVERRIDE
+    {
+      if(fs.isOpened())
+      {
+        fs << "name" << getDefaultName();
+        fs << "nfeatures" << nfeatures;
+        fs << "qualityLevel" << qualityLevel;
+        fs << "minDistance" << minDistance;
+        fs << "blockSize" << blockSize;
+        fs << "gradSize" << gradSize;
+        fs << "useHarrisDetector" << useHarrisDetector;
+        fs << "k" << k;
+      }
+    }
+
+    void setMaxFeatures(int maxFeatures) CV_OVERRIDE { nfeatures = maxFeatures; }
+    int getMaxFeatures() const CV_OVERRIDE { return nfeatures; }
+
+    void setQualityLevel(double qlevel) CV_OVERRIDE { qualityLevel = qlevel; }
+    double getQualityLevel() const CV_OVERRIDE { return qualityLevel; }
+
+    void setMinDistance(double minDistance_) CV_OVERRIDE { minDistance = minDistance_; }
+    double getMinDistance() const CV_OVERRIDE { return minDistance; }
+
+    void setBlockSize(int blockSize_) CV_OVERRIDE { blockSize = blockSize_; }
+    int getBlockSize() const CV_OVERRIDE { return blockSize; }
+
+    void setGradientSize(int gradientSize_) CV_OVERRIDE { gradSize = gradientSize_; }
+    int getGradientSize() CV_OVERRIDE { return gradSize; }
+
+    void setHarrisDetector(bool val) CV_OVERRIDE { useHarrisDetector = val; }
+    bool getHarrisDetector() const CV_OVERRIDE { return useHarrisDetector; }
+
+    void setK(double k_) CV_OVERRIDE { k = k_; }
+    double getK() const CV_OVERRIDE { return k; }
+
+    void detect( InputArray _image, std::vector<KeyPoint>& keypoints, InputArray _mask ) CV_OVERRIDE
+    {
+        CV_INSTRUMENT_REGION();
+
+        if(_image.empty())
+        {
+            keypoints.clear();
+            return;
+        }
+
+        std::vector<Point2f> corners;
+        std::vector<float> cornersQuality;
+
+        if (_image.isUMat())
+        {
+            UMat ugrayImage;
+            if( _image.type() != CV_8U )
+                cvtColor( _image, ugrayImage, COLOR_BGR2GRAY );
+            else
+                ugrayImage = _image.getUMat();
+
+            goodFeaturesToTrack( ugrayImage, corners, nfeatures, qualityLevel, minDistance, _mask,
+                                 cornersQuality, blockSize, gradSize, useHarrisDetector, k );
+        }
+        else
+        {
+            Mat image = _image.getMat(), grayImage = image;
+            if( image.type() != CV_8U )
+                cvtColor( image, grayImage, COLOR_BGR2GRAY );
+
+            goodFeaturesToTrack( grayImage, corners, nfeatures, qualityLevel, minDistance, _mask,
+                                 cornersQuality, blockSize, gradSize, useHarrisDetector, k );
+        }
+
+        CV_Assert(corners.size() == cornersQuality.size());
+
+        keypoints.resize(corners.size());
+        for (size_t i = 0; i < corners.size(); i++)
+            keypoints[i] = KeyPoint(corners[i], (float)blockSize, -1, cornersQuality[i]);
+
+    }
+
+    int nfeatures;
+    double qualityLevel;
+    double minDistance;
+    int blockSize;
+    int gradSize;
+    bool useHarrisDetector;
+    double k;
+};
+
+
+Ptr<GFTTDetector> GFTTDetector::create( int _nfeatures, double _qualityLevel,
+                         double _minDistance, int _blockSize, int _gradientSize,
+                         bool _useHarrisDetector, double _k )
+{
+    return makePtr<GFTTDetector_Impl>(_nfeatures, _qualityLevel,
+                                      _minDistance, _blockSize, _gradientSize, _useHarrisDetector, _k);
+}
+
+Ptr<GFTTDetector> GFTTDetector::create( int _nfeatures, double _qualityLevel,
+                         double _minDistance, int _blockSize,
+                         bool _useHarrisDetector, double _k )
+{
+    return makePtr<GFTTDetector_Impl>(_nfeatures, _qualityLevel,
+                                      _minDistance, _blockSize, 3, _useHarrisDetector, _k);
+}
+
+String GFTTDetector::getDefaultName() const
+{
+    return (Feature2D::getDefaultName() + ".GFTTDetector");
+}
+
+}

modules/features/src/hal_replacement.hpp

@@ -0,0 +1,145 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2017, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef OPENCV_FEATURES_HAL_REPLACEMENT_HPP
+#define OPENCV_FEATURES_HAL_REPLACEMENT_HPP
+
+#include "opencv2/core/hal/interface.h"
+
+#if defined(__clang__)  // clang or MSVC clang
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wunused-parameter"
+#elif defined(_MSC_VER)
+#pragma warning(push)
+#pragma warning(disable : 4100)
+#elif defined(__GNUC__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+#endif
+
+//! @addtogroup features_hal_interface
+//! @note Define your functions to override default implementations:
+//! @code
+//! #undef hal_add8u
+//! #define hal_add8u my_add8u
+//! @endcode
+//! @{
+/**
+   @brief Detects corners using the FAST algorithm, returns mask.
+   @param src_data Source image data
+   @param src_step Source image step
+   @param dst_data Destination mask data
+   @param dst_step Destination mask step
+   @param width Source image width
+   @param height Source image height
+   @param type FAST type
+*/
+inline int hal_ni_FAST_dense(const uchar* src_data, size_t src_step, uchar* dst_data, size_t dst_step, int width, int height, cv::FastFeatureDetector::DetectorType type) { return CV_HAL_ERROR_NOT_IMPLEMENTED; }
+
+//! @cond IGNORED
+#define cv_hal_FAST_dense hal_ni_FAST_dense
+//! @endcond
+
+/**
+   @brief Non-maximum suppression for FAST_9_16.
+   @param src_data,src_step Source mask
+   @param dst_data,dst_step Destination mask after NMS
+   @param width,height Source mask dimensions
+*/
+inline int hal_ni_FAST_NMS(const uchar* src_data, size_t src_step, uchar* dst_data, size_t dst_step, int width, int height) { return CV_HAL_ERROR_NOT_IMPLEMENTED; }
+
+//! @cond IGNORED
+#define cv_hal_FAST_NMS hal_ni_FAST_NMS
+//! @endcond
+
+/**
+   @brief Detects corners using the FAST algorithm.
+   @param src_data Source image data
+   @param src_step Source image step
+   @param width Source image width
+   @param height Source image height
+   @param keypoints_data Pointer to keypoints
+   @param keypoints_count Count of keypoints
+   @param threshold Threshold for keypoint
+   @param nonmax_suppression Indicates if make nonmaxima suppression or not.
+   @param type FAST type
+*/
+inline int hal_ni_FAST(const uchar* src_data, size_t src_step, int width, int height, uchar* keypoints_data, size_t* keypoints_count, int threshold, bool nonmax_suppression, int /*cv::FastFeatureDetector::DetectorType*/ type) { return CV_HAL_ERROR_NOT_IMPLEMENTED; }
+
+//! @cond IGNORED
+#define cv_hal_FAST hal_ni_FAST
+//! @endcond
+
+//! @}
+
+
+#if defined(__clang__)
+#pragma clang diagnostic pop
+#elif defined(_MSC_VER)
+#pragma warning(pop)
+#elif defined(__GNUC__)
+#pragma GCC diagnostic pop
+#endif
+
+#include "custom_hal.hpp"
+
+//! @cond IGNORED
+#define CALL_HAL_RET(name, fun, retval, ...) \
+    int res = __CV_EXPAND(fun(__VA_ARGS__, &retval)); \
+    if (res == CV_HAL_ERROR_OK) \
+        return retval; \
+    else if (res != CV_HAL_ERROR_NOT_IMPLEMENTED) \
+        CV_Error_(cv::Error::StsInternal, \
+            ("HAL implementation " CVAUX_STR(name) " ==> " CVAUX_STR(fun) " returned %d (0x%08x)", res, res));
+
+
+#define CALL_HAL(name, fun, ...) \
+{                                           \
+    int res = __CV_EXPAND(fun(__VA_ARGS__)); \
+    if (res == CV_HAL_ERROR_OK) \
+        return; \
+    else if (res != CV_HAL_ERROR_NOT_IMPLEMENTED) \
+        CV_Error_(cv::Error::StsInternal, \
+            ("HAL implementation " CVAUX_STR(name) " ==> " CVAUX_STR(fun) " returned %d (0x%08x)", res, res)); \
+}
+//! @endcond
+
+#endif

modules/features/src/keypoint.cpp

@@ -0,0 +1,293 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2008, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+struct KeypointResponseGreaterThanOrEqualToThreshold
+{
+    KeypointResponseGreaterThanOrEqualToThreshold(float _value) :
+    value(_value)
+    {
+    }
+    inline bool operator()(const KeyPoint& kpt) const
+    {
+        return kpt.response >= value;
+    }
+    float value;
+};
+
+struct KeypointResponseGreater
+{
+    inline bool operator()(const KeyPoint& kp1, const KeyPoint& kp2) const
+    {
+        return kp1.response > kp2.response;
+    }
+};
+
+// takes keypoints and culls them by the response
+void KeyPointsFilter::retainBest(std::vector<KeyPoint>& keypoints, int n_points)
+{
+    //this is only necessary if the keypoints size is greater than the number of desired points.
+    if( n_points >= 0 && keypoints.size() > (size_t)n_points )
+    {
+        if (n_points==0)
+        {
+            keypoints.clear();
+            return;
+        }
+        //first use nth element to partition the keypoints into the best and worst.
+        std::nth_element(keypoints.begin(), keypoints.begin() + n_points - 1, keypoints.end(), KeypointResponseGreater());
+        //this is the boundary response, and in the case of FAST may be ambiguous
+        float ambiguous_response = keypoints[n_points - 1].response;
+        //use std::partition to grab all of the keypoints with the boundary response.
+        std::vector<KeyPoint>::const_iterator new_end =
+        std::partition(keypoints.begin() + n_points, keypoints.end(),
+                       KeypointResponseGreaterThanOrEqualToThreshold(ambiguous_response));
+        //resize the keypoints, given this new end point. nth_element and partition reordered the points inplace
+        keypoints.resize(new_end - keypoints.begin());
+    }
+}
+
+struct RoiPredicate
+{
+    RoiPredicate( const Rect& _r ) : r(_r)
+    {}
+
+    bool operator()( const KeyPoint& keyPt ) const
+    {
+        // workaround for https://github.com/opencv/opencv/issues/26016
+        // To keep its behaviour, keyPt.pt casts to Point_<int>.
+        return !r.contains( Point_<int>(keyPt.pt) );
+    }
+
+    Rect r;
+};
+
+void KeyPointsFilter::runByImageBorder( std::vector<KeyPoint>& keypoints, Size imageSize, int borderSize )
+{
+    if( borderSize > 0)
+    {
+        if (imageSize.height <= borderSize * 2 || imageSize.width <= borderSize * 2)
+            keypoints.clear();
+        else
+            keypoints.erase( std::remove_if(keypoints.begin(), keypoints.end(),
+                                       RoiPredicate(Rect(Point(borderSize, borderSize),
+                                                         Point(imageSize.width - borderSize, imageSize.height - borderSize)))),
+                             keypoints.end() );
+    }
+}
+
+struct SizePredicate
+{
+    SizePredicate( float _minSize, float _maxSize ) : minSize(_minSize), maxSize(_maxSize)
+    {}
+
+    bool operator()( const KeyPoint& keyPt ) const
+    {
+        float size = keyPt.size;
+        return (size < minSize) || (size > maxSize);
+    }
+
+    float minSize, maxSize;
+};
+
+void KeyPointsFilter::runByKeypointSize( std::vector<KeyPoint>& keypoints, float minSize, float maxSize )
+{
+    CV_Assert( minSize >= 0 );
+    CV_Assert( maxSize >= 0);
+    CV_Assert( minSize <= maxSize );
+
+    keypoints.erase( std::remove_if(keypoints.begin(), keypoints.end(), SizePredicate(minSize, maxSize)),
+                     keypoints.end() );
+}
+
+class MaskPredicate
+{
+public:
+    MaskPredicate( const Mat& _mask ) : mask(_mask) {}
+    bool operator() (const KeyPoint& key_pt) const
+    {
+        return mask.at<uchar>( (int)(key_pt.pt.y + 0.5f), (int)(key_pt.pt.x + 0.5f) ) == 0;
+    }
+    MaskPredicate& operator=(const MaskPredicate&) = delete;
+    // To avoid -Wdeprecated-copy warning, copy constructor is needed.
+    MaskPredicate(const MaskPredicate&) = default;
+
+private:
+    const Mat mask;
+};
+
+void KeyPointsFilter::runByPixelsMask( std::vector<KeyPoint>& keypoints, const Mat& mask )
+{
+    CV_INSTRUMENT_REGION();
+
+    if( mask.empty() )
+        return;
+
+    keypoints.erase(std::remove_if(keypoints.begin(), keypoints.end(), MaskPredicate(mask)), keypoints.end());
+}
+/*
+ * Remove objects from some image and a vector by mask for pixels of this image
+ */
+template <typename T>
+void runByPixelsMask2(std::vector<KeyPoint> &keypoints, std::vector<T> &removeFrom, const Mat &mask)
+{
+    if (mask.empty())
+        return;
+
+    MaskPredicate maskPredicate(mask);
+    removeFrom.erase(std::remove_if(removeFrom.begin(), removeFrom.end(),
+                                    [&](const T &x)
+                                    {
+                                        auto index = &x - &removeFrom.front();
+                                        return maskPredicate(keypoints[index]);
+                                    }),
+                    removeFrom.end());
+    keypoints.erase(std::remove_if(keypoints.begin(), keypoints.end(), maskPredicate), keypoints.end());
+}
+void KeyPointsFilter::runByPixelsMask2VectorPoint(std::vector<KeyPoint> &keypoints, std::vector<std::vector<Point> > &removeFrom, const Mat &mask)
+{
+    runByPixelsMask2(keypoints, removeFrom, mask);
+}
+
+struct KeyPoint_LessThan
+{
+    KeyPoint_LessThan(const std::vector<KeyPoint>& _kp) : kp(&_kp) {}
+    bool operator()(int i, int j) const
+    {
+        const KeyPoint& kp1 = (*kp)[i];
+        const KeyPoint& kp2 = (*kp)[j];
+        if( kp1.pt.x != kp2.pt.x )
+            return kp1.pt.x < kp2.pt.x;
+        if( kp1.pt.y != kp2.pt.y )
+            return kp1.pt.y < kp2.pt.y;
+        if( kp1.size != kp2.size )
+            return kp1.size > kp2.size;
+        if( kp1.angle != kp2.angle )
+            return kp1.angle < kp2.angle;
+        if( kp1.response != kp2.response )
+            return kp1.response > kp2.response;
+        if( kp1.octave != kp2.octave )
+            return kp1.octave > kp2.octave;
+        if( kp1.class_id != kp2.class_id )
+            return kp1.class_id > kp2.class_id;
+
+        return i < j;
+    }
+    const std::vector<KeyPoint>* kp;
+};
+
+void KeyPointsFilter::removeDuplicated( std::vector<KeyPoint>& keypoints )
+{
+    int i, j, n = (int)keypoints.size();
+    std::vector<int> kpidx(n);
+    std::vector<uchar> mask(n, (uchar)1);
+
+    for( i = 0; i < n; i++ )
+        kpidx[i] = i;
+    std::sort(kpidx.begin(), kpidx.end(), KeyPoint_LessThan(keypoints));
+    for( i = 1, j = 0; i < n; i++ )
+    {
+        KeyPoint& kp1 = keypoints[kpidx[i]];
+        KeyPoint& kp2 = keypoints[kpidx[j]];
+        if( kp1.pt.x != kp2.pt.x || kp1.pt.y != kp2.pt.y ||
+            kp1.size != kp2.size || kp1.angle != kp2.angle )
+            j = i;
+        else
+            mask[kpidx[i]] = 0;
+    }
+
+    for( i = j = 0; i < n; i++ )
+    {
+        if( mask[i] )
+        {
+            if( i != j )
+                keypoints[j] = keypoints[i];
+            j++;
+        }
+    }
+    keypoints.resize(j);
+}
+
+struct KeyPoint12_LessThan
+{
+    bool operator()(const KeyPoint &kp1, const KeyPoint &kp2) const
+    {
+        if( kp1.pt.x != kp2.pt.x )
+            return kp1.pt.x < kp2.pt.x;
+        if( kp1.pt.y != kp2.pt.y )
+            return kp1.pt.y < kp2.pt.y;
+        if( kp1.size != kp2.size )
+            return kp1.size > kp2.size;
+        if( kp1.angle != kp2.angle )
+            return kp1.angle < kp2.angle;
+        if( kp1.response != kp2.response )
+            return kp1.response > kp2.response;
+        if( kp1.octave != kp2.octave )
+            return kp1.octave > kp2.octave;
+        return kp1.class_id > kp2.class_id;
+    }
+};
+
+void KeyPointsFilter::removeDuplicatedSorted( std::vector<KeyPoint>& keypoints )
+{
+    int i, j, n = (int)keypoints.size();
+
+    if (n < 2) return;
+
+    std::sort(keypoints.begin(), keypoints.end(), KeyPoint12_LessThan());
+
+    for( i = 0, j = 1; j < n; ++j )
+    {
+        const KeyPoint& kp1 = keypoints[i];
+        const KeyPoint& kp2 = keypoints[j];
+        if( kp1.pt.x != kp2.pt.x || kp1.pt.y != kp2.pt.y ||
+            kp1.size != kp2.size || kp1.angle != kp2.angle ) {
+            keypoints[++i] = keypoints[j];
+        }
+    }
+    keypoints.resize(i + 1);
+}
+
+}

modules/features/src/main.cpp

@@ -0,0 +1,52 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Copyright (C) 2015, Itseez Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+//
+// Library initialization file
+//
+
+#include "precomp.hpp"
+
+IPP_INITIALIZER_AUTO
+
+/* End of file. */

modules/features/src/opencl/brute_force_match.cl

@@ -0,0 +1,560 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2010-2012, Multicoreware, Inc., all rights reserved.
+// Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// @Authors
+//    Nathan, liujun@multicorewareinc.com
+//    Peng Xiao, pengxiao@outlook.com
+//    Baichuan Su, baichuan@multicorewareinc.com
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics:enable
+#define MAX_FLOAT 3.40282e+038f
+
+#ifndef T
+#define T float
+#endif
+
+#ifndef BLOCK_SIZE
+#define BLOCK_SIZE 16
+#endif
+#ifndef MAX_DESC_LEN
+#define MAX_DESC_LEN 64
+#endif
+
+#define BLOCK_SIZE_ODD          (BLOCK_SIZE + 1)
+#ifndef SHARED_MEM_SZ
+#  if (BLOCK_SIZE < MAX_DESC_LEN)
+#    define SHARED_MEM_SZ      (kercn * (BLOCK_SIZE * MAX_DESC_LEN + BLOCK_SIZE * BLOCK_SIZE))
+#  else
+#    define SHARED_MEM_SZ      (kercn * 2 * BLOCK_SIZE_ODD * BLOCK_SIZE)
+#  endif
+#endif
+
+#ifndef DIST_TYPE
+#define DIST_TYPE 2
+#endif
+
+// dirty fix for non-template support
+#if (DIST_TYPE == 2) // L1Dist
+#   ifdef T_FLOAT
+        typedef float result_type;
+#       if (8 == kercn)
+            typedef float8 value_type;
+#           define DIST(x, y) {value_type d = fabs((x) - (y)); result += d.s0 + d.s1 + d.s2 + d.s3 + d.s4 + d.s5 + d.s6 + d.s7;}
+#       elif (4 == kercn)
+            typedef float4 value_type;
+#           define DIST(x, y) {value_type d = fabs((x) - (y)); result += d.s0 + d.s1 + d.s2 + d.s3;}
+#       else
+            typedef float value_type;
+#           define DIST(x, y) result += fabs((x) - (y))
+#       endif
+#   else
+        typedef int result_type;
+#       if (8 == kercn)
+            typedef int8 value_type;
+#           define DIST(x, y) {value_type d = abs((x) - (y)); result += d.s0 + d.s1 + d.s2 + d.s3 + d.s4 + d.s5 + d.s6 + d.s7;}
+#       elif (4 == kercn)
+            typedef int4 value_type;
+#           define DIST(x, y) {value_type d = abs((x) - (y)); result += d.s0 + d.s1 + d.s2 + d.s3;}
+#       else
+            typedef int  value_type;
+#           define DIST(x, y) result += abs((x) - (y))
+#       endif
+#   endif
+#   define DIST_RES(x) (x)
+#elif (DIST_TYPE == 4) // L2Dist
+    typedef float result_type;
+#   if (8 == kercn)
+        typedef float8 value_type;
+#       define DIST(x, y)   {value_type d = ((x) - (y)); result += dot(d.s0123, d.s0123) + dot(d.s4567, d.s4567);}
+#   elif (4 == kercn)
+        typedef float4      value_type;
+#       define DIST(x, y)   {value_type d = ((x) - (y)); result += dot(d, d);}
+#   else
+        typedef float       value_type;
+#       define DIST(x, y)   {value_type d = ((x) - (y)); result = mad(d, d, result);}
+#   endif
+#   define DIST_RES(x) sqrt(x)
+#elif (DIST_TYPE == 6) // Hamming
+#   if (8 == kercn)
+        typedef int8 value_type;
+#   elif (4 == kercn)
+        typedef int4 value_type;
+#   else
+        typedef int value_type;
+#   endif
+    typedef int result_type;
+#   define DIST(x, y) result += popcount( (x) ^ (y) )
+#   define DIST_RES(x) (x)
+#endif
+
+inline result_type reduce_block(
+    __local value_type *s_query,
+    __local value_type *s_train,
+    int lidx,
+    int lidy
+    )
+{
+    result_type result = 0;
+    #pragma unroll
+    for (int j = 0 ; j < BLOCK_SIZE ; j++)
+    {
+        DIST(s_query[lidy * BLOCK_SIZE_ODD + j], s_train[j * BLOCK_SIZE_ODD + lidx]);
+    }
+    return DIST_RES(result);
+}
+
+inline result_type reduce_block_match(
+    __local value_type *s_query,
+    __local value_type *s_train,
+    int lidx,
+    int lidy
+    )
+{
+    result_type result = 0;
+    #pragma unroll
+    for (int j = 0 ; j < BLOCK_SIZE ; j++)
+    {
+        DIST(s_query[lidy * BLOCK_SIZE_ODD + j], s_train[j * BLOCK_SIZE_ODD + lidx]);
+    }
+    return result;
+}
+
+inline result_type reduce_multi_block(
+    __local value_type *s_query,
+    __local value_type *s_train,
+    int block_index,
+    int lidx,
+    int lidy
+    )
+{
+    result_type result = 0;
+    #pragma unroll
+    for (int j = 0 ; j < BLOCK_SIZE ; j++)
+    {
+        DIST(s_query[lidy * MAX_DESC_LEN + block_index * BLOCK_SIZE + j], s_train[j * BLOCK_SIZE + lidx]);
+    }
+    return result;
+}
+
+__kernel void BruteForceMatch_Match(
+    __global T *query,
+    __global T *train,
+    __global int *bestTrainIdx,
+    __global float *bestDistance,
+    int query_rows,
+    int query_cols,
+    int train_rows,
+    int train_cols,
+    int step
+)
+{
+    const int lidx = get_local_id(0);
+    const int lidy = get_local_id(1);
+    const int groupidx = get_group_id(0);
+
+    const int queryIdx = mad24(BLOCK_SIZE, groupidx, lidy);
+    const int queryOffset = min(queryIdx, query_rows - 1) * step;
+    __global TN *query_vec = (__global TN *)(query + queryOffset);
+    query_cols /= kercn;
+
+    __local float sharebuffer[SHARED_MEM_SZ];
+    __local value_type *s_query = (__local value_type *)sharebuffer;
+
+#if 0 < MAX_DESC_LEN
+    __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE * MAX_DESC_LEN;
+    // load the query into local memory.
+    #pragma unroll
+    for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; i++)
+    {
+        const int loadx = mad24(BLOCK_SIZE, i, lidx);
+        s_query[mad24(MAX_DESC_LEN, lidy, loadx)] = loadx < query_cols ? query_vec[loadx] : 0;
+    }
+#else
+    __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE;
+    const int s_query_i = mad24(BLOCK_SIZE_ODD, lidy, lidx);
+    const int s_train_i = mad24(BLOCK_SIZE_ODD, lidx, lidy);
+#endif
+
+    float myBestDistance = MAX_FLOAT;
+    int myBestTrainIdx = -1;
+
+    // loopUnrolledCached to find the best trainIdx and best distance.
+    for (int t = 0, endt = (train_rows + BLOCK_SIZE - 1) / BLOCK_SIZE; t < endt; t++)
+    {
+        result_type result = 0;
+
+        const int trainOffset = min(mad24(BLOCK_SIZE, t, lidy), train_rows - 1) * step;
+        __global TN *train_vec = (__global TN *)(train + trainOffset);
+#if 0 < MAX_DESC_LEN
+        #pragma unroll
+        for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; i++)
+        {
+            //load a BLOCK_SIZE * BLOCK_SIZE block into local train.
+            const int loadx = mad24(BLOCK_SIZE, i, lidx);
+            s_train[mad24(BLOCK_SIZE, lidx, lidy)] = loadx < train_cols ? train_vec[loadx] : 0;
+
+            //synchronize to make sure each elem for reduceIteration in share memory is written already.
+            barrier(CLK_LOCAL_MEM_FENCE);
+
+            result += reduce_multi_block(s_query, s_train, i, lidx, lidy);
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+#else
+        for (int i = 0, endq = (query_cols + BLOCK_SIZE - 1) / BLOCK_SIZE; i < endq; i++)
+        {
+            const int loadx = mad24(i, BLOCK_SIZE, lidx);
+            //load query and train into local memory
+            if (loadx < query_cols)
+            {
+                s_query[s_query_i] = query_vec[loadx];
+                s_train[s_train_i] = train_vec[loadx];
+            }
+            else
+            {
+                s_query[s_query_i] = 0;
+                s_train[s_train_i] = 0;
+            }
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+
+            result += reduce_block_match(s_query, s_train, lidx, lidy);
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+#endif
+        result = DIST_RES(result);
+
+        const int trainIdx = mad24(BLOCK_SIZE, t, lidx);
+
+        if (queryIdx < query_rows && trainIdx < train_rows && result < myBestDistance /*&& mask(queryIdx, trainIdx)*/)
+        {
+            myBestDistance = result;
+            myBestTrainIdx = trainIdx;
+        }
+    }
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    __local float *s_distance = (__local float *)sharebuffer;
+    __local int *s_trainIdx = (__local int *)(sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE);
+
+    //findBestMatch
+    s_distance += lidy * BLOCK_SIZE_ODD;
+    s_trainIdx += lidy * BLOCK_SIZE_ODD;
+    s_distance[lidx] = myBestDistance;
+    s_trainIdx[lidx] = myBestTrainIdx;
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    //reduce -- now all reduce implement in each threads.
+    #pragma unroll
+    for (int k = 0 ; k < BLOCK_SIZE; k++)
+    {
+        if (myBestDistance > s_distance[k])
+        {
+            myBestDistance = s_distance[k];
+            myBestTrainIdx = s_trainIdx[k];
+        }
+    }
+
+    if (queryIdx < query_rows && lidx == 0)
+    {
+        bestTrainIdx[queryIdx] = myBestTrainIdx;
+        bestDistance[queryIdx] = myBestDistance;
+    }
+}
+
+//radius_match
+__kernel void BruteForceMatch_RadiusMatch(
+    __global T *query,
+    __global T *train,
+    float maxDistance,
+    __global int *bestTrainIdx,
+    __global float *bestDistance,
+    __global int *nMatches,
+    int query_rows,
+    int query_cols,
+    int train_rows,
+    int train_cols,
+    int bestTrainIdx_cols,
+    int step,
+    int ostep
+)
+{
+    const int lidx = get_local_id(0);
+    const int lidy = get_local_id(1);
+    const int groupidx = get_group_id(0);
+    const int groupidy = get_group_id(1);
+
+    const int queryIdx = mad24(BLOCK_SIZE, groupidy, lidy);
+    const int queryOffset = min(queryIdx, query_rows - 1) * step;
+    __global TN *query_vec = (__global TN *)(query + queryOffset);
+
+    const int trainIdx = mad24(BLOCK_SIZE, groupidx, lidx);
+    const int trainOffset = min(mad24(BLOCK_SIZE, groupidx, lidy), train_rows - 1) * step;
+    __global TN *train_vec = (__global TN *)(train + trainOffset);
+
+    query_cols /= kercn;
+
+    __local float sharebuffer[SHARED_MEM_SZ];
+    __local value_type *s_query = (__local value_type *)sharebuffer;
+    __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE;
+
+    result_type result = 0;
+    const int s_query_i = mad24(BLOCK_SIZE_ODD, lidy, lidx);
+    const int s_train_i = mad24(BLOCK_SIZE_ODD, lidx, lidy);
+    for (int i = 0 ; i < (query_cols + BLOCK_SIZE - 1) / BLOCK_SIZE ; ++i)
+    {
+        //load a BLOCK_SIZE * BLOCK_SIZE block into local train.
+        const int loadx = mad24(BLOCK_SIZE, i, lidx);
+
+        if (loadx < query_cols)
+        {
+            s_query[s_query_i] = query_vec[loadx];
+            s_train[s_train_i] = train_vec[loadx];
+        }
+        else
+        {
+            s_query[s_query_i] = 0;
+            s_train[s_train_i] = 0;
+        }
+
+        //synchronize to make sure each elem for reduceIteration in share memory is written already.
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        result += reduce_block(s_query, s_train, lidx, lidy);
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+    }
+    if (queryIdx < query_rows && trainIdx < train_rows && convert_float(result) < maxDistance)
+    {
+        int ind = atom_inc(nMatches + queryIdx);
+
+        if(ind < bestTrainIdx_cols)
+        {
+            bestTrainIdx[mad24(queryIdx, ostep, ind)] = trainIdx;
+            bestDistance[mad24(queryIdx, ostep, ind)] = result;
+        }
+    }
+}
+
+__kernel void BruteForceMatch_knnMatch(
+    __global T *query,
+    __global T *train,
+    __global int2 *bestTrainIdx,
+    __global float2 *bestDistance,
+    int query_rows,
+    int query_cols,
+    int train_rows,
+    int train_cols,
+    int step
+)
+{
+    const int lidx = get_local_id(0);
+    const int lidy = get_local_id(1);
+    const int groupidx = get_group_id(0);
+
+    const int queryIdx = mad24(BLOCK_SIZE, groupidx, lidy);
+    const int queryOffset = min(queryIdx, query_rows - 1) * step;
+    __global TN *query_vec = (__global TN *)(query + queryOffset);
+    query_cols /= kercn;
+
+    __local float sharebuffer[SHARED_MEM_SZ];
+    __local value_type *s_query = (__local value_type *)sharebuffer;
+
+#if 0 < MAX_DESC_LEN
+    __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE * MAX_DESC_LEN;
+    // load the query into local memory.
+    #pragma unroll
+    for (int i = 0 ;  i <  MAX_DESC_LEN / BLOCK_SIZE; i ++)
+    {
+        int loadx = mad24(BLOCK_SIZE, i, lidx);
+        s_query[mad24(MAX_DESC_LEN, lidy, loadx)] = loadx < query_cols ? query_vec[loadx] : 0;
+    }
+#else
+    __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE;
+    const int s_query_i = mad24(BLOCK_SIZE_ODD, lidy, lidx);
+    const int s_train_i = mad24(BLOCK_SIZE_ODD, lidx, lidy);
+#endif
+
+    float myBestDistance1 = MAX_FLOAT;
+    float myBestDistance2 = MAX_FLOAT;
+    int myBestTrainIdx1 = -1;
+    int myBestTrainIdx2 = -1;
+
+    for (int t = 0, endt = (train_rows + BLOCK_SIZE - 1) / BLOCK_SIZE; t < endt ; t++)
+    {
+        result_type result = 0;
+
+        int trainOffset = min(mad24(BLOCK_SIZE, t, lidy), train_rows - 1) * step;
+        __global TN *train_vec = (__global TN *)(train + trainOffset);
+#if 0 < MAX_DESC_LEN
+        #pragma unroll
+        for (int i = 0 ; i < MAX_DESC_LEN / BLOCK_SIZE ; i++)
+        {
+            //load a BLOCK_SIZE * BLOCK_SIZE block into local train.
+            const int loadx = mad24(BLOCK_SIZE, i, lidx);
+            s_train[mad24(BLOCK_SIZE, lidx, lidy)] = loadx < train_cols ? train_vec[loadx] : 0;
+
+            //synchronize to make sure each elem for reduceIteration in share memory is written already.
+            barrier(CLK_LOCAL_MEM_FENCE);
+
+            result += reduce_multi_block(s_query, s_train, i, lidx, lidy);
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+#else
+        for (int i = 0, endq = (query_cols + BLOCK_SIZE -1) / BLOCK_SIZE; i < endq ; i++)
+        {
+            const int loadx = mad24(BLOCK_SIZE, i, lidx);
+            //load query and train into local memory
+            if (loadx < query_cols)
+            {
+                s_query[s_query_i] = query_vec[loadx];
+                s_train[s_train_i] = train_vec[loadx];
+            }
+            else
+            {
+                s_query[s_query_i] = 0;
+                s_train[s_train_i] = 0;
+            }
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+
+            result += reduce_block_match(s_query, s_train, lidx, lidy);
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+#endif
+        result = DIST_RES(result);
+
+        const int trainIdx = mad24(BLOCK_SIZE, t, lidx);
+
+        if (queryIdx < query_rows && trainIdx < train_rows)
+        {
+            if (result < myBestDistance1)
+            {
+                myBestDistance2 = myBestDistance1;
+                myBestTrainIdx2 = myBestTrainIdx1;
+                myBestDistance1 = result;
+                myBestTrainIdx1 = trainIdx;
+            }
+            else if (result < myBestDistance2)
+            {
+                myBestDistance2 = result;
+                myBestTrainIdx2 = trainIdx;
+            }
+        }
+    }
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    __local float *s_distance = (__local float *)sharebuffer;
+    __local int *s_trainIdx = (__local int *)(sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE);
+
+    // find BestMatch
+    s_distance += lidy * BLOCK_SIZE_ODD;
+    s_trainIdx += lidy * BLOCK_SIZE_ODD;
+    s_distance[lidx] = myBestDistance1;
+    s_trainIdx[lidx] = myBestTrainIdx1;
+
+    float bestDistance1 = MAX_FLOAT;
+    float bestDistance2 = MAX_FLOAT;
+    int bestTrainIdx1 = -1;
+    int bestTrainIdx2 = -1;
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    if (lidx == 0)
+    {
+        for (int i = 0 ; i < BLOCK_SIZE ; i++)
+        {
+            float val = s_distance[i];
+            if (val < bestDistance1)
+            {
+                bestDistance2 = bestDistance1;
+                bestTrainIdx2 = bestTrainIdx1;
+
+                bestDistance1 = val;
+                bestTrainIdx1 = s_trainIdx[i];
+            }
+            else if (val < bestDistance2)
+            {
+                bestDistance2 = val;
+                bestTrainIdx2 = s_trainIdx[i];
+            }
+        }
+    }
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    s_distance[lidx] = myBestDistance2;
+    s_trainIdx[lidx] = myBestTrainIdx2;
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    if (lidx == 0)
+    {
+        for (int i = 0 ; i < BLOCK_SIZE ; i++)
+        {
+            float val = s_distance[i];
+
+            if (val < bestDistance2)
+            {
+                bestDistance2 = val;
+                bestTrainIdx2 = s_trainIdx[i];
+            }
+        }
+    }
+
+    myBestDistance1 = bestDistance1;
+    myBestDistance2 = bestDistance2;
+
+    myBestTrainIdx1 = bestTrainIdx1;
+    myBestTrainIdx2 = bestTrainIdx2;
+
+    if (queryIdx < query_rows && lidx == 0)
+    {
+        bestTrainIdx[queryIdx] = (int2)(myBestTrainIdx1, myBestTrainIdx2);
+        bestDistance[queryIdx] = (float2)(myBestDistance1, myBestDistance2);
+    }
+}

modules/features/src/opencl/fast.cl

@@ -0,0 +1,162 @@
+// OpenCL port of the FAST corner detector.
+// Copyright (C) 2014, Itseez Inc. See the license at http://opencv.org
+
+inline int cornerScore(__global const uchar* img, int step)
+{
+    int k, tofs, v = img[0], a0 = 0, b0;
+    int d[16];
+    #define LOAD2(idx, ofs) \
+        tofs = ofs; d[idx] = (short)(v - img[tofs]); d[idx+8] = (short)(v - img[-tofs])
+    LOAD2(0, 3);
+    LOAD2(1, -step+3);
+    LOAD2(2, -step*2+2);
+    LOAD2(3, -step*3+1);
+    LOAD2(4, -step*3);
+    LOAD2(5, -step*3-1);
+    LOAD2(6, -step*2-2);
+    LOAD2(7, -step-3);
+
+    #pragma unroll
+    for( k = 0; k < 16; k += 2 )
+    {
+        int a = min((int)d[(k+1)&15], (int)d[(k+2)&15]);
+        a = min(a, (int)d[(k+3)&15]);
+        a = min(a, (int)d[(k+4)&15]);
+        a = min(a, (int)d[(k+5)&15]);
+        a = min(a, (int)d[(k+6)&15]);
+        a = min(a, (int)d[(k+7)&15]);
+        a = min(a, (int)d[(k+8)&15]);
+        a0 = max(a0, min(a, (int)d[k&15]));
+        a0 = max(a0, min(a, (int)d[(k+9)&15]));
+    }
+
+    b0 = -a0;
+    #pragma unroll
+    for( k = 0; k < 16; k += 2 )
+    {
+        int b = max((int)d[(k+1)&15], (int)d[(k+2)&15]);
+        b = max(b, (int)d[(k+3)&15]);
+        b = max(b, (int)d[(k+4)&15]);
+        b = max(b, (int)d[(k+5)&15]);
+        b = max(b, (int)d[(k+6)&15]);
+        b = max(b, (int)d[(k+7)&15]);
+        b = max(b, (int)d[(k+8)&15]);
+
+        b0 = min(b0, max(b, (int)d[k]));
+        b0 = min(b0, max(b, (int)d[(k+9)&15]));
+    }
+
+    return -b0-1;
+}
+
+__kernel
+void FAST_findKeypoints(
+    __global const uchar * _img, int step, int img_offset,
+    int img_rows, int img_cols,
+    volatile __global int* kp_loc,
+    int max_keypoints, int threshold )
+{
+    int j = get_global_id(0) + 3;
+    int i = get_global_id(1) + 3;
+
+    if (i < img_rows - 3 && j < img_cols - 3)
+    {
+        __global const uchar* img = _img + mad24(i, step, j + img_offset);
+        int v = img[0], t0 = v - threshold, t1 = v + threshold;
+        int k, tofs, v0, v1;
+        int m0 = 0, m1 = 0;
+
+        #define UPDATE_MASK(idx, ofs) \
+            tofs = ofs; v0 = img[tofs]; v1 = img[-tofs]; \
+            m0 |= ((v0 < t0) << idx) | ((v1 < t0) << (8 + idx)); \
+            m1 |= ((v0 > t1) << idx) | ((v1 > t1) << (8 + idx))
+
+        UPDATE_MASK(0, 3);
+        if( (m0 | m1) == 0 )
+            return;
+
+        UPDATE_MASK(2, -step*2+2);
+        UPDATE_MASK(4, -step*3);
+        UPDATE_MASK(6, -step*2-2);
+
+        #define EVEN_MASK (1+4+16+64)
+
+        if( ((m0 | (m0 >> 8)) & EVEN_MASK) != EVEN_MASK &&
+            ((m1 | (m1 >> 8)) & EVEN_MASK) != EVEN_MASK )
+            return;
+
+        UPDATE_MASK(1, -step+3);
+        UPDATE_MASK(3, -step*3+1);
+        UPDATE_MASK(5, -step*3-1);
+        UPDATE_MASK(7, -step-3);
+        if( ((m0 | (m0 >> 8)) & 255) != 255 &&
+            ((m1 | (m1 >> 8)) & 255) != 255 )
+            return;
+
+        m0 |= m0 << 16;
+        m1 |= m1 << 16;
+
+        #define CHECK0(i) ((m0 & (511 << i)) == (511 << i))
+        #define CHECK1(i) ((m1 & (511 << i)) == (511 << i))
+
+        if( CHECK0(0) + CHECK0(1) + CHECK0(2) + CHECK0(3) +
+            CHECK0(4) + CHECK0(5) + CHECK0(6) + CHECK0(7) +
+            CHECK0(8) + CHECK0(9) + CHECK0(10) + CHECK0(11) +
+            CHECK0(12) + CHECK0(13) + CHECK0(14) + CHECK0(15) +
+
+            CHECK1(0) + CHECK1(1) + CHECK1(2) + CHECK1(3) +
+            CHECK1(4) + CHECK1(5) + CHECK1(6) + CHECK1(7) +
+            CHECK1(8) + CHECK1(9) + CHECK1(10) + CHECK1(11) +
+            CHECK1(12) + CHECK1(13) + CHECK1(14) + CHECK1(15) == 0 )
+            return;
+
+        {
+            int idx = atomic_inc(kp_loc);
+            if( idx < max_keypoints )
+            {
+                kp_loc[1 + 2*idx] = j;
+                kp_loc[2 + 2*idx] = i;
+            }
+        }
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////
+// nonmaxSupression
+
+__kernel
+void FAST_nonmaxSupression(
+    __global const int* kp_in, volatile __global int* kp_out,
+    __global const uchar * _img, int step, int img_offset,
+    int rows, int cols, int counter, int max_keypoints)
+{
+    const int idx = get_global_id(0);
+
+    if (idx < counter)
+    {
+        int x = kp_in[1 + 2*idx];
+        int y = kp_in[2 + 2*idx];
+        __global const uchar* img = _img + mad24(y, step, x + img_offset);
+
+        int s = cornerScore(img, step);
+
+        if( (x < 4 || s > cornerScore(img-1, step)) +
+            (y < 4 || s > cornerScore(img-step, step)) != 2 )
+            return;
+        if( (x >= cols - 4 || s > cornerScore(img+1, step)) +
+            (y >= rows - 4 || s > cornerScore(img+step, step)) +
+            (x < 4 || y < 4 || s > cornerScore(img-step-1, step)) +
+            (x >= cols - 4 || y < 4 || s > cornerScore(img-step+1, step)) +
+            (x < 4 || y >= rows - 4 || s > cornerScore(img+step-1, step)) +
+            (x >= cols - 4 || y >= rows - 4 || s > cornerScore(img+step+1, step)) == 6)
+        {
+            int new_idx = atomic_inc(kp_out);
+            if( new_idx < max_keypoints )
+            {
+                kp_out[1 + 3*new_idx] = x;
+                kp_out[2 + 3*new_idx] = y;
+                kp_out[3 + 3*new_idx] = s;
+            }
+        }
+    }
+}

modules/features/src/opencl/orb.cl

@@ -0,0 +1,254 @@
+// OpenCL port of the ORB feature detector and descriptor extractor
+// Copyright (C) 2014, Itseez Inc. See the license at http://opencv.org
+//
+// The original code has been contributed by Peter Andreas Entschev, peter@entschev.com
+
+#define LAYERINFO_SIZE 1
+#define LAYERINFO_OFS 0
+#define KEYPOINT_SIZE 3
+#define ORIENTED_KEYPOINT_SIZE 4
+#define KEYPOINT_X 0
+#define KEYPOINT_Y 1
+#define KEYPOINT_Z 2
+#define KEYPOINT_ANGLE 3
+
+/////////////////////////////////////////////////////////////
+
+#ifdef ORB_RESPONSES
+
+__kernel void
+ORB_HarrisResponses(__global const uchar* imgbuf, int imgstep, int imgoffset0,
+                    __global const int* layerinfo, __global const int* keypoints,
+                    __global float* responses, int nkeypoints )
+{
+    int idx = get_global_id(0);
+    if( idx < nkeypoints )
+    {
+        __global const int* kpt = keypoints + idx*KEYPOINT_SIZE;
+        __global const int* layer = layerinfo + kpt[KEYPOINT_Z]*LAYERINFO_SIZE;
+        __global const uchar* img = imgbuf + imgoffset0 + layer[LAYERINFO_OFS] +
+            (kpt[KEYPOINT_Y] - blockSize/2)*imgstep + (kpt[KEYPOINT_X] - blockSize/2);
+
+        int i, j;
+        int a = 0, b = 0, c = 0;
+        for( i = 0; i < blockSize; i++, img += imgstep-blockSize )
+        {
+            for( j = 0; j < blockSize; j++, img++ )
+            {
+                int Ix = (img[1] - img[-1])*2 + img[-imgstep+1] - img[-imgstep-1] + img[imgstep+1] - img[imgstep-1];
+                int Iy = (img[imgstep] - img[-imgstep])*2 + img[imgstep-1] - img[-imgstep-1] + img[imgstep+1] - img[-imgstep+1];
+                a += Ix*Ix;
+                b += Iy*Iy;
+                c += Ix*Iy;
+            }
+        }
+        responses[idx] = ((float)a * b - (float)c * c - HARRIS_K * (float)(a + b) * (a + b))*scale_sq_sq;
+    }
+}
+
+#endif
+
+/////////////////////////////////////////////////////////////
+
+#ifdef ORB_ANGLES
+
+#define _DBL_EPSILON 2.2204460492503131e-16f
+#define atan2_p1 (0.9997878412794807f*57.29577951308232f)
+#define atan2_p3 (-0.3258083974640975f*57.29577951308232f)
+#define atan2_p5 (0.1555786518463281f*57.29577951308232f)
+#define atan2_p7 (-0.04432655554792128f*57.29577951308232f)
+
+inline float fastAtan2( float y, float x )
+{
+    float ax = fabs(x), ay = fabs(y);
+    float a, c, c2;
+    if( ax >= ay )
+    {
+        c = ay/(ax + _DBL_EPSILON);
+        c2 = c*c;
+        a = (((atan2_p7*c2 + atan2_p5)*c2 + atan2_p3)*c2 + atan2_p1)*c;
+    }
+    else
+    {
+        c = ax/(ay + _DBL_EPSILON);
+        c2 = c*c;
+        a = 90.f - (((atan2_p7*c2 + atan2_p5)*c2 + atan2_p3)*c2 + atan2_p1)*c;
+    }
+    if( x < 0 )
+        a = 180.f - a;
+    if( y < 0 )
+        a = 360.f - a;
+    return a;
+}
+
+
+__kernel void
+ORB_ICAngle(__global const uchar* imgbuf, int imgstep, int imgoffset0,
+            __global const int* layerinfo, __global const int* keypoints,
+            __global float* responses, const __global int* u_max,
+            int nkeypoints, int half_k )
+{
+    int idx = get_global_id(0);
+    if( idx < nkeypoints )
+    {
+        __global const int* kpt = keypoints + idx*KEYPOINT_SIZE;
+
+        __global const int* layer = layerinfo + kpt[KEYPOINT_Z]*LAYERINFO_SIZE;
+        __global const uchar* center = imgbuf + imgoffset0 + layer[LAYERINFO_OFS] +
+            kpt[KEYPOINT_Y]*imgstep + kpt[KEYPOINT_X];
+
+        int u, v, m_01 = 0, m_10 = 0;
+
+        // Treat the center line differently, v=0
+        for( u = -half_k; u <= half_k; u++ )
+            m_10 += u * center[u];
+
+        // Go line by line in the circular patch
+        for( v = 1; v <= half_k; v++ )
+        {
+            // Proceed over the two lines
+            int v_sum = 0;
+            int d = u_max[v];
+            for( u = -d; u <= d; u++ )
+            {
+                int val_plus = center[u + v*imgstep], val_minus = center[u - v*imgstep];
+                v_sum += (val_plus - val_minus);
+                m_10 += u * (val_plus + val_minus);
+            }
+            m_01 += v * v_sum;
+        }
+
+        // we do not use OpenCL's atan2 intrinsic,
+        // because we want to get _exactly_ the same results as the CPU version
+        responses[idx] = fastAtan2((float)m_01, (float)m_10);
+    }
+}
+
+#endif
+
+/////////////////////////////////////////////////////////////
+
+#ifdef ORB_DESCRIPTORS
+
+__kernel void
+ORB_computeDescriptor(__global const uchar* imgbuf, int imgstep, int imgoffset0,
+                      __global const int* layerinfo, __global const int* keypoints,
+                      __global uchar* _desc, const __global int* pattern,
+                      int nkeypoints, int dsize )
+{
+    int idx = get_global_id(0);
+    if( idx < nkeypoints )
+    {
+        int i;
+        __global const int* kpt = keypoints + idx*ORIENTED_KEYPOINT_SIZE;
+
+        __global const int* layer = layerinfo + kpt[KEYPOINT_Z]*LAYERINFO_SIZE;
+        __global const uchar* center = imgbuf + imgoffset0 + layer[LAYERINFO_OFS] +
+                                kpt[KEYPOINT_Y]*imgstep + kpt[KEYPOINT_X];
+        float angle = as_float(kpt[KEYPOINT_ANGLE]);
+        angle *= 0.01745329251994329547f;
+
+        float cosa;
+        float sina = sincos(angle, &cosa);
+
+        __global uchar* desc = _desc + idx*dsize;
+
+        #define GET_VALUE(idx) \
+            center[mad24(convert_int_rte(pattern[(idx)*2] * sina + pattern[(idx)*2+1] * cosa), imgstep, \
+                        convert_int_rte(pattern[(idx)*2] * cosa - pattern[(idx)*2+1] * sina))]
+
+        for( i = 0; i < dsize; i++ )
+        {
+            int val;
+        #if WTA_K == 2
+            int t0, t1;
+
+            t0 = GET_VALUE(0); t1 = GET_VALUE(1);
+            val = t0 < t1;
+
+            t0 = GET_VALUE(2); t1 = GET_VALUE(3);
+            val |= (t0 < t1) << 1;
+
+            t0 = GET_VALUE(4); t1 = GET_VALUE(5);
+            val |= (t0 < t1) << 2;
+
+            t0 = GET_VALUE(6); t1 = GET_VALUE(7);
+            val |= (t0 < t1) << 3;
+
+            t0 = GET_VALUE(8); t1 = GET_VALUE(9);
+            val |= (t0 < t1) << 4;
+
+            t0 = GET_VALUE(10); t1 = GET_VALUE(11);
+            val |= (t0 < t1) << 5;
+
+            t0 = GET_VALUE(12); t1 = GET_VALUE(13);
+            val |= (t0 < t1) << 6;
+
+            t0 = GET_VALUE(14); t1 = GET_VALUE(15);
+            val |= (t0 < t1) << 7;
+
+            pattern += 16*2;
+
+        #elif WTA_K == 3
+            int t0, t1, t2;
+
+            t0 = GET_VALUE(0); t1 = GET_VALUE(1); t2 = GET_VALUE(2);
+            val = t2 > t1 ? (t2 > t0 ? 2 : 0) : (t1 > t0);
+
+            t0 = GET_VALUE(3); t1 = GET_VALUE(4); t2 = GET_VALUE(5);
+            val |= (t2 > t1 ? (t2 > t0 ? 2 : 0) : (t1 > t0)) << 2;
+
+            t0 = GET_VALUE(6); t1 = GET_VALUE(7); t2 = GET_VALUE(8);
+            val |= (t2 > t1 ? (t2 > t0 ? 2 : 0) : (t1 > t0)) << 4;
+
+            t0 = GET_VALUE(9); t1 = GET_VALUE(10); t2 = GET_VALUE(11);
+            val |= (t2 > t1 ? (t2 > t0 ? 2 : 0) : (t1 > t0)) << 6;
+
+            pattern += 12*2;
+
+        #elif WTA_K == 4
+            int t0, t1, t2, t3, k;
+            int a, b;
+
+            t0 = GET_VALUE(0); t1 = GET_VALUE(1);
+            t2 = GET_VALUE(2); t3 = GET_VALUE(3);
+            a = 0, b = 2;
+            if( t1 > t0 ) t0 = t1, a = 1;
+            if( t3 > t2 ) t2 = t3, b = 3;
+            k = t0 > t2 ? a : b;
+            val = k;
+
+            t0 = GET_VALUE(4); t1 = GET_VALUE(5);
+            t2 = GET_VALUE(6); t3 = GET_VALUE(7);
+            a = 0, b = 2;
+            if( t1 > t0 ) t0 = t1, a = 1;
+            if( t3 > t2 ) t2 = t3, b = 3;
+            k = t0 > t2 ? a : b;
+            val |= k << 2;
+
+            t0 = GET_VALUE(8); t1 = GET_VALUE(9);
+            t2 = GET_VALUE(10); t3 = GET_VALUE(11);
+            a = 0, b = 2;
+            if( t1 > t0 ) t0 = t1, a = 1;
+            if( t3 > t2 ) t2 = t3, b = 3;
+            k = t0 > t2 ? a : b;
+            val |= k << 4;
+
+            t0 = GET_VALUE(12); t1 = GET_VALUE(13);
+            t2 = GET_VALUE(14); t3 = GET_VALUE(15);
+            a = 0, b = 2;
+            if( t1 > t0 ) t0 = t1, a = 1;
+            if( t3 > t2 ) t2 = t3, b = 3;
+            k = t0 > t2 ? a : b;
+            val |= k << 6;
+
+            pattern += 16*2;
+        #else
+            #error "unknown/undefined WTA_K value; should be 2, 3 or 4"
+        #endif
+            desc[i] = (uchar)val;
+        }
+    }
+}
+
+#endif

modules/features/src/precomp.hpp

@@ -0,0 +1,56 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                          License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_PRECOMP_H__
+#define __OPENCV_PRECOMP_H__
+
+#include "opencv2/features.hpp"
+#include "opencv2/imgproc.hpp"
+
+#include "opencv2/core/utility.hpp"
+#include "opencv2/core/private.hpp"
+#include "opencv2/core/ocl.hpp"
+#include "opencv2/core/hal/hal.hpp"
+
+#include <algorithm>
+
+#endif

modules/features/src/sift.dispatch.cpp

@@ -0,0 +1,622 @@
+// 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.
+//
+// Copyright (c) 2006-2010, Rob Hess <hess@eecs.oregonstate.edu>
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Copyright (C) 2020, Intel Corporation, all rights reserved.
+
+/**********************************************************************************************\
+ Implementation of SIFT is based on the code from http://blogs.oregonstate.edu/hess/code/sift/
+ Below is the original copyright.
+ Patent US6711293 expired in March 2020.
+
+//    Copyright (c) 2006-2010, Rob Hess <hess@eecs.oregonstate.edu>
+//    All rights reserved.
+
+//    The following patent has been issued for methods embodied in this
+//    software: "Method and apparatus for identifying scale invariant features
+//    in an image and use of same for locating an object in an image," David
+//    G. Lowe, US Patent 6,711,293 (March 23, 2004). Provisional application
+//    filed March 8, 1999. Asignee: The University of British Columbia. For
+//    further details, contact David Lowe (lowe@cs.ubc.ca) or the
+//    University-Industry Liaison Office of the University of British
+//    Columbia.
+
+//    Note that restrictions imposed by this patent (and possibly others)
+//    exist independently of and may be in conflict with the freedoms granted
+//    in this license, which refers to copyright of the program, not patents
+//    for any methods that it implements.  Both copyright and patent law must
+//    be obeyed to legally use and redistribute this program and it is not the
+//    purpose of this license to induce you to infringe any patents or other
+//    property right claims or to contest validity of any such claims.  If you
+//    redistribute or use the program, then this license merely protects you
+//    from committing copyright infringement.  It does not protect you from
+//    committing patent infringement.  So, before you do anything with this
+//    program, make sure that you have permission to do so not merely in terms
+//    of copyright, but also in terms of patent law.
+
+//    Please note that this license is not to be understood as a guarantee
+//    either.  If you use the program according to this license, but in
+//    conflict with patent law, it does not mean that the licensor will refund
+//    you for any losses that you incur if you are sued for your patent
+//    infringement.
+
+//    Redistribution and use in source and binary forms, with or without
+//    modification, are permitted provided that the following conditions are
+//    met:
+//        * Redistributions of source code must retain the above copyright and
+//          patent notices, this list of conditions and the following
+//          disclaimer.
+//        * Redistributions in binary form must reproduce the above copyright
+//          notice, this list of conditions and the following disclaimer in
+//          the documentation and/or other materials provided with the
+//          distribution.
+//        * Neither the name of Oregon State University nor the names of its
+//          contributors may be used to endorse or promote products derived
+//          from this software without specific prior written permission.
+
+//    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+//    IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+//    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+//    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+//    HOLDER BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+//    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+//    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+//    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+//    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+//    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+//    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+\**********************************************************************************************/
+
+#include "precomp.hpp"
+#include <opencv2/core/hal/hal.hpp>
+#include <opencv2/core/utils/tls.hpp>
+#include <opencv2/core/utils/logger.hpp>
+
+#include "sift.simd.hpp"
+#include "sift.simd_declarations.hpp" // defines CV_CPU_DISPATCH_MODES_ALL=AVX2,...,BASELINE based on CMakeLists.txt content
+
+namespace cv {
+
+/*!
+ SIFT implementation.
+
+ The class implements SIFT algorithm by D. Lowe.
+ */
+class SIFT_Impl : public SIFT
+{
+public:
+    explicit SIFT_Impl( int nfeatures = 0, int nOctaveLayers = 3,
+                          double contrastThreshold = 0.04, double edgeThreshold = 10,
+                          double sigma = 1.6, int descriptorType = CV_32F,
+                          bool enable_precise_upscale = true );
+
+    //! returns the descriptor size in floats (128)
+    int descriptorSize() const CV_OVERRIDE;
+
+    //! returns the descriptor type
+    int descriptorType() const CV_OVERRIDE;
+
+    //! returns the default norm type
+    int defaultNorm() const CV_OVERRIDE;
+
+    //! finds the keypoints and computes descriptors for them using SIFT algorithm.
+    //! Optionally it can compute descriptors for the user-provided keypoints
+    void detectAndCompute(InputArray img, InputArray mask,
+                    std::vector<KeyPoint>& keypoints,
+                    OutputArray descriptors,
+                    bool useProvidedKeypoints = false) CV_OVERRIDE;
+
+    void buildGaussianPyramid( const Mat& base, std::vector<Mat>& pyr, int nOctaves ) const;
+    void buildDoGPyramid( const std::vector<Mat>& pyr, std::vector<Mat>& dogpyr ) const;
+    void findScaleSpaceExtrema( const std::vector<Mat>& gauss_pyr, const std::vector<Mat>& dog_pyr,
+                               std::vector<KeyPoint>& keypoints ) const;
+
+    void read( const FileNode& fn) CV_OVERRIDE;
+    void write( FileStorage& fs) const CV_OVERRIDE;
+
+    void setNFeatures(int maxFeatures) CV_OVERRIDE { nfeatures = maxFeatures; }
+    int getNFeatures() const CV_OVERRIDE { return nfeatures; }
+
+    void setNOctaveLayers(int nOctaveLayers_) CV_OVERRIDE { nOctaveLayers = nOctaveLayers_; }
+    int getNOctaveLayers() const CV_OVERRIDE { return nOctaveLayers; }
+
+    void setContrastThreshold(double contrastThreshold_) CV_OVERRIDE  { contrastThreshold = contrastThreshold_; }
+    double getContrastThreshold() const CV_OVERRIDE { return contrastThreshold; }
+
+    void setEdgeThreshold(double edgeThreshold_) CV_OVERRIDE  { edgeThreshold = edgeThreshold_; }
+    double getEdgeThreshold() const CV_OVERRIDE { return edgeThreshold; }
+
+    void setSigma(double sigma_) CV_OVERRIDE  { sigma = sigma_; }
+    double getSigma() const CV_OVERRIDE { return sigma; }
+
+protected:
+    CV_PROP_RW int nfeatures;
+    CV_PROP_RW int nOctaveLayers;
+    CV_PROP_RW double contrastThreshold;
+    CV_PROP_RW double edgeThreshold;
+    CV_PROP_RW double sigma;
+    CV_PROP_RW int descriptor_type;
+    CV_PROP_RW bool enable_precise_upscale;
+};
+
+Ptr<SIFT> SIFT::create( int _nfeatures, int _nOctaveLayers,
+                     double _contrastThreshold, double _edgeThreshold, double _sigma, bool enable_precise_upscale )
+{
+    CV_TRACE_FUNCTION();
+
+    return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, CV_32F, enable_precise_upscale);
+}
+
+Ptr<SIFT> SIFT::create( int _nfeatures, int _nOctaveLayers,
+                     double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType, bool enable_precise_upscale )
+{
+    CV_TRACE_FUNCTION();
+
+    // SIFT descriptor supports 32bit floating point and 8bit unsigned int.
+    CV_Assert(_descriptorType == CV_32F || _descriptorType == CV_8U);
+    return makePtr<SIFT_Impl>(_nfeatures, _nOctaveLayers, _contrastThreshold, _edgeThreshold, _sigma, _descriptorType, enable_precise_upscale);
+}
+
+String SIFT::getDefaultName() const
+{
+    return (Feature2D::getDefaultName() + ".SIFT");
+}
+
+static inline void
+unpackOctave(const KeyPoint& kpt, int& octave, int& layer, float& scale)
+{
+    octave = kpt.octave & 255;
+    layer = (kpt.octave >> 8) & 255;
+    octave = octave < 128 ? octave : (-128 | octave);
+    scale = octave >= 0 ? 1.f/(1 << octave) : (float)(1 << -octave);
+}
+
+static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma, bool enable_precise_upscale )
+{
+    CV_TRACE_FUNCTION();
+
+    Mat gray, gray_fpt;
+    if( img.channels() == 3 || img.channels() == 4 )
+    {
+        cvtColor(img, gray, COLOR_BGR2GRAY);
+        gray.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);
+    }
+    else
+        img.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);
+
+    float sig_diff;
+
+    if( doubleImageSize )
+    {
+        sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f) );
+
+        Mat dbl;
+        if (enable_precise_upscale) {
+            dbl.create(Size(gray_fpt.cols*2, gray_fpt.rows*2), gray_fpt.type());
+            Mat H = Mat::zeros(2, 3, CV_32F);
+            H.at<float>(0, 0) = 0.5f;
+            H.at<float>(1, 1) = 0.5f;
+
+            cv::warpAffine(gray_fpt, dbl, H, dbl.size(), INTER_LINEAR | WARP_INVERSE_MAP, BORDER_REFLECT);
+        } else {
+#if DoG_TYPE_SHORT
+            resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR_EXACT);
+#else
+            resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR);
+#endif
+        }
+        Mat result;
+        GaussianBlur(dbl, result, Size(), sig_diff, sig_diff);
+        return result;
+    }
+    else
+    {
+        sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA, 0.01f) );
+        Mat result;
+        GaussianBlur(gray_fpt, result, Size(), sig_diff, sig_diff);
+        return result;
+    }
+}
+
+
+void SIFT_Impl::buildGaussianPyramid( const Mat& base, std::vector<Mat>& pyr, int nOctaves ) const
+{
+    CV_TRACE_FUNCTION();
+
+    std::vector<double> sig(nOctaveLayers + 3);
+    pyr.resize(nOctaves*(nOctaveLayers + 3));
+
+    // precompute Gaussian sigmas using the following formula:
+    //  \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2
+    sig[0] = sigma;
+    double k = std::pow( 2., 1. / nOctaveLayers );
+    for( int i = 1; i < nOctaveLayers + 3; i++ )
+    {
+        double sig_prev = std::pow(k, (double)(i-1))*sigma;
+        double sig_total = sig_prev*k;
+        sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev);
+    }
+
+    for( int o = 0; o < nOctaves; o++ )
+    {
+        for( int i = 0; i < nOctaveLayers + 3; i++ )
+        {
+            Mat& dst = pyr[o*(nOctaveLayers + 3) + i];
+            if( o == 0  &&  i == 0 )
+                dst = base;
+            // base of new octave is halved image from end of previous octave
+            else if( i == 0 )
+            {
+                const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers];
+                resize(src, dst, Size(src.cols/2, src.rows/2),
+                       0, 0, INTER_NEAREST);
+            }
+            else
+            {
+                const Mat& src = pyr[o*(nOctaveLayers + 3) + i-1];
+                GaussianBlur(src, dst, Size(), sig[i], sig[i]);
+            }
+        }
+    }
+}
+
+
+class buildDoGPyramidComputer : public ParallelLoopBody
+{
+public:
+    buildDoGPyramidComputer(
+        int _nOctaveLayers,
+        const std::vector<Mat>& _gpyr,
+        std::vector<Mat>& _dogpyr)
+        : nOctaveLayers(_nOctaveLayers),
+          gpyr(_gpyr),
+          dogpyr(_dogpyr) { }
+
+    void operator()( const cv::Range& range ) const CV_OVERRIDE
+    {
+        CV_TRACE_FUNCTION();
+
+        const int begin = range.start;
+        const int end = range.end;
+
+        for( int a = begin; a < end; a++ )
+        {
+            const int o = a / (nOctaveLayers + 2);
+            const int i = a % (nOctaveLayers + 2);
+
+            const Mat& src1 = gpyr[o*(nOctaveLayers + 3) + i];
+            const Mat& src2 = gpyr[o*(nOctaveLayers + 3) + i + 1];
+            Mat& dst = dogpyr[o*(nOctaveLayers + 2) + i];
+            subtract(src2, src1, dst, noArray(), DataType<sift_wt>::type);
+        }
+    }
+
+private:
+    int nOctaveLayers;
+    const std::vector<Mat>& gpyr;
+    std::vector<Mat>& dogpyr;
+};
+
+void SIFT_Impl::buildDoGPyramid( const std::vector<Mat>& gpyr, std::vector<Mat>& dogpyr ) const
+{
+    CV_TRACE_FUNCTION();
+
+    int nOctaves = (int)gpyr.size()/(nOctaveLayers + 3);
+    dogpyr.resize( nOctaves*(nOctaveLayers + 2) );
+
+    parallel_for_(Range(0, nOctaves * (nOctaveLayers + 2)), buildDoGPyramidComputer(nOctaveLayers, gpyr, dogpyr));
+}
+
+class findScaleSpaceExtremaComputer : public ParallelLoopBody
+{
+public:
+    findScaleSpaceExtremaComputer(
+        int _o,
+        int _i,
+        int _threshold,
+        int _idx,
+        int _step,
+        int _cols,
+        int _nOctaveLayers,
+        double _contrastThreshold,
+        double _edgeThreshold,
+        double _sigma,
+        const std::vector<Mat>& _gauss_pyr,
+        const std::vector<Mat>& _dog_pyr,
+        TLSData<std::vector<KeyPoint> > &_tls_kpts_struct)
+
+        : o(_o),
+          i(_i),
+          threshold(_threshold),
+          idx(_idx),
+          step(_step),
+          cols(_cols),
+          nOctaveLayers(_nOctaveLayers),
+          contrastThreshold(_contrastThreshold),
+          edgeThreshold(_edgeThreshold),
+          sigma(_sigma),
+          gauss_pyr(_gauss_pyr),
+          dog_pyr(_dog_pyr),
+          tls_kpts_struct(_tls_kpts_struct) { }
+    void operator()( const cv::Range& range ) const CV_OVERRIDE
+    {
+        CV_TRACE_FUNCTION();
+
+        std::vector<KeyPoint>& kpts = tls_kpts_struct.getRef();
+
+        CV_CPU_DISPATCH(findScaleSpaceExtrema, (o, i, threshold, idx, step, cols, nOctaveLayers, contrastThreshold, edgeThreshold, sigma, gauss_pyr, dog_pyr, kpts, range),
+            CV_CPU_DISPATCH_MODES_ALL);
+    }
+private:
+    int o, i;
+    int threshold;
+    int idx, step, cols;
+    int nOctaveLayers;
+    double contrastThreshold;
+    double edgeThreshold;
+    double sigma;
+    const std::vector<Mat>& gauss_pyr;
+    const std::vector<Mat>& dog_pyr;
+    TLSData<std::vector<KeyPoint> > &tls_kpts_struct;
+};
+
+//
+// Detects features at extrema in DoG scale space.  Bad features are discarded
+// based on contrast and ratio of principal curvatures.
+void SIFT_Impl::findScaleSpaceExtrema( const std::vector<Mat>& gauss_pyr, const std::vector<Mat>& dog_pyr,
+                                  std::vector<KeyPoint>& keypoints ) const
+{
+    CV_TRACE_FUNCTION();
+
+    const int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3);
+    const int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);
+
+    keypoints.clear();
+    TLSDataAccumulator<std::vector<KeyPoint> > tls_kpts_struct;
+
+    for( int o = 0; o < nOctaves; o++ )
+        for( int i = 1; i <= nOctaveLayers; i++ )
+        {
+            const int idx = o*(nOctaveLayers+2)+i;
+            const Mat& img = dog_pyr[idx];
+            const int step = (int)img.step1();
+            const int rows = img.rows, cols = img.cols;
+
+            parallel_for_(Range(SIFT_IMG_BORDER, rows-SIFT_IMG_BORDER),
+                findScaleSpaceExtremaComputer(
+                    o, i, threshold, idx, step, cols,
+                    nOctaveLayers,
+                    contrastThreshold,
+                    edgeThreshold,
+                    sigma,
+                    gauss_pyr, dog_pyr, tls_kpts_struct));
+        }
+
+    std::vector<std::vector<KeyPoint>*> kpt_vecs;
+    tls_kpts_struct.gather(kpt_vecs);
+    for (size_t i = 0; i < kpt_vecs.size(); ++i) {
+        keypoints.insert(keypoints.end(), kpt_vecs[i]->begin(), kpt_vecs[i]->end());
+    }
+}
+
+
+static
+void calcSIFTDescriptor(
+        const Mat& img, Point2f ptf, float ori, float scl,
+        int d, int n, Mat& dst, int row
+)
+{
+    CV_TRACE_FUNCTION();
+
+    CV_CPU_DISPATCH(calcSIFTDescriptor, (img, ptf, ori, scl, d, n, dst, row),
+        CV_CPU_DISPATCH_MODES_ALL);
+}
+
+class calcDescriptorsComputer : public ParallelLoopBody
+{
+public:
+    calcDescriptorsComputer(const std::vector<Mat>& _gpyr,
+                            const std::vector<KeyPoint>& _keypoints,
+                            Mat& _descriptors,
+                            int _nOctaveLayers,
+                            int _firstOctave)
+        : gpyr(_gpyr),
+          keypoints(_keypoints),
+          descriptors(_descriptors),
+          nOctaveLayers(_nOctaveLayers),
+          firstOctave(_firstOctave) { }
+
+    void operator()( const cv::Range& range ) const CV_OVERRIDE
+    {
+        CV_TRACE_FUNCTION();
+
+        const int begin = range.start;
+        const int end = range.end;
+
+        static const int d = SIFT_DESCR_WIDTH, n = SIFT_DESCR_HIST_BINS;
+
+        for ( int i = begin; i<end; i++ )
+        {
+            KeyPoint kpt = keypoints[i];
+            int octave, layer;
+            float scale;
+            unpackOctave(kpt, octave, layer, scale);
+            CV_Assert(octave >= firstOctave && layer <= nOctaveLayers+2);
+            float size=kpt.size*scale;
+            Point2f ptf(kpt.pt.x*scale, kpt.pt.y*scale);
+            const Mat& img = gpyr[(octave - firstOctave)*(nOctaveLayers + 3) + layer];
+
+            float angle = 360.f - kpt.angle;
+            if(std::abs(angle - 360.f) < FLT_EPSILON)
+                angle = 0.f;
+            calcSIFTDescriptor(img, ptf, angle, size*0.5f, d, n, descriptors, i);
+        }
+    }
+private:
+    const std::vector<Mat>& gpyr;
+    const std::vector<KeyPoint>& keypoints;
+    Mat& descriptors;
+    int nOctaveLayers;
+    int firstOctave;
+};
+
+static void calcDescriptors(const std::vector<Mat>& gpyr, const std::vector<KeyPoint>& keypoints,
+                            Mat& descriptors, int nOctaveLayers, int firstOctave )
+{
+    CV_TRACE_FUNCTION();
+    parallel_for_(Range(0, static_cast<int>(keypoints.size())), calcDescriptorsComputer(gpyr, keypoints, descriptors, nOctaveLayers, firstOctave));
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+SIFT_Impl::SIFT_Impl( int _nfeatures, int _nOctaveLayers,
+           double _contrastThreshold, double _edgeThreshold, double _sigma, int _descriptorType, bool _enable_precise_upscale)
+    : nfeatures(_nfeatures), nOctaveLayers(_nOctaveLayers),
+    contrastThreshold(_contrastThreshold), edgeThreshold(_edgeThreshold), sigma(_sigma), descriptor_type(_descriptorType),
+    enable_precise_upscale(_enable_precise_upscale)
+{
+    if (!enable_precise_upscale) {
+        CV_LOG_ONCE_INFO(NULL, "precise upscale disabled, this is now deprecated as it was found to induce a location bias");
+    }
+}
+
+int SIFT_Impl::descriptorSize() const
+{
+    return SIFT_DESCR_WIDTH*SIFT_DESCR_WIDTH*SIFT_DESCR_HIST_BINS;
+}
+
+int SIFT_Impl::descriptorType() const
+{
+    return descriptor_type;
+}
+
+int SIFT_Impl::defaultNorm() const
+{
+    return NORM_L2;
+}
+
+
+void SIFT_Impl::detectAndCompute(InputArray _image, InputArray _mask,
+                      std::vector<KeyPoint>& keypoints,
+                      OutputArray _descriptors,
+                      bool useProvidedKeypoints)
+{
+    CV_TRACE_FUNCTION();
+
+    int firstOctave = -1, actualNOctaves = 0, actualNLayers = 0;
+    Mat image = _image.getMat(), mask = _mask.getMat();
+
+    if( image.empty() || image.depth() != CV_8U )
+        CV_Error( Error::StsBadArg, "image is empty or has incorrect depth (!=CV_8U)" );
+
+    if( !mask.empty() && mask.type() != CV_8UC1 )
+        CV_Error( Error::StsBadArg, "mask has incorrect type (!=CV_8UC1)" );
+
+    if( useProvidedKeypoints )
+    {
+        firstOctave = 0;
+        int maxOctave = INT_MIN;
+        for( size_t i = 0; i < keypoints.size(); i++ )
+        {
+            int octave, layer;
+            float scale;
+            unpackOctave(keypoints[i], octave, layer, scale);
+            firstOctave = std::min(firstOctave, octave);
+            maxOctave = std::max(maxOctave, octave);
+            actualNLayers = std::max(actualNLayers, layer-2);
+        }
+
+        firstOctave = std::min(firstOctave, 0);
+        CV_Assert( firstOctave >= -1 && actualNLayers <= nOctaveLayers );
+        actualNOctaves = maxOctave - firstOctave + 1;
+    }
+
+    Mat base = createInitialImage(image, firstOctave < 0, (float)sigma, enable_precise_upscale);
+    std::vector<Mat> gpyr;
+    int nOctaves = actualNOctaves > 0 ? actualNOctaves : cvRound(std::log( (double)std::min( base.cols, base.rows ) ) / std::log(2.) - 2) - firstOctave;
+
+    //double t, tf = getTickFrequency();
+    //t = (double)getTickCount();
+    buildGaussianPyramid(base, gpyr, nOctaves);
+
+    //t = (double)getTickCount() - t;
+    //printf("pyramid construction time: %g\n", t*1000./tf);
+
+    if( !useProvidedKeypoints )
+    {
+        std::vector<Mat> dogpyr;
+        buildDoGPyramid(gpyr, dogpyr);
+        //t = (double)getTickCount();
+        findScaleSpaceExtrema(gpyr, dogpyr, keypoints);
+        KeyPointsFilter::removeDuplicatedSorted( keypoints );
+
+        if( nfeatures > 0 )
+            KeyPointsFilter::retainBest(keypoints, nfeatures);
+        //t = (double)getTickCount() - t;
+        //printf("keypoint detection time: %g\n", t*1000./tf);
+
+        if( firstOctave < 0 )
+            for( size_t i = 0; i < keypoints.size(); i++ )
+            {
+                KeyPoint& kpt = keypoints[i];
+                float scale = 1.f/(float)(1 << -firstOctave);
+                kpt.octave = (kpt.octave & ~255) | ((kpt.octave + firstOctave) & 255);
+                kpt.pt *= scale;
+                kpt.size *= scale;
+            }
+
+        if( !mask.empty() )
+            KeyPointsFilter::runByPixelsMask( keypoints, mask );
+    }
+    else
+    {
+        // filter keypoints by mask
+        //KeyPointsFilter::runByPixelsMask( keypoints, mask );
+    }
+
+    if( _descriptors.needed() )
+    {
+        //t = (double)getTickCount();
+        int dsize = descriptorSize();
+        _descriptors.create((int)keypoints.size(), dsize, descriptor_type);
+
+        Mat descriptors = _descriptors.getMat();
+        calcDescriptors(gpyr, keypoints, descriptors, nOctaveLayers, firstOctave);
+        //t = (double)getTickCount() - t;
+        //printf("descriptor extraction time: %g\n", t*1000./tf);
+    }
+}
+
+void SIFT_Impl::read( const FileNode& fn)
+{
+  // if node is empty, keep previous value
+  if (!fn["nfeatures"].empty())
+    fn["nfeatures"] >> nfeatures;
+  if (!fn["nOctaveLayers"].empty())
+    fn["nOctaveLayers"] >> nOctaveLayers;
+  if (!fn["contrastThreshold"].empty())
+    fn["contrastThreshold"] >> contrastThreshold;
+  if (!fn["edgeThreshold"].empty())
+    fn["edgeThreshold"] >> edgeThreshold;
+  if (!fn["sigma"].empty())
+    fn["sigma"] >> sigma;
+  if (!fn["descriptorType"].empty())
+    fn["descriptorType"] >> descriptor_type;
+}
+void SIFT_Impl::write( FileStorage& fs) const
+{
+  if(fs.isOpened())
+  {
+    fs << "name" << getDefaultName();
+    fs << "nfeatures" << nfeatures;
+    fs << "nOctaveLayers" << nOctaveLayers;
+    fs << "contrastThreshold" << contrastThreshold;
+    fs << "edgeThreshold" << edgeThreshold;
+    fs << "sigma" << sigma;
+    fs << "descriptorType" << descriptor_type;
+  }
+}
+
+}