Merge pull request #847 from GregorKovalcik:PCTSignatures

modules/xfeatures2d/doc/xfeatures2d.bib

@@ -92,3 +92,22 @@
   booktitle = "Computer Vision and Pattern Recognition",
   year = {2013}
 }
+
+@article{KrulisLS16,
+  author    = {Martin Krulis and Jakub Lokoc and Tomas Skopal},
+  title     = {Efficient extraction of clustering-based feature signatures using {GPU} architectures},
+  journal   = {Multimedia Tools Appl.},
+  volume    = {75},
+  number    = {13},
+  pages     = {8071--8103},
+  year      = {2016}
+}
+
+@inproceedings{BeecksUS10,
+  author    = {Christian Beecks and Merih Seran Uysal and Thomas Seidl},
+  title     = {Signature Quadratic Form Distance},
+  booktitle = {{CIVR}},
+  pages     = {438--445},
+  publisher = {{ACM}},
+  year      = {2010}
+}

modules/xfeatures2d/include/opencv2/xfeatures2d.hpp

@@ -359,6 +359,479 @@ public:
                     bool use_scale_orientation = true, float scale_factor = 6.25f );
 };
 
+
+/*
+* Position-Color-Texture signatures
+*/
+
+/**
+* @brief Class implementing PCT (position-color-texture) signature extraction
+*       as described in @cite KrulisLS16.
+*       The algorithm is divided to a feature sampler and a clusterizer.
+*       Feature sampler produces samples at given set of coordinates.
+*       Clusterizer then produces clusters of these samples using k-means algorithm.
+*       Resulting set of clusters is the signature of the input image.
+*
+*       A signature is an array of SIGNATURE_DIMENSION-dimensional points.
+*       Used dimensions are:
+*       weight, x, y position; lab color, contrast, entropy.
+* @cite KrulisLS16
+* @cite BeecksUS10
+*/
+class CV_EXPORTS_W PCTSignatures : public Algorithm
+{
+public:
+    /**
+    * @brief Lp distance function selector.
+    */
+    enum DistanceFunction
+    {
+        L0_25, L0_5, L1, L2, L2SQUARED, L5, L_INFINITY
+    };
+
+    /**
+    * @brief Point distributions supported by random point generator.
+    */
+    enum PointDistribution
+    {
+        UNIFORM,    //!< Generate numbers uniformly.
+        REGULAR,    //!< Generate points in a regular grid.
+        NORMAL      //!< Generate points with normal (gaussian) distribution.
+    };
+
+    /**
+    * @brief Similarity function selector.
+    * @see
+    *       Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    *       Signature quadratic form distance.
+    *       In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    *       ACM, 2010.
+    * @cite BeecksUS10
+    * @note For selected distance function: \f[ d(c_i, c_j) \f]  and parameter: \f[ \alpha \f]
+    */
+    enum SimilarityFunction
+    {
+        MINUS,      //!< \f[ -d(c_i, c_j) \f]
+        GAUSSIAN,   //!< \f[ e^{ -\alpha * d^2(c_i, c_j)} \f]
+        HEURISTIC   //!< \f[ \frac{1}{\alpha + d(c_i, c_j)} \f]
+    };
+
+
+    /**
+    * @brief Creates PCTSignatures algorithm using sample and seed count.
+    *       It generates its own sets of sampling points and clusterization seed indexes.
+    * @param initSampleCount Number of points used for image sampling.
+    * @param initSeedCount Number of initial clusterization seeds.
+    *       Must be lower or equal to initSampleCount
+    * @param pointDistribution Distribution of generated points. Default: UNIFORM.
+    *       Available: UNIFORM, REGULAR, NORMAL.
+    * @return Created algorithm.
+    */
+    CV_WRAP static Ptr<PCTSignatures> create(
+        const int initSampleCount = 2000,
+        const int initSeedCount = 400,
+        const int pointDistribution = 0);
+
+    /**
+    * @brief Creates PCTSignatures algorithm using pre-generated sampling points
+    *       and number of clusterization seeds. It uses the provided
+    *       sampling points and generates its own clusterization seed indexes.
+    * @param initSamplingPoints Sampling points used in image sampling.
+    * @param initSeedCount Number of initial clusterization seeds.
+    *       Must be lower or equal to initSamplingPoints.size().
+    * @return Created algorithm.
+    */
+    CV_WRAP static Ptr<PCTSignatures> create(
+        const std::vector<Point2f>& initSamplingPoints,
+        const int initSeedCount);
+
+    /**
+    * @brief Creates PCTSignatures algorithm using pre-generated sampling points
+    *       and clusterization seeds indexes.
+    * @param initSamplingPoints Sampling points used in image sampling.
+    * @param initClusterSeedIndexes Indexes of initial clusterization seeds.
+    *       Its size must be lower or equal to initSamplingPoints.size().
+    * @return Created algorithm.
+    */
+    CV_WRAP static Ptr<PCTSignatures> create(
+        const std::vector<Point2f>& initSamplingPoints,
+        const std::vector<int>& initClusterSeedIndexes);
+
+
+
+    /**
+    * @brief Computes signature of given image.
+    * @param image Input image of CV_8U type.
+    * @param signature Output computed signature.
+    */
+    CV_WRAP virtual void computeSignature(
+        InputArray image,
+        OutputArray signature) const = 0;
+
+    /**
+    * @brief Computes signatures for multiple images in parallel.
+    * @param images Vector of input images of CV_8U type.
+    * @param signatures Vector of computed signatures.
+    */
+    CV_WRAP virtual void computeSignatures(
+        const std::vector<Mat>& images,
+        std::vector<Mat>& signatures) const = 0;
+
+    /**
+    * @brief Draws signature in the source image and outputs the result.
+    *       Signatures are visualized as a circle
+    *       with radius based on signature weight
+    *       and color based on signature color.
+    *       Contrast and entropy are not visualized.
+    * @param source Source image.
+    * @param signature Image signature.
+    * @param result Output result.
+    * @param radiusToShorterSideRatio Determines maximal radius of signature in the output image.
+    * @param borderThickness Border thickness of the visualized signature.
+    */
+    CV_WRAP static void drawSignature(
+        InputArray source,
+        InputArray signature,
+        OutputArray result,
+        float radiusToShorterSideRatio = 1.0 / 8,
+        int borderThickness = 1);
+
+    /**
+    * @brief Generates initial sampling points according to selected point distribution.
+    * @param initPoints Output vector where the generated points will be saved.
+    * @param count Number of points to generate.
+    * @param pointDistribution Point distribution selector.
+    *       Available: UNIFORM, REGULAR, NORMAL.
+    * @note Generated coordinates are in range [0..1)
+    */
+    CV_WRAP static void generateInitPoints(
+        std::vector<Point2f>& initPoints,
+        const int count,
+        int pointDistribution);
+
+
+    /**** sampler ****/
+
+    /**
+    * @brief Number of initial samples taken from the image.
+    */
+    CV_WRAP virtual int getSampleCount() const = 0;
+
+    /**
+    * @brief Color resolution of the greyscale bitmap represented in allocated bits
+    *       (i.e., value 4 means that 16 shades of grey are used).
+    *       The greyscale bitmap is used for computing contrast and entropy values.
+    */
+    CV_WRAP virtual int getGrayscaleBits() const = 0;
+    /**
+    * @brief Color resolution of the greyscale bitmap represented in allocated bits
+    *       (i.e., value 4 means that 16 shades of grey are used).
+    *       The greyscale bitmap is used for computing contrast and entropy values.
+    */
+    CV_WRAP virtual void setGrayscaleBits(int grayscaleBits) = 0;
+
+    /**
+    * @brief Size of the texture sampling window used to compute contrast and entropy
+    *       (center of the window is always in the pixel selected by x,y coordinates
+    *       of the corresponding feature sample).
+    */
+    CV_WRAP virtual int getWindowRadius() const = 0;
+    /**
+    * @brief Size of the texture sampling window used to compute contrast and entropy
+    *       (center of the window is always in the pixel selected by x,y coordinates
+    *       of the corresponding feature sample).
+    */
+    CV_WRAP virtual void setWindowRadius(int radius) = 0;
+
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual float getWeightX() const = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual void setWeightX(float weight) = 0;
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual float getWeightY() const = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual void setWeightY(float weight) = 0;
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual float getWeightL() const = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual void setWeightL(float weight) = 0;
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual float getWeightA() const = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual void setWeightA(float weight) = 0;
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual float getWeightB() const = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual void setWeightB(float weight) = 0;
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual float getWeightConstrast() const = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual void setWeightContrast(float weight) = 0;
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual float getWeightEntropy() const = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space
+    *       (x,y = position; L,a,b = color in CIE Lab space; c = contrast. e = entropy)
+    */
+    CV_WRAP virtual void setWeightEntropy(float weight) = 0;
+
+    /**
+    * @brief Initial samples taken from the image.
+    *       These sampled features become the input for clustering.
+    */
+    CV_WRAP virtual std::vector<Point2f> getSamplingPoints() const = 0;
+
+
+
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space.
+    * @param idx ID of the weight
+    * @param value Value of the weight
+    * @note
+    *       WEIGHT_IDX = 0;
+    *       X_IDX = 1;
+    *       Y_IDX = 2;
+    *       L_IDX = 3;
+    *       A_IDX = 4;
+    *       B_IDX = 5;
+    *       CONTRAST_IDX = 6;
+    *       ENTROPY_IDX = 7;
+    */
+    CV_WRAP virtual void setWeight(int idx, float value) = 0;
+    /**
+    * @brief Weights (multiplicative constants) that linearly stretch individual axes of the feature space.
+    * @param weights Values of all weights.
+    * @note
+    *       WEIGHT_IDX = 0;
+    *       X_IDX = 1;
+    *       Y_IDX = 2;
+    *       L_IDX = 3;
+    *       A_IDX = 4;
+    *       B_IDX = 5;
+    *       CONTRAST_IDX = 6;
+    *       ENTROPY_IDX = 7;
+    */
+    CV_WRAP virtual void setWeights(const std::vector<float>& weights) = 0;
+
+    /**
+    * @brief Translations of the individual axes of the feature space.
+    * @param idx ID of the translation
+    * @param value Value of the translation
+    * @note
+    *       WEIGHT_IDX = 0;
+    *       X_IDX = 1;
+    *       Y_IDX = 2;
+    *       L_IDX = 3;
+    *       A_IDX = 4;
+    *       B_IDX = 5;
+    *       CONTRAST_IDX = 6;
+    *       ENTROPY_IDX = 7;
+    */
+    CV_WRAP virtual void setTranslation(int idx, float value) = 0;
+    /**
+    * @brief Translations of the individual axes of the feature space.
+    * @param translations Values of all translations.
+    * @note
+    *       WEIGHT_IDX = 0;
+    *       X_IDX = 1;
+    *       Y_IDX = 2;
+    *       L_IDX = 3;
+    *       A_IDX = 4;
+    *       B_IDX = 5;
+    *       CONTRAST_IDX = 6;
+    *       ENTROPY_IDX = 7;
+    */
+    CV_WRAP virtual void setTranslations(const std::vector<float>& translations) = 0;
+
+    /**
+    * @brief Sets sampling points used to sample the input image.
+    * @param samplingPoints Vector of sampling points in range [0..1)
+    * @note Number of sampling points must be greater or equal to clusterization seed count.
+    */
+    CV_WRAP virtual void setSamplingPoints(std::vector<Point2f> samplingPoints) = 0;
+
+
+
+    /**** clusterizer ****/
+    /**
+    * @brief Initial seeds (initial number of clusters) for the k-means algorithm.
+    */
+    CV_WRAP virtual std::vector<int> getInitSeedIndexes() const = 0;
+    /**
+    * @brief Initial seed indexes for the k-means algorithm.
+    */
+    CV_WRAP virtual void setInitSeedIndexes(std::vector<int> initSeedIndexes) = 0;
+    /**
+    * @brief Number of initial seeds (initial number of clusters) for the k-means algorithm.
+    */
+    CV_WRAP virtual int getInitSeedCount() const = 0;
+
+    /**
+    * @brief Number of iterations of the k-means clustering.
+    *       We use fixed number of iterations, since the modified clustering is pruning clusters
+    *       (not iteratively refining k clusters).
+    */
+    CV_WRAP virtual int getIterationCount() const = 0;
+    /**
+    * @brief Number of iterations of the k-means clustering.
+    *       We use fixed number of iterations, since the modified clustering is pruning clusters
+    *       (not iteratively refining k clusters).
+    */
+    CV_WRAP virtual void setIterationCount(int iterationCount) = 0;
+
+    /**
+    * @brief Maximal number of generated clusters. If the number is exceeded,
+    *       the clusters are sorted by their weights and the smallest clusters are cropped.
+    */
+    CV_WRAP virtual int getMaxClustersCount() const = 0;
+    /**
+    * @brief Maximal number of generated clusters. If the number is exceeded,
+    *       the clusters are sorted by their weights and the smallest clusters are cropped.
+    */
+    CV_WRAP virtual void setMaxClustersCount(int maxClustersCount) = 0;
+
+    /**
+    * @brief This parameter multiplied by the index of iteration gives lower limit for cluster size.
+    *       Clusters containing fewer points than specified by the limit have their centroid dismissed
+    *       and points are reassigned.
+    */
+    CV_WRAP virtual int getClusterMinSize() const = 0;
+    /**
+    * @brief This parameter multiplied by the index of iteration gives lower limit for cluster size.
+    *       Clusters containing fewer points than specified by the limit have their centroid dismissed
+    *       and points are reassigned.
+    */
+    CV_WRAP virtual void setClusterMinSize(int clusterMinSize) = 0;
+
+    /**
+    * @brief Threshold euclidean distance between two centroids.
+    *       If two cluster centers are closer than this distance,
+    *       one of the centroid is dismissed and points are reassigned.
+    */
+    CV_WRAP virtual float getJoiningDistance() const = 0;
+    /**
+    * @brief Threshold euclidean distance between two centroids.
+    *       If two cluster centers are closer than this distance,
+    *       one of the centroid is dismissed and points are reassigned.
+    */
+    CV_WRAP virtual void setJoiningDistance(float joiningDistance) = 0;
+
+    /**
+    * @brief Remove centroids in k-means whose weight is lesser or equal to given threshold.
+    */
+    CV_WRAP virtual float getDropThreshold() const = 0;
+    /**
+    * @brief Remove centroids in k-means whose weight is lesser or equal to given threshold.
+    */
+    CV_WRAP virtual void setDropThreshold(float dropThreshold) = 0;
+
+    /**
+    * @brief Distance function selector used for measuring distance between two points in k-means.
+    */
+    CV_WRAP virtual int getDistanceFunction() const = 0;
+    /**
+    * @brief Distance function selector used for measuring distance between two points in k-means.
+    *       Available: L0_25, L0_5, L1, L2, L2SQUARED, L5, L_INFINITY.
+    */
+    CV_WRAP virtual void setDistanceFunction(int distanceFunction) = 0;
+
+};
+
+/**
+* @brief Class implementing Signature Quadratic Form Distance (SQFD).
+* @see Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+*   Signature quadratic form distance.
+*   In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+*   ACM, 2010.
+* @cite BeecksUS10
+*/
+class CV_EXPORTS_W PCTSignaturesSQFD : public Algorithm
+{
+public:
+
+    /**
+    * @brief Creates the algorithm instance using selected distance function,
+    *       similarity function and similarity function parameter.
+    * @param distanceFunction Distance function selector. Default: L2
+    *       Available: L0_25, L0_5, L1, L2, L2SQUARED, L5, L_INFINITY
+    * @param similarityFunction Similarity function selector. Default: HEURISTIC
+    *       Available: MINUS, GAUSSIAN, HEURISTIC
+    * @param similarityParameter Parameter of the similarity function.
+    */
+    CV_WRAP static Ptr<PCTSignaturesSQFD> create(
+        const int distanceFunction = 3,
+        const int similarityFunction = 2,
+        const float similarityParameter = 1.0f);
+
+    /**
+    * @brief Computes Signature Quadratic Form Distance of two signatures.
+    * @param _signature0 The first signature.
+    * @param _signature1 The second signature.
+    */
+    CV_WRAP virtual float computeQuadraticFormDistance(
+        InputArray _signature0,
+        InputArray _signature1) const = 0;
+
+    /**
+    * @brief Computes Signature Quadratic Form Distance between the reference signature
+    *       and each of the other image signatures.
+    * @param sourceSignature The signature to measure distance of other signatures from.
+    * @param imageSignatures Vector of signatures to measure distance from the source signature.
+    * @param distances Output vector of measured distances.
+    */
+    CV_WRAP virtual void computeQuadraticFormDistances(
+        const Mat& sourceSignature,
+        const std::vector<Mat>& imageSignatures,
+        std::vector<float>& distances) const = 0;
+
+};
+
 //! @}
 
 }

modules/xfeatures2d/samples/pct_signatures.cpp

@@ -0,0 +1,172 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+#include <opencv2/core.hpp>
+#include <opencv2/highgui.hpp>
+#include <opencv2/xfeatures2d.hpp>
+
+#include <iostream>
+#include <string>
+
+using namespace std;
+using namespace cv;
+using namespace xfeatures2d;
+
+
+void printHelpMessage(void);
+void printHelpMessage(void)
+{
+    cout << "Example of the PCTSignatures algorithm computing and visualizing\n"
+        "image signature for one image, or comparing multiple images with the first\n"
+        "image using the signature quadratic form distance.\n\n"
+        "Usage: pct_signatures ImageToProcessAndDisplay\n"
+        "or:    pct_signatures ReferenceImage [ImagesToCompareWithTheReferenceImage]\n\n"
+        "The program has 2 modes:\n"
+        "- single argument: program computes and visualizes the image signature\n"
+        "- multiple arguments: program compares the first image to the others\n"
+        "  using pct signatures and signature quadratic form distance (SQFD)";
+}
+
+/** @brief
+
+Example of the PCTSignatures algorithm.
+
+The program has 2 modes:
+- single argument mode, where the program computes and visualizes the image signature
+- multiple argument mode, where the program compares the first image to the others
+using signatures and signature quadratic form distance (SQFD)
+
+*/
+int main(int argc, char** argv)
+{
+    if (argc < 2)                               // Check arguments
+    {
+        printHelpMessage();
+        return 1;
+    }
+
+    Mat source;
+    source = imread(argv[1]);                   // Read the file
+
+    if (!source.data)                           // Check for invalid input
+    {
+        cerr << "Could not open or find the image: " << argv[1];
+        return -1;
+    }
+
+    Mat signature, result;                      // define variables
+    int initSampleCount = 2000;
+    int initSeedCount = 400;
+    int grayscaleBitsPerPixel = 4;
+    vector<Point2f> initPoints;
+
+    namedWindow("Source", WINDOW_AUTOSIZE);     // Create windows for display.
+    namedWindow("Result", WINDOW_AUTOSIZE);
+
+                                                // create the algorithm
+    PCTSignatures::generateInitPoints(initPoints, initSampleCount, PCTSignatures::UNIFORM);
+    Ptr<PCTSignatures> pctSignatures = PCTSignatures::create(initPoints, initSeedCount);
+    pctSignatures->setGrayscaleBits(grayscaleBitsPerPixel);
+
+                                                // compute and visualize the first image
+    double start = (double)getTickCount();
+    pctSignatures->computeSignature(source, signature);
+    double end = (double)getTickCount();
+    cout << "Signature of the reference image computed in " << (end - start) / (getTickFrequency() * 1.0f) << " seconds." << endl;
+    PCTSignatures::drawSignature(source, signature, result);
+
+    imshow("Source", source);                   // show the result
+    imshow("Result", result);
+
+    if (argc == 2)          // single image -> finish right after the visualization
+    {
+        waitKey(0);         // Wait for user input
+        return 0;
+    }
+    // multiple images -> compare to the first one
+    else
+    {
+        vector<Mat> images;
+        vector<Mat> signatures;
+        vector<float> distances;
+
+        for (int i = 2; i < argc; i++)
+        {
+            Mat image = imread(argv[i]);
+            if (!source.data)                               // Check for invalid input
+            {
+                cerr << "Could not open or find the image: " << argv[i] << std::endl;
+                return 1;
+            }
+            images.push_back(image);
+        }
+
+        pctSignatures->computeSignatures(images, signatures);
+        Ptr<PCTSignaturesSQFD> pctSQFD = PCTSignaturesSQFD::create();
+        pctSQFD->computeQuadraticFormDistances(signature, signatures, distances);
+
+        for (int i = 0; i < (int)(distances.size()); i++)
+        {
+            cout << "Image: " << argv[i + 2] << ", similarity: " << distances[i] << endl;
+        }
+        waitKey(0); // Wait for user input
+    }
+
+    return 0;
+}

modules/xfeatures2d/samples/pct_webcam.cpp

@@ -0,0 +1,170 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+#include <opencv2/core.hpp>
+#include <opencv2/highgui.hpp>
+#include <opencv2/xfeatures2d.hpp>
+
+#include <iostream>
+#include <string>
+
+using namespace std;
+using namespace cv;
+using namespace xfeatures2d;
+
+
+void printHelpMessage(void);
+void printHelpMessage(void)
+{
+    cout << "Example of the PCTSignatures algorithm.\n\n"
+        "This program computes and visualizes position-color-texture signatures\n"
+        "using images from webcam if available.\n\n"
+        "Usage:\n"
+        "pct_webcam [sample_count] [seed_count]\n"
+        "Note: sample_count must be greater or equal to seed_count.";
+}
+
+
+/** @brief
+
+Example of the PCTSignatures algorithm.
+
+This program computes and visualizes position-color-texture signatures
+of images taken from webcam if available.
+*/
+int main(int argc, char** argv)
+{
+    // define variables
+    Mat frame, signature, result;
+    int initSampleCount = 2000;
+    int initSeedCount = 400;
+    int grayscaleBitsPerPixel = 4;
+
+    // parse for help argument
+    {
+        for (int i = 1; i < argc; i++)
+        {
+            if ((string)argv[i] == "-h" || (string)argv[i] == "--help")
+            {
+                printHelpMessage();
+                return 0;
+            }
+        }
+    }
+
+    // parse optional arguments
+    if (argc > 1)               // sample count
+    {
+        initSampleCount = atoi(argv[1]);
+        if (initSampleCount <= 0)
+        {
+            cerr << "Sample count have to be a positive integer: " << argv[1] << endl;
+            return 1;
+        }
+        initSeedCount = (int)floor(initSampleCount / 4);
+        initSeedCount = std::max(1, initSeedCount);     // fallback if sample count == 1
+    }
+    if (argc > 2)               // seed count
+    {
+        initSeedCount = atoi(argv[2]);
+        if (initSeedCount <= 0)
+        {
+            cerr << "Seed count have to be a positive integer: " << argv[2] << endl;
+            return 1;
+        }
+        if (initSeedCount > initSampleCount)
+        {
+            cerr << "Seed count have to be lower or equal to sample count!" << endl;
+            return 1;
+        }
+    }
+
+    // create algorithm
+    Ptr<PCTSignatures> pctSignatures = PCTSignatures::create(initSampleCount, initSeedCount, PCTSignatures::UNIFORM);
+    pctSignatures->setGrayscaleBits(grayscaleBitsPerPixel);
+
+    // open video capture device
+    VideoCapture videoCapture;
+    if (!videoCapture.open(0))
+    {
+        cerr << "Unable to open the first video capture device with ID = 0!" << endl;
+        return 1;
+    }
+
+    // Create windows for display.
+    namedWindow("Source", WINDOW_AUTOSIZE);
+    namedWindow("Result", WINDOW_AUTOSIZE);
+
+    // run drawing loop
+    for (;;)
+    {
+        videoCapture >> frame;
+        if (frame.empty()) break; // end of video stream
+
+        pctSignatures->computeSignature(frame, signature);
+        PCTSignatures::drawSignature(Mat::zeros(frame.size(), frame.type()), signature, result);
+
+        imshow("Source", frame);                // Show our images inside the windows.
+        imshow("Result", result);
+
+        if (waitKey(1) == 27) break;            // stop videocapturing by pressing ESC
+    }
+
+    return 0;
+}

modules/xfeatures2d/src/pct_signatures.cpp

@@ -0,0 +1,443 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+#include "precomp.hpp"
+
+#include "pct_signatures/constants.hpp"
+#include "pct_signatures/pct_sampler.hpp"
+#include "pct_signatures/pct_clusterizer.hpp"
+
+using namespace cv::xfeatures2d::pct_signatures;
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            class PCTSignatures_Impl : public PCTSignatures
+            {
+            public:
+                PCTSignatures_Impl(const std::vector<Point2f>& initSamplingPoints, int initSeedCount)
+                {
+                    if (initSamplingPoints.size() == 0)
+                    {
+                        CV_Error(Error::StsBadArg, "No sampling points provided!");
+                    }
+                    if (initSeedCount <= 0)
+                    {
+                        CV_Error(Error::StsBadArg, "Not enough initial seeds, at least 1 required.");
+                    }
+
+                    mSampler = PCTSampler::create(initSamplingPoints);
+
+                    initSeedCount = std::min(initSeedCount, (int)initSamplingPoints.size());
+                    std::vector<int> initClusterSeedIndexes = pickRandomClusterSeedIndexes(initSeedCount);
+                    mClusterizer = PCTClusterizer::create(initClusterSeedIndexes);
+                }
+
+                PCTSignatures_Impl(
+                    const std::vector<Point2f>& initSamplingPoints,
+                    const std::vector<int>& initClusterSeedIndexes)
+                {
+                    if (initSamplingPoints.size() == 0)
+                    {
+                        CV_Error(Error::StsBadArg, "No sampling points provided!");
+                    }
+                    if (initClusterSeedIndexes.size() == 0)
+                    {
+                        CV_Error(Error::StsBadArg, "Not enough initial seeds, at least 1 required.");
+                    }
+                    if (initClusterSeedIndexes.size() > initSamplingPoints.size())
+                    {
+                        CV_Error(Error::StsBadArg, "Too much cluster seeds or not enough sampling points.");
+                    }
+                    for (int iCluster = 0; iCluster < (int)(initClusterSeedIndexes.size()); iCluster++)
+                    {
+                        if (initClusterSeedIndexes[iCluster] < 0
+                            || initClusterSeedIndexes[iCluster] >= (int)(initSamplingPoints.size()))
+                        {
+                            CV_Error(Error::StsBadArg,
+                                "Initial cluster seed indexes contain an index outside the range of the sampling point list.");
+                        }
+                    }
+
+                    mSampler = PCTSampler::create(initSamplingPoints);
+                    mClusterizer = PCTClusterizer::create(initClusterSeedIndexes);
+                }
+
+                void computeSignature(InputArray image, OutputArray signature) const;
+
+                void computeSignatures(const std::vector<Mat>& images, std::vector<Mat>& signatures) const;
+
+                void getGrayscaleBitmap(OutputArray _grayscaleBitmap, bool normalize) const;
+
+                /**** sampler ****/
+
+                int getSampleCount() const                      { return mSampler->getGrayscaleBits(); }
+                int getGrayscaleBits() const                    { return mSampler->getGrayscaleBits(); }
+                int getWindowRadius() const                     { return mSampler->getWindowRadius(); }
+                float getWeightX() const                        { return mSampler->getWeightX(); }
+                float getWeightY() const                        { return mSampler->getWeightY(); }
+                float getWeightL() const                        { return mSampler->getWeightL(); }
+                float getWeightA() const                        { return mSampler->getWeightA(); }
+                float getWeightB() const                        { return mSampler->getWeightB(); }
+                float getWeightConstrast() const                { return mSampler->getWeightConstrast(); }
+                float getWeightEntropy() const                  { return mSampler->getWeightEntropy(); }
+
+                std::vector<Point2f> getSamplingPoints() const  { return mSampler->getSamplingPoints(); }
+
+
+                void setGrayscaleBits(int grayscaleBits)        { mSampler->setGrayscaleBits(grayscaleBits); }
+                void setWindowRadius(int windowRadius)          { mSampler->setWindowRadius(windowRadius); }
+                void setWeightX(float weight)                   { mSampler->setWeightX(weight); }
+                void setWeightY(float weight)                   { mSampler->setWeightY(weight); }
+                void setWeightL(float weight)                   { mSampler->setWeightL(weight); }
+                void setWeightA(float weight)                   { mSampler->setWeightA(weight); }
+                void setWeightB(float weight)                   { mSampler->setWeightB(weight); }
+                void setWeightContrast(float weight)            { mSampler->setWeightContrast(weight); }
+                void setWeightEntropy(float weight)             { mSampler->setWeightEntropy(weight); }
+
+                void setWeight(int idx, float value)                           { mSampler->setWeight(idx, value); }
+                void setWeights(const std::vector<float>& weights)             { mSampler->setWeights(weights); }
+                void setTranslation(int idx, float value)                      { mSampler->setTranslation(idx, value); }
+                void setTranslations(const std::vector<float>& translations)   { mSampler->setTranslations(translations); }
+
+                void setSamplingPoints(std::vector<Point2f> samplingPoints)    { mSampler->setSamplingPoints(samplingPoints); }
+
+
+                /**** clusterizer ****/
+
+                int getIterationCount() const                       { return mClusterizer->getIterationCount(); }
+                std::vector<int> getInitSeedIndexes() const         { return mClusterizer->getInitSeedIndexes(); }
+                int getInitSeedCount() const                        { return (int)mClusterizer->getInitSeedIndexes().size(); }
+                int getMaxClustersCount() const                     { return mClusterizer->getMaxClustersCount(); }
+                int getClusterMinSize() const                       { return mClusterizer->getClusterMinSize(); }
+                float getJoiningDistance() const                    { return mClusterizer->getJoiningDistance(); }
+                float getDropThreshold() const                      { return mClusterizer->getDropThreshold(); }
+                int getDistanceFunction() const
+                                                                    { return mClusterizer->getDistanceFunction(); }
+
+                void setIterationCount(int iterations)              { mClusterizer->setIterationCount(iterations); }
+                void setInitSeedIndexes(std::vector<int> initSeedIndexes)
+                                                                    { mClusterizer->setInitSeedIndexes(initSeedIndexes); }
+                void setMaxClustersCount(int maxClusters)           { mClusterizer->setMaxClustersCount(maxClusters); }
+                void setClusterMinSize(int clusterMinSize)          { mClusterizer->setClusterMinSize(clusterMinSize); }
+                void setJoiningDistance(float joiningDistance)      { mClusterizer->setJoiningDistance(joiningDistance); }
+                void setDropThreshold(float dropThreshold)          { mClusterizer->setDropThreshold(dropThreshold); }
+                void setDistanceFunction(int distanceFunction)
+                                                                    { mClusterizer->setDistanceFunction(distanceFunction); }
+
+            private:
+                /**
+                * @brief Samples used for sampling the input image and producing list of samples.
+                */
+                Ptr<PCTSampler> mSampler;
+
+                /**
+                * @brief Clusterizer using k-means algorithm to produce list of centroids from sampled points - the image signature.
+                */
+                Ptr<PCTClusterizer> mClusterizer;
+
+                /**
+                * @brief Creates vector of random indexes of sampling points
+                *       which will be used as initial centroids for k-means clusterization.
+                * @param initSeedCount Number of indexes of initial centroids to be produced.
+                * @return The generated vector of random indexes.
+                */
+                static std::vector<int> pickRandomClusterSeedIndexes(int initSeedCount)
+                {
+                    std::vector<int> seedIndexes;
+                    for (int iSeed = 0; iSeed < initSeedCount; iSeed++)
+                    {
+                        seedIndexes.push_back(iSeed);
+                    }
+
+                    randShuffle(seedIndexes);
+                    return seedIndexes;
+                }
+            };
+
+
+            /**
+            * @brief Class implementing parallel computing of signatures for multiple images.
+            */
+            class Parallel_computeSignatures : public ParallelLoopBody
+            {
+            private:
+                const PCTSignatures* mPctSignaturesAlgorithm;
+                const std::vector<Mat>* mImages;
+                std::vector<Mat>* mSignatures;
+
+            public:
+                Parallel_computeSignatures(
+                    const PCTSignatures* pctSignaturesAlgorithm,
+                    const std::vector<Mat>* images,
+                    std::vector<Mat>* signatures)
+                    : mPctSignaturesAlgorithm(pctSignaturesAlgorithm),
+                    mImages(images),
+                    mSignatures(signatures)
+                {
+                    mSignatures->resize(images->size());
+                }
+
+                void operator()(const Range& range) const
+                {
+                    for (int i = range.start; i < range.end; i++)
+                    {
+                        mPctSignaturesAlgorithm->computeSignature((*mImages)[i], (*mSignatures)[i]);
+                    }
+                }
+            };
+
+
+            /**
+            * @brief Computes signature for one image.
+            */
+            void PCTSignatures_Impl::computeSignature(InputArray _image, OutputArray _signature) const
+            {
+                if (_image.empty())
+                {
+                    _signature.create(_image.size(), CV_32FC1);
+                    return;
+                }
+
+                Mat image = _image.getMat();
+                CV_Assert(image.depth() == CV_8U);
+
+                // TODO: OpenCL
+                //if (ocl::useOpenCL())
+                //{
+                //}
+
+                // sample features
+                Mat samples;
+                mSampler->sample(image, samples);               // HOT PATH: 40%
+
+                // kmeans clusterize, use feature samples, produce signature clusters
+                Mat signature;
+                mClusterizer->clusterize(samples, signature);   // HOT PATH: 60%
+
+
+                // set result
+                _signature.create(signature.size(), signature.type());
+                Mat result = _signature.getMat();
+                signature.copyTo(result);
+            }
+
+            /**
+            * @brief Parallel function computing signatures for multiple images.
+            */
+            void PCTSignatures_Impl::computeSignatures(const std::vector<Mat>& images, std::vector<Mat>& signatures) const
+            {
+                parallel_for_(Range(0, (int)images.size()), Parallel_computeSignatures(this, &images, &signatures));
+            }
+
+        } // end of namespace pct_signatures
+
+
+
+        Ptr<PCTSignatures> PCTSignatures::create(
+            const int initSampleCount,
+            const int initSeedCount,
+            const int pointDistribution)
+        {
+            std::vector<Point2f> initPoints;
+            generateInitPoints(initPoints, initSampleCount, pointDistribution);
+            return create(initPoints, initSeedCount);
+        }
+
+        Ptr<PCTSignatures> PCTSignatures::create(
+            const std::vector<Point2f>& initPoints,
+            const int initSeedCount)
+        {
+            return makePtr<PCTSignatures_Impl>(initPoints, initSeedCount);
+        }
+
+        Ptr<PCTSignatures> PCTSignatures::create(
+            const std::vector<Point2f>& initPoints,
+            const std::vector<int>& initClusterSeedIndexes)
+        {
+            return makePtr<PCTSignatures_Impl>(initPoints, initClusterSeedIndexes);
+        }
+
+
+        void PCTSignatures::drawSignature(
+            InputArray _source,
+            InputArray _signature,
+            OutputArray _result,
+            float radiusToShorterSideRatio,
+            int borderThickness)
+        {
+            // check source
+            if (_source.empty())
+            {
+                return;
+            }
+            Mat source = _source.getMat();
+
+            // create result
+            _result.create(source.size(), source.type());
+            Mat result = _result.getMat();
+            source.copyTo(result);
+
+            // check signature
+            if (_signature.empty())
+            {
+                return;
+            }
+            Mat signature = _signature.getMat();
+            if (signature.type() != CV_32F || signature.cols != SIGNATURE_DIMENSION)
+            {
+                CV_Error_(Error::StsBadArg, ("Invalid signature format. Type must be CV_32F and signature.cols must be %d.", SIGNATURE_DIMENSION));
+            }
+
+            // compute max radius using given ratio of shorter image side
+            float maxRadius = ((source.rows < source.cols) ? source.rows : source.cols) * radiusToShorterSideRatio;
+
+            // draw signature
+            for (int i = 0; i < signature.rows; i++)
+            {
+                Vec3f labColor(
+                    signature.at<float>(i, L_IDX) * L_COLOR_RANGE,      // convert Lab pixel to BGR
+                    signature.at<float>(i, A_IDX) * A_COLOR_RANGE,
+                    signature.at<float>(i, B_IDX) * B_COLOR_RANGE);
+                Mat labPixel(1, 1, CV_32FC3);
+                labPixel.at<Vec3f>(0, 0) = labColor;
+                Mat rgbPixel;
+                cvtColor(labPixel, rgbPixel, COLOR_Lab2BGR);
+                rgbPixel.convertTo(rgbPixel, CV_8UC3, 255);
+                Vec3b rgbColor = rgbPixel.at<Vec3b>(0, 0);
+
+                // precompute variables
+                Point center((int)(signature.at<float>(i, X_IDX) * source.cols), (int)(signature.at<float>(i, Y_IDX) * source.rows));
+                int radius = (int)(maxRadius * signature.at<float>(i, WEIGHT_IDX));
+                Vec3b borderColor(0, 0, 0);
+
+                // draw filled circle
+                circle(result, center, radius, rgbColor, -1);
+                // draw circle outline
+                circle(result, center, radius, borderColor, borderThickness);
+            }
+        }
+
+
+        void PCTSignatures::generateInitPoints(
+            std::vector<Point2f>& initPoints,
+            const int count,
+            const int pointDistribution)
+        {
+            RNG random;
+            random.state = getTickCount();
+            initPoints.resize(count);
+
+            switch (pointDistribution)
+            {
+            case UNIFORM:
+                for (int i = 0; i < count; i++)
+                {
+                    // returns uniformly distributed float random number from [0, 1) range
+                    initPoints[i] = (Point2f(random.uniform((float)0.0, (float)1.0), random.uniform((float)0.0, (float)1.0)));
+                }
+                break;
+            case REGULAR:
+            {
+                int gridSize = (int)ceil(sqrt(count));
+                const float step = 1.0f / gridSize;
+                const float halfStep = step / 2;
+                float x = halfStep;
+                float y = halfStep;
+                for (int i = 0; i < count; i++)
+                {
+                    // returns regular grid
+                    initPoints[i] = Point2f(x, y);
+                    if ((i + 1) % gridSize == 0)
+                    {
+                        x = halfStep;
+                        y += step;
+                    }
+                    else
+                    {
+                        x += step;
+                    }
+                }
+                break;
+            }
+            case NORMAL:
+                for (int i = 0; i < count; i++)
+                {
+                    // returns normally distributed float random number from (0, 1) range with mean 0.5
+                    float sigma = 0.2f;
+                    float x = (float)random.gaussian(sigma);
+                    float y = (float)random.gaussian(sigma);
+                    while (x <= -0.5f || x >= 0.5f)
+                        x = (float)random.gaussian(sigma);
+                    while (y <= -0.5f || y >= 0.5f)
+                        y = (float)random.gaussian(sigma);
+                    initPoints[i] = Point2f(x, y) + Point2f(0.5, 0.5);
+                }
+                break;
+            default:
+                CV_Error(Error::StsNotImplemented, "Generation of this init point distribution is not implemented!");
+                break;
+            }
+        }
+
+
+    } // end of namespace xfeatures2d
+} // end of namespace cv

modules/xfeatures2d/src/pct_signatures/constants.hpp

@@ -0,0 +1,96 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+
+#ifndef _OPENCV_XFEATURES_2D_PCT_SIGNATURES_CONSTANTS_HPP_
+#define _OPENCV_XFEATURES_2D_PCT_SIGNATURES_CONSTANTS_HPP_
+
+#ifdef __cplusplus
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            const int SIGNATURE_DIMENSION = 8;
+
+            const int WEIGHT_IDX = 0;
+            const int X_IDX = 1;
+            const int Y_IDX = 2;
+            const int L_IDX = 3;
+            const int A_IDX = 4;
+            const int B_IDX = 5;
+            const int CONTRAST_IDX = 6;
+            const int ENTROPY_IDX = 7;
+
+            const float L_COLOR_RANGE = 100;
+            const float A_COLOR_RANGE = 127;
+            const float B_COLOR_RANGE = 127;
+
+            const float SAMPLER_CONTRAST_NORMALIZER = 25.0;   // both determined empirically
+            const float SAMPLER_ENTROPY_NORMALIZER = 4.0;
+
+        }
+    }
+}
+
+#endif
+
+#endif

modules/xfeatures2d/src/pct_signatures/distance.hpp

@@ -0,0 +1,221 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+
+#ifndef _OPENCV_XFEATURES_2D_PCT_SIGNATURES_DISTANCE_HPP_
+#define _OPENCV_XFEATURES_2D_PCT_SIGNATURES_DISTANCE_HPP_
+
+#ifdef __cplusplus
+#include "precomp.hpp"
+
+#include "constants.hpp"
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+
+            static inline float distanceL0_25(
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float result = (float)0.0;
+                for (int d = 1; d < SIGNATURE_DIMENSION; ++d)
+                {
+                    float difference = points1.at<float>(idx1, d) - points2.at<float>(idx2, d);
+                    result += std::sqrt(std::sqrt(std::abs(difference)));
+                }
+                result *= result;
+                return result * result;
+            }
+
+
+            static inline float distanceL0_5(
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float result = (float)0.0;
+                for (int d = 1; d < SIGNATURE_DIMENSION; ++d)
+                {
+                    float difference = points1.at<float>(idx1, d) - points2.at<float>(idx2, d);
+                    result += std::sqrt(std::abs(difference));
+                }
+                return result * result;
+            }
+
+
+            static inline float distanceL1(
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float result = (float)0.0;
+                for (int d = 1; d < SIGNATURE_DIMENSION; ++d)
+                {
+                    float difference = points1.at<float>(idx1, d) - points2.at<float>(idx2, d);
+                    result += std::abs(difference);
+                }
+                return result;
+            }
+
+
+            static inline float distanceL2(
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float result = (float)0.0;
+                for (int d = 1; d < SIGNATURE_DIMENSION; ++d)
+                {
+                    float difference = points1.at<float>(idx1, d) - points2.at<float>(idx2, d);
+                    result += difference * difference;
+                }
+                return (float)std::sqrt(result);
+            }
+
+
+            static inline float distanceL2Squared(
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float result = (float)0.0;
+                for (int d = 1; d < SIGNATURE_DIMENSION; ++d)
+                {
+                    float difference = points1.at<float>(idx1, d) - points2.at<float>(idx2, d);
+                    result += difference * difference;
+                }
+                return result;
+            }
+
+
+            static inline float distanceL5(
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float result = (float)0.0;
+                for (int d = 1; d < SIGNATURE_DIMENSION; ++d)
+                {
+                    float difference = points1.at<float>(idx1, d) - points2.at<float>(idx2, d);
+                    result += std::abs(difference) * difference * difference * difference * difference;
+                }
+                return std::pow(result, (float)0.2);
+            }
+
+
+            static inline float distanceLInfinity(
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float result = (float)0.0;
+                for (int d = 1; d < SIGNATURE_DIMENSION; ++d)
+                {
+                    float difference = points1.at<float>(idx1, d) - points2.at<float>(idx2, d);
+                    if (difference > result)
+                    {
+                        result = difference;
+                    }
+                }
+                return result;
+            }
+
+
+            /**
+            * @brief Computed distance between two centroids using given distance function.
+            * @param distanceFunction Distance function selector.
+            * @param points1 The first signature matrix - one centroid in each row.
+            * @param idx1 ID of centroid in the first signature
+            * @param points2 The second signature matrix - one centroid in each row.
+            * @param idx2 ID of centroid in the first signature
+            * @note The first column of a signature contains weights,
+            *       so only rows 1 to SIGNATURE_DIMENSION are used.
+            */
+            static inline float computeDistance(
+                const int distanceFunction,
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                switch (distanceFunction)
+                {
+                case PCTSignatures::L0_25:
+                    return distanceL0_25(points1, idx1, points2, idx2);
+                case PCTSignatures::L0_5:
+                    return distanceL0_5(points1, idx1, points2, idx2);
+                case PCTSignatures::L1:
+                    return distanceL1(points1, idx1, points2, idx2);
+                case PCTSignatures::L2:
+                    return distanceL2(points1, idx1, points2, idx2);
+                case PCTSignatures::L2SQUARED:
+                    return distanceL2Squared(points1, idx1, points2, idx2);
+                case PCTSignatures::L5:
+                    return distanceL5(points1, idx1, points2, idx2);
+                case PCTSignatures::L_INFINITY:
+                    return distanceLInfinity(points1, idx1, points2, idx2);
+                default:
+                    CV_Error(Error::StsBadArg, "Distance function not implemented!");
+                    return -1;
+                }
+            }
+        }
+    }
+}
+
+
+#endif
+
+#endif

modules/xfeatures2d/src/pct_signatures/grayscale_bitmap.cpp

@@ -0,0 +1,181 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+#include "precomp.hpp"
+
+#include "grayscale_bitmap.hpp"
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            GrayscaleBitmap::GrayscaleBitmap(InputArray _bitmap, int bitsPerPixel)
+                : mBitsPerPixel(bitsPerPixel)
+            {
+                Mat bitmap = _bitmap.getMat();
+                if (bitmap.empty())
+                {
+                    CV_Error(Error::StsBadArg, "Input bitmap is empty");
+                }
+                if (bitmap.depth() != CV_8U && bitmap.depth() != CV_16U)
+                {
+                    CV_Error(Error::StsUnsupportedFormat, "Input bitmap depth must be CV_8U or CV_16U");
+                }
+                if (bitmap.depth() == CV_8U)
+                {
+                    bitmap.convertTo(bitmap, CV_16U, 257);
+                }
+
+                Mat grayscaleBitmap;
+                cvtColor(bitmap, grayscaleBitmap, COLOR_BGR2GRAY);
+
+                mWidth = bitmap.cols;
+                mHeight = bitmap.rows;
+
+                if (bitsPerPixel <= 0 || bitsPerPixel > 8)
+                {
+                    CV_Error_(Error::StsBadArg, ("Invalid number of bits per pixel %d. Only values in range [1..8] are accepted.", bitsPerPixel));
+                }
+
+                // Allocate space for pixel data.
+                int pixelsPerItem = 32 / mBitsPerPixel;
+                mData.resize((mWidth*mHeight + pixelsPerItem - 1) / pixelsPerItem);
+
+                // Convert the bitmap to grayscale and fill the pixel data.
+                CV_Assert(grayscaleBitmap.depth() == CV_16U);
+                for (int y = 0; y < mHeight; y++)
+                {
+                    for (int x = 0; x < mWidth; x++)
+                    {
+                        uint grayVal = ((uint)grayscaleBitmap.at<ushort>((int)y, (int)x)) >> (16 - mBitsPerPixel);
+                        setPixel(x, y, grayVal);
+                    }
+                }
+                // Prepare the preallocated contrast matrix for contrast-entropy computations
+                mCoOccurrenceMatrix.resize(1 << (mBitsPerPixel * 2));   // mCoOccurrenceMatrix size = maxPixelValue^2
+            }
+
+
+            // HOT PATH 30%
+            void GrayscaleBitmap::getContrastEntropy(int x, int y, float& contrast, float& entropy, int radius)
+            {
+                int fromX = (x > radius) ? x - radius : 0;
+                int fromY = (y > radius) ? y - radius : 0;
+                int toX = std::min<int>(mWidth - 1, x + radius + 1);
+                int toY = std::min<int>(mHeight - 1, y + radius + 1);
+                for (int j = fromY; j < toY; ++j)
+                {
+                    for (int i = fromX; i < toX; ++i)                               // for each pixel in the window
+                    {
+                        updateCoOccurrenceMatrix(getPixel(i, j), getPixel(i, j + 1));        // match every pixel with all 8 its neighbours
+                        updateCoOccurrenceMatrix(getPixel(i, j), getPixel(i + 1, j));
+                        updateCoOccurrenceMatrix(getPixel(i, j), getPixel(i + 1, j + 1));
+                        updateCoOccurrenceMatrix(getPixel(i + 1, j), getPixel(i, j + 1));    // 4 updates per pixel in the window
+                    }
+                }
+
+                contrast = 0.0;
+                entropy = 0.0;
+
+                uint pixelsScale = 1 << mBitsPerPixel;
+                float normalizer = (float)((toX - fromX) * (toY - fromY) * 4);            // four increments per pixel in the window (see above)
+                for (int j = 0; j < (int)pixelsScale; ++j)                              // iterate row in a 2D histogram
+                {
+                    for (int i = 0; i <= j; ++i)                                        // iterate column up to the diagonal in 2D histogram
+                    {
+                        if (mCoOccurrenceMatrix[j*pixelsScale + i] != 0)                // consider only non-zero values
+                        {
+                            float value = (float)mCoOccurrenceMatrix[j*pixelsScale + i] / normalizer; // normalize value by number of histogram updates
+                            contrast += (i - j) * (i - j) * value;          // compute contrast
+                            entropy -= value * std::log(value);             // compute entropy
+                            mCoOccurrenceMatrix[j*pixelsScale + i] = 0;     // clear the histogram array for the next computation
+                        }
+                    }
+                }
+            }
+
+
+            void GrayscaleBitmap::convertToMat(OutputArray _bitmap, bool normalize) const
+            {
+                _bitmap.create((int)getHeight(), (int)getWidth(), CV_8U);
+                Mat bitmap = _bitmap.getMat();
+
+                for (int y = 0; y < getHeight(); ++y)
+                {
+                    for (int x = 0; x < getWidth(); ++x)
+                    {
+                        uint pixel = getPixel(x, y);
+                        if (normalize && mBitsPerPixel < 8)
+                        {
+                            pixel <<= 8 - mBitsPerPixel;
+                        }
+                        else if (normalize && mBitsPerPixel > 8)
+                        {
+                            pixel >>= mBitsPerPixel - 8;
+                        }
+                        bitmap.at<uchar>((int)y, (int)x) = (uchar)pixel;
+                    }
+                }
+            }
+
+        }
+    }
+}

modules/xfeatures2d/src/pct_signatures/grayscale_bitmap.hpp

@@ -0,0 +1,219 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+
+#ifndef _OPENCV_XFEATURES_2D_PCT_SIGNATURES_GRAYSCALE_BITMAP_HPP_
+#define _OPENCV_XFEATURES_2D_PCT_SIGNATURES_GRAYSCALE_BITMAP_HPP_
+
+#ifdef __cplusplus
+
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            /**
+            * @brief Specific implementation of grayscale bitmap. This bitmap is used
+            *       to compute contrast and entropy features from surroundings of a pixel.
+            */
+            class GrayscaleBitmap
+            {
+            public:
+                /**
+                * @brief Initialize the grayscale bitmap from regular bitmap.
+                * @param bitmap Bitmap used as source of data.
+                * @param bitsPerPixel How many bits occupy one pixel in grayscale (e.g., 8 ~ 256 grayscale values).
+                *       Must be within [1..8] range.
+                */
+                GrayscaleBitmap(InputArray bitmap, int bitsPerPixel = 4);
+
+                /**
+                * @brief Return the width of the image in pixels.
+                */
+                int getWidth() const
+                {
+                    return mWidth;
+                }
+
+                /**
+                * @brief Return the height of the image in pixels.
+                */
+                int getHeight() const
+                {
+                    return mHeight;
+                }
+
+                /**
+                * @brief Return the height of the image in pixels.
+                */
+                int getBitsPerPixel() const
+                {
+                    return mBitsPerPixel;
+                }
+
+
+                /**
+                * @brief Compute contrast and entropy at selected coordinates.
+                * @param x The horizontal coordinate of the pixel (0..width-1).
+                * @param y The vertical coordinate of the pixel (0..height-1).
+                * @param contrast Output variable where contrast value is saved.
+                * @param entropy Output variable where entropy value is saved.
+                * @param radius Radius of the rectangular window around selected pixel used for computing
+                *       contrast and entropy. Size of the window side is (2*radius + 1).
+                *       The window is cropped if [x,y] is too near the image border.
+                */
+                void getContrastEntropy(
+                    int x,
+                    int y,
+                    float& contrast,
+                    float& entropy,
+                    int windowRadius = 3);
+
+                /**
+                * @brief Converts to OpenCV CV_8U Mat for debug and visualization purposes.
+                * @param bitmap OutputArray proxy where Mat will be written.
+                * @param normalize Whether to normalize the bitmap to the whole range of 255 values
+                *       to improve contrast.
+                */
+                void convertToMat(OutputArray bitmap, bool normalize = true) const;
+
+
+            private:
+                /**
+                * @brief Width of the image.
+                */
+                int mWidth;
+
+                /**
+                * @brief Height of the image.
+                */
+                int mHeight;
+
+                /**
+                * @brief Number of bits per pixel.
+                */
+                int mBitsPerPixel;
+
+                /**
+                * @brief Pixel data packed in 32-bit uints.
+                */
+                std::vector<uint> mData;
+
+                /**
+                * @brief Tmp matrix used for computing contrast and entropy.
+                */
+                std::vector<uint> mCoOccurrenceMatrix;
+
+
+                /**
+                * @brief Get pixel from packed data vector.
+                * @param x The horizontal coordinate of the pixel (0..width-1).
+                * @param y The vertical coordinate of the pixel (0..height-1).
+                * @return Grayscale value (0 ~ black, 2^bitPerPixel - 1 ~ white).
+                */
+                uint inline getPixel(int x, int y) const
+                {
+                    int pixelsPerItem = 32 / mBitsPerPixel;
+                    int offset = y*mWidth + x;
+                    int shift = (offset % pixelsPerItem) * mBitsPerPixel;
+                    uint mask = (1 << mBitsPerPixel) - 1;
+                    return (mData[offset / pixelsPerItem] >> shift) & mask;
+                }
+
+                /**
+                * @brief Set pixel in the packed data vector.
+                * @param x The horizontal coordinate of the pixel (0..width-1).
+                * @param y The vertical coordinate of the pixel (0..height-1).
+                * @param val Grayscale value (0 ~ black, 2^bitPerPixel - 1 ~ white).
+                */
+                void inline setPixel(int x, int y, uint val)
+                {
+                    int pixelsPerItem = 32 / mBitsPerPixel;
+                    int offset = y*mWidth + x;
+                    int shift = (offset % pixelsPerItem) * mBitsPerPixel;
+                    uint mask = (1 << mBitsPerPixel) - 1;
+                    val &= mask;
+                    mData[offset / pixelsPerItem] &= ~(mask << shift);
+                    mData[offset / pixelsPerItem] |= val << shift;
+                }
+
+                /**
+                * @brief Perform an update of contrast matrix.
+                */
+                void inline updateCoOccurrenceMatrix(uint a, uint b)
+                {
+                    // co-occurrence matrix is symmetric
+                    // merge to a variable with greater higher bits
+                    // to accumulate just in upper triangle in co-occurrence matrix for efficiency
+                    int offset = (int)((a > b) ? (a << mBitsPerPixel) + b : a + (b << mBitsPerPixel));
+                    mCoOccurrenceMatrix[offset]++;
+                }
+
+
+            };
+        }
+    }
+}
+
+
+#endif
+
+#endif

modules/xfeatures2d/src/pct_signatures/pct_clusterizer.cpp

@@ -0,0 +1,453 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+#include "precomp.hpp"
+
+#include "opencv2/core/core_c.h"    // <- because CV_REDUCE_SUM was undeclared without it
+#include "pct_clusterizer.hpp"
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            class PCTClusterizer_Impl : public PCTClusterizer
+            {
+            public:
+
+                PCTClusterizer_Impl(
+                    const std::vector<int>& initSeedIndexes,
+                    int iterationCount,
+                    int maxClustersCount,
+                    int clusterMinSize,
+                    float joiningDistance,
+                    float dropThreshold,
+                    int distanceFunction)
+                    : mInitSeedIndexes(initSeedIndexes),
+                    mIterationCount(iterationCount),
+                    mMaxClustersCount(maxClustersCount),
+                    mClusterMinSize(clusterMinSize),
+                    mJoiningDistance(joiningDistance),
+                    mDropThreshold(dropThreshold),
+                    mDistanceFunction(distanceFunction)
+                {
+                }
+
+                ~PCTClusterizer_Impl()
+                {
+                }
+
+
+                int getIterationCount() const                       { return mIterationCount; }
+                std::vector<int> getInitSeedIndexes() const         { return mInitSeedIndexes; }
+                int getMaxClustersCount() const                     { return mMaxClustersCount; }
+                int getClusterMinSize() const                       { return mClusterMinSize; }
+                float getJoiningDistance() const                    { return mJoiningDistance; }
+                float getDropThreshold() const                      { return mDropThreshold; }
+                int getDistanceFunction() const                     { return mDistanceFunction; }
+
+                void setIterationCount(int iterationCount)                  { mIterationCount = iterationCount; }
+                void setInitSeedIndexes(std::vector<int> initSeedIndexes)   { mInitSeedIndexes = initSeedIndexes; }
+                void setMaxClustersCount(int maxClustersCount)              { mMaxClustersCount = maxClustersCount; }
+                void setClusterMinSize(int clusterMinSize)                  { mClusterMinSize = clusterMinSize; }
+                void setJoiningDistance(float joiningDistance)              { mJoiningDistance = joiningDistance; }
+                void setDropThreshold(float dropThreshold)                  { mDropThreshold = dropThreshold; }
+                void setDistanceFunction(int distanceFunction)              { mDistanceFunction = distanceFunction; }
+
+
+                /**
+                * @brief K-means algorithm over the sampled points producing centroids as signatures.
+                * @param samples List of sampled points.
+                * @param signature Output list of computed centroids - the signature of the image.
+                */
+                void clusterize(InputArray _samples, OutputArray _signature)
+                {
+                    CV_Assert(!_samples.empty());
+
+                    // prepare matrices
+                    Mat samples = _samples.getMat();
+
+                    if ((int)(this->mInitSeedIndexes.size()) > samples.rows)
+                    {
+                        CV_Error_(Error::StsBadArg, ("Number of seeds %d must be less or equal to the number of samples %d.",
+                            mInitSeedIndexes.size(), samples.rows));
+                    }
+
+                    // Prepare initial centroids.
+                    Mat clusters;
+                    pickRandomClusters(samples, clusters);
+                    clusters(Rect(WEIGHT_IDX, 0, 1, clusters.rows)) = 1;    // set initial weight to 1
+
+
+                    // prepare for iterating
+                    joinCloseClusters(clusters);
+                    dropLightPoints(clusters);
+
+
+                    // Main iterations cycle. Our implementation has fixed number of iterations.
+                    for (int iteration = 0; iteration < mIterationCount; iteration++)
+                    {
+                        // Prepare space for new centroid values.
+                        Mat tmpCentroids(clusters.size(), clusters.type());
+                        tmpCentroids = 0;
+
+                        // Clear weights for new iteration.
+                        clusters(Rect(WEIGHT_IDX, 0, 1, clusters.rows)) = 0;
+
+                        // Compute affiliation of points and sum new coordinates for centroids.
+                        for (int iSample = 0; iSample < samples.rows; iSample++)
+                        {
+                            int iClosest = findClosestCluster(clusters, samples, iSample);
+                            for (int iDimension = 1; iDimension < SIGNATURE_DIMENSION; iDimension++)
+                            {
+                                tmpCentroids.at<float>(iClosest, iDimension) += samples.at<float>(iSample, iDimension);
+                            }
+                            clusters.at<float>(iClosest, WEIGHT_IDX)++;
+                        }
+
+                        // Compute average from tmp coordinates and throw away too small clusters.
+                        int lastIdx = 0;
+                        for (int i = 0; (int)i < tmpCentroids.rows; ++i)
+                        {
+                            // Skip clusters that are too small (large-enough clusters are packed right away)
+                            if (clusters.at<float>(i, WEIGHT_IDX) >(iteration + 1) * this->mClusterMinSize)
+                            {
+                                for (int d = 1; d < SIGNATURE_DIMENSION; d++)
+                                {
+                                    clusters.at<float>(lastIdx, d) = tmpCentroids.at<float>(i, d) / clusters.at<float>(i, WEIGHT_IDX);
+                                }
+                                // weights must be compacted as well
+                                clusters.at<float>(lastIdx, WEIGHT_IDX) = clusters.at<float>(i, WEIGHT_IDX);
+                                lastIdx++;
+                            }
+                        }
+
+                        // Crop the arrays if some centroids were dropped.
+                        clusters.resize(lastIdx);
+                        if (clusters.rows == 0)
+                        {
+                            break;
+                        }
+
+                        // Finally join clusters with too close centroids.
+                        joinCloseClusters(clusters);
+                        dropLightPoints(clusters);
+                    }
+
+                    // The result must not be empty!
+                    if (clusters.rows == 0)
+                    {
+                        singleClusterFallback(samples, clusters);
+                    }
+
+                    cropClusters(clusters);
+
+                    normalizeWeights(clusters);
+
+                    // save the result
+                    _signature.create(clusters.rows, SIGNATURE_DIMENSION, clusters.type());
+                    Mat signature = _signature.getMat();
+                    clusters.copyTo(signature);
+                }
+
+
+            private:
+
+                /**
+                * @brief Join clusters that are closer than joining distance.
+                *       If two clusters are joined one of them gets its weight set to 0.
+                * @param clusters List of clusters to be scaned and joined.
+                */
+                void joinCloseClusters(Mat& clusters)
+                {
+                    for (int i = 0; i < clusters.rows - 1; i++)
+                    {
+                        if (clusters.at<float>(i, WEIGHT_IDX) == 0)
+                        {
+                            continue;
+                        }
+
+                        for (int j = i + 1; j < clusters.rows; j++)
+                        {
+                            if (clusters.at<float>(j, WEIGHT_IDX) > 0
+                                && computeDistance(mDistanceFunction, clusters, i, clusters, j) <= mJoiningDistance)
+                            {
+                                clusters.at<float>(i, WEIGHT_IDX) = 0;
+                                break;
+                            }
+                        }
+                    }
+                }
+
+
+                /**
+                * @brief Remove points whose weight is lesser or equal to given threshold.
+                *       The point list is compacted and relative order of points is maintained.
+                * @param clusters List of clusters to be scaned and dropped.
+                */
+                void dropLightPoints(Mat& clusters)
+                {
+                    int frontIdx = 0;
+
+                    // Skip leading continuous part of weighted-enough points.
+                    while (frontIdx < clusters.rows && clusters.at<float>(frontIdx, WEIGHT_IDX) > mDropThreshold)
+                    {
+                        ++frontIdx;
+                    }
+
+                    // Mark first emptied position and advance front index.
+                    int tailIdx = frontIdx++;
+
+                    while (frontIdx < clusters.rows)
+                    {
+                        if (clusters.at<float>(frontIdx, WEIGHT_IDX) > mDropThreshold)
+                        {
+                            // Current (front) item is not dropped -> copy it to the tail.
+                            clusters.row(frontIdx).copyTo(clusters.row(tailIdx));
+                            ++tailIdx;  // grow the tail
+                        }
+                        ++frontIdx;
+                    }
+                    clusters.resize(tailIdx);
+                }
+
+
+                /**
+                * @brief Find closest cluster to selected point.
+                * @param clusters List of cluster centroids.
+                * @param points List of points.
+                * @param pointIdx Index to the list of points (the point for which the closest cluster is being found).
+                * @return Index to clusters list pointing at the closest cluster.
+                */
+                int findClosestCluster(const Mat& clusters, const Mat& points, const int pointIdx) const    // HOT PATH: 35%
+                {
+                    int iClosest = 0;
+                    float minDistance = computeDistance(mDistanceFunction, clusters, 0, points, pointIdx);
+
+                    for (int iCluster = 1; iCluster < clusters.rows; iCluster++)
+                    {
+                        float distance = computeDistance(mDistanceFunction, clusters, iCluster, points, pointIdx);
+                        if (distance < minDistance)
+                        {
+                            iClosest = iCluster;
+                            minDistance = distance;
+                        }
+                    }
+                    return iClosest;
+                }
+
+
+                /**
+                * @brief Make sure that the number of clusters does not exceed maxClusters parameter.
+                *       If it does, the clusters are sorted by their weights and the smallest clusters
+                *       are cropped.
+                * @param clusters List of clusters being scanned and cropped.
+                * @note This method does nothing iff the number of clusters is within the range.
+                */
+                void cropClusters(Mat& clusters) const
+                {
+                    if (clusters.rows > mMaxClustersCount)
+                    {
+                        Mat duplicate(clusters);                // save original clusters
+                        Mat sortedIdx;                          // sort using weight column
+                        sortIdx(clusters(Rect(WEIGHT_IDX, 0, 1, clusters.rows)), sortedIdx, SORT_EVERY_COLUMN + SORT_DESCENDING);
+
+                        clusters.resize(mMaxClustersCount);     // crop to max clusters
+                        for (int i = 0; i < mMaxClustersCount; i++)
+                        {                                       // copy sorted rows
+                            duplicate.row(sortedIdx.at<int>(i, 0)).copyTo(clusters.row(i));
+                        }
+                    }
+                }
+
+
+                /**
+                * @brief Fallback procedure invoked when the clustering eradicates all clusters.
+                *       Single cluster consisting of all points is then created.
+                * @param points List of sampled points.
+                * @param clusters Output list of clusters where the single cluster will be written.
+                */
+                void singleClusterFallback(const Mat& points, Mat& clusters) const
+                {
+                    if (clusters.rows != 0)
+                    {
+                        return;
+                    }
+                    clusters.create(1, points.cols, CV_32FC1);
+
+                    // Sum all points.
+                    reduce(points, clusters, 0, CV_REDUCE_SUM, CV_32FC1);
+
+                    // Sum all weights, all points have the same weight -> sum is the point count
+                    clusters.at<float>(0, WEIGHT_IDX) = static_cast<float>(points.rows);
+
+                    // Divide centroid by number of points -> compute average in each dimension.
+                    clusters /= clusters.at<float>(0, WEIGHT_IDX);
+                }
+
+
+                /**
+                * @brief Normalizes cluster weights, so they fall into range [0..1]
+                *       The other values are preserved.
+                * @param clusters List of clusters which will be normalized.
+                * @note Only weight of the clusters will be changed (the first column);
+                */
+                void normalizeWeights(Mat& clusters)
+                {
+                    // get max weight
+                    float maxWeight = clusters.at<float>(0, WEIGHT_IDX);
+                    for (int i = 1; i < clusters.rows; i++)
+                    {
+                        if (clusters.at<float>(i, WEIGHT_IDX) > maxWeight)
+                        {
+                            maxWeight = clusters.at<float>(i, WEIGHT_IDX);
+                        }
+                    }
+
+                    // normalize weight
+                    float weightNormalizer = 1 / maxWeight;
+                    for (int i = 0; i < clusters.rows; i++)
+                    {
+                        clusters.at<float>(i, WEIGHT_IDX) = clusters.at<float>(i, WEIGHT_IDX) * weightNormalizer;
+                    }
+                }
+
+
+                /**
+                * @brief Chooses which sample points will be used as cluster centers,
+                *       using list of random indexes that is stored in mInitSeedIndexes.
+                * @param samples List of sampled points.
+                * @param clusters Output list of random sampled points.
+                */
+                void pickRandomClusters(Mat& samples, Mat& clusters)
+                {
+                    clusters.create(0, SIGNATURE_DIMENSION, samples.type());
+                    clusters.reserve(mInitSeedIndexes.size());
+
+                    for (int iSeed = 0; iSeed < (int)(mInitSeedIndexes.size()); iSeed++)
+                    {
+                        clusters.push_back(samples.row(mInitSeedIndexes[iSeed]));
+                    }
+                }
+
+
+
+                /**
+                * @brief Indexes of initial clusters from sampling points.
+                */
+                std::vector<int> mInitSeedIndexes;
+
+                /**
+                * @brief Number of iterations performed (the number is fixed).
+                */
+                int mIterationCount;
+
+                /**
+                * @brief Maximal number of clusters. If the number of clusters is exceeded after
+                *       all iterations are completed, the smallest clusters are thrown away.
+                */
+                int mMaxClustersCount;
+
+                /**
+                * @brief Minimal size of the cluster. After each (i-th) iteration, the cluster sizes are checked.
+                *       Cluster smaller than i*mClusterMinSize are disposed of.
+                */
+                int mClusterMinSize;
+
+                /**
+                * @brief Distance threshold between two centroids. If two centroids become closer
+                *       than this distance, the clusters are joined.
+                */
+                float mJoiningDistance;
+
+                /**
+                * @brief Weight threshold of centroids. If weight of a centroid is below this value at the end of the iteration, it is thrown away.
+                */
+                float mDropThreshold;
+
+                /**
+                * @brief L_p metric is used for computing distances.
+                */
+                int mDistanceFunction;
+
+            };
+
+
+            Ptr<PCTClusterizer> PCTClusterizer::create(
+                const std::vector<int>& initSeedIndexes,
+                int iterationCount,
+                int maxClustersCount,
+                int clusterMinSize,
+                float joiningDistance,
+                float dropThreshold,
+                int distanceFunction)
+            {
+                return makePtr<PCTClusterizer_Impl>(
+                    initSeedIndexes,
+                    iterationCount,
+                    maxClustersCount,
+                    clusterMinSize,
+                    joiningDistance,
+                    dropThreshold,
+                    distanceFunction);
+            }
+        }
+    }
+}

modules/xfeatures2d/src/pct_signatures/pct_clusterizer.hpp

@@ -0,0 +1,124 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+
+#ifndef _OPENCV_XFEATURES_2D_PCT_SIGNATURES_CLUSTERIZER_HPP_
+#define _OPENCV_XFEATURES_2D_PCT_SIGNATURES_CLUSTERIZER_HPP_
+
+#ifdef __cplusplus
+
+
+#include "constants.hpp"
+#include "distance.hpp"
+
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            class PCTClusterizer : public Algorithm
+            {
+            public:
+
+                static Ptr<PCTClusterizer> create(
+                    const std::vector<int>& initSeedIndexes,
+                    int iterations = 10,
+                    int maxClusters = 768,  // max for Fermi GPU architecture
+                    int clusterMinSize = 2,
+                    float joiningDistance = 0.2,
+                    float dropThreshold = 0,
+                    int distanceFunction = PCTSignatures::L2);
+
+
+                /**** Accessors ****/
+
+                virtual int getIterationCount() const = 0;
+                virtual std::vector<int> getInitSeedIndexes() const = 0;
+                virtual int getMaxClustersCount() const = 0;
+                virtual int getClusterMinSize() const = 0;
+                virtual float getJoiningDistance() const = 0;
+                virtual float getDropThreshold() const = 0;
+                virtual int getDistanceFunction() const = 0;
+
+
+                virtual void setIterationCount(int iterationCount) = 0;
+                virtual void setInitSeedIndexes(std::vector<int> initSeedIndexes) = 0;
+                virtual void setMaxClustersCount(int maxClustersCount) = 0;
+                virtual void setClusterMinSize(int clusterMinSize) = 0;
+                virtual void setJoiningDistance(float joiningDistance) = 0;
+                virtual void setDropThreshold(float dropThreshold) = 0;
+                virtual void setDistanceFunction(int distanceFunction) = 0;
+
+                /**
+                * @brief K-means algorithm over the sampled points producing centroids as signatures.
+                * @param samples List of sampled points.
+                * @param signature Output list of computed centroids - the signature of the image.
+                */
+                virtual void clusterize(InputArray samples, OutputArray signature) = 0;
+            };
+        }
+    }
+}
+
+
+#endif
+
+#endif

modules/xfeatures2d/src/pct_signatures/pct_sampler.cpp

@@ -0,0 +1,248 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+#include "precomp.hpp"
+
+#include "pct_sampler.hpp"
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            class PCTSampler_Impl : public PCTSampler
+            {
+            private:
+                /**
+                * @brief Initial sampling point coordinates.
+                */
+                std::vector<Point2f> mInitSamplingPoints;
+
+                /**
+                * @brief Number of bits per pixel in grayscale image used for computing contrast and entropy.
+                */
+                int mGrayscaleBits;
+
+                /**
+                * @brief Radius of scanning window around a sampled point used for computing contrast and entropy.
+                */
+                int mWindowRadius;
+
+                /**
+                * @brief Weights of different feauture dimensions.
+                *       Default values are 1;
+                */
+                std::vector<float> mWeights;
+
+                /**
+                * @brief Translation of different feauture dimensions.
+                *       Default values are 0;
+                */
+                std::vector<float> mTranslations;
+
+
+            public:
+
+                PCTSampler_Impl(
+                    const std::vector<Point2f>& initSamplingPoints,
+                    int grayscaleBits,
+                    int windowRadius)
+                    : mInitSamplingPoints(initSamplingPoints),
+                    mGrayscaleBits(grayscaleBits),
+                    mWindowRadius(windowRadius)
+                {
+                    // Initialize weights and translation vectors to neutral values.
+                    for (int i = 0; i < SIGNATURE_DIMENSION; i++)
+                    {
+                        mWeights.push_back(1.0);
+                        mTranslations.push_back(0.0);
+                    }
+                }
+
+
+                /**** Acessors ****/
+
+                int getSampleCount() const      { return (int)mInitSamplingPoints.size(); }
+                int getGrayscaleBits() const    { return mGrayscaleBits; }
+                int getWindowRadius() const     { return mWindowRadius; }
+
+                float getWeightX() const               { return mWeights[X_IDX]; }
+                float getWeightY() const               { return mWeights[Y_IDX]; }
+                float getWeightL() const               { return mWeights[L_IDX]; }
+                float getWeightA() const               { return mWeights[A_IDX]; }
+                float getWeightB() const               { return mWeights[B_IDX]; }
+                float getWeightConstrast() const       { return mWeights[CONTRAST_IDX]; }
+                float getWeightEntropy() const         { return mWeights[ENTROPY_IDX]; }
+
+                std::vector<Point2f> getSamplingPoints() const
+                                                        { return mInitSamplingPoints; }
+
+
+                void setGrayscaleBits(int grayscaleBits)    { mGrayscaleBits = grayscaleBits; }
+                void setWindowRadius(int windowRadius)      { mWindowRadius = windowRadius; }
+
+                void setWeightX(float weight)          { mWeights[X_IDX] = weight; }
+                void setWeightY(float weight)          { mWeights[Y_IDX] = weight; }
+                void setWeightL(float weight)          { mWeights[L_IDX] = weight; }
+                void setWeightA(float weight)          { mWeights[A_IDX] = weight; }
+                void setWeightB(float weight)          { mWeights[B_IDX] = weight; }
+                void setWeightContrast(float weight)   { mWeights[CONTRAST_IDX] = weight; }
+                void setWeightEntropy(float weight)    { mWeights[ENTROPY_IDX] = weight; }
+
+
+                void setWeight(int idx, float value)
+                {
+                    mWeights[idx] = value;
+                }
+
+                void setWeights(const std::vector<float>& weights)
+                {
+                    if (weights.size() != mWeights.size())
+                    {
+                        CV_Error_(Error::StsUnmatchedSizes,
+                            ("Invalid weights dimension %d (max %d)", weights.size(), mWeights.size()));
+                    }
+                    else
+                    {
+                        for (int i = 0; i < (int)(mWeights.size()); ++i)
+                        {
+                            mWeights[i] = weights[i];
+                        }
+                    }
+                }
+
+                void setTranslation(int idx, float value)
+                {
+                    mTranslations[idx] = value;
+                }
+
+                void setTranslations(const std::vector<float>& translations)
+                {
+                    if (translations.size() != mTranslations.size())
+                    {
+                        CV_Error_(Error::StsUnmatchedSizes,
+                            ("Invalid translations dimension %d (max %d)", translations.size(), mTranslations.size()));
+                    }
+                    else
+                    {
+                        for (int i = 0; i < (int)(mTranslations.size()); ++i)
+                        {
+                            mTranslations[i] = translations[i];
+                        }
+                    }
+                }
+
+                void setSamplingPoints(std::vector<Point2f> samplingPoints) { mInitSamplingPoints = samplingPoints; }
+
+
+                void sample(InputArray _image, OutputArray _samples) const
+                {
+                    // prepare matrices
+                    Mat image = _image.getMat();
+                    _samples.create((int)(mInitSamplingPoints.size()), SIGNATURE_DIMENSION, CV_32F);
+                    Mat samples = _samples.getMat();
+                    GrayscaleBitmap grayscaleBitmap(image, mGrayscaleBits);
+
+                    // sample each sample point
+                    for (int iSample = 0; iSample < (int)(mInitSamplingPoints.size()); iSample++)
+                    {
+                        // sampling points are in range [0..1)
+                        int x = (int)(mInitSamplingPoints[iSample].x * (image.cols));
+                        int y = (int)(mInitSamplingPoints[iSample].y * (image.rows));
+
+                        // x, y normalized
+                        samples.at<float>(iSample, X_IDX) = (float)((float)x / (float)image.cols * mWeights[X_IDX] + mTranslations[X_IDX]);
+                        samples.at<float>(iSample, Y_IDX) = (float)((float)y / (float)image.rows * mWeights[Y_IDX] + mTranslations[Y_IDX]);
+
+                        // get Lab pixel color
+                        Mat rgbPixel(image, Rect(x, y, 1, 1));
+                        Mat labPixel;
+                        rgbPixel.convertTo(rgbPixel, CV_32FC3, 1.0 / 255);
+                        cvtColor(rgbPixel, labPixel, COLOR_BGR2Lab);
+                        Vec3f labColor = labPixel.at<Vec3f>(0, 0);
+
+                        // Lab color normalized
+                        samples.at<float>(iSample, L_IDX) = (float)(std::floor(labColor[0] + 0.5) / L_COLOR_RANGE * mWeights[L_IDX] + mTranslations[L_IDX]);
+                        samples.at<float>(iSample, A_IDX) = (float)(std::floor(labColor[1] + 0.5) / A_COLOR_RANGE * mWeights[A_IDX] + mTranslations[A_IDX]);
+                        samples.at<float>(iSample, B_IDX) = (float)(std::floor(labColor[2] + 0.5) / B_COLOR_RANGE * mWeights[B_IDX] + mTranslations[B_IDX]);
+
+                        // contrast and entropy
+                        float contrast = 0.0, entropy = 0.0;
+                        grayscaleBitmap.getContrastEntropy(x, y, contrast, entropy, mWindowRadius);     // HOT PATH: 30%
+                        samples.at<float>(iSample, CONTRAST_IDX)
+                            = (float)(contrast / SAMPLER_CONTRAST_NORMALIZER * mWeights[CONTRAST_IDX] + mTranslations[CONTRAST_IDX]);
+                        samples.at<float>(iSample, ENTROPY_IDX)
+                            = (float)(entropy / SAMPLER_ENTROPY_NORMALIZER * mWeights[ENTROPY_IDX] + mTranslations[ENTROPY_IDX]);
+                    }
+                }
+
+            };
+
+
+            Ptr<PCTSampler> PCTSampler::create(
+                const std::vector<Point2f>& initPoints,
+                int                         grayscaleBits,
+                int                         windowRadius)
+            {
+                return makePtr<PCTSampler_Impl>(initPoints, grayscaleBits, windowRadius);
+            }
+        }
+    }
+}

modules/xfeatures2d/src/pct_signatures/pct_sampler.hpp

@@ -0,0 +1,138 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+
+#ifndef _OPENCV_XFEATURES_2D_PCT_SIGNATURES_SAMPLER_HPP_
+#define _OPENCV_XFEATURES_2D_PCT_SIGNATURES_SAMPLER_HPP_
+
+#ifdef __cplusplus
+
+
+#include "constants.hpp"
+#include "grayscale_bitmap.hpp"
+
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            class PCTSampler : public Algorithm
+            {
+            public:
+
+                static Ptr<PCTSampler> create(
+                    const std::vector<Point2f>& initSamplingPoints,
+                    int                     grayscaleBits = 4,
+                    int                     windowRadius = 3);
+
+
+                /**
+                * @brief Sampling algorithm that produces samples of the input image
+                *       using the stored sampling point locations.
+                * @param image Input image.
+                * @param signature Output list of computed image samples.
+                */
+                virtual void sample(InputArray image, OutputArray samples) const = 0;
+
+
+                /**** accessors ****/
+
+                virtual int getSampleCount() const = 0;
+                virtual int getGrayscaleBits() const = 0;
+                virtual int getWindowRadius() const = 0;
+                virtual float getWeightX() const = 0;
+                virtual float getWeightY() const = 0;
+                virtual float getWeightL() const = 0;
+                virtual float getWeightA() const = 0;
+                virtual float getWeightB() const = 0;
+                virtual float getWeightConstrast() const = 0;
+                virtual float getWeightEntropy() const = 0;
+
+                virtual std::vector<Point2f> getSamplingPoints() const = 0;
+
+
+                virtual void setGrayscaleBits(int grayscaleBits) = 0;
+                virtual void setWindowRadius(int radius) = 0;
+                virtual void setWeightX(float weight) = 0;
+                virtual void setWeightY(float weight) = 0;
+                virtual void setWeightL(float weight) = 0;
+                virtual void setWeightA(float weight) = 0;
+                virtual void setWeightB(float weight) = 0;
+                virtual void setWeightContrast(float weight) = 0;
+                virtual void setWeightEntropy(float weight) = 0;
+
+                virtual void setWeight(int idx, float value) = 0;
+                virtual void setWeights(const std::vector<float>& weights) = 0;
+                virtual void setTranslation(int idx, float value) = 0;
+                virtual void setTranslations(const std::vector<float>& translations) = 0;
+
+                virtual void setSamplingPoints(std::vector<Point2f> samplingPoints) = 0;
+
+            };
+
+        }
+    }
+}
+
+
+#endif
+
+#endif

modules/xfeatures2d/src/pct_signatures/similarity.hpp

@@ -0,0 +1,133 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+
+#ifndef _OPENCV_XFEATURES_2D_PCT_SIGNATURES_SIMILARITY_HPP_
+#define _OPENCV_XFEATURES_2D_PCT_SIGNATURES_SIMILARITY_HPP_
+
+#ifdef __cplusplus
+
+#include "precomp.hpp"
+
+#include "distance.hpp"
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+
+            static inline float minusSimilarity(
+                const int distancefunction,
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                return -computeDistance(distancefunction, points1, idx1, points2, idx2);
+            }
+
+
+            static inline float gaussianSimilarity(
+                const int distancefunction,
+                const float alpha,
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                float distance = computeDistance(distancefunction, points1, idx1, points2, idx2);
+                return exp(-alpha + distance * distance);
+            }
+
+
+            static inline float heuristicSimilarity(
+                const int distancefunction,
+                const float alpha,
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                return 1 / (alpha + computeDistance(distancefunction, points1, idx1, points2, idx2));
+            }
+
+
+
+            static inline float computeSimilarity(
+                const int distancefunction,
+                const int similarity,
+                const float similarityParameter,
+                const Mat& points1, int idx1,
+                const Mat& points2, int idx2)
+            {
+                switch (similarity)
+                {
+                case PCTSignatures::MINUS:
+                    return minusSimilarity(distancefunction, points1, idx1, points2, idx2);
+                case PCTSignatures::GAUSSIAN:
+                    return gaussianSimilarity(distancefunction, similarityParameter, points1, idx1, points2, idx2);
+                case PCTSignatures::HEURISTIC:
+                    return heuristicSimilarity(distancefunction, similarityParameter, points1, idx1, points2, idx2);
+                default:
+                    CV_Error(Error::StsNotImplemented, "Similarity function not implemented!");
+                    return -1;
+                }
+            }
+        }
+    }
+}
+
+#endif
+
+#endif

modules/xfeatures2d/src/pct_signatures_sqfd.cpp

@@ -0,0 +1,224 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2000-2016, Intel Corporation, all rights reserved.
+Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+Copyright (C) 2009-2016, NVIDIA Corporation, all rights reserved.
+Copyright (C) 2010-2013, Advanced Micro Devices, Inc., all rights reserved.
+Copyright (C) 2015-2016, OpenCV Foundation, all rights reserved.
+Copyright (C) 2015-2016, 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:
+
+  * 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 names of the copyright holders nor the names of the 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 copyright holders 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.
+*/
+
+/*
+Contributed by Gregor Kovalcik <gregor dot kovalcik at gmail dot com>
+    based on code provided by Martin Krulis, Jakub Lokoc and Tomas Skopal.
+
+References:
+    Martin Krulis, Jakub Lokoc, Tomas Skopal.
+    Efficient Extraction of Clustering-Based Feature Signatures Using GPU Architectures.
+    Multimedia tools and applications, 75(13), pp.: 8071–8103, Springer, ISSN: 1380-7501, 2016
+
+    Christian Beecks, Merih Seran Uysal, Thomas Seidl.
+    Signature quadratic form distance.
+    In Proceedings of the ACM International Conference on Image and Video Retrieval, pages 438-445.
+    ACM, 2010.
+*/
+#include "precomp.hpp"
+
+#include "pct_signatures/constants.hpp"
+#include "pct_signatures/similarity.hpp"
+
+namespace cv
+{
+    namespace xfeatures2d
+    {
+        namespace pct_signatures
+        {
+            class PCTSignaturesSQFD_Impl : public PCTSignaturesSQFD
+            {
+            public:
+                PCTSignaturesSQFD_Impl(
+                    const int distanceFunction,
+                    const int similarityFunction,
+                    const float similarityParameter)
+                    : mDistanceFunction(distanceFunction),
+                    mSimilarityFunction(similarityFunction),
+                    mSimilarityParameter(similarityParameter)
+                {
+
+                }
+
+
+                float computeQuadraticFormDistance(
+                    InputArray _signature0,
+                    InputArray _signature1) const;
+
+                void computeQuadraticFormDistances(
+                    const Mat& sourceSignature,
+                    const std::vector<Mat>& imageSignatures,
+                    std::vector<float>& distances) const;
+
+
+            private:
+                int mDistanceFunction;
+                int mSimilarityFunction;
+                float mSimilarityParameter;
+
+                float computePartialSQFD(
+                    const Mat& signature0,
+                    const Mat& signature1) const;
+
+            };
+
+
+            /**
+            * @brief Class implementing parallel computing of SQFD distance for multiple images.
+            */
+            class Parallel_computeSQFDs : public ParallelLoopBody
+            {
+            private:
+                const PCTSignaturesSQFD* mPctSignaturesSQFDAlgorithm;
+                const Mat* mSourceSignature;
+                const std::vector<Mat>* mImageSignatures;
+                std::vector<float>* mDistances;
+
+            public:
+                Parallel_computeSQFDs(
+                    const PCTSignaturesSQFD* pctSignaturesSQFDAlgorithm,
+                    const Mat* sourceSignature,
+                    const std::vector<Mat>* imageSignatures,
+                    std::vector<float>* distances)
+                    : mPctSignaturesSQFDAlgorithm(pctSignaturesSQFDAlgorithm),
+                    mSourceSignature(sourceSignature),
+                    mImageSignatures(imageSignatures),
+                    mDistances(distances)
+                {
+                    mDistances->resize(imageSignatures->size());
+                }
+
+                void operator()(const Range& range) const
+                {
+                    if (mSourceSignature->empty())
+                    {
+                        CV_Error(Error::StsBadArg, "Source signature is empty!");
+                    }
+
+                    for (int i = range.start; i < range.end; i++)
+                    {
+                        if (mImageSignatures[i].empty())
+                        {
+                            CV_Error_(Error::StsBadArg, ("Signature ID: %d is empty!", i));
+                        }
+
+                        (*mDistances)[i] = mPctSignaturesSQFDAlgorithm->computeQuadraticFormDistance(
+                            *mSourceSignature, (*mImageSignatures)[i]);
+                    }
+                }
+            };
+
+
+            float PCTSignaturesSQFD_Impl::computeQuadraticFormDistance(
+                      InputArray _signature0,
+                      InputArray _signature1) const
+            {
+                // check input
+                if (_signature0.empty() || _signature1.empty())
+                {
+                    CV_Error(Error::StsBadArg, "Empty signature!");
+                }
+
+                Mat signature0 = _signature0.getMat();
+                Mat signature1 = _signature1.getMat();
+
+                if (signature0.cols != SIGNATURE_DIMENSION || signature1.cols != SIGNATURE_DIMENSION)
+                {
+                    CV_Error_(Error::StsBadArg, ("Signature dimension must be %d!", SIGNATURE_DIMENSION));
+                }
+
+                if (signature0.rows <= 0 || signature1.rows <= 0)
+                {
+                    CV_Error(Error::StsBadArg, "Signature count must be greater than 0!");
+                }
+
+                // compute sqfd
+                float result = 0;
+                result += computePartialSQFD(signature0, signature0);
+                result += computePartialSQFD(signature1, signature1);
+                result -= computePartialSQFD(signature0, signature1) * 2;
+
+                return sqrt(result);
+            }
+
+            void PCTSignaturesSQFD_Impl::computeQuadraticFormDistances(
+                      const Mat& sourceSignature,
+                      const std::vector<Mat>& imageSignatures,
+                      std::vector<float>& distances) const
+            {
+                parallel_for_(Range(0, (int)imageSignatures.size()),
+                    Parallel_computeSQFDs(this, &sourceSignature, &imageSignatures, &distances));
+            }
+
+            float PCTSignaturesSQFD_Impl::computePartialSQFD(
+                      const Mat& signature0,
+                      const Mat& signature1) const
+            {
+                float result = 0;
+                for (int i = 0; i < signature0.rows; i++)
+                {
+                    for (int j = 0; j < signature1.rows; j++)
+                    {
+                        result += signature0.at<float>(i, WEIGHT_IDX) * signature1.at<float>(j, WEIGHT_IDX)
+                            * computeSimilarity(mDistanceFunction, mSimilarityFunction, mSimilarityParameter, signature0, i, signature1, j);
+                    }
+                }
+                return result;
+            }
+
+
+
+        }// end of namespace pct_signatures
+
+
+        Ptr<PCTSignaturesSQFD> PCTSignaturesSQFD::create(
+            const int distanceFunction,
+            const int similarityFunction,
+            const float similarityParameter)
+        {
+            return makePtr<pct_signatures::PCTSignaturesSQFD_Impl>(distanceFunction, similarityFunction, similarityParameter);
+        }
+
+    }
+}