Merge pull request #1527 from vpisarev:shape_module

modules/shape/CMakeLists.txt

@@ -0,0 +1,2 @@
+set(the_description "Shape descriptors and matchers.")
+ocv_define_module(shape opencv_core opencv_imgproc opencv_video)

modules/shape/doc/emdL1.rst

@@ -0,0 +1,11 @@
+EMD-L1
+======
+Computes the "minimal work" distance between two weighted point configurations base on the papers "EMD-L1: An efficient and Robust Algorithm
+for comparing histogram-based descriptors", by Haibin Ling and Kazunori Okuda; and "The Earth Mover's Distance is the Mallows Distance:
+Some Insights from Statistics", by Elizaveta Levina and Peter Bickel.
+
+.. ocv:function:: float EMDL1( InputArray signature1, InputArray signature2 )
+
+    :param signature1: First signature, a single column floating-point matrix. Each row is the value of the histogram in each bin.
+
+    :param signature2: Second signature of the same format and size as  ``signature1``.

modules/shape/doc/histogram_cost_matrix.rst

@@ -0,0 +1,82 @@
+Cost Matrix for Histograms Common Interface
+===========================================
+
+.. highlight:: cpp
+
+A common interface is defined to ease the implementation of some algorithms pipelines, such
+as the Shape Context Matching Algorithm. A common class is defined, so any object that implements
+a Cost Matrix builder inherits the
+:ocv:class:`HistogramCostExtractor` interface.
+
+HistogramCostExtractor
+----------------------
+.. ocv:class:: HistogramCostExtractor : public Algorithm
+
+Abstract base class for histogram cost algorithms. ::
+
+    class CV_EXPORTS_W HistogramCostExtractor : public Algorithm
+    {
+    public:
+        CV_WRAP virtual void buildCostMatrix(InputArray descriptors1, InputArray descriptors2, OutputArray costMatrix) = 0;
+
+        CV_WRAP virtual void setNDummies(int nDummies) = 0;
+        CV_WRAP virtual int getNDummies() const = 0;
+
+        CV_WRAP virtual void setDefaultCost(float defaultCost) = 0;
+        CV_WRAP virtual float getDefaultCost() const = 0;
+    };
+
+NormHistogramCostExtractor
+--------------------------
+.. ocv:class:: NormHistogramCostExtractor : public HistogramCostExtractor
+
+A norm based cost extraction. ::
+
+    class CV_EXPORTS_W NormHistogramCostExtractor : public HistogramCostExtractor
+    {
+    public:
+        CV_WRAP virtual void setNormFlag(int flag) = 0;
+        CV_WRAP virtual int getNormFlag() const = 0;
+    };
+
+    CV_EXPORTS_W Ptr<HistogramCostExtractor>
+        createNormHistogramCostExtractor(int flag=cv::DIST_L2, int nDummies=25, float defaultCost=0.2);
+
+EMDHistogramCostExtractor
+-------------------------
+.. ocv:class:: EMDHistogramCostExtractor : public HistogramCostExtractor
+
+An EMD based cost extraction. ::
+
+    class CV_EXPORTS_W EMDHistogramCostExtractor : public HistogramCostExtractor
+    {
+    public:
+        CV_WRAP virtual void setNormFlag(int flag) = 0;
+        CV_WRAP virtual int getNormFlag() const = 0;
+    };
+
+    CV_EXPORTS_W Ptr<HistogramCostExtractor>
+        createEMDHistogramCostExtractor(int flag=cv::DIST_L2, int nDummies=25, float defaultCost=0.2);
+
+ChiHistogramCostExtractor
+-------------------------
+.. ocv:class:: ChiHistogramCostExtractor : public HistogramCostExtractor
+
+An Chi based cost extraction. ::
+
+    class CV_EXPORTS_W ChiHistogramCostExtractor : public HistogramCostExtractor
+    {};
+
+    CV_EXPORTS_W Ptr<HistogramCostExtractor> createChiHistogramCostExtractor(int nDummies=25, float defaultCost=0.2);
+
+EMDL1HistogramCostExtractor
+---------------------------
+.. ocv:class:: EMDL1HistogramCostExtractor : public HistogramCostExtractor
+
+An EMD-L1 based cost extraction. ::
+
+    class CV_EXPORTS_W EMDL1HistogramCostExtractor : public HistogramCostExtractor
+    {};
+
+    CV_EXPORTS_W Ptr<HistogramCostExtractor>
+        createEMDL1HistogramCostExtractor(int nDummies=25, float defaultCost=0.2);

modules/shape/doc/shape.rst

@@ -0,0 +1,15 @@
+**********************************
+shape. Shape Distance and Matching
+**********************************
+
+The module contains algorithms that embed a notion of shape distance.
+These algorithms may be used for shape matching and retrieval, or shape
+comparison.
+
+.. toctree::
+    :maxdepth: 2
+
+    shape_distances
+    shape_transformers
+    histogram_cost_matrix
+    emdL1

modules/shape/doc/shape_distances.rst

@@ -0,0 +1,231 @@
+Shape Distance and Common Interfaces
+====================================
+
+.. highlight:: cpp
+
+Shape Distance algorithms in OpenCV are derivated from a common interface that allows you to
+switch between them in a practical way for solving the same problem with different methods.
+Thus, all objects that implement shape distance measures inherit the
+:ocv:class:`ShapeDistanceExtractor` interface.
+
+
+ShapeDistanceExtractor
+----------------------
+.. ocv:class:: ShapeDistanceExtractor : public Algorithm
+
+Abstract base class for shape distance algorithms. ::
+
+    class CV_EXPORTS_W ShapeDistanceExtractor : public Algorithm
+    {
+    public:
+        CV_WRAP virtual float computeDistance(InputArray contour1, InputArray contour2) = 0;
+    };
+
+ShapeDistanceExtractor::computeDistance
+---------------------------------------
+Compute the shape distance between two shapes defined by its contours.
+
+.. ocv:function:: float ShapeDistanceExtractor::computeDistance( InputArray contour1, InputArray contour2 )
+
+    :param contour1: Contour defining first shape.
+
+    :param contour2: Contour defining second shape.
+
+ShapeContextDistanceExtractor
+-----------------------------
+.. ocv:class:: ShapeContextDistanceExtractor : public ShapeDistanceExtractor
+
+Implementation of the Shape Context descriptor and matching algorithm proposed by Belongie et al. in
+"Shape Matching and Object Recognition Using Shape Contexts" (PAMI 2002).
+This implementation is packaged in a generic scheme, in order to allow you the implementation of the
+common variations of the original pipeline. ::
+
+    class CV_EXPORTS_W ShapeContextDistanceExtractor : public ShapeDistanceExtractor
+    {
+    public:
+        CV_WRAP virtual void setAngularBins(int nAngularBins) = 0;
+        CV_WRAP virtual int getAngularBins() const = 0;
+
+        CV_WRAP virtual void setRadialBins(int nRadialBins) = 0;
+        CV_WRAP virtual int getRadialBins() const = 0;
+
+        CV_WRAP virtual void setInnerRadius(float innerRadius) = 0;
+        CV_WRAP virtual float getInnerRadius() const = 0;
+
+        CV_WRAP virtual void setOuterRadius(float outerRadius) = 0;
+        CV_WRAP virtual float getOuterRadius() const = 0;
+
+        CV_WRAP virtual void setRotationInvariant(bool rotationInvariant) = 0;
+        CV_WRAP virtual bool getRotationInvariant() const = 0;
+
+        CV_WRAP virtual void setShapeContextWeight(float shapeContextWeight) = 0;
+        CV_WRAP virtual float getShapeContextWeight() const = 0;
+
+        CV_WRAP virtual void setImageAppearanceWeight(float imageAppearanceWeight) = 0;
+        CV_WRAP virtual float getImageAppearanceWeight() const = 0;
+
+        CV_WRAP virtual void setBendingEnergyWeight(float bendingEnergyWeight) = 0;
+        CV_WRAP virtual float getBendingEnergyWeight() const = 0;
+
+        CV_WRAP virtual void setImages(InputArray image1, InputArray image2) = 0;
+        CV_WRAP virtual void getImages(OutputArray image1, OutputArray image2) const = 0;
+
+        CV_WRAP virtual void setIterations(int iterations) = 0;
+        CV_WRAP virtual int getIterations() const = 0;
+
+        CV_WRAP virtual void setCostExtractor(Ptr<HistogramCostExtractor> comparer) = 0;
+        CV_WRAP virtual Ptr<HistogramCostExtractor> getCostExtractor() const = 0;
+
+        CV_WRAP virtual void setTransformAlgorithm(Ptr<ShapeTransformer> transformer) = 0;
+        CV_WRAP virtual Ptr<ShapeTransformer> getTransformAlgorithm() const = 0;
+    };
+
+    /* Complete constructor */
+    CV_EXPORTS_W Ptr<ShapeContextDistanceExtractor>
+        createShapeContextDistanceExtractor(int nAngularBins=12, int nRadialBins=4,
+                                            float innerRadius=0.2, float outerRadius=2, int iterations=3,
+                                            const Ptr<HistogramCostExtractor> &comparer = createChiHistogramCostExtractor(),
+                                            const Ptr<ShapeTransformer> &transformer = createThinPlateSplineShapeTransformer());
+
+ShapeContextDistanceExtractor::setAngularBins
+---------------------------------------------
+Establish the number of angular bins for the Shape Context Descriptor used in the shape matching pipeline.
+
+.. ocv:function:: void setAngularBins( int nAngularBins )
+
+    :param nAngularBins: The number of angular bins in the shape context descriptor.
+
+ShapeContextDistanceExtractor::setRadialBins
+--------------------------------------------
+Establish the number of radial bins for the Shape Context Descriptor used in the shape matching pipeline.
+
+.. ocv:function:: void setRadialBins( int nRadialBins )
+
+    :param nRadialBins: The number of radial bins in the shape context descriptor.
+
+ShapeContextDistanceExtractor::setInnerRadius
+---------------------------------------------
+Set the inner radius of the shape context descriptor.
+
+.. ocv:function:: void setInnerRadius(float innerRadius)
+
+    :param innerRadius: The value of the inner radius.
+
+ShapeContextDistanceExtractor::setOuterRadius
+---------------------------------------------
+Set the outer radius of the shape context descriptor.
+
+.. ocv:function:: void setOuterRadius(float outerRadius)
+
+    :param outerRadius: The value of the outer radius.
+
+ShapeContextDistanceExtractor::setShapeContextWeight
+----------------------------------------------------
+Set the weight of the shape context distance in the final value of the shape distance.
+The shape context distance between two shapes is defined as the symmetric sum of shape
+context matching costs over best matching points.
+The final value of the shape distance is a user-defined linear combination of the shape
+context distance, an image appearance distance, and a bending energy.
+
+.. ocv:function:: void setShapeContextWeight( float shapeContextWeight )
+
+    :param shapeContextWeight: The weight of the shape context distance in the final distance value.
+
+ShapeContextDistanceExtractor::setImageAppearanceWeight
+-------------------------------------------------------
+Set the weight of the Image Appearance cost in the final value of the shape distance.
+The image appearance cost is defined as the sum of squared brightness differences in
+Gaussian windows around corresponding image points.
+The final value of the shape distance is a user-defined linear combination of the shape
+context distance, an image appearance distance, and a bending energy.
+If this value is set to a number different from 0, is mandatory to set the images that
+correspond to each shape.
+
+.. ocv:function:: void setImageAppearanceWeight( float imageAppearanceWeight )
+
+    :param imageAppearanceWeight: The weight of the appearance cost in the final distance value.
+
+ShapeContextDistanceExtractor::setBendingEnergyWeight
+-----------------------------------------------------
+Set the weight of the Bending Energy in the final value of the shape distance.
+The bending energy definition depends on what transformation is being used to align the
+shapes.
+The final value of the shape distance is a user-defined linear combination of the shape
+context distance, an image appearance distance, and a bending energy.
+
+.. ocv:function:: void setBendingEnergyWeight( float bendingEnergyWeight )
+
+    :param bendingEnergyWeight: The weight of the Bending Energy in the final distance value.
+
+ShapeContextDistanceExtractor::setImages
+----------------------------------------
+Set the images that correspond to each shape. This images are used in the calculation of the
+Image Appearance cost.
+
+.. ocv:function:: void setImages( InputArray image1, InputArray image2 )
+
+    :param image1: Image corresponding to the shape defined by ``contours1``.
+
+    :param image2: Image corresponding to the shape defined by ``contours2``.
+
+ShapeContextDistanceExtractor::setCostExtractor
+-----------------------------------------------
+Set the algorithm used for building the shape context descriptor cost matrix.
+
+.. ocv:function:: void setCostExtractor( Ptr<HistogramCostExtractor> comparer )
+
+    :param comparer: Smart pointer to a HistogramCostExtractor, an algorithm that defines the cost matrix between descriptors.
+
+ShapeContextDistanceExtractor::setStdDev
+----------------------------------------
+Set the value of the standard deviation for the Gaussian window for the image appearance cost.
+
+.. ocv:function:: void setStdDev( float sigma )
+
+    :param sigma: Standard Deviation.
+
+ShapeContextDistanceExtractor::setTransformAlgorithm
+----------------------------------------------------
+Set the algorithm used for aligning the shapes.
+
+.. ocv:function:: void setTransformAlgorithm( Ptr<ShapeTransformer> transformer )
+
+    :param comparer: Smart pointer to a ShapeTransformer, an algorithm that defines the aligning transformation.
+
+HausdorffDistanceExtractor
+--------------------------
+.. ocv:class:: HausdorffDistanceExtractor : public ShapeDistanceExtractor
+
+A simple Hausdorff distance measure between shapes defined by contours,
+according to the paper "Comparing Images using the Hausdorff distance." by
+D.P. Huttenlocher, G.A. Klanderman, and W.J. Rucklidge. (PAMI 1993). ::
+
+    class CV_EXPORTS_W HausdorffDistanceExtractor : public ShapeDistanceExtractor
+    {
+    public:
+        CV_WRAP virtual void setDistanceFlag(int distanceFlag) = 0;
+        CV_WRAP virtual int getDistanceFlag() const = 0;
+
+        CV_WRAP virtual void setRankProportion(float rankProportion) = 0;
+        CV_WRAP virtual float getRankProportion() const = 0;
+    };
+
+    /* Constructor */
+    CV_EXPORTS_W Ptr<HausdorffDistanceExtractor> createHausdorffDistanceExtractor(int distanceFlag=cv::NORM_L2, float rankProp=0.6);
+
+HausdorffDistanceExtractor::setDistanceFlag
+-------------------------------------------
+Set the norm used to compute the Hausdorff value between two shapes. It can be L1 or L2 norm.
+
+.. ocv:function:: void setDistanceFlag( int distanceFlag )
+
+    :param distanceFlag: Flag indicating which norm is used to compute the Hausdorff distance (NORM_L1, NORM_L2).
+
+HausdorffDistanceExtractor::setRankProportion
+---------------------------------------------
+This method sets the rank proportion (or fractional value) that establish the Kth ranked value of the
+partial Hausdorff distance. Experimentally had been shown that 0.6 is a good value to compare shapes.
+
+.. ocv:function:: void setRankProportion( float rankProportion )
+
+    :param rankProportion: fractional value (between 0 and 1).

modules/shape/doc/shape_transformers.rst

@@ -0,0 +1,108 @@
+Shape Transformers and Interfaces
+=================================
+
+.. highlight:: cpp
+
+A virtual interface that ease the use of transforming algorithms in some pipelines, such as
+the Shape Context Matching Algorithm. Thus, all objects that implement shape transformation
+techniques inherit the
+:ocv:class:`ShapeTransformer` interface.
+
+ShapeTransformer
+----------------
+.. ocv:class:: ShapeTransformer : public Algorithm
+
+Abstract base class for shape transformation algorithms. ::
+
+    class CV_EXPORTS_W ShapeTransformer : public Algorithm
+    {
+    public:
+        CV_WRAP virtual void estimateTransformation(InputArray transformingShape, InputArray targetShape,
+                                                     std::vector<DMatch>& matches) = 0;
+
+        CV_WRAP virtual float applyTransformation(InputArray input, OutputArray output=noArray()) = 0;
+
+        CV_WRAP virtual void warpImage(InputArray transformingImage, OutputArray output,
+                                       int flags=INTER_LINEAR, int borderMode=BORDER_CONSTANT,
+                                       const Scalar& borderValue=Scalar()) const = 0;
+    };
+
+ShapeTransformer::estimateTransformation
+----------------------------------------
+Estimate the transformation parameters of the current transformer algorithm, based on point matches.
+
+.. ocv:function:: void estimateTransformation( InputArray transformingShape, InputArray targetShape, std::vector<DMatch>& matches )
+
+    :param transformingShape: Contour defining first shape.
+
+    :param targetShape: Contour defining second shape (Target).
+
+    :param matches: Standard vector of Matches between points.
+
+ShapeTransformer::applyTransformation
+-------------------------------------
+Apply a transformation, given a pre-estimated transformation parameters.
+
+.. ocv:function:: float applyTransformation( InputArray input, OutputArray output=noArray() )
+
+    :param input: Contour (set of points) to apply the transformation.
+
+    :param output: Output contour.
+
+ShapeTransformer::warpImage
+---------------------------
+Apply a transformation, given a pre-estimated transformation parameters, to an Image.
+
+.. ocv:function:: void warpImage( InputArray transformingImage, OutputArray output, int flags=INTER_LINEAR, int borderMode=BORDER_CONSTANT, const Scalar& borderValue=Scalar() )
+
+    :param transformingImage: Input image.
+
+    :param output: Output image.
+
+    :param flags: Image interpolation method.
+
+    :param borderMode: border style.
+
+    :param borderValue: border value.
+
+ThinPlateSplineShapeTransformer
+-------------------------------
+.. ocv:class:: ThinPlateSplineShapeTransformer : public Algorithm
+
+Definition of the transformation ocupied in the paper "Principal Warps: Thin-Plate Splines and Decomposition
+of Deformations", by F.L. Bookstein (PAMI 1989). ::
+
+    class CV_EXPORTS_W ThinPlateSplineShapeTransformer : public ShapeTransformer
+    {
+    public:
+        CV_WRAP virtual void setRegularizationParameter(double beta) = 0;
+        CV_WRAP virtual double getRegularizationParameter() const = 0;
+    };
+
+    /* Complete constructor */
+    CV_EXPORTS_W Ptr<ThinPlateSplineShapeTransformer>
+        createThinPlateSplineShapeTransformer(double regularizationParameter=0);
+
+ThinPlateSplineShapeTransformer::setRegularizationParameter
+-----------------------------------------------------------
+Set the regularization parameter for relaxing the exact interpolation requirements of the TPS algorithm.
+
+.. ocv:function:: void setRegularizationParameter( double beta )
+
+    :param beta: value of the regularization parameter.
+
+AffineTransformer
+-----------------
+.. ocv:class:: AffineTransformer : public Algorithm
+
+Wrapper class for the OpenCV Affine Transformation algorithm. ::
+
+    class CV_EXPORTS_W AffineTransformer : public ShapeTransformer
+    {
+    public:
+        CV_WRAP virtual void setFullAffine(bool fullAffine) = 0;
+        CV_WRAP virtual bool getFullAffine() const = 0;
+    };
+
+    /* Complete constructor */
+    CV_EXPORTS_W Ptr<AffineTransformer> createAffineTransformer(bool fullAffine);

modules/shape/include/opencv2/shape.hpp

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

modules/shape/include/opencv2/shape/emdL1.hpp

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

modules/shape/include/opencv2/shape/hist_cost.hpp

@@ -0,0 +1,103 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                          License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_HIST_COST_HPP__
+#define __OPENCV_HIST_COST_HPP__
+
+#include "opencv2/imgproc.hpp"
+
+namespace cv
+{
+
+/*!
+ * The base class for HistogramCostExtractor.
+ */
+class CV_EXPORTS_W HistogramCostExtractor : public Algorithm
+{
+public:
+    CV_WRAP virtual void buildCostMatrix(InputArray descriptors1, InputArray descriptors2, OutputArray costMatrix) = 0;
+
+    CV_WRAP virtual void setNDummies(int nDummies) = 0;
+    CV_WRAP virtual int getNDummies() const = 0;
+
+    CV_WRAP virtual void setDefaultCost(float defaultCost) = 0;
+    CV_WRAP virtual float getDefaultCost() const = 0;
+};
+
+/*!  */
+class CV_EXPORTS_W NormHistogramCostExtractor : public HistogramCostExtractor
+{
+public:
+    CV_WRAP virtual void setNormFlag(int flag) = 0;
+    CV_WRAP virtual int getNormFlag() const = 0;
+};
+
+CV_EXPORTS_W Ptr<HistogramCostExtractor>
+    createNormHistogramCostExtractor(int flag=DIST_L2, int nDummies=25, float defaultCost=0.2);
+
+/*!  */
+class CV_EXPORTS_W EMDHistogramCostExtractor : public HistogramCostExtractor
+{
+public:
+    CV_WRAP virtual void setNormFlag(int flag) = 0;
+    CV_WRAP virtual int getNormFlag() const = 0;
+};
+
+CV_EXPORTS_W Ptr<HistogramCostExtractor>
+    createEMDHistogramCostExtractor(int flag=DIST_L2, int nDummies=25, float defaultCost=0.2);
+
+/*!  */
+class CV_EXPORTS_W ChiHistogramCostExtractor : public HistogramCostExtractor
+{};
+
+CV_EXPORTS_W Ptr<HistogramCostExtractor> createChiHistogramCostExtractor(int nDummies=25, float defaultCost=0.2);
+
+/*!  */
+class CV_EXPORTS_W EMDL1HistogramCostExtractor : public HistogramCostExtractor
+{};
+
+CV_EXPORTS_W Ptr<HistogramCostExtractor>
+    createEMDL1HistogramCostExtractor(int nDummies=25, float defaultCost=0.2);
+
+} // cv
+#endif

modules/shape/include/opencv2/shape/shape.hpp

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

modules/shape/include/opencv2/shape/shape_distance.hpp

@@ -0,0 +1,143 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                          License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_SHAPE_SHAPE_DISTANCE_HPP__
+#define __OPENCV_SHAPE_SHAPE_DISTANCE_HPP__
+#include "opencv2/core.hpp"
+#include "opencv2/shape/hist_cost.hpp"
+#include "opencv2/shape/shape_transformer.hpp"
+
+namespace cv
+{
+
+/*!
+ * The base class for ShapeDistanceExtractor.
+ * This is just to define the common interface for
+ * shape comparisson techniques.
+ */
+class CV_EXPORTS_W ShapeDistanceExtractor : public Algorithm
+{
+public:
+    CV_WRAP virtual float computeDistance(InputArray contour1, InputArray contour2) = 0;
+};
+
+/***********************************************************************************/
+/***********************************************************************************/
+/***********************************************************************************/
+/*!
+ * Shape Context implementation.
+ * The SCD class implements SCD algorithm proposed by Belongie et al.in
+ * "Shape Matching and Object Recognition Using Shape Contexts".
+ * Implemented by Juan M. Perez for the GSOC 2013.
+ */
+class CV_EXPORTS_W ShapeContextDistanceExtractor : public ShapeDistanceExtractor
+{
+public:
+    CV_WRAP virtual void setAngularBins(int nAngularBins) = 0;
+    CV_WRAP virtual int getAngularBins() const = 0;
+
+    CV_WRAP virtual void setRadialBins(int nRadialBins) = 0;
+    CV_WRAP virtual int getRadialBins() const = 0;
+
+    CV_WRAP virtual void setInnerRadius(float innerRadius) = 0;
+    CV_WRAP virtual float getInnerRadius() const = 0;
+
+    CV_WRAP virtual void setOuterRadius(float outerRadius) = 0;
+    CV_WRAP virtual float getOuterRadius() const = 0;
+
+    CV_WRAP virtual void setRotationInvariant(bool rotationInvariant) = 0;
+    CV_WRAP virtual bool getRotationInvariant() const = 0;
+
+    CV_WRAP virtual void setShapeContextWeight(float shapeContextWeight) = 0;
+    CV_WRAP virtual float getShapeContextWeight() const = 0;
+
+    CV_WRAP virtual void setImageAppearanceWeight(float imageAppearanceWeight) = 0;
+    CV_WRAP virtual float getImageAppearanceWeight() const = 0;
+
+    CV_WRAP virtual void setBendingEnergyWeight(float bendingEnergyWeight) = 0;
+    CV_WRAP virtual float getBendingEnergyWeight() const = 0;
+
+    CV_WRAP virtual void setImages(InputArray image1, InputArray image2) = 0;
+    CV_WRAP virtual void getImages(OutputArray image1, OutputArray image2) const = 0;
+
+    CV_WRAP virtual void setIterations(int iterations) = 0;
+    CV_WRAP virtual int getIterations() const = 0;
+
+    CV_WRAP virtual void setCostExtractor(Ptr<HistogramCostExtractor> comparer) = 0;
+    CV_WRAP virtual Ptr<HistogramCostExtractor> getCostExtractor() const = 0;
+
+    CV_WRAP virtual void setStdDev(float sigma) = 0;
+    CV_WRAP virtual float getStdDev() const = 0;
+
+    CV_WRAP virtual void setTransformAlgorithm(Ptr<ShapeTransformer> transformer) = 0;
+    CV_WRAP virtual Ptr<ShapeTransformer> getTransformAlgorithm() const = 0;
+};
+
+/* Complete constructor */
+CV_EXPORTS_W Ptr<ShapeContextDistanceExtractor>
+    createShapeContextDistanceExtractor(int nAngularBins=12, int nRadialBins=4,
+                                        float innerRadius=0.2, float outerRadius=2, int iterations=3,
+                                        const Ptr<HistogramCostExtractor> &comparer = createChiHistogramCostExtractor(),
+                                        const Ptr<ShapeTransformer> &transformer = createThinPlateSplineShapeTransformer());
+
+/***********************************************************************************/
+/***********************************************************************************/
+/***********************************************************************************/
+/*!
+ * Hausdorff distace implementation based on
+ */
+class CV_EXPORTS_W HausdorffDistanceExtractor : public ShapeDistanceExtractor
+{
+public:
+    CV_WRAP virtual void setDistanceFlag(int distanceFlag) = 0;
+    CV_WRAP virtual int getDistanceFlag() const = 0;
+
+    CV_WRAP virtual void setRankProportion(float rankProportion) = 0;
+    CV_WRAP virtual float getRankProportion() const = 0;
+};
+
+/* Constructor */
+CV_EXPORTS_W Ptr<HausdorffDistanceExtractor> createHausdorffDistanceExtractor(int distanceFlag=cv::NORM_L2, float rankProp=0.6);
+
+} // cv
+#endif

modules/shape/include/opencv2/shape/shape_transformer.hpp

@@ -0,0 +1,109 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                          License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_SHAPE_SHAPE_TRANSFORM_HPP__
+#define __OPENCV_SHAPE_SHAPE_TRANSFORM_HPP__
+#include <vector>
+#include "opencv2/core.hpp"
+#include "opencv2/imgproc.hpp"
+
+namespace cv
+{
+
+/*!
+ * The base class for ShapeTransformer.
+ * This is just to define the common interface for
+ * shape transformation techniques.
+ */
+class CV_EXPORTS_W ShapeTransformer : public Algorithm
+{
+public:
+    /* Estimate, Apply Transformation and return Transforming cost*/
+    CV_WRAP virtual void estimateTransformation(InputArray transformingShape, InputArray targetShape,
+                                                 std::vector<DMatch>& matches) = 0;
+
+    CV_WRAP virtual float applyTransformation(InputArray input, OutputArray output=noArray()) = 0;
+
+    CV_WRAP virtual void warpImage(InputArray transformingImage, OutputArray output,
+                                   int flags=INTER_LINEAR, int borderMode=BORDER_CONSTANT,
+                                   const Scalar& borderValue=Scalar()) const = 0;
+};
+
+/***********************************************************************************/
+/***********************************************************************************/
+/*!
+ * Thin Plate Spline Transformation
+ * Implementation of the TPS transformation
+ * according to "Principal Warps: Thin-Plate Splines and the
+ * Decomposition of Deformations" by Juan Manuel Perez for the GSOC 2013
+ */
+
+class CV_EXPORTS_W ThinPlateSplineShapeTransformer : public ShapeTransformer
+{
+public:
+    CV_WRAP virtual void setRegularizationParameter(double beta) = 0;
+    CV_WRAP virtual double getRegularizationParameter() const = 0;
+};
+
+/* Complete constructor */
+CV_EXPORTS_W Ptr<ThinPlateSplineShapeTransformer>
+    createThinPlateSplineShapeTransformer(double regularizationParameter=0);
+
+/***********************************************************************************/
+/***********************************************************************************/
+/*!
+ * Affine Transformation as a derivated from ShapeTransformer
+ */
+
+class CV_EXPORTS_W AffineTransformer : public ShapeTransformer
+{
+public:
+    CV_WRAP virtual void setFullAffine(bool fullAffine) = 0;
+    CV_WRAP virtual bool getFullAffine() const = 0;
+};
+
+/* Complete constructor */
+CV_EXPORTS_W Ptr<AffineTransformer> createAffineTransformer(bool fullAffine);
+
+} // cv
+#endif

modules/shape/src/aff_trans.cpp

@@ -0,0 +1,266 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+class AffineTransformerImpl : public AffineTransformer
+{
+public:
+    /* Constructors */
+    AffineTransformerImpl()
+    {
+        fullAffine = true;
+        name_ = "ShapeTransformer.AFF";
+    }
+
+    AffineTransformerImpl(bool _fullAffine)
+    {
+        fullAffine = _fullAffine;
+        name_ = "ShapeTransformer.AFF";
+    }
+
+    /* Destructor */
+    ~AffineTransformerImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operator
+    virtual void estimateTransformation(InputArray transformingShape, InputArray targetShape, std::vector<DMatch> &matches);
+    virtual float applyTransformation(InputArray input, OutputArray output=noArray());
+    virtual void warpImage(InputArray transformingImage, OutputArray output,
+                           int flags, int borderMode, const Scalar& borderValue) const;
+
+    //! Setters/Getters
+    virtual void setFullAffine(bool _fullAffine) {fullAffine=_fullAffine;}
+    virtual bool getFullAffine() const {return fullAffine;}
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "affine_type" << int(fullAffine);
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        fullAffine = int(fn["affine_type"])?true:false;
+    }
+
+private:
+    bool fullAffine;
+    Mat affineMat;
+    float transformCost;
+
+protected:
+    String name_;
+};
+
+void AffineTransformerImpl::warpImage(InputArray transformingImage, OutputArray output,
+                                      int flags, int borderMode, const Scalar& borderValue) const
+{
+    CV_Assert(!affineMat.empty());
+    warpAffine(transformingImage, output, affineMat, transformingImage.getMat().size(), flags, borderMode, borderValue);
+}
+
+
+static Mat _localAffineEstimate(const std::vector<Point2f>& shape1, const std::vector<Point2f>& shape2,
+                                bool fullAfine)
+{
+    Mat out(2,3,CV_32F);
+    int siz=2*shape1.size();
+
+    if (fullAfine)
+    {
+        Mat matM(siz, 6, CV_32F);
+        Mat matP(siz,1,CV_32F);
+        int contPt=0;
+        for (int ii=0; ii<siz; ii++)
+        {
+            Mat therow = Mat::zeros(1,6,CV_32F);
+            if (ii%2==0)
+            {
+                therow.at<float>(0,0)=shape1[contPt].x;
+                therow.at<float>(0,1)=shape1[contPt].y;
+                therow.at<float>(0,2)=1;
+                therow.row(0).copyTo(matM.row(ii));
+                matP.at<float>(ii,0) = shape2[contPt].x;
+            }
+            else
+            {
+                therow.at<float>(0,3)=shape1[contPt].x;
+                therow.at<float>(0,4)=shape1[contPt].y;
+                therow.at<float>(0,5)=1;
+                therow.row(0).copyTo(matM.row(ii));
+                matP.at<float>(ii,0) = shape2[contPt].y;
+                contPt++;
+            }
+        }
+        Mat sol;
+        solve(matM, matP, sol, DECOMP_SVD);
+        out = sol.reshape(0,2);
+    }
+    else
+    {
+        Mat matM(siz, 4, CV_32F);
+        Mat matP(siz,1,CV_32F);
+        int contPt=0;
+        for (int ii=0; ii<siz; ii++)
+        {
+            Mat therow = Mat::zeros(1,4,CV_32F);
+            if (ii%2==0)
+            {
+                therow.at<float>(0,0)=shape1[contPt].x;
+                therow.at<float>(0,1)=shape1[contPt].y;
+                therow.at<float>(0,2)=1;
+                therow.row(0).copyTo(matM.row(ii));
+                matP.at<float>(ii,0) = shape2[contPt].x;
+            }
+            else
+            {
+                therow.at<float>(0,0)=-shape1[contPt].y;
+                therow.at<float>(0,1)=shape1[contPt].x;
+                therow.at<float>(0,3)=1;
+                therow.row(0).copyTo(matM.row(ii));
+                matP.at<float>(ii,0) = shape2[contPt].y;
+                contPt++;
+            }
+        }
+        Mat sol;
+        solve(matM, matP, sol, DECOMP_SVD);
+        out.at<float>(0,0)=sol.at<float>(0,0);
+        out.at<float>(0,1)=sol.at<float>(1,0);
+        out.at<float>(0,2)=sol.at<float>(2,0);
+        out.at<float>(1,0)=-sol.at<float>(1,0);
+        out.at<float>(1,1)=sol.at<float>(0,0);
+        out.at<float>(1,2)=sol.at<float>(3,0);
+    }
+    return out;
+}
+
+void AffineTransformerImpl::estimateTransformation(InputArray _pts1, InputArray _pts2, std::vector<DMatch>& _matches)
+{
+    Mat pts1 = _pts1.getMat();
+    Mat pts2 = _pts2.getMat();
+    CV_Assert((pts1.channels()==2) && (pts1.cols>0) && (pts2.channels()==2) && (pts2.cols>0));
+    CV_Assert(_matches.size()>1);
+
+    if (pts1.type() != CV_32F)
+        pts1.convertTo(pts1, CV_32F);
+    if (pts2.type() != CV_32F)
+        pts2.convertTo(pts2, CV_32F);
+
+    // Use only valid matchings //
+    std::vector<DMatch> matches;
+    for (size_t i=0; i<_matches.size(); i++)
+    {
+        if (_matches[i].queryIdx<pts1.cols &&
+            _matches[i].trainIdx<pts2.cols)
+        {
+            matches.push_back(_matches[i]);
+        }
+    }
+
+    // Organizing the correspondent points in vector style //
+    std::vector<Point2f> shape1; // transforming shape
+    std::vector<Point2f> shape2; // target shape
+    for (size_t i=0; i<matches.size(); i++)
+    {
+        Point2f pt1=pts1.at<Point2f>(0,matches[i].queryIdx);
+        shape1.push_back(pt1);
+
+        Point2f pt2=pts2.at<Point2f>(0,matches[i].trainIdx);
+        shape2.push_back(pt2);
+    }
+
+    // estimateRigidTransform //
+    Mat affine;
+    estimateRigidTransform(shape1, shape2, fullAffine).convertTo(affine, CV_32F);
+
+    if (affine.empty())
+        affine=_localAffineEstimate(shape1, shape2, fullAffine); //In case there is not good solution, just give a LLS based one
+
+    affineMat = affine;
+}
+
+float AffineTransformerImpl::applyTransformation(InputArray inPts, OutputArray outPts)
+{
+    Mat pts1 = inPts.getMat();
+    CV_Assert((pts1.channels()==2) && (pts1.cols>0));
+
+    //Apply transformation in the complete set of points
+    Mat fAffine;
+    transform(pts1, fAffine, affineMat);
+
+    // Ensambling output //
+    if (outPts.needed())
+    {
+        outPts.create(1,fAffine.cols, CV_32FC2);
+        Mat outMat = outPts.getMat();
+        for (int i=0; i<fAffine.cols; i++)
+            outMat.at<Point2f>(0,i)=fAffine.at<Point2f>(0,i);
+    }
+
+    // Updating Transform Cost //
+    Mat Af(2, 2, CV_32F);
+    Af.at<float>(0,0)=affineMat.at<float>(0,0);
+    Af.at<float>(0,1)=affineMat.at<float>(1,0);
+    Af.at<float>(1,0)=affineMat.at<float>(0,1);
+    Af.at<float>(1,1)=affineMat.at<float>(1,1);
+    SVD mysvd(Af, SVD::NO_UV);
+    Mat singVals=mysvd.w;
+    transformCost=std::log((singVals.at<float>(0,0)+FLT_MIN)/(singVals.at<float>(1,0)+FLT_MIN));
+
+    return transformCost;
+}
+
+Ptr <AffineTransformer> createAffineTransformer(bool fullAffine)
+{
+    return Ptr<AffineTransformer>( new AffineTransformerImpl(fullAffine) );
+}
+
+} // cv

modules/shape/src/emdL1.cpp

@@ -0,0 +1,793 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+/*
+ * Implementation of an optimized EMD for histograms based in
+ * the papers "EMD-L1: An efficient and Robust Algorithm
+ * for comparing histogram-based descriptors", by Haibin Ling and
+ * Kazunori Okuda; and "The Earth Mover's Distance is the Mallows
+ * Distance: Some Insights from Statistics", by Elizaveta Levina and
+ * Peter Bickel, based on HAIBIN LING AND KAZUNORI OKADA implementation.
+ */
+
+#include "precomp.hpp"
+#include "emdL1_def.hpp"
+#include <limits>
+
+/****************************************************************************************\
+*                                   EMDL1 Class                                         *
+\****************************************************************************************/
+
+float EmdL1::getEMDL1(cv::Mat &sig1, cv::Mat &sig2)
+{
+    // Initialization
+    CV_Assert((sig1.rows==sig2.rows) && (sig1.cols==sig2.cols) && (!sig1.empty()) && (!sig2.empty()));
+    if(!initBaseTrees(sig1.rows, 1))
+        return -1;
+
+    float *H1=new float[sig1.rows], *H2 = new float[sig2.rows];
+    for (int ii=0; ii<sig1.rows; ii++)
+    {
+        H1[ii]=sig1.at<float>(ii,0);
+        H2[ii]=sig2.at<float>(ii,0);
+    }
+
+    fillBaseTrees(H1,H2); // Initialize histograms
+    greedySolution(); // Construct an initial Basic Feasible solution
+    initBVTree(); // Initialize BVTree
+
+    // Iteration
+    bool bOptimal = false;
+    m_nItr = 0;
+    while(!bOptimal && m_nItr<nMaxIt)
+    {
+        // Derive U=(u_ij) for row i and column j
+        if(m_nItr==0) updateSubtree(m_pRoot);
+        else updateSubtree(m_pEnter->pChild);
+
+        // Optimality test
+        bOptimal = isOptimal();
+
+        // Find new solution
+        if(!bOptimal)
+            findNewSolution();
+        ++m_nItr;
+    }
+    delete [] H1;
+    delete [] H2;
+    // Output the total flow
+    return compuTotalFlow();
+}
+
+void EmdL1::setMaxIteration(int _nMaxIt)
+{
+    nMaxIt=_nMaxIt;
+}
+
+//-- SubFunctions called in the EMD algorithm
+bool EmdL1::initBaseTrees(int n1, int n2, int n3)
+{
+    if(binsDim1==n1 && binsDim2==n2 && binsDim3==n3)
+        return true;
+    binsDim1 = n1;
+    binsDim2 = n2;
+    binsDim3 = n3;
+    if(binsDim1==0 || binsDim2==0) dimension = 0;
+    else dimension	= (binsDim3==0)?2:3;
+
+    if(dimension==2)
+    {
+        m_Nodes.resize(binsDim1);
+        m_EdgesUp.resize(binsDim1);
+        m_EdgesRight.resize(binsDim1);
+        for(int i1=0; i1<binsDim1; i1++)
+        {
+            m_Nodes[i1].resize(binsDim2);
+            m_EdgesUp[i1].resize(binsDim2);
+            m_EdgesRight[i1].resize(binsDim2);
+        }
+        m_NBVEdges.resize(binsDim1*binsDim2*4+2);
+        m_auxQueue.resize(binsDim1*binsDim2+2);
+        m_fromLoop.resize(binsDim1*binsDim2+2);
+        m_toLoop.resize(binsDim1*binsDim2+2);
+    }
+    else if(dimension==3)
+    {
+        m_3dNodes.resize(binsDim1);
+        m_3dEdgesUp.resize(binsDim1);
+        m_3dEdgesRight.resize(binsDim1);
+        m_3dEdgesDeep.resize(binsDim1);
+        for(int i1=0; i1<binsDim1; i1++)
+        {
+            m_3dNodes[i1].resize(binsDim2);
+            m_3dEdgesUp[i1].resize(binsDim2);
+            m_3dEdgesRight[i1].resize(binsDim2);
+            m_3dEdgesDeep[i1].resize(binsDim2);
+            for(int i2=0; i2<binsDim2; i2++)
+            {
+                m_3dNodes[i1][i2].resize(binsDim3);
+                m_3dEdgesUp[i1][i2].resize(binsDim3);
+                m_3dEdgesRight[i1][i2].resize(binsDim3);
+                m_3dEdgesDeep[i1][i2].resize(binsDim3);
+            }
+        }
+        m_NBVEdges.resize(binsDim1*binsDim2*binsDim3*6+4);
+        m_auxQueue.resize(binsDim1*binsDim2*binsDim3+4);
+        m_fromLoop.resize(binsDim1*binsDim2*binsDim3+4);
+        m_toLoop.resize(binsDim1*binsDim2*binsDim3+2);
+    }
+    else
+        return false;
+
+    return true;
+}
+
+bool EmdL1::fillBaseTrees(float *H1, float *H2)
+{
+    //- Set global counters
+    m_pRoot	= NULL;
+    // Graph initialization
+    float *p1 = H1;
+    float *p2 = H2;
+    if(dimension==2)
+    {
+        for(int c=0; c<binsDim2; c++)
+        {
+            for(int r=0; r<binsDim1; r++)
+            {
+                //- initialize nodes and links
+                m_Nodes[r][c].pos[0] = r;
+                m_Nodes[r][c].pos[1] = c;
+                m_Nodes[r][c].d = *(p1++)-*(p2++);
+                m_Nodes[r][c].pParent = NULL;
+                m_Nodes[r][c].pChild = NULL;
+                m_Nodes[r][c].iLevel = -1;
+
+                //- initialize edges
+                // to the right
+                m_EdgesRight[r][c].pParent = &(m_Nodes[r][c]);
+                m_EdgesRight[r][c].pChild = &(m_Nodes[r][(c+1)%binsDim2]);
+                m_EdgesRight[r][c].flow	= 0;
+                m_EdgesRight[r][c].iDir	= 1;
+                m_EdgesRight[r][c].pNxt	= NULL;
+
+                // to the upward
+                m_EdgesUp[r][c].pParent	= &(m_Nodes[r][c]);
+                m_EdgesUp[r][c].pChild	= &(m_Nodes[(r+1)%binsDim1][c]);
+                m_EdgesUp[r][c].flow = 0;
+                m_EdgesUp[r][c].iDir = 1;
+                m_EdgesUp[r][c].pNxt = NULL;
+            }
+        }
+    }
+    else if(dimension==3)
+    {
+        for(int z=0; z<binsDim3; z++)
+        {
+            for(int c=0; c<binsDim2; c++)
+            {
+                for(int r=0; r<binsDim1; r++)
+                {
+                    //- initialize nodes and edges
+                    m_3dNodes[r][c][z].pos[0] = r;
+                    m_3dNodes[r][c][z].pos[1] = c;
+                    m_3dNodes[r][c][z].pos[2] = z;
+                    m_3dNodes[r][c][z].d = *(p1++)-*(p2++);
+                    m_3dNodes[r][c][z].pParent = NULL;
+                    m_3dNodes[r][c][z].pChild = NULL;
+                    m_3dNodes[r][c][z].iLevel = -1;
+
+                    //- initialize edges
+                    // to the upward
+                    m_3dEdgesUp[r][c][z].pParent= &(m_3dNodes[r][c][z]);
+                    m_3dEdgesUp[r][c][z].pChild	= &(m_3dNodes[(r+1)%binsDim1][c][z]);
+                    m_3dEdgesUp[r][c][z].flow = 0;
+                    m_3dEdgesUp[r][c][z].iDir = 1;
+                    m_3dEdgesUp[r][c][z].pNxt = NULL;
+
+                    // to the right
+                    m_3dEdgesRight[r][c][z].pParent	= &(m_3dNodes[r][c][z]);
+                    m_3dEdgesRight[r][c][z].pChild	= &(m_3dNodes[r][(c+1)%binsDim2][z]);
+                    m_3dEdgesRight[r][c][z].flow	= 0;
+                    m_3dEdgesRight[r][c][z].iDir	= 1;
+                    m_3dEdgesRight[r][c][z].pNxt	= NULL;
+
+                    // to the deep
+                    m_3dEdgesDeep[r][c][z].pParent	= &(m_3dNodes[r][c][z]);
+                    m_3dEdgesDeep[r][c][z].pChild	= &(m_3dNodes[r][c])[(z+1)%binsDim3];
+                    m_3dEdgesDeep[r][c][z].flow = 0;
+                    m_3dEdgesDeep[r][c][z].iDir = 1;
+                    m_3dEdgesDeep[r][c][z].pNxt = NULL;
+                }
+            }
+        }
+    }
+    return true;
+}
+
+bool EmdL1::greedySolution()
+{
+    return dimension==2?greedySolution2():greedySolution3();
+}
+
+bool EmdL1::greedySolution2()
+{
+    //- Prepare auxiliary array, D=H1-H2
+    int	c,r;
+    floatArray2D D(binsDim1);
+    for(r=0; r<binsDim1; r++)
+    {
+        D[r].resize(binsDim2);
+        for(c=0; c<binsDim2; c++) D[r][c] = m_Nodes[r][c].d;
+    }
+    // compute integrated values along each dimension
+    std::vector<float> d2s(binsDim2);
+    d2s[0] = 0;
+    for(c=0; c<binsDim2-1; c++)
+    {
+        d2s[c+1] = d2s[c];
+        for(r=0; r<binsDim1; r++) d2s[c+1]-= D[r][c];
+    }
+
+    std::vector<float> d1s(binsDim1);
+    d1s[0] = 0;
+    for(r=0; r<binsDim1-1; r++)
+    {
+        d1s[r+1] = d1s[r];
+        for(c=0; c<binsDim2; c++) d1s[r+1]-= D[r][c];
+    }
+
+    //- Greedy algorithm for initial solution
+    cvPEmdEdge pBV;
+    float dFlow;
+    bool bUpward = false;
+    nNBV = 0; // number of NON-BV edges
+
+    for(c=0; c<binsDim2-1; c++)
+    for(r=0; r<binsDim1; r++)
+    {
+        dFlow = D[r][c];
+        bUpward = (r<binsDim1-1) && (fabs(dFlow+d2s[c+1]) > fabs(dFlow+d1s[r+1]));	// Move upward or right
+
+        // modify basic variables, record BV and related values
+        if(bUpward)
+        {
+            // move to up
+            pBV	= &(m_EdgesUp[r][c]);
+            m_NBVEdges[nNBV++]	= &(m_EdgesRight[r][c]);
+            D[r+1][c] += dFlow;		// auxilary matrix maintanence
+            d1s[r+1] += dFlow;		// auxilary matrix maintanence
+        }
+        else
+        {
+            // move to right, no other choice
+            pBV	= &(m_EdgesRight[r][c]);
+            if(r<binsDim1-1)
+                m_NBVEdges[nNBV++]	= &(m_EdgesUp[r][c]);
+
+            D[r][c+1] += dFlow;		// auxilary matrix maintanence
+            d2s[c+1] += dFlow;		// auxilary matrix maintanence
+        }
+        pBV->pParent->pChild = pBV;
+        pBV->flow = fabs(dFlow);
+        pBV->iDir = dFlow>0;		// 1:outward, 0:inward
+    }
+
+    //- rightmost column, no choice but move upward
+    c = binsDim2-1;
+    for(r=0; r<binsDim1-1; r++)
+    {
+        dFlow = D[r][c];
+        pBV = &(m_EdgesUp[r][c]);
+        D[r+1][c] += dFlow;		// auxilary matrix maintanence
+        pBV->pParent->pChild= pBV;
+        pBV->flow = fabs(dFlow);
+        pBV->iDir = dFlow>0;		// 1:outward, 0:inward
+    }
+    return true;
+}
+
+bool EmdL1::greedySolution3()
+{
+    //- Prepare auxiliary array, D=H1-H2
+    int i1,i2,i3;
+    std::vector<floatArray2D> D(binsDim1);
+    for(i1=0; i1<binsDim1; i1++)
+    {
+        D[i1].resize(binsDim2);
+        for(i2=0; i2<binsDim2; i2++)
+        {
+            D[i1][i2].resize(binsDim3);
+            for(i3=0; i3<binsDim3; i3++)
+                D[i1][i2][i3] = m_3dNodes[i1][i2][i3].d;
+        }
+    }
+
+    // compute integrated values along each dimension
+    std::vector<float> d1s(binsDim1);
+    d1s[0] = 0;
+    for(i1=0; i1<binsDim1-1; i1++)
+    {
+        d1s[i1+1] = d1s[i1];
+        for(i2=0; i2<binsDim2; i2++)
+        {
+            for(i3=0; i3<binsDim3; i3++)
+                d1s[i1+1] -= D[i1][i2][i3];
+        }
+    }
+
+    std::vector<float> d2s(binsDim2);
+    d2s[0] = 0;
+    for(i2=0; i2<binsDim2-1; i2++)
+    {
+        d2s[i2+1] = d2s[i2];
+        for(i1=0; i1<binsDim1; i1++)
+        {
+            for(i3=0; i3<binsDim3; i3++)
+                d2s[i2+1] -= D[i1][i2][i3];
+        }
+    }
+
+    std::vector<float> d3s(binsDim3);
+    d3s[0] = 0;
+    for(i3=0; i3<binsDim3-1; i3++)
+    {
+        d3s[i3+1]	= d3s[i3];
+        for(i1=0; i1<binsDim1; i1++)
+        {
+            for(i2=0; i2<binsDim2; i2++)
+                d3s[i3+1] -= D[i1][i2][i3];
+        }
+    }
+
+    //- Greedy algorithm for initial solution
+    cvPEmdEdge pBV;
+    float dFlow, f1,f2,f3;
+    nNBV = 0; // number of NON-BV edges
+    for(i3=0; i3<binsDim3; i3++)
+    {
+        for(i2=0; i2<binsDim2; i2++)
+        {
+            for(i1=0; i1<binsDim1; i1++)
+            {
+                if(i3==binsDim3-1 && i2==binsDim2-1 && i1==binsDim1-1) break;
+
+                //- determine which direction to move, either right or upward
+                dFlow = D[i1][i2][i3];
+                f1 = (i1<(binsDim1-1))?fabs(dFlow+d1s[i1+1]):std::numeric_limits<float>::max();
+                f2 = (i2<(binsDim2-1))?fabs(dFlow+d2s[i2+1]):std::numeric_limits<float>::max();
+                f3 = (i3<(binsDim3-1))?fabs(dFlow+d3s[i3+1]):std::numeric_limits<float>::max();
+
+                if(f1<f2 && f1<f3)
+                {
+                    pBV	= &(m_3dEdgesUp[i1][i2][i3]); // up
+                    if(i2<binsDim2-1) m_NBVEdges[nNBV++] = &(m_3dEdgesRight[i1][i2][i3]);	// right
+                    if(i3<binsDim3-1) m_NBVEdges[nNBV++] = &(m_3dEdgesDeep[i1][i2][i3]); // deep
+                    D[i1+1][i2][i3]	+= dFlow; // maintain auxilary matrix
+                    d1s[i1+1] += dFlow;
+                }
+                else if(f2<f3)
+                {
+                    pBV	= &(m_3dEdgesRight[i1][i2][i3]); // right
+                    if(i1<binsDim1-1) m_NBVEdges[nNBV++] = &(m_3dEdgesUp[i1][i2][i3]); // up
+                    if(i3<binsDim3-1) m_NBVEdges[nNBV++] = &(m_3dEdgesDeep[i1][i2][i3]); // deep
+                    D[i1][i2+1][i3]	+= dFlow; // maintain auxilary matrix
+                    d2s[i2+1] += dFlow;
+                }
+                else
+                {
+                    pBV	= &(m_3dEdgesDeep[i1][i2][i3]); // deep
+                    if(i2<binsDim2-1) m_NBVEdges[nNBV++] = &(m_3dEdgesRight[i1][i2][i3]);	// right
+                    if(i1<binsDim1-1) m_NBVEdges[nNBV++] = &(m_3dEdgesUp[i1][i2][i3]); // up
+                    D[i1][i2][i3+1]	+= dFlow; // maintain auxilary matrix
+                    d3s[i3+1] += dFlow;
+                }
+
+                pBV->flow = fabs(dFlow);
+                pBV->iDir = dFlow>0; // 1:outward, 0:inward
+                pBV->pParent->pChild= pBV;
+            }
+        }
+    }
+    return true;
+}
+
+void EmdL1::initBVTree()
+{
+    // initialize BVTree from the initial BF solution
+    //- Using the center of the graph as the root
+    int r = (int)(0.5*binsDim1-.5);
+    int c = (int)(0.5*binsDim2-.5);
+    int z = (int)(0.5*binsDim3-.5);
+    m_pRoot	= dimension==2 ? &(m_Nodes[r][c]) : &(m_3dNodes[r][c][z]);
+    m_pRoot->u = 0;
+    m_pRoot->iLevel	= 0;
+    m_pRoot->pParent= NULL;
+    m_pRoot->pPEdge	= NULL;
+
+    //- Prepare a queue
+    m_auxQueue[0] = m_pRoot;
+    int nQueue = 1; // length of queue
+    int iQHead = 0; // head of queue
+
+    //- Recursively build subtrees
+    cvPEmdEdge pCurE=NULL, pNxtE=NULL;
+    cvPEmdNode pCurN=NULL, pNxtN=NULL;
+    int	nBin = binsDim1*binsDim2*std::max(binsDim3,1);
+    while(iQHead<nQueue && nQueue<nBin)
+    {
+        pCurN = m_auxQueue[iQHead++];	// pop out from queue
+        r = pCurN->pos[0];
+        c = pCurN->pos[1];
+        z = pCurN->pos[2];
+
+        // check connection from itself
+        pCurE = pCurN->pChild;	// the initial child from initial solution
+        if(pCurE)
+        {
+            pNxtN = pCurE->pChild;
+            pNxtN->pParent = pCurN;
+            pNxtN->pPEdge = pCurE;
+            m_auxQueue[nQueue++] = pNxtN;
+        }
+
+        // check four neighbor nodes
+        int	nNB	= dimension==2?4:6;
+        for(int k=0;k<nNB;k++)
+        {
+            if(dimension==2)
+            {
+                if(k==0 && c>0) pNxtN = &(m_Nodes[r][c-1]);		// left
+                else if(k==1 && r>0) pNxtN	= &(m_Nodes[r-1][c]);		// down
+                else if(k==2 && c<binsDim2-1) pNxtN	= &(m_Nodes[r][c+1]);		// right
+                else if(k==3 && r<binsDim1-1) pNxtN	= &(m_Nodes[r+1][c]);		// up
+                else continue;
+            }
+            else if(dimension==3)
+            {
+                if(k==0 && c>0) pNxtN = &(m_3dNodes[r][c-1][z]); // left
+                else if(k==1 && c<binsDim2-1) pNxtN	= &(m_3dNodes[r][c+1][z]); // right
+                else if(k==2 && r>0) pNxtN	= &(m_3dNodes[r-1][c][z]); // down
+                else if(k==3 && r<binsDim1-1) pNxtN	= &(m_3dNodes[r+1][c][z]); // up
+                else if(k==4 && z>0) pNxtN = &(m_3dNodes[r][c][z-1]); // shallow
+                else if(k==5 && z<binsDim3-1) pNxtN	= &(m_3dNodes[r][c][z+1]); // deep
+                else continue;
+            }
+            if(pNxtN != pCurN->pParent)
+            {
+                pNxtE = pNxtN->pChild;
+                if(pNxtE && pNxtE->pChild==pCurN) // has connection
+                {
+                    pNxtN->pParent = pCurN;
+                    pNxtN->pPEdge = pNxtE;
+                    pNxtN->pChild = NULL;
+                    m_auxQueue[nQueue++] = pNxtN;
+
+                    pNxtE->pParent = pCurN; // reverse direction
+                    pNxtE->pChild = pNxtN;
+                    pNxtE->iDir = !pNxtE->iDir;
+
+                    if(pCurE) pCurE->pNxt = pNxtE;	// add to edge list
+                    else pCurN->pChild = pNxtE;
+                    pCurE = pNxtE;
+                }
+            }
+        }
+    }
+}
+
+void EmdL1::updateSubtree(cvPEmdNode pRoot)
+{
+    // Initialize auxiliary queue
+    m_auxQueue[0] = pRoot;
+    int nQueue = 1; // queue length
+    int iQHead = 0; // head of queue
+
+    // BFS browing
+    cvPEmdNode pCurN=NULL,pNxtN=NULL;
+    cvPEmdEdge pCurE=NULL;
+    while(iQHead<nQueue)
+    {
+        pCurN = m_auxQueue[iQHead++];	// pop out from queue
+        pCurE = pCurN->pChild;
+
+        // browsing all children
+        while(pCurE)
+        {
+            pNxtN = pCurE->pChild;
+            pNxtN->iLevel = pCurN->iLevel+1;
+            pNxtN->u = pCurE->iDir ? (pCurN->u - 1) : (pCurN->u + 1);
+            pCurE = pCurE->pNxt;
+            m_auxQueue[nQueue++] = pNxtN;
+        }
+    }
+}
+
+bool EmdL1::isOptimal()
+{
+    int iC, iMinC = 0;
+    cvPEmdEdge pE;
+    m_pEnter = NULL;
+    m_iEnter = -1;
+
+    // test each NON-BV edges
+    for(int k=0; k<nNBV; ++k)
+    {
+        pE = m_NBVEdges[k];
+        iC = 1 - pE->pParent->u + pE->pChild->u;
+        if(iC<iMinC)
+        {
+            iMinC = iC;
+            m_iEnter= k;
+        }
+        else
+        {
+            // Try reversing the direction
+            iC	= 1 + pE->pParent->u - pE->pChild->u;
+            if(iC<iMinC)
+            {
+                iMinC = iC;
+                m_iEnter= k;
+            }
+        }
+    }
+
+    if(m_iEnter>=0)
+    {
+        m_pEnter = m_NBVEdges[m_iEnter];
+        if(iMinC == (1 - m_pEnter->pChild->u + m_pEnter->pParent->u))	{
+            // reverse direction
+            cvPEmdNode pN = m_pEnter->pParent;
+            m_pEnter->pParent = m_pEnter->pChild;
+            m_pEnter->pChild = pN;
+        }
+
+        m_pEnter->iDir = 1;
+    }
+    return m_iEnter==-1;
+}
+
+void EmdL1::findNewSolution()
+{
+    // Find loop formed by adding the Enter BV edge.
+    findLoopFromEnterBV();
+    // Modify flow values along the loop
+    cvPEmdEdge pE = NULL;
+    float	minFlow = m_pLeave->flow;
+    int k;
+    for(k=0; k<m_iFrom; k++)
+    {
+        pE = m_fromLoop[k];
+        if(pE->iDir) pE->flow += minFlow; // outward
+        else pE->flow -= minFlow; // inward
+    }
+    for(k=0; k<m_iTo; k++)
+    {
+        pE = m_toLoop[k];
+        if(pE->iDir) pE->flow -= minFlow; // outward
+        else pE->flow += minFlow; // inward
+    }
+
+    // Update BV Tree, removing the Leaving-BV edge
+    cvPEmdNode pLParentN = m_pLeave->pParent;
+    cvPEmdNode pLChildN = m_pLeave->pChild;
+    cvPEmdEdge pPreE = pLParentN->pChild;
+    if(pPreE==m_pLeave)
+    {
+        pLParentN->pChild = m_pLeave->pNxt; // Leaving-BV is the first child
+    }
+    else
+    {
+        while(pPreE->pNxt != m_pLeave)
+            pPreE	= pPreE->pNxt;
+        pPreE->pNxt	= m_pLeave->pNxt; // remove Leaving-BV from child list
+    }
+    pLChildN->pParent = NULL;
+    pLChildN->pPEdge = NULL;
+
+    m_NBVEdges[m_iEnter]= m_pLeave; // put the leaving-BV into the NBV array
+
+    // Add the Enter BV edge
+    cvPEmdNode pEParentN = m_pEnter->pParent;
+    cvPEmdNode pEChildN = m_pEnter->pChild;
+    m_pEnter->flow = minFlow;
+    m_pEnter->pNxt = pEParentN->pChild;		// insert the Enter BV as the first child
+    pEParentN->pChild = m_pEnter;					//		of its parent
+
+    // Recursively update the tree start from pEChildN
+    cvPEmdNode pPreN = pEParentN;
+    cvPEmdNode pCurN = pEChildN;
+    cvPEmdNode pNxtN;
+    cvPEmdEdge pNxtE, pPreE0;
+    pPreE = m_pEnter;
+    while(pCurN)
+    {
+        pNxtN = pCurN->pParent;
+        pNxtE = pCurN->pPEdge;
+        pCurN->pParent = pPreN;
+        pCurN->pPEdge = pPreE;
+        if(pNxtN)
+        {
+            // remove the edge from pNxtN's child list
+            if(pNxtN->pChild==pNxtE)
+            {
+                pNxtN->pChild	= pNxtE->pNxt;			// first child
+            }
+            else
+            {
+                pPreE0	= pNxtN->pChild;
+                while(pPreE0->pNxt != pNxtE)
+                    pPreE0	= pPreE0->pNxt;
+                pPreE0->pNxt	= pNxtE->pNxt;			// remove Leaving-BV from child list
+            }
+            // reverse the parent-child direction
+            pNxtE->pParent = pCurN;
+            pNxtE->pChild = pNxtN;
+            pNxtE->iDir = !pNxtE->iDir;
+            pNxtE->pNxt = pCurN->pChild;
+            pCurN->pChild = pNxtE;
+            pPreE = pNxtE;
+            pPreN = pCurN;
+        }
+        pCurN = pNxtN;
+    }
+
+    // Update U at the child of the Enter BV
+    pEChildN->u = m_pEnter->iDir?(pEParentN->u-1):(pEParentN->u + 1);
+    pEChildN->iLevel = pEParentN->iLevel+1;
+}
+
+void EmdL1::findLoopFromEnterBV()
+{
+    // Initialize Leaving-BV edge
+    float minFlow	= std::numeric_limits<float>::max();
+    cvPEmdEdge pE = NULL;
+    int iLFlag = 0;	// 0: in the FROM list, 1: in the TO list
+
+    // Using two loop list to store the loop nodes
+    cvPEmdNode pFrom = m_pEnter->pParent;
+    cvPEmdNode pTo = m_pEnter->pChild;
+    m_iFrom	= 0;
+    m_iTo = 0;
+    m_pLeave = NULL;
+
+    // Trace back to make pFrom and pTo at the same level
+    while(pFrom->iLevel > pTo->iLevel)
+    {
+        pE = pFrom->pPEdge;
+        m_fromLoop[m_iFrom++] = pE;
+        if(!pE->iDir && pE->flow<minFlow)
+        {
+            minFlow = pE->flow;
+            m_pLeave = pE;
+            iLFlag = 0;	// 0: in the FROM list
+        }
+        pFrom = pFrom->pParent;
+    }
+
+    while(pTo->iLevel > pFrom->iLevel)
+    {
+        pE = pTo->pPEdge;
+        m_toLoop[m_iTo++] = pE;
+        if(pE->iDir && pE->flow<minFlow)
+        {
+            minFlow = pE->flow;
+            m_pLeave = pE;
+            iLFlag = 1;	// 1: in the TO list
+        }
+        pTo	= pTo->pParent;
+    }
+
+    // Trace pTo and pFrom simultaneously till find their common ancester
+    while(pTo!=pFrom)
+    {
+        pE = pFrom->pPEdge;
+        m_fromLoop[m_iFrom++] = pE;
+        if(!pE->iDir && pE->flow<minFlow)
+        {
+            minFlow = pE->flow;
+            m_pLeave = pE;
+            iLFlag = 0;	// 0: in the FROM list, 1: in the TO list
+        }
+        pFrom = pFrom->pParent;
+
+        pE = pTo->pPEdge;
+        m_toLoop[m_iTo++] = pE;
+        if(pE->iDir && pE->flow<minFlow)
+        {
+            minFlow = pE->flow;
+            m_pLeave = pE;
+            iLFlag = 1;	// 0: in the FROM list, 1: in the TO list
+        }
+        pTo	= pTo->pParent;
+    }
+
+    // Reverse the direction of the Enter BV edge if necessary
+    if(iLFlag==0)
+    {
+        cvPEmdNode pN = m_pEnter->pParent;
+        m_pEnter->pParent = m_pEnter->pChild;
+        m_pEnter->pChild = pN;
+        m_pEnter->iDir = !m_pEnter->iDir;
+    }
+}
+
+float EmdL1::compuTotalFlow()
+{
+    // Computing the total flow as the final distance
+    float f = 0;
+
+    // Initialize auxiliary queue
+    m_auxQueue[0] = m_pRoot;
+    int nQueue = 1; // length of queue
+    int iQHead = 0; // head of queue
+
+    // BFS browing the tree
+    cvPEmdNode pCurN=NULL,pNxtN=NULL;
+    cvPEmdEdge pCurE=NULL;
+    while(iQHead<nQueue)
+    {
+        pCurN = m_auxQueue[iQHead++];	// pop out from queue
+        pCurE = pCurN->pChild;
+
+        // browsing all children
+        while(pCurE)
+        {
+            f += pCurE->flow;
+            pNxtN = pCurE->pChild;
+            pCurE = pCurE->pNxt;
+            m_auxQueue[nQueue++] = pNxtN;
+        }
+    }
+    return f;
+}
+
+/****************************************************************************************\
+*                                   EMDL1 Function                                      *
+\****************************************************************************************/
+
+float cv::EMDL1(InputArray _signature1, InputArray _signature2)
+{
+    Mat signature1 = _signature1.getMat(), signature2 = _signature2.getMat();
+    EmdL1 emdl1;
+    return emdl1.getEMDL1(signature1, signature2);
+}

modules/shape/src/emdL1_def.hpp

@@ -0,0 +1,140 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include <stdlib.h>
+#include <math.h>
+#include <vector>
+
+/****************************************************************************************\
+*                                   For EMDL1 Framework                                 *
+\****************************************************************************************/
+typedef struct cvEMDEdge* cvPEmdEdge;
+typedef struct cvEMDNode* cvPEmdNode;
+struct cvEMDNode
+{
+    int pos[3]; // grid position
+    float d; // initial value
+    int u;
+    // tree maintainance
+    int iLevel; // level in the tree, 0 means root
+    cvPEmdNode pParent; // pointer to its parent
+    cvPEmdEdge pChild;
+    cvPEmdEdge pPEdge; // point to the edge coming out from its parent
+};
+struct cvEMDEdge
+{
+    float flow; // initial value
+    int iDir; // 1:outward, 0:inward
+    // tree maintainance
+    cvPEmdNode pParent; // point to its parent
+    cvPEmdNode pChild; // the child node
+    cvPEmdEdge pNxt; // next child/edge
+};
+typedef std::vector<cvEMDNode> cvEMDNodeArray;
+typedef std::vector<cvEMDEdge> cvEMDEdgeArray;
+typedef std::vector<cvEMDNodeArray> cvEMDNodeArray2D;
+typedef std::vector<cvEMDEdgeArray> cvEMDEdgeArray2D;
+typedef std::vector<float> floatArray;
+typedef std::vector<floatArray> floatArray2D;
+
+/****************************************************************************************\
+*                                   EMDL1 Class                                         *
+\****************************************************************************************/
+class EmdL1
+{
+public:
+    EmdL1()
+    {
+        m_pRoot	= NULL;
+        binsDim1 = 0;
+        binsDim2 = 0;
+        binsDim3 = 0;
+        dimension = 0;
+        nMaxIt = 500;
+    }
+
+    ~EmdL1()
+    {
+    }
+
+    float getEMDL1(cv::Mat &sig1, cv::Mat &sig2);
+    void setMaxIteration(int _nMaxIt);
+
+private:
+    //-- SubFunctions called in the EMD algorithm
+    bool initBaseTrees(int n1=0, int n2=0, int n3=0);
+    bool fillBaseTrees(float *H1, float *H2);
+    bool greedySolution();
+    bool greedySolution2();
+    bool greedySolution3();
+    void initBVTree();
+    void updateSubtree(cvPEmdNode pRoot);
+    bool isOptimal();
+    void findNewSolution();
+    void findLoopFromEnterBV();
+    float compuTotalFlow();
+
+private:
+    int dimension;
+    int binsDim1, binsDim2, binsDim3; // the hitogram contains m_n1 rows and m_n2 columns
+    int nNBV; // number of Non-Basic Variables (NBV)
+    int nMaxIt;
+    cvEMDNodeArray2D m_Nodes; // all nodes
+    cvEMDEdgeArray2D m_EdgesRight; // all edges to right
+    cvEMDEdgeArray2D m_EdgesUp; // all edges to upward
+    std::vector<cvEMDNodeArray2D>	m_3dNodes; // all nodes for 3D
+    std::vector<cvEMDEdgeArray2D>	m_3dEdgesRight; // all edges to right, 3D
+    std::vector<cvEMDEdgeArray2D>	m_3dEdgesUp; // all edges to upward, 3D
+    std::vector<cvEMDEdgeArray2D>	m_3dEdgesDeep; // all edges to deep, 3D
+    std::vector<cvPEmdEdge> m_NBVEdges; // pointers to all NON-BV edges
+    std::vector<cvPEmdNode> m_auxQueue; // auxiliary node queue
+    cvPEmdNode m_pRoot; // root of the BV Tree
+    cvPEmdEdge m_pEnter; // Enter BV edge
+    int m_iEnter; // Enter BV edge, index in m_NBVEdges
+    cvPEmdEdge m_pLeave; // Leave BV edge
+    int m_nItr; // number of iteration
+    // auxiliary variables for searching a new loop
+    std::vector<cvPEmdEdge> m_fromLoop;
+    std::vector<cvPEmdEdge> m_toLoop;
+    int	m_iFrom;
+    int m_iTo;
+};

modules/shape/src/haus_dis.cpp

@@ -0,0 +1,149 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+class HausdorffDistanceExtractorImpl : public HausdorffDistanceExtractor
+{
+public:
+    /* Constructor */
+    HausdorffDistanceExtractorImpl(int _distanceFlag = NORM_L1, float _rankProportion=0.6)
+    {
+        distanceFlag = _distanceFlag;
+        rankProportion = _rankProportion;
+        name_ = "ShapeDistanceExtractor.HAU";
+    }
+
+    /* Destructor */
+    ~HausdorffDistanceExtractorImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operator
+    virtual float computeDistance(InputArray contour1, InputArray contour2);
+
+    //! Setters/Getters
+    virtual void setDistanceFlag(int _distanceFlag) {distanceFlag=_distanceFlag;}
+    virtual int getDistanceFlag() const {return distanceFlag;}
+
+    virtual void setRankProportion(float _rankProportion)
+    {
+        CV_Assert((_rankProportion>0) && (_rankProportion<=1));
+        rankProportion=_rankProportion;
+    }
+    virtual float getRankProportion() const {return rankProportion;}
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "distance" << distanceFlag
+           << "rank" << rankProportion;
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        distanceFlag = (int)fn["distance"];
+        rankProportion = (float)fn["rank"];
+    }
+
+private:
+    int distanceFlag;
+    float rankProportion;
+
+protected:
+    String name_;
+};
+
+//! Hausdorff distance for a pair of set of points
+static float _apply(const Mat &set1, const Mat &set2, int distType, double propRank)
+{
+    // Building distance matrix //
+    Mat disMat(set1.cols, set2.cols, CV_32F);
+    int K = int(propRank*(disMat.rows-1));
+
+    for (int r=0; r<disMat.rows; r++)
+    {
+        for (int c=0; c<disMat.cols; c++)
+        {
+            Point2f diff = set1.at<Point2f>(0,r)-set2.at<Point2f>(0,c);
+            disMat.at<float>(r,c) = (float)norm(Mat(diff), distType);
+        }
+    }
+
+    Mat shortest(disMat.rows,1,CV_32F);
+    for (int ii=0; ii<disMat.rows; ii++)
+    {
+        Mat therow = disMat.row(ii);
+        double mindis;
+        minMaxIdx(therow, &mindis);
+        shortest.at<float>(ii,0) = float(mindis);
+    }
+    Mat sorted;
+    cv::sort(shortest, sorted, SORT_EVERY_ROW | SORT_DESCENDING);
+    return sorted.at<float>(K,0);
+}
+
+float HausdorffDistanceExtractorImpl::computeDistance(InputArray contour1, InputArray contour2)
+{
+    Mat set1=contour1.getMat(), set2=contour2.getMat();
+    if (set1.type() != CV_32F)
+        set1.convertTo(set1, CV_32F);
+    if (set2.type() != CV_32F)
+        set2.convertTo(set2, CV_32F);
+    CV_Assert((set1.channels()==2) && (set1.cols>0));
+    CV_Assert((set2.channels()==2) && (set2.cols>0));
+    return std::max( _apply(set1, set2, distanceFlag, rankProportion),
+                     _apply(set2, set1, distanceFlag, rankProportion) );
+}
+
+Ptr <HausdorffDistanceExtractor> createHausdorffDistanceExtractor(int distanceFlag, float rankProp)
+{
+    return Ptr<HausdorffDistanceExtractor>(new HausdorffDistanceExtractorImpl(distanceFlag, rankProp));
+}
+
+} // cv

modules/shape/src/hist_cost.cpp

@@ -0,0 +1,545 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+/*!  */
+class NormHistogramCostExtractorImpl : public NormHistogramCostExtractor
+{
+public:
+    /* Constructors */
+    NormHistogramCostExtractorImpl(int _flag, int _nDummies, float _defaultCost)
+    {
+        flag=_flag;
+        nDummies=_nDummies;
+        defaultCost=_defaultCost;
+        name_ = "HistogramCostExtractor.NOR";
+    }
+
+    /* Destructor */
+    ~NormHistogramCostExtractorImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operator
+    virtual void buildCostMatrix(InputArray descriptors1, InputArray descriptors2, OutputArray costMatrix);
+
+    //! Setters/Getters
+    void setNDummies(int _nDummies)
+    {
+        nDummies=_nDummies;
+    }
+
+    int getNDummies() const
+    {
+        return nDummies;
+    }
+
+    void setDefaultCost(float _defaultCost)
+    {
+        defaultCost=_defaultCost;
+    }
+
+    float getDefaultCost() const
+    {
+        return defaultCost;
+    }
+
+    virtual void setNormFlag(int _flag)
+    {
+        flag=_flag;
+    }
+
+    virtual int getNormFlag() const
+    {
+        return flag;
+    }
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "flag" << flag
+           << "dummies" << nDummies
+           << "default" << defaultCost;
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        flag = (int)fn["flag"];
+        nDummies = (int)fn["dummies"];
+        defaultCost = (float)fn["default"];
+    }
+
+private:
+    int flag;
+    int nDummies;
+    float defaultCost;
+
+protected:
+    String name_;
+};
+
+void NormHistogramCostExtractorImpl::buildCostMatrix(InputArray _descriptors1, InputArray _descriptors2, OutputArray _costMatrix)
+{
+    // size of the costMatrix with dummies //
+    Mat descriptors1=_descriptors1.getMat();
+    Mat descriptors2=_descriptors2.getMat();
+    int costrows = std::max(descriptors1.rows, descriptors2.rows)+nDummies;
+    _costMatrix.create(costrows, costrows, CV_32F);
+    Mat costMatrix=_costMatrix.getMat();
+
+
+    // Obtain copies of the descriptors //
+    cv::Mat scd1 = descriptors1.clone();
+    cv::Mat scd2 = descriptors2.clone();
+
+    // row normalization //
+    for(int i=0; i<scd1.rows; i++)
+    {
+        scd1.row(i)/=(sum(scd1.row(i))[0]+FLT_EPSILON);
+    }
+    for(int i=0; i<scd2.rows; i++)
+    {
+        scd2.row(i)/=(sum(scd2.row(i))[0]+FLT_EPSILON);
+    }
+
+    // Compute the Cost Matrix //
+    for(int i=0; i<costrows; i++)
+    {
+        for(int j=0; j<costrows; j++)
+        {
+            if (i<scd1.rows && j<scd2.rows)
+            {
+                Mat columnDiff = scd1.row(i)-scd2.row(j);
+                costMatrix.at<float>(i,j)=(float)norm(columnDiff, flag);
+            }
+            else
+            {
+                costMatrix.at<float>(i,j)=defaultCost;
+            }
+        }
+    }
+}
+
+Ptr <HistogramCostExtractor> createNormHistogramCostExtractor(int flag, int nDummies, float defaultCost)
+{
+    return Ptr <HistogramCostExtractor>( new NormHistogramCostExtractorImpl(flag, nDummies, defaultCost) );
+}
+
+/*!  */
+class EMDHistogramCostExtractorImpl : public EMDHistogramCostExtractor
+{
+public:
+    /* Constructors */
+    EMDHistogramCostExtractorImpl(int _flag, int _nDummies, float _defaultCost)
+    {
+        flag=_flag;
+        nDummies=_nDummies;
+        defaultCost=_defaultCost;
+        name_ = "HistogramCostExtractor.EMD";
+    }
+
+    /* Destructor */
+    ~EMDHistogramCostExtractorImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operator
+    virtual void buildCostMatrix(InputArray descriptors1, InputArray descriptors2, OutputArray costMatrix);
+
+    //! Setters/Getters
+    void setNDummies(int _nDummies)
+    {
+        nDummies=_nDummies;
+    }
+
+    int getNDummies() const
+    {
+        return nDummies;
+    }
+
+    void setDefaultCost(float _defaultCost)
+    {
+        defaultCost=_defaultCost;
+    }
+
+    float getDefaultCost() const
+    {
+        return defaultCost;
+    }
+
+    virtual void setNormFlag(int _flag)
+    {
+        flag=_flag;
+    }
+
+    virtual int getNormFlag() const
+    {
+        return flag;
+    }
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "flag" << flag
+           << "dummies" << nDummies
+           << "default" << defaultCost;
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        flag = (int)fn["flag"];
+        nDummies = (int)fn["dummies"];
+        defaultCost = (float)fn["default"];
+    }
+
+private:
+    int flag;
+    int nDummies;
+    float defaultCost;
+
+protected:
+    String name_;
+};
+
+void EMDHistogramCostExtractorImpl::buildCostMatrix(InputArray _descriptors1, InputArray _descriptors2, OutputArray _costMatrix)
+{
+    // size of the costMatrix with dummies //
+    Mat descriptors1=_descriptors1.getMat();
+    Mat descriptors2=_descriptors2.getMat();
+    int costrows = std::max(descriptors1.rows, descriptors2.rows)+nDummies;
+    _costMatrix.create(costrows, costrows, CV_32F);
+    Mat costMatrix=_costMatrix.getMat();
+
+    // Obtain copies of the descriptors //
+    cv::Mat scd1=descriptors1.clone();
+    cv::Mat scd2=descriptors2.clone();
+
+    // row normalization //
+    for(int i=0; i<scd1.rows; i++)
+    {
+        cv::Mat row = scd1.row(i);
+        scd1.row(i)/=(sum(row)[0]+FLT_EPSILON);
+    }
+    for(int i=0; i<scd2.rows; i++)
+    {
+        cv::Mat row = scd2.row(i);
+        scd2.row(i)/=(sum(row)[0]+FLT_EPSILON);
+    }
+
+    // Compute the Cost Matrix //
+    for(int i=0; i<costrows; i++)
+    {
+        for(int j=0; j<costrows; j++)
+        {
+            if (i<scd1.rows && j<scd2.rows)
+            {
+                cv::Mat sig1(scd1.cols,2,CV_32F), sig2(scd2.cols,2,CV_32F);
+                sig1.col(0)=scd1.row(i).t();
+                sig2.col(0)=scd2.row(j).t();
+                for (int k=0; k<sig1.rows; k++)
+                {
+                    sig1.at<float>(k,1)=float(k);
+                }
+                for (int k=0; k<sig2.rows; k++)
+                {
+                    sig2.at<float>(k,1)=float(k);
+                }
+
+                costMatrix.at<float>(i,j) = cv::EMD(sig1, sig2, flag);
+            }
+            else
+            {
+                costMatrix.at<float>(i,j) = defaultCost;
+            }
+        }
+    }
+}
+
+Ptr <HistogramCostExtractor> createEMDHistogramCostExtractor(int flag, int nDummies, float defaultCost)
+{
+    return Ptr <HistogramCostExtractor>( new EMDHistogramCostExtractorImpl(flag, nDummies, defaultCost) );
+}
+
+/*!  */
+class ChiHistogramCostExtractorImpl : public ChiHistogramCostExtractor
+{
+public:
+    /* Constructors */
+    ChiHistogramCostExtractorImpl(int _nDummies, float _defaultCost)
+    {
+        name_ = "HistogramCostExtractor.CHI";
+        nDummies=_nDummies;
+        defaultCost=_defaultCost;
+    }
+
+    /* Destructor */
+    ~ChiHistogramCostExtractorImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operator
+    virtual void buildCostMatrix(InputArray descriptors1, InputArray descriptors2, OutputArray costMatrix);
+
+    //! setters / getters
+    void setNDummies(int _nDummies)
+    {
+        nDummies=_nDummies;
+    }
+
+    int getNDummies() const
+    {
+        return nDummies;
+    }
+
+    void setDefaultCost(float _defaultCost)
+    {
+        defaultCost=_defaultCost;
+    }
+
+    float getDefaultCost() const
+    {
+        return defaultCost;
+    }
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "dummies" << nDummies
+           << "default" << defaultCost;
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        nDummies = (int)fn["dummies"];
+        defaultCost = (float)fn["default"];
+    }
+
+protected:
+    String name_;
+    int nDummies;
+    float defaultCost;
+};
+
+void ChiHistogramCostExtractorImpl::buildCostMatrix(InputArray _descriptors1, InputArray _descriptors2, OutputArray _costMatrix)
+{
+    // size of the costMatrix with dummies //
+    Mat descriptors1=_descriptors1.getMat();
+    Mat descriptors2=_descriptors2.getMat();
+    int costrows = std::max(descriptors1.rows, descriptors2.rows)+nDummies;
+    _costMatrix.create(costrows, costrows, CV_32FC1);
+    Mat costMatrix=_costMatrix.getMat();
+
+    // Obtain copies of the descriptors //
+    cv::Mat scd1=descriptors1.clone();
+    cv::Mat scd2=descriptors2.clone();
+
+    // row normalization //
+    for(int i=0; i<scd1.rows; i++)
+    {
+        cv::Mat row = scd1.row(i);
+        scd1.row(i)/=(sum(row)[0]+FLT_EPSILON);
+    }
+    for(int i=0; i<scd2.rows; i++)
+    {
+       cv::Mat row = scd2.row(i);
+        scd2.row(i)/=(sum(row)[0]+FLT_EPSILON);
+    }
+
+    // Compute the Cost Matrix //
+    for(int i=0; i<costrows; i++)
+    {
+        for(int j=0; j<costrows; j++)
+        {
+            if (i<scd1.rows && j<scd2.rows)
+            {
+                float csum = 0;
+                for(int k=0; k<scd2.cols; k++)
+                {
+                    float resta=scd1.at<float>(i,k)-scd2.at<float>(j,k);
+                    float suma=scd1.at<float>(i,k)+scd2.at<float>(j,k);
+                    csum += resta*resta/(FLT_EPSILON+suma);
+                }
+                costMatrix.at<float>(i,j)=csum/2;
+            }
+            else
+            {
+                costMatrix.at<float>(i,j)=defaultCost;
+            }
+        }
+    }
+}
+
+Ptr <HistogramCostExtractor> createChiHistogramCostExtractor(int nDummies, float defaultCost)
+{
+    return Ptr <HistogramCostExtractor>( new ChiHistogramCostExtractorImpl(nDummies, defaultCost) );
+}
+
+/*!  */
+class EMDL1HistogramCostExtractorImpl : public EMDL1HistogramCostExtractor
+{
+public:
+    /* Constructors */
+    EMDL1HistogramCostExtractorImpl(int _nDummies, float _defaultCost)
+    {
+        name_ = "HistogramCostExtractor.CHI";
+        nDummies=_nDummies;
+        defaultCost=_defaultCost;
+    }
+
+    /* Destructor */
+    ~EMDL1HistogramCostExtractorImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operator
+    virtual void buildCostMatrix(InputArray descriptors1, InputArray descriptors2, OutputArray costMatrix);
+
+    //! setters / getters
+    void setNDummies(int _nDummies)
+    {
+        nDummies=_nDummies;
+    }
+
+    int getNDummies() const
+    {
+        return nDummies;
+    }
+
+    void setDefaultCost(float _defaultCost)
+    {
+        defaultCost=_defaultCost;
+    }
+
+    float getDefaultCost() const
+    {
+        return defaultCost;
+    }
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "dummies" << nDummies
+           << "default" << defaultCost;
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        nDummies = (int)fn["dummies"];
+        defaultCost = (float)fn["default"];
+    }
+
+protected:
+    String name_;
+    int nDummies;
+    float defaultCost;
+};
+
+void EMDL1HistogramCostExtractorImpl::buildCostMatrix(InputArray _descriptors1, InputArray _descriptors2, OutputArray _costMatrix)
+{
+    // size of the costMatrix with dummies //
+    Mat descriptors1=_descriptors1.getMat();
+    Mat descriptors2=_descriptors2.getMat();
+    int costrows = std::max(descriptors1.rows, descriptors2.rows)+nDummies;
+    _costMatrix.create(costrows, costrows, CV_32F);
+    Mat costMatrix=_costMatrix.getMat();
+
+    // Obtain copies of the descriptors //
+    cv::Mat scd1=descriptors1.clone();
+    cv::Mat scd2=descriptors2.clone();
+
+    // row normalization //
+    for(int i=0; i<scd1.rows; i++)
+    {
+        cv::Mat row = scd1.row(i);
+        scd1.row(i)/=(sum(row)[0]+FLT_EPSILON);
+    }
+    for(int i=0; i<scd2.rows; i++)
+    {
+        cv::Mat row = scd2.row(i);
+        scd2.row(i)/=(sum(row)[0]+FLT_EPSILON);
+    }
+
+    // Compute the Cost Matrix //
+    for(int i=0; i<costrows; i++)
+    {
+        for(int j=0; j<costrows; j++)
+        {
+            if (i<scd1.rows && j<scd2.rows)
+            {
+                cv::Mat sig1(scd1.cols,1,CV_32F), sig2(scd2.cols,1,CV_32F);
+                sig1.col(0)=scd1.row(i).t();
+                sig2.col(0)=scd2.row(j).t();
+                costMatrix.at<float>(i,j) = cv::EMDL1(sig1, sig2);
+            }
+            else
+            {
+                costMatrix.at<float>(i,j) = defaultCost;
+            }
+        }
+    }
+}
+
+Ptr <HistogramCostExtractor> createEMDL1HistogramCostExtractor(int nDummies, float defaultCost)
+{
+    return Ptr <HistogramCostExtractor>( new EMDL1HistogramCostExtractorImpl(nDummies, defaultCost) );
+}
+
+} // cv

modules/shape/src/precomp.cpp

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

modules/shape/src/precomp.hpp

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

modules/shape/src/sc_dis.cpp

@@ -0,0 +1,781 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+/*
+ * Implementation of the paper Shape Matching and Object Recognition Using Shape Contexts
+ * Belongie et al., 2002 by Juan Manuel Perez for GSoC 2013.
+ */
+#include "precomp.hpp"
+#include "opencv2/core.hpp"
+#include "scd_def.hpp"
+
+namespace cv
+{
+class ShapeContextDistanceExtractorImpl : public ShapeContextDistanceExtractor
+{
+public:
+    /* Constructors */
+    ShapeContextDistanceExtractorImpl(int _nAngularBins, int _nRadialBins, float _innerRadius, float _outerRadius, int _iterations,
+                                      const Ptr<HistogramCostExtractor> &_comparer, const Ptr<ShapeTransformer> &_transformer)
+    {
+        nAngularBins=_nAngularBins;
+        nRadialBins=_nRadialBins;
+        innerRadius=_innerRadius;
+        outerRadius=_outerRadius;
+        rotationInvariant=false;
+        comparer=_comparer;
+        iterations=_iterations;
+        transformer=_transformer;
+        bendingEnergyWeight=0.3;
+        imageAppearanceWeight=0.0;
+        shapeContextWeight=1.0;
+        sigma=10;
+        name_ = "ShapeDistanceExtractor.SCD";
+    }
+
+    /* Destructor */
+    ~ShapeContextDistanceExtractorImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operator
+    virtual float computeDistance(InputArray contour1, InputArray contour2);
+
+    //! Setters/Getters
+    virtual void setAngularBins(int _nAngularBins){CV_Assert(_nAngularBins>0); nAngularBins=_nAngularBins;}
+    virtual int getAngularBins() const {return nAngularBins;}
+
+    virtual void setRadialBins(int _nRadialBins){CV_Assert(_nRadialBins>0); nRadialBins=_nRadialBins;}
+    virtual int getRadialBins() const {return nRadialBins;}
+
+    virtual void setInnerRadius(float _innerRadius) {CV_Assert(_innerRadius>0); innerRadius=_innerRadius;}
+    virtual float getInnerRadius() const {return innerRadius;}
+
+    virtual void setOuterRadius(float _outerRadius) {CV_Assert(_outerRadius>0); outerRadius=_outerRadius;}
+    virtual float getOuterRadius() const {return outerRadius;}
+
+    virtual void setRotationInvariant(bool _rotationInvariant) {rotationInvariant=_rotationInvariant;}
+    virtual bool getRotationInvariant() const {return rotationInvariant;}
+
+    virtual void setCostExtractor(Ptr<HistogramCostExtractor> _comparer) { comparer = _comparer; }
+    virtual Ptr<HistogramCostExtractor> getCostExtractor() const { return comparer; }
+
+    virtual void setShapeContextWeight(float _shapeContextWeight) {shapeContextWeight=_shapeContextWeight;}
+    virtual float getShapeContextWeight() const {return shapeContextWeight;}
+
+    virtual void setImageAppearanceWeight(float _imageAppearanceWeight) {imageAppearanceWeight=_imageAppearanceWeight;}
+    virtual float getImageAppearanceWeight() const {return imageAppearanceWeight;}
+
+    virtual void setBendingEnergyWeight(float _bendingEnergyWeight) {bendingEnergyWeight=_bendingEnergyWeight;}
+    virtual float getBendingEnergyWeight() const {return bendingEnergyWeight;}
+
+    virtual void setStdDev(float _sigma) {sigma=_sigma;}
+    virtual float getStdDev() const {return sigma;}
+
+    virtual void setImages(InputArray _image1, InputArray _image2)
+    {
+        Mat image1_=_image1.getMat(), image2_=_image2.getMat();
+        CV_Assert((image1_.depth()==0) && (image2_.depth()==0));
+        image1=image1_;
+        image2=image2_;
+    }
+
+    virtual void getImages(OutputArray _image1, OutputArray _image2) const
+    {
+        CV_Assert((!image1.empty()) && (!image2.empty()));
+        _image1.create(image1.size(), image1.type());
+        _image2.create(image2.size(), image2.type());
+        _image1.getMat()=image1;
+        _image2.getMat()=image2;
+    }
+
+    virtual void setIterations(int _iterations) {CV_Assert(_iterations>0); iterations=_iterations;}
+    virtual int getIterations() const {return iterations;}
+
+    virtual void setTransformAlgorithm(Ptr<ShapeTransformer> _transformer) {transformer=_transformer;}
+    virtual Ptr<ShapeTransformer> getTransformAlgorithm() const {return transformer;}
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "nRads" << nRadialBins
+           << "nAngs" << nAngularBins
+           << "iters" << iterations
+           << "img_1" << image1
+           << "img_2" << image2
+           << "beWei" << bendingEnergyWeight
+           << "scWei" << shapeContextWeight
+           << "iaWei" << imageAppearanceWeight
+           << "costF" << costFlag
+           << "rotIn" << rotationInvariant
+           << "sigma" << sigma;
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        nRadialBins = (int)fn["nRads"];
+        nAngularBins = (int)fn["nAngs"];
+        iterations = (int)fn["iters"];
+        bendingEnergyWeight = (float)fn["beWei"];
+        shapeContextWeight = (float)fn["scWei"];
+        imageAppearanceWeight = (float)fn["iaWei"];
+        costFlag = (int)fn["costF"];
+        sigma = (float)fn["sigma"];
+    }
+
+protected:
+    int nAngularBins;
+    int nRadialBins;
+    float innerRadius;
+    float outerRadius;
+    bool rotationInvariant;
+    int costFlag;
+    int iterations;
+    Ptr<ShapeTransformer> transformer;
+    Ptr<HistogramCostExtractor> comparer;
+    Mat image1;
+    Mat image2;
+    float bendingEnergyWeight;
+    float imageAppearanceWeight;
+    float shapeContextWeight;
+    float sigma;
+    String name_;
+};
+
+float ShapeContextDistanceExtractorImpl::computeDistance(InputArray contour1, InputArray contour2)
+{
+    // Checking //
+    Mat sset1=contour1.getMat(), sset2=contour2.getMat(), set1, set2;
+    if (set1.type() != CV_32F)
+        sset1.convertTo(set1, CV_32F);
+    else
+        sset1.copyTo(set1);
+
+    if (set2.type() != CV_32F)
+        sset2.convertTo(set2, CV_32F);
+    else
+        sset1.copyTo(set2);
+
+    CV_Assert((set1.channels()==2) && (set1.cols>0));
+    CV_Assert((set2.channels()==2) && (set2.cols>0));
+    if (imageAppearanceWeight!=0)
+    {
+        CV_Assert((!image1.empty()) && (!image2.empty()));
+    }
+
+    // Initializing Extractor, Descriptor structures and Matcher //
+    SCD set1SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, rotationInvariant);
+    Mat set1SCD;
+    SCD set2SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, rotationInvariant);
+    Mat set2SCD;
+    SCDMatcher matcher;
+    std::vector<DMatch> matches;
+
+    // Distance components (The output is a linear combination of these 3) //
+    float sDistance=0, bEnergy=0, iAppearance=0;
+    float beta;
+
+    // Initializing some variables //
+    std::vector<int> inliers1, inliers2;
+
+    Ptr<ThinPlateSplineShapeTransformer> transDown = transformer.dynamicCast<ThinPlateSplineShapeTransformer>();
+
+    Mat warpedImage;
+    int ii, jj, pt;
+
+    for (ii=0; ii<iterations; ii++)
+    {
+        // Extract SCD descriptor in the set1 //
+        set1SCE.extractSCD(set1, set1SCD, inliers1);
+
+        // Extract SCD descriptor of the set2 (TARGET) //
+        set2SCE.extractSCD(set2, set2SCD, inliers2, set1SCE.getMeanDistance());
+
+        // regularization parameter with annealing rate annRate //
+        beta=set1SCE.getMeanDistance();
+        beta *= beta;
+
+        // match //
+        matcher.matchDescriptors(set1SCD, set2SCD, matches, comparer, inliers1, inliers2);
+
+        // apply TPS transform //
+        if ( !transDown.empty() )
+            transDown->setRegularizationParameter(beta);
+        transformer->estimateTransformation(set1, set2, matches);
+        bEnergy += transformer->applyTransformation(set1, set1);
+
+        // Image appearance //
+        if (imageAppearanceWeight!=0)
+        {
+            // Have to accumulate the transformation along all the iterations
+            if (ii==0)
+            {
+                if ( !transDown.empty() )
+                {
+                    image2.copyTo(warpedImage);
+                }
+                else
+                {
+                    image1.copyTo(warpedImage);
+                }
+            }
+            transformer->warpImage(warpedImage, warpedImage);
+        }
+    }
+
+    Mat gaussWindow, diffIm;
+    if (imageAppearanceWeight!=0)
+    {
+        // compute appearance cost
+        if ( !transDown.empty() )
+        {
+            resize(warpedImage, warpedImage, image1.size());
+            Mat temp=(warpedImage-image1);
+            multiply(temp, temp, diffIm);
+        }
+        else
+        {
+            resize(warpedImage, warpedImage, image2.size());
+            Mat temp=(warpedImage-image2);
+            multiply(temp, temp, diffIm);
+        }
+        gaussWindow = Mat::zeros(warpedImage.rows, warpedImage.cols, CV_32F);
+        for (pt=0; pt<sset1.cols; pt++)
+        {
+            Point2f p = sset1.at<Point2f>(0,pt);
+            for (ii=0; ii<diffIm.rows; ii++)
+            {
+                for (jj=0; jj<diffIm.cols; jj++)
+                {
+                    float val = float(std::exp( -float( (p.x-jj)*(p.x-jj) + (p.y-ii)*(p.y-ii) )/(2*sigma*sigma) ) / (sigma*sigma*2*CV_PI));
+                    gaussWindow.at<float>(ii,jj) += val;
+                }
+            }
+        }
+
+        Mat appIm(diffIm.rows, diffIm.cols, CV_32F);
+        for (ii=0; ii<diffIm.rows; ii++)
+        {
+            for (jj=0; jj<diffIm.cols; jj++)
+            {
+                float elema=float( diffIm.at<uchar>(ii,jj) )/255;
+                float elemb=gaussWindow.at<float>(ii,jj);
+                appIm.at<float>(ii,jj) = elema*elemb;
+            }
+        }
+        iAppearance = float(cv::sum(appIm)[0]/sset1.cols);
+    }
+    sDistance = matcher.getMatchingCost();
+
+    return (sDistance*shapeContextWeight+bEnergy*bendingEnergyWeight+iAppearance*imageAppearanceWeight);
+}
+
+Ptr <ShapeContextDistanceExtractor> createShapeContextDistanceExtractor(int nAngularBins, int nRadialBins, float innerRadius, float outerRadius, int iterations,
+                                                                        const Ptr<HistogramCostExtractor> &comparer, const Ptr<ShapeTransformer> &transformer)
+{
+    return Ptr <ShapeContextDistanceExtractor> ( new ShapeContextDistanceExtractorImpl(nAngularBins, nRadialBins, innerRadius,
+                                                                                       outerRadius, iterations, comparer, transformer) );
+}
+
+//! SCD
+void SCD::extractSCD(cv::Mat &contour, cv::Mat &descriptors, const std::vector<int> &queryInliers, const float _meanDistance)
+{
+    cv::Mat contourMat = contour;
+    cv::Mat disMatrix = cv::Mat::zeros(contourMat.cols, contourMat.cols, CV_32F);
+    cv::Mat angleMatrix = cv::Mat::zeros(contourMat.cols, contourMat.cols, CV_32F);
+
+    std::vector<double> logspaces, angspaces;
+    logarithmicSpaces(logspaces);
+    angularSpaces(angspaces);
+    buildNormalizedDistanceMatrix(contourMat, disMatrix, queryInliers, _meanDistance);
+    buildAngleMatrix(contourMat, angleMatrix);
+
+    // Now, build the descriptor matrix (each row is a point) //
+    descriptors = cv::Mat::zeros(contourMat.cols, descriptorSize(), CV_32F);
+
+    for (int ptidx=0; ptidx<contourMat.cols; ptidx++)
+    {
+        for (int cmp=0; cmp<contourMat.cols; cmp++)
+        {
+            if (ptidx==cmp) continue;
+            if ((int)queryInliers.size()>0)
+            {
+                if (queryInliers[ptidx]==0 || queryInliers[cmp]==0) continue; //avoid outliers
+            }
+
+            int angidx=-1, radidx=-1;
+            for (int i=0; i<nRadialBins; i++)
+            {
+                if (disMatrix.at<float>(ptidx, cmp)<logspaces[i])
+                {
+                    radidx=i;
+                    break;
+                }
+            }
+            for (int i=0; i<nAngularBins; i++)
+            {
+                if (angleMatrix.at<float>(ptidx, cmp)<angspaces[i])
+                {
+                    angidx=i;
+                    break;
+                }
+            }
+            if (angidx!=-1 && radidx!=-1)
+            {
+                int idx = angidx+radidx*nAngularBins;
+                descriptors.at<float>(ptidx, idx)++;
+            }
+        }
+    }
+}
+
+void SCD::logarithmicSpaces(std::vector<double> &vecSpaces) const
+{
+    double logmin=log10(innerRadius);
+    double logmax=log10(outerRadius);
+    double delta=(logmax-logmin)/(nRadialBins-1);
+    double accdelta=0;
+
+    for (int i=0; i<nRadialBins; i++)
+    {
+        double val = std::pow(10,logmin+accdelta);
+        vecSpaces.push_back(val);
+        accdelta += delta;
+    }
+}
+
+void SCD::angularSpaces(std::vector<double> &vecSpaces) const
+{
+    double delta=2*CV_PI/nAngularBins;
+    double val=0;
+
+    for (int i=0; i<nAngularBins; i++)
+    {
+        val += delta;
+        vecSpaces.push_back(val);
+    }
+}
+
+void SCD::buildNormalizedDistanceMatrix(cv::Mat &contour, cv::Mat &disMatrix, const std::vector<int> &queryInliers, const float _meanDistance)
+{
+    cv::Mat contourMat = contour;
+    cv::Mat mask(disMatrix.rows, disMatrix.cols, CV_8U);
+
+    for (int i=0; i<contourMat.cols; i++)
+    {
+      for (int j=0; j<contourMat.cols; j++)
+      {
+          disMatrix.at<float>(i,j) = (float)norm( cv::Mat(contourMat.at<cv::Point2f>(0,i)-contourMat.at<cv::Point2f>(0,j)), cv::NORM_L2 );
+          if (_meanDistance<0)
+          {
+              if (queryInliers.size()>0)
+              {
+                  mask.at<char>(i,j)=char(queryInliers[j] && queryInliers[i]);
+              }
+              else
+              {
+                  mask.at<char>(i,j)=1;
+              }
+          }
+      }
+    }
+
+    if (_meanDistance<0)
+    {
+      meanDistance=(float)mean(disMatrix, mask)[0];
+    }
+    else
+    {
+      meanDistance=_meanDistance;
+    }
+    disMatrix/=meanDistance+FLT_EPSILON;
+}
+
+void SCD::buildAngleMatrix(cv::Mat &contour, cv::Mat &angleMatrix) const
+{
+    cv::Mat contourMat = contour;
+
+    // if descriptor is rotationInvariant compute massCenter //
+    cv::Point2f massCenter(0,0);
+    if (rotationInvariant)
+    {
+        for (int i=0; i<contourMat.cols; i++)
+        {
+            massCenter+=contourMat.at<cv::Point2f>(0,i);
+        }
+        massCenter.x=massCenter.x/(float)contourMat.cols;
+        massCenter.y=massCenter.y/(float)contourMat.cols;
+    }
+
+
+    for (int i=0; i<contourMat.cols; i++)
+    {
+        for (int j=0; j<contourMat.cols; j++)
+        {
+            if (i==j)
+            {
+                angleMatrix.at<float>(i,j)=0.0;
+            }
+            else
+            {
+                cv::Point2f dif = contourMat.at<cv::Point2f>(0,i) - contourMat.at<cv::Point2f>(0,j);
+                angleMatrix.at<float>(i,j) = std::atan2(dif.y, dif.x);
+
+                if (rotationInvariant)
+                {
+                    cv::Point2f refPt = contourMat.at<cv::Point2f>(0,i) - massCenter;
+                    float refAngle = atan2(refPt.y, refPt.x);
+                    angleMatrix.at<float>(i,j) -= refAngle;
+                }
+                angleMatrix.at<float>(i,j) = float(fmod(double(angleMatrix.at<float>(i,j)+(double)FLT_EPSILON),2*CV_PI)+CV_PI);
+            }
+        }
+    }
+}
+
+//! SCDMatcher
+void SCDMatcher::matchDescriptors(cv::Mat &descriptors1, cv::Mat &descriptors2, std::vector<cv::DMatch> &matches,
+                                  cv::Ptr<cv::HistogramCostExtractor> &comparer, std::vector<int> &inliers1, std::vector<int> &inliers2)
+{
+    matches.clear();
+
+    // Build the cost Matrix between descriptors //
+    cv::Mat costMat;
+    buildCostMatrix(descriptors1, descriptors2, costMat, comparer);
+
+    // Solve the matching problem using the hungarian method //
+    hungarian(costMat, matches, inliers1, inliers2, descriptors1.rows, descriptors2.rows);
+}
+
+void SCDMatcher::buildCostMatrix(const cv::Mat &descriptors1, const cv::Mat &descriptors2,
+                                 cv::Mat &costMatrix, cv::Ptr<cv::HistogramCostExtractor> &comparer) const
+{
+    comparer->buildCostMatrix(descriptors1, descriptors2, costMatrix);
+}
+
+void SCDMatcher::hungarian(cv::Mat &costMatrix, std::vector<cv::DMatch> &outMatches, std::vector<int> &inliers1,
+                           std::vector<int> &inliers2, int sizeScd1, int sizeScd2)
+{
+    std::vector<int> free(costMatrix.rows, 0), collist(costMatrix.rows, 0);
+    std::vector<int> matches(costMatrix.rows, 0), colsol(costMatrix.rows), rowsol(costMatrix.rows);
+    std::vector<float> d(costMatrix.rows), pred(costMatrix.rows), v(costMatrix.rows);
+
+    const float LOWV=1e-10;
+    bool unassignedfound;
+    int  i=0, imin=0, numfree=0, prvnumfree=0, f=0, i0=0, k=0, freerow=0;
+    int  j=0, j1=0, j2=0, endofpath=0, last=0, low=0, up=0;
+    float min=0, h=0, umin=0, usubmin=0, v2=0;
+
+    // COLUMN REDUCTION //
+    for (j = costMatrix.rows-1; j >= 0; j--)
+    {
+        // find minimum cost over rows.
+        min = costMatrix.at<float>(0,j);
+        imin = 0;
+        for (i = 1; i < costMatrix.rows; i++)
+        if (costMatrix.at<float>(i,j) < min)
+        {
+            min = costMatrix.at<float>(i,j);
+            imin = i;
+        }
+        v[j] = min;
+
+        if (++matches[imin] == 1)
+        {
+            rowsol[imin] = j;
+            colsol[j] = imin;
+        }
+        else
+        {
+            colsol[j]=-1;
+        }
+    }
+
+    // REDUCTION TRANSFER //
+    for (i=0; i<costMatrix.rows; i++)
+    {
+        if (matches[i] == 0)
+        {
+            free[numfree++] = i;
+        }
+        else
+        {
+            if (matches[i] == 1)
+            {
+                j1=rowsol[i];
+                min=std::numeric_limits<float>::max();
+                for (j=0; j<costMatrix.rows; j++)
+                {
+                    if (j!=j1)
+                    {
+                        if (costMatrix.at<float>(i,j)-v[j] < min)
+                        {
+                            min=costMatrix.at<float>(i,j)-v[j];
+                        }
+                    }
+                }
+                v[j1] = v[j1]-min;
+            }
+        }
+    }
+    // AUGMENTING ROW REDUCTION //
+    int loopcnt = 0;
+    do
+    {
+        loopcnt++;
+        k=0;
+        prvnumfree=numfree;
+        numfree=0;
+        while (k < prvnumfree)
+        {
+            i=free[k];
+            k++;
+            umin = costMatrix.at<float>(i,0)-v[0];
+            j1=0;
+            usubmin = std::numeric_limits<float>::max();
+            for (j=1; j<costMatrix.rows; j++)
+            {
+                h = costMatrix.at<float>(i,j)-v[j];
+                if (h < usubmin)
+                {
+                    if (h >= umin)
+                    {
+                        usubmin = h;
+                        j2 = j;
+                    }
+                    else
+                    {
+                        usubmin = umin;
+                        umin = h;
+                        j2 = j1;
+                        j1 = j;
+                    }
+                }
+            }
+            i0 = colsol[j1];
+
+            if (fabs(umin-usubmin) > LOWV) //if( umin < usubmin )
+            {
+                v[j1] = v[j1] - (usubmin - umin);
+            }
+            else // minimum and subminimum equal.
+            {
+                if (i0 >= 0) // minimum column j1 is assigned.
+                {
+                    j1 = j2;
+                    i0 = colsol[j2];
+                }
+            }
+            // (re-)assign i to j1, possibly de-assigning an i0.
+            rowsol[i]=j1;
+            colsol[j1]=i;
+
+            if (i0 >= 0)
+            {
+                //if( umin < usubmin )
+                if (fabs(umin-usubmin) > LOWV)
+                {
+                    free[--k] = i0;
+                }
+                else
+                {
+                    free[numfree++] = i0;
+                }
+            }
+        }
+    }while (loopcnt<2); // repeat once.
+
+    // AUGMENT SOLUTION for each free row //
+    for (f = 0; f<numfree; f++)
+    {
+        freerow = free[f];       // start row of augmenting path.
+        // Dijkstra shortest path algorithm.
+        // runs until unassigned column added to shortest path tree.
+        for (j = 0; j < costMatrix.rows; j++)
+        {
+            d[j] = costMatrix.at<float>(freerow,j) - v[j];
+            pred[j] = float(freerow);
+            collist[j] = j;        // init column list.
+        }
+
+        low=0; // columns in 0..low-1 are ready, now none.
+        up=0;  // columns in low..up-1 are to be scanned for current minimum, now none.
+        unassignedfound = false;
+        do
+        {
+            if (up == low)
+            {
+                last=low-1;
+                min = d[collist[up++]];
+                for (k = up; k < costMatrix.rows; k++)
+                {
+                    j = collist[k];
+                    h = d[j];
+                    if (h <= min)
+                    {
+                        if (h < min) // new minimum.
+                        {
+                            up = low; // restart list at index low.
+                            min = h;
+                        }
+                        collist[k] = collist[up];
+                        collist[up++] = j;
+                    }
+                }
+                for (k=low; k<up; k++)
+                {
+                    if (colsol[collist[k]] < 0)
+                    {
+                        endofpath = collist[k];
+                        unassignedfound = true;
+                        break;
+                    }
+                }
+            }
+
+            if (!unassignedfound)
+            {
+                // update 'distances' between freerow and all unscanned columns, via next scanned column.
+                j1 = collist[low];
+                low++;
+                i = colsol[j1];
+                h = costMatrix.at<float>(i,j1)-v[j1]-min;
+
+                for (k = up; k < costMatrix.rows; k++)
+                {
+                    j = collist[k];
+                    v2 = costMatrix.at<float>(i,j) - v[j] - h;
+                    if (v2 < d[j])
+                    {
+                        pred[j] = float(i);
+                        if (v2 == min)
+                        {
+                            if (colsol[j] < 0)
+                            {
+                                // if unassigned, shortest augmenting path is complete.
+                                endofpath = j;
+                                unassignedfound = true;
+                                break;
+                            }
+                            else
+                            {
+                                collist[k] = collist[up];
+                                collist[up++] = j;
+                            }
+                        }
+                        d[j] = v2;
+                    }
+                }
+            }
+        }while (!unassignedfound);
+
+        // update column prices.
+        for (k = 0; k <= last; k++)
+        {
+            j1 = collist[k];
+            v[j1] = v[j1] + d[j1] - min;
+        }
+
+        // reset row and column assignments along the alternating path.
+        do
+        {
+            i = int(pred[endofpath]);
+            colsol[endofpath] = i;
+            j1 = endofpath;
+            endofpath = rowsol[i];
+            rowsol[i] = j1;
+        }while (i != freerow);
+    }
+
+    // calculate symmetric shape context cost
+    cv::Mat trueCostMatrix(costMatrix, cv::Rect(0,0,sizeScd1, sizeScd2));
+    float leftcost = 0;
+    for (int nrow=0; nrow<trueCostMatrix.rows; nrow++)
+    {
+        double minval;
+        minMaxIdx(trueCostMatrix.row(nrow), &minval);
+        leftcost+=float(minval);
+    }
+    leftcost /= trueCostMatrix.rows;
+
+    float rightcost = 0;
+    for (int ncol=0; ncol<trueCostMatrix.cols; ncol++)
+    {
+        double minval;
+        minMaxIdx(trueCostMatrix.col(ncol), &minval);
+        rightcost+=float(minval);
+    }
+    rightcost /= trueCostMatrix.cols;
+
+    minMatchCost = std::max(leftcost,rightcost);
+
+    // Save in a DMatch vector
+    for (i=0;i<costMatrix.cols;i++)
+    {
+        cv::DMatch singleMatch(colsol[i],i,costMatrix.at<float>(colsol[i],i));//queryIdx,trainIdx,distance
+        outMatches.push_back(singleMatch);
+    }
+
+    // Update inliers
+    inliers1.reserve(sizeScd1);
+    for (size_t kc = 0; kc<inliers1.size(); kc++)
+    {
+        if (rowsol[kc]<sizeScd1) // if a real match
+            inliers1[kc]=1;
+        else
+            inliers1[kc]=0;
+    }
+    inliers2.reserve(sizeScd2);
+    for (size_t kc = 0; kc<inliers2.size(); kc++)
+    {
+        if (colsol[kc]<sizeScd2) // if a real match
+            inliers2[kc]=1;
+        else
+            inliers2[kc]=0;
+    }
+}
+
+}

modules/shape/src/scd_def.hpp

@@ -0,0 +1,132 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include <stdlib.h>
+#include <math.h>
+#include <vector>
+
+namespace cv
+{
+/*
+ * ShapeContextDescriptor class
+ */
+class SCD
+{
+public:
+    //! the full constructor taking all the necessary parameters
+    explicit SCD(int _nAngularBins=12, int _nRadialBins=5,
+                 double _innerRadius=0.1, double _outerRadius=1, bool _rotationInvariant=false)
+    {
+        setAngularBins(_nAngularBins);
+        setRadialBins(_nRadialBins);
+        setInnerRadius(_innerRadius);
+        setOuterRadius(_outerRadius);
+        setRotationInvariant(_rotationInvariant);
+    }
+
+    void extractSCD(cv::Mat& contour, cv::Mat& descriptors,
+                    const std::vector<int>& queryInliers=std::vector<int>(),
+                    const float _meanDistance=-1);
+
+    int descriptorSize() {return nAngularBins*nRadialBins;}
+    void setAngularBins(int angularBins) { nAngularBins=angularBins; }
+    void setRadialBins(int radialBins) { nRadialBins=radialBins; }
+    void setInnerRadius(double _innerRadius) { innerRadius=_innerRadius; }
+    void setOuterRadius(double _outerRadius) { outerRadius=_outerRadius; }
+    void setRotationInvariant(bool _rotationInvariant) { rotationInvariant=_rotationInvariant; }
+    int getAngularBins() const { return nAngularBins; }
+    int getRadialBins() const { return nRadialBins; }
+    double getInnerRadius() const { return innerRadius; }
+    double getOuterRadius() const { return outerRadius; }
+    bool getRotationInvariant() const { return rotationInvariant; }
+    float getMeanDistance() const { return meanDistance; }
+
+private:
+    int nAngularBins;
+    int nRadialBins;
+    double innerRadius;
+    double outerRadius;
+    bool rotationInvariant;
+    float meanDistance;
+
+protected:
+    void logarithmicSpaces(std::vector<double>& vecSpaces) const;
+    void angularSpaces(std::vector<double>& vecSpaces) const;
+
+    void buildNormalizedDistanceMatrix(cv::Mat& contour,
+                          cv::Mat& disMatrix, const std::vector<int> &queryInliers,
+                          const float _meanDistance=-1);
+
+    void buildAngleMatrix(cv::Mat& contour,
+                              cv::Mat& angleMatrix) const;
+};
+
+/*
+ * Matcher
+ */
+class SCDMatcher
+{
+public:
+    // the full constructor
+    SCDMatcher()
+    {
+    }
+
+    // the matcher function using Hungarian method
+    void matchDescriptors(cv::Mat& descriptors1,  cv::Mat& descriptors2, std::vector<cv::DMatch>& matches, cv::Ptr<cv::HistogramCostExtractor>& comparer,
+                                      std::vector<int>& inliers1, std::vector<int> &inliers2);
+
+    // matching cost
+    float getMatchingCost() const {return minMatchCost;}
+
+private:
+    float minMatchCost;
+    float betaAdditional;
+protected:
+    void buildCostMatrix(const cv::Mat& descriptors1, const cv::Mat& descriptors2,
+                                     cv::Mat& costMatrix, cv::Ptr<cv::HistogramCostExtractor>& comparer) const;
+    void hungarian(cv::Mat& costMatrix, std::vector<cv::DMatch>& outMatches, std::vector<int> &inliers1,
+                   std::vector<int> &inliers2, int sizeScd1=0, int sizeScd2=0);
+
+};
+
+}

modules/shape/src/tps_trans.cpp

@@ -0,0 +1,288 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+namespace cv
+{
+
+class ThinPlateSplineShapeTransformerImpl : public ThinPlateSplineShapeTransformer
+{
+public:
+    /* Constructors */
+    ThinPlateSplineShapeTransformerImpl()
+    {
+        regularizationParameter=0;
+        name_ = "ShapeTransformer.TPS";
+        tpsComputed=false;
+    }
+
+    ThinPlateSplineShapeTransformerImpl(double _regularizationParameter)
+    {
+        regularizationParameter=_regularizationParameter;
+        name_ = "ShapeTransformer.TPS";
+        tpsComputed=false;
+    }
+
+    /* Destructor */
+    ~ThinPlateSplineShapeTransformerImpl()
+    {
+    }
+
+    virtual AlgorithmInfo* info() const { return 0; }
+
+    //! the main operators
+    virtual void estimateTransformation(InputArray transformingShape, InputArray targetShape, std::vector<DMatch> &matches);
+    virtual float applyTransformation(InputArray inPts, OutputArray output=noArray());
+    virtual void warpImage(InputArray transformingImage, OutputArray output,
+                           int flags, int borderMode, const Scalar& borderValue) const;
+
+    //! Setters/Getters
+    virtual void setRegularizationParameter(double _regularizationParameter) {regularizationParameter=_regularizationParameter;}
+    virtual double getRegularizationParameter() const {return regularizationParameter;}
+
+    //! write/read
+    virtual void write(FileStorage& fs) const
+    {
+        fs << "name" << name_
+           << "regularization" << regularizationParameter;
+    }
+
+    virtual void read(const FileNode& fn)
+    {
+        CV_Assert( (String)fn["name"] == name_ );
+        regularizationParameter = (int)fn["regularization"];
+    }
+
+private:
+    bool tpsComputed;
+    double regularizationParameter;
+    float transformCost;
+    Mat tpsParameters;
+    Mat shapeReference;
+
+protected:
+    String name_;
+};
+
+static float distance(Point2f p, Point2f q)
+{
+    Point2f diff = p - q;
+    float norma = diff.x*diff.x + diff.y*diff.y;// - 2*diff.x*diff.y;
+    if (norma<0) norma=0;
+    //else norma = std::sqrt(norma);
+    norma = norma*std::log(norma+FLT_EPSILON);
+    return norma;
+}
+
+static Point2f _applyTransformation(const Mat &shapeRef, const Point2f point, const Mat &tpsParameters)
+{
+    Point2f out;
+    for (int i=0; i<2; i++)
+    {
+        float a1=tpsParameters.at<float>(tpsParameters.rows-3,i);
+        float ax=tpsParameters.at<float>(tpsParameters.rows-2,i);
+        float ay=tpsParameters.at<float>(tpsParameters.rows-1,i);
+
+        float affine=a1+ax*point.x+ay*point.y;
+        float nonrigid=0;
+        for (int j=0; j<shapeRef.rows; j++)
+        {
+            nonrigid+=tpsParameters.at<float>(j,i)*
+                    distance(Point2f(shapeRef.at<float>(j,0),shapeRef.at<float>(j,1)),
+                            point);
+        }
+        if (i==0)
+        {
+            out.x=affine+nonrigid;
+        }
+        if (i==1)
+        {
+            out.y=affine+nonrigid;
+        }
+    }
+    return out;
+}
+
+/* public methods */
+void ThinPlateSplineShapeTransformerImpl::warpImage(InputArray transformingImage, OutputArray output,
+                                      int flags, int borderMode, const Scalar& borderValue) const
+{
+    CV_Assert(tpsComputed==true);
+
+    Mat theinput = transformingImage.getMat();
+    Mat mapX(theinput.rows, theinput.cols, CV_32FC1);
+    Mat mapY(theinput.rows, theinput.cols, CV_32FC1);
+
+    for (int row = 0; row < theinput.rows; row++)
+    {
+        for (int col = 0; col < theinput.cols; col++)
+        {
+            Point2f pt = _applyTransformation(shapeReference, Point2f(float(col), float(row)), tpsParameters);
+            mapX.at<float>(row, col) = pt.x;
+            mapY.at<float>(row, col) = pt.y;
+        }
+    }
+    remap(transformingImage, output, mapX, mapY, flags, borderMode, borderValue);
+}
+
+float ThinPlateSplineShapeTransformerImpl::applyTransformation(InputArray inPts, OutputArray outPts)
+{
+    CV_Assert(tpsComputed);
+    Mat pts1 = inPts.getMat();
+    CV_Assert((pts1.channels()==2) && (pts1.cols>0));
+
+    //Apply transformation in the complete set of points
+    // Ensambling output //
+    if (outPts.needed())
+    {
+        outPts.create(1,pts1.cols, CV_32FC2);
+        Mat outMat = outPts.getMat();
+        for (int i=0; i<pts1.cols; i++)
+        {
+            Point2f pt=pts1.at<Point2f>(0,i);
+            outMat.at<Point2f>(0,i)=_applyTransformation(shapeReference, pt, tpsParameters);
+        }
+    }
+
+    return transformCost;
+}
+
+void ThinPlateSplineShapeTransformerImpl::estimateTransformation(InputArray _pts1, InputArray _pts2,
+                                                               std::vector<DMatch>& _matches )
+{
+    Mat pts1 = _pts1.getMat();
+    Mat pts2 = _pts2.getMat();
+    CV_Assert((pts1.channels()==2) && (pts1.cols>0) && (pts2.channels()==2) && (pts2.cols>0));
+    CV_Assert(_matches.size()>1);
+
+    if (pts1.type() != CV_32F)
+        pts1.convertTo(pts1, CV_32F);
+    if (pts2.type() != CV_32F)
+        pts2.convertTo(pts2, CV_32F);
+
+    // Use only valid matchings //
+    std::vector<DMatch> matches;
+    for (size_t i=0; i<_matches.size(); i++)
+    {
+        if (_matches[i].queryIdx<pts1.cols &&
+            _matches[i].trainIdx<pts2.cols)
+        {
+            matches.push_back(_matches[i]);
+        }
+    }
+
+    // Organizing the correspondent points in matrix style //
+    Mat shape1(matches.size(),2,CV_32F); // transforming shape
+    Mat shape2(matches.size(),2,CV_32F); // target shape
+    for (size_t i=0; i<matches.size(); i++)
+    {
+        Point2f pt1=pts1.at<Point2f>(0,matches[i].queryIdx);
+        shape1.at<float>(i,0) = pt1.x;
+        shape1.at<float>(i,1) = pt1.y;
+
+        Point2f pt2=pts2.at<Point2f>(0,matches[i].trainIdx);
+        shape2.at<float>(i,0) = pt2.x;
+        shape2.at<float>(i,1) = pt2.y;
+    }
+    shape1.copyTo(shapeReference);
+
+    // Building the matrices for solving the L*(w|a)=(v|0) problem with L={[K|P];[P'|0]}
+
+    //Building K and P (Neede to buil L)
+    Mat matK(matches.size(),matches.size(),CV_32F);
+    Mat matP(matches.size(),3,CV_32F);
+    for (size_t i=0; i<matches.size(); i++)
+    {
+        for (size_t j=0; j<matches.size(); j++)
+        {
+            if (i==j)
+            {
+                matK.at<float>(i,j)=float(regularizationParameter);
+            }
+            else
+            {
+                matK.at<float>(i,j) = distance(Point2f(shape1.at<float>(i,0),shape1.at<float>(i,1)),
+                                               Point2f(shape1.at<float>(j,0),shape1.at<float>(j,1)));
+            }
+        }
+        matP.at<float>(i,0) = 1;
+        matP.at<float>(i,1) = shape1.at<float>(i,0);
+        matP.at<float>(i,2) = shape1.at<float>(i,1);
+    }
+
+    //Building L
+    Mat matL=Mat::zeros(matches.size()+3,matches.size()+3,CV_32F);
+    Mat matLroi(matL, Rect(0,0,matches.size(),matches.size())); //roi for K
+    matK.copyTo(matLroi);
+    matLroi = Mat(matL,Rect(matches.size(),0,3,matches.size())); //roi for P
+    matP.copyTo(matLroi);
+    Mat matPt;
+    transpose(matP,matPt);
+    matLroi = Mat(matL,Rect(0,matches.size(),matches.size(),3)); //roi for P'
+    matPt.copyTo(matLroi);
+
+    //Building B (v|0)
+    Mat matB = Mat::zeros(matches.size()+3,2,CV_32F);
+    for (size_t i=0; i<matches.size(); i++)
+    {
+        matB.at<float>(i,0) = shape2.at<float>(i,0); //x's
+        matB.at<float>(i,1) = shape2.at<float>(i,1); //y's
+    }
+
+    //Obtaining transformation params (w|a)
+    solve(matL, matB, tpsParameters, DECOMP_LU);
+    //tpsParameters = matL.inv()*matB;
+
+    //Setting transform Cost and Shape reference
+    Mat w(tpsParameters, Rect(0,0,2,tpsParameters.rows-3));
+    Mat Q=w.t()*matK*w;
+    transformCost=fabs(Q.at<float>(0,0)*Q.at<float>(1,1));//fabs(mean(Q.diag(0))[0]);//std::max(Q.at<float>(0,0),Q.at<float>(1,1));
+    tpsComputed=true;
+}
+
+Ptr <ThinPlateSplineShapeTransformer> createThinPlateSplineShapeTransformer(double regularizationParameter)
+{
+    return Ptr<ThinPlateSplineShapeTransformer>( new ThinPlateSplineShapeTransformerImpl(regularizationParameter) );
+}
+
+} // cv

modules/shape/test/test_emdl1.cpp

@@ -0,0 +1,266 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "test_precomp.hpp"
+
+using namespace cv;
+using namespace std;
+
+const int angularBins=12;
+const int radialBins=4;
+const float minRad=0.2;
+const float maxRad=2;
+const int NSN=5;//10;//20; //number of shapes per class
+const int NP=100; //number of points sympliying the contour
+const float outlierWeight=0.1;
+const int numOutliers=20;
+const float CURRENT_MAX_ACCUR=95; //98% and 99% reached in several tests, 95 is fixed as minimum boundary
+
+class CV_ShapeEMDTest : public cvtest::BaseTest
+{
+public:
+    CV_ShapeEMDTest();
+    ~CV_ShapeEMDTest();
+protected:
+    void run(int);
+
+private:
+    void mpegTest();
+    void listShapeNames(vector<string> &listHeaders);
+    vector<Point2f> convertContourType(const Mat &, int n=0 );
+    float computeShapeDistance(vector <Point2f>& queryNormal,
+                               vector <Point2f>& queryFlipped1,
+                               vector <Point2f>& queryFlipped2,
+                               vector<Point2f>& testq);
+    void displayMPEGResults();
+};
+
+CV_ShapeEMDTest::CV_ShapeEMDTest()
+{
+}
+CV_ShapeEMDTest::~CV_ShapeEMDTest()
+{
+}
+
+vector <Point2f> CV_ShapeEMDTest::convertContourType(const Mat& currentQuery, int n)
+{
+    vector<vector<Point> > _contoursQuery;
+    vector <Point2f> contoursQuery;
+    findContours(currentQuery, _contoursQuery, RETR_LIST, CHAIN_APPROX_NONE);
+    for (size_t border=0; border<_contoursQuery.size(); border++)
+    {
+        for (size_t p=0; p<_contoursQuery[border].size(); p++)
+        {
+            contoursQuery.push_back(Point2f((float)_contoursQuery[border][p].x,
+                                            (float)_contoursQuery[border][p].y));
+        }
+    }
+
+    // In case actual number of points is less than n
+    int dum=0;
+    for (int add=contoursQuery.size()-1; add<n; add++)
+    {
+        contoursQuery.push_back(contoursQuery[dum++]); //adding dummy values
+    }
+
+    // Uniformly sampling
+    random_shuffle(contoursQuery.begin(), contoursQuery.end());
+    int nStart=n;
+    vector<Point2f> cont;
+    for (int i=0; i<nStart; i++)
+    {
+        cont.push_back(contoursQuery[i]);
+    }
+    return cont;
+}
+
+void CV_ShapeEMDTest::listShapeNames( vector<string> &listHeaders)
+{
+    listHeaders.push_back("apple"); //ok
+    listHeaders.push_back("children"); // ok
+    listHeaders.push_back("device7"); // ok
+    listHeaders.push_back("Heart"); // ok
+    listHeaders.push_back("teddy"); // ok
+}
+float CV_ShapeEMDTest::computeShapeDistance(vector <Point2f>& query1, vector <Point2f>& query2,
+                                         vector <Point2f>& query3, vector <Point2f>& testq)
+{
+    //waitKey(0);
+    Ptr <ShapeContextDistanceExtractor> mysc = createShapeContextDistanceExtractor(angularBins, radialBins, minRad, maxRad);
+    //Ptr <HistogramCostExtractor> cost = createNormHistogramCostExtractor(cv::DIST_L1);
+    //Ptr <HistogramCostExtractor> cost = createChiHistogramCostExtractor(30,0.15);
+    //Ptr <HistogramCostExtractor> cost = createEMDHistogramCostExtractor();
+    // Ptr <HistogramCostExtractor> cost = createEMDL1HistogramCostExtractor();
+    mysc->setIterations(1); //(3)
+    mysc->setCostExtractor( createEMDL1HistogramCostExtractor() );
+    //mysc->setTransformAlgorithm(createAffineTransformer(true));
+    mysc->setTransformAlgorithm( createThinPlateSplineShapeTransformer() );
+    //mysc->setImageAppearanceWeight(1.6);
+    //mysc->setImageAppearanceWeight(0.0);
+    //mysc->setImages(im1,imtest);
+    return ( std::min( mysc->computeDistance(query1, testq),
+                       std::min(mysc->computeDistance(query2, testq), mysc->computeDistance(query3, testq) )));
+}
+
+void CV_ShapeEMDTest::mpegTest()
+{
+    string baseTestFolder="shape/mpeg_test/";
+    string path = cvtest::TS::ptr()->get_data_path() + baseTestFolder;
+    vector<string> namesHeaders;
+    listShapeNames(namesHeaders);
+
+    // distance matrix //
+    Mat distanceMat=Mat::zeros(NSN*namesHeaders.size(), NSN*namesHeaders.size(), CV_32F);
+
+    // query contours (normal v flipped, h flipped) and testing contour //
+    vector<Point2f> contoursQuery1, contoursQuery2, contoursQuery3, contoursTesting;
+
+    // reading query and computing its properties //
+    int counter=0;
+    const int loops=NSN*namesHeaders.size()*NSN*namesHeaders.size();
+    for (size_t n=0; n<namesHeaders.size(); n++)
+    {
+        for (int i=1; i<=NSN; i++)
+        {
+            // read current image //
+            stringstream thepathandname;
+            thepathandname<<path+namesHeaders[n]<<"-"<<i<<".png";
+            Mat currentQuery, flippedHQuery, flippedVQuery;
+            currentQuery=imread(thepathandname.str(), IMREAD_GRAYSCALE);
+            Mat currentQueryBuf=currentQuery.clone();
+            flip(currentQuery, flippedHQuery, 0);
+            flip(currentQuery, flippedVQuery, 1);
+            // compute border of the query and its flipped versions //
+            vector<Point2f> origContour;
+            contoursQuery1=convertContourType(currentQuery, NP);
+            origContour=contoursQuery1;
+            contoursQuery2=convertContourType(flippedHQuery, NP);
+            contoursQuery3=convertContourType(flippedVQuery, NP);
+
+            // compare with all the rest of the images: testing //
+            for (size_t nt=0; nt<namesHeaders.size(); nt++)
+            {
+                for (int it=1; it<=NSN; it++)
+                {
+                    // skip self-comparisson //
+                    counter++;
+                    if (nt==n && it==i)
+                    {
+                        distanceMat.at<float>(NSN*n+i-1,
+                                              NSN*nt+it-1)=0;
+                        continue;
+                    }
+                    // read testing image //
+                    stringstream thetestpathandname;
+                    thetestpathandname<<path+namesHeaders[nt]<<"-"<<it<<".png";
+                    Mat currentTest;
+                    currentTest=imread(thetestpathandname.str().c_str(), 0);
+                    // compute border of the testing //
+                    contoursTesting=convertContourType(currentTest, NP);
+
+                    // compute shape distance //
+                    std::cout<<std::endl<<"Progress: "<<counter<<"/"<<loops<<": "<<100*double(counter)/loops<<"% *******"<<std::endl;
+                    std::cout<<"Computing shape distance between "<<namesHeaders[n]<<i<<
+                               " and "<<namesHeaders[nt]<<it<<": ";
+                    distanceMat.at<float>(NSN*n+i-1, NSN*nt+it-1)=
+                            computeShapeDistance(contoursQuery1, contoursQuery2, contoursQuery3, contoursTesting);
+                    std::cout<<distanceMat.at<float>(NSN*n+i-1, NSN*nt+it-1)<<std::endl;
+                }
+            }
+        }
+    }
+    // save distance matrix //
+    FileStorage fs(cvtest::TS::ptr()->get_data_path() + baseTestFolder + "distanceMatrixMPEGTest.yml", FileStorage::WRITE);
+    fs << "distanceMat" << distanceMat;
+}
+
+const int FIRST_MANY=2*NSN;
+void CV_ShapeEMDTest::displayMPEGResults()
+{
+    string baseTestFolder="shape/mpeg_test/";
+    Mat distanceMat;
+    FileStorage fs(cvtest::TS::ptr()->get_data_path() + baseTestFolder + "distanceMatrixMPEGTest.yml", FileStorage::READ);
+    vector<string> namesHeaders;
+    listShapeNames(namesHeaders);
+
+    // Read generated MAT //
+    fs["distanceMat"]>>distanceMat;
+
+    int corrects=0;
+    int divi=0;
+    for (int row=0; row<distanceMat.rows; row++)
+    {
+        if (row%NSN==0) //another group
+        {
+            divi+=NSN;
+        }
+        for (int col=divi-NSN; col<divi; col++)
+        {
+            int nsmall=0;
+            for (int i=0; i<distanceMat.cols; i++)
+            {
+                if (distanceMat.at<float>(row,col)>distanceMat.at<float>(row,i))
+                {
+                    nsmall++;
+                }
+            }
+            if (nsmall<=FIRST_MANY)
+            {
+                corrects++;
+            }
+        }
+    }
+    float porc = 100*float(corrects)/(NSN*distanceMat.rows);
+    std::cout<<"%="<<porc<<std::endl;
+    if (porc >= CURRENT_MAX_ACCUR)
+        ts->set_failed_test_info(cvtest::TS::OK);
+    else
+        ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
+
+}
+
+void CV_ShapeEMDTest::run( int /*start_from*/ )
+{
+    mpegTest();
+    displayMPEGResults();
+}
+
+TEST(ShapeEMD_SCD, regression) { CV_ShapeEMDTest test; test.safe_run(); }

modules/shape/test/test_hausdorff.cpp

@@ -0,0 +1,280 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "test_precomp.hpp"
+#include <stdlib.h>
+
+using namespace cv;
+using namespace std;
+
+const int NSN=5;//10;//20; //number of shapes per class
+const float CURRENT_MAX_ACCUR=85; //90% and 91% reached in several tests, 85 is fixed as minimum boundary
+
+class CV_HaussTest : public cvtest::BaseTest
+{
+public:
+    CV_HaussTest();
+    ~CV_HaussTest();
+protected:
+    void run(int);
+private:
+    float computeShapeDistance(vector<Point> &query1, vector<Point> &query2,
+                               vector<Point> &query3, vector<Point> &testq);
+    vector <Point> convertContourType(const Mat& currentQuery, int n=180);
+    vector<Point2f> normalizeContour(const vector <Point>& contour);
+    void listShapeNames( vector<string> &listHeaders);
+    void mpegTest();
+    void displayMPEGResults();
+};
+
+CV_HaussTest::CV_HaussTest()
+{
+}
+CV_HaussTest::~CV_HaussTest()
+{
+}
+
+vector<Point2f> CV_HaussTest::normalizeContour(const vector<Point> &contour)
+{
+    vector<Point2f> output(contour.size());
+    Mat disMat(contour.size(),contour.size(),CV_32F);
+    Point2f meanpt(0,0);
+    float meanVal=1;
+
+    for (size_t ii=0; ii<contour.size(); ii++)
+    {
+        for (size_t jj=0; jj<contour.size(); jj++)
+        {
+            if (ii==jj) disMat.at<float>(ii,jj)=0;
+            else
+            {
+                disMat.at<float>(ii,jj)=
+                    float(fabs(double(contour[ii].x*contour[jj].x)))+float(fabs(double(contour[ii].y*contour[jj].y)));
+            }
+        }
+        meanpt.x+=contour[ii].x;
+        meanpt.y+=contour[ii].y;
+    }
+    meanpt.x/=contour.size();
+    meanpt.y/=contour.size();
+    meanVal=float(cv::mean(disMat)[0]);
+    for (size_t ii=0; ii<contour.size(); ii++)
+    {
+        output[ii].x = (contour[ii].x-meanpt.x)/meanVal;
+        output[ii].y = (contour[ii].y-meanpt.y)/meanVal;
+    }
+    return output;
+}
+
+void CV_HaussTest::listShapeNames( vector<string> &listHeaders)
+{
+    listHeaders.push_back("apple"); //ok
+    listHeaders.push_back("children"); // ok
+    listHeaders.push_back("device7"); // ok
+    listHeaders.push_back("Heart"); // ok
+    listHeaders.push_back("teddy"); // ok
+}
+
+
+vector <Point> CV_HaussTest::convertContourType(const Mat& currentQuery, int n)
+{
+    vector<vector<Point> > _contoursQuery;
+    vector <Point> contoursQuery;
+    findContours(currentQuery, _contoursQuery, RETR_LIST, CHAIN_APPROX_NONE);
+    for (size_t border=0; border<_contoursQuery.size(); border++)
+    {
+        for (size_t p=0; p<_contoursQuery[border].size(); p++)
+        {
+            contoursQuery.push_back(_contoursQuery[border][p]);
+        }
+    }
+
+    // In case actual number of points is less than n
+    for (int add=contoursQuery.size()-1; add<n; add++)
+    {
+        contoursQuery.push_back(contoursQuery[contoursQuery.size()-add+1]); //adding dummy values
+    }
+
+    // Uniformly sampling
+    random_shuffle(contoursQuery.begin(), contoursQuery.end());
+    int nStart=n;
+    vector<Point> cont;
+    for (int i=0; i<nStart; i++)
+    {
+        cont.push_back(contoursQuery[i]);
+    }
+    return cont;
+}
+
+float CV_HaussTest::computeShapeDistance(vector <Point>& query1, vector <Point>& query2,
+                                         vector <Point>& query3, vector <Point>& testq)
+{
+    Ptr <HausdorffDistanceExtractor> haus = createHausdorffDistanceExtractor();
+    return std::min(haus->computeDistance(query1,testq), std::min(haus->computeDistance(query2,testq),
+                             haus->computeDistance(query3,testq)));
+}
+
+void CV_HaussTest::mpegTest()
+{
+    string baseTestFolder="shape/mpeg_test/";
+    string path = cvtest::TS::ptr()->get_data_path() + baseTestFolder;
+    vector<string> namesHeaders;
+    listShapeNames(namesHeaders);
+
+    // distance matrix //
+    Mat distanceMat=Mat::zeros(NSN*namesHeaders.size(), NSN*namesHeaders.size(), CV_32F);
+
+    // query contours (normal v flipped, h flipped) and testing contour //
+    vector<Point> contoursQuery1, contoursQuery2, contoursQuery3, contoursTesting;
+
+    // reading query and computing its properties //
+    int counter=0;
+    const int loops=NSN*namesHeaders.size()*NSN*namesHeaders.size();
+    for (size_t n=0; n<namesHeaders.size(); n++)
+    {
+        for (int i=1; i<=NSN; i++)
+        {
+            // read current image //
+            stringstream thepathandname;
+            thepathandname<<path+namesHeaders[n]<<"-"<<i<<".png";
+            Mat currentQuery, flippedHQuery, flippedVQuery;
+            currentQuery=imread(thepathandname.str(), IMREAD_GRAYSCALE);
+            flip(currentQuery, flippedHQuery, 0);
+            flip(currentQuery, flippedVQuery, 1);
+            // compute border of the query and its flipped versions //
+            vector<Point> origContour;
+            contoursQuery1=convertContourType(currentQuery);
+            origContour=contoursQuery1;
+            contoursQuery2=convertContourType(flippedHQuery);
+            contoursQuery3=convertContourType(flippedVQuery);
+
+            // compare with all the rest of the images: testing //
+            for (size_t nt=0; nt<namesHeaders.size(); nt++)
+            {
+                for (int it=1; it<=NSN; it++)
+                {
+                    /* skip self-comparisson */
+                    counter++;
+                    if (nt==n && it==i)
+                    {
+                        distanceMat.at<float>(NSN*n+i-1,
+                                              NSN*nt+it-1)=0;
+                        continue;
+                    }
+                    // read testing image //
+                    stringstream thetestpathandname;
+                    thetestpathandname<<path+namesHeaders[nt]<<"-"<<it<<".png";
+                    Mat currentTest;
+                    currentTest=imread(thetestpathandname.str().c_str(), 0);
+
+                    // compute border of the testing //
+                    contoursTesting=convertContourType(currentTest);
+
+                    // compute shape distance //
+                    std::cout<<std::endl<<"Progress: "<<counter<<"/"<<loops<<": "<<100*double(counter)/loops<<"% *******"<<std::endl;
+                    std::cout<<"Computing shape distance between "<<namesHeaders[n]<<i<<
+                               " and "<<namesHeaders[nt]<<it<<": ";
+                    distanceMat.at<float>(NSN*n+i-1, NSN*nt+it-1)=
+                            computeShapeDistance(contoursQuery1, contoursQuery2, contoursQuery3, contoursTesting);
+                    std::cout<<distanceMat.at<float>(NSN*n+i-1, NSN*nt+it-1)<<std::endl;
+                }
+            }
+        }
+    }
+    // save distance matrix //
+    FileStorage fs(cvtest::TS::ptr()->get_data_path() + baseTestFolder + "distanceMatrixMPEGTest.yml", FileStorage::WRITE);
+    fs << "distanceMat" << distanceMat;
+}
+
+const int FIRST_MANY=2*NSN;
+void CV_HaussTest::displayMPEGResults()
+{
+    string baseTestFolder="shape/mpeg_test/";
+    Mat distanceMat;
+    FileStorage fs(cvtest::TS::ptr()->get_data_path() + baseTestFolder + "distanceMatrixMPEGTest.yml", FileStorage::READ);
+    vector<string> namesHeaders;
+    listShapeNames(namesHeaders);
+
+    // Read generated MAT //
+    fs["distanceMat"]>>distanceMat;
+
+    int corrects=0;
+    int divi=0;
+    for (int row=0; row<distanceMat.rows; row++)
+    {
+        if (row%NSN==0) //another group
+        {
+            divi+=NSN;
+        }
+        for (int col=divi-NSN; col<divi; col++)
+        {
+            int nsmall=0;
+            for (int i=0; i<distanceMat.cols; i++)
+            {
+                if (distanceMat.at<float>(row,col)>distanceMat.at<float>(row,i))
+                {
+                    nsmall++;
+                }
+            }
+            if (nsmall<=FIRST_MANY)
+            {
+                corrects++;
+            }
+        }
+    }
+    float porc = 100*float(corrects)/(NSN*distanceMat.rows);
+    std::cout<<"%="<<porc<<std::endl;
+    if (porc >= CURRENT_MAX_ACCUR)
+        ts->set_failed_test_info(cvtest::TS::OK);
+    else
+        ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
+
+}
+
+
+void CV_HaussTest::run(int /* */)
+{
+    mpegTest();
+    displayMPEGResults();
+    ts->set_failed_test_info(cvtest::TS::OK);
+}
+
+TEST(Hauss, regression) { CV_HaussTest test; test.safe_run(); }

modules/shape/test/test_main.cpp

@@ -0,0 +1,3 @@
+#include "test_precomp.hpp"
+
+CV_TEST_MAIN("cv")

modules/shape/test/test_precomp.cpp

@@ -0,0 +1 @@
+#include "test_precomp.hpp"

modules/shape/test/test_precomp.hpp

@@ -0,0 +1,21 @@
+#ifdef __GNUC__
+#  pragma GCC diagnostic ignored "-Wmissing-declarations"
+#  if defined __clang__ || defined __APPLE__
+#    pragma GCC diagnostic ignored "-Wmissing-prototypes"
+#    pragma GCC diagnostic ignored "-Wextra"
+#  endif
+#endif
+
+#ifndef __OPENCV_TEST_PRECOMP_HPP__
+#define __OPENCV_TEST_PRECOMP_HPP__
+
+#include <iostream>
+#include "opencv2/ts.hpp"
+#include "opencv2/core.hpp"
+#include "opencv2/imgproc.hpp"
+#include "opencv2/highgui.hpp"
+#include "opencv2/shape.hpp"
+
+#include "opencv2/opencv_modules.hpp"
+
+#endif

modules/shape/test/test_shape.cpp

@@ -0,0 +1,267 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "test_precomp.hpp"
+
+using namespace cv;
+using namespace std;
+
+const int angularBins=12;
+const int radialBins=4;
+const float minRad=0.2;
+const float maxRad=2;
+const int NSN=5;//10;//20; //number of shapes per class
+const int NP=120; //number of points sympliying the contour
+const float outlierWeight=0.1;
+const int numOutliers=20;
+const float CURRENT_MAX_ACCUR=95.0; //99% and 100% reached in several tests, 95 is fixed as minimum boundary
+
+class CV_ShapeTest : public cvtest::BaseTest
+{
+public:
+    CV_ShapeTest();
+    ~CV_ShapeTest();
+protected:
+    void run(int);
+
+private:
+    void mpegTest();
+    void listShapeNames(vector<string> &listHeaders);
+    vector<Point2f> convertContourType(const Mat &, int n=0 );
+    float computeShapeDistance(vector <Point2f>& queryNormal,
+                               vector <Point2f>& queryFlipped1,
+                               vector <Point2f>& queryFlipped2,
+                               vector<Point2f>& testq);
+    void displayMPEGResults();
+};
+
+CV_ShapeTest::CV_ShapeTest()
+{
+}
+CV_ShapeTest::~CV_ShapeTest()
+{
+}
+
+vector <Point2f> CV_ShapeTest::convertContourType(const Mat& currentQuery, int n)
+{
+    vector<vector<Point> > _contoursQuery;
+    vector <Point2f> contoursQuery;
+    findContours(currentQuery, _contoursQuery, RETR_LIST, CHAIN_APPROX_NONE);
+    for (size_t border=0; border<_contoursQuery.size(); border++)
+    {
+        for (size_t p=0; p<_contoursQuery[border].size(); p++)
+        {
+            contoursQuery.push_back(Point2f((float)_contoursQuery[border][p].x,
+                                            (float)_contoursQuery[border][p].y));
+        }
+    }
+
+    // In case actual number of points is less than n
+    for (int add=contoursQuery.size()-1; add<n; add++)
+    {
+        contoursQuery.push_back(contoursQuery[contoursQuery.size()-add+1]); //adding dummy values
+    }
+
+    // Uniformly sampling
+    random_shuffle(contoursQuery.begin(), contoursQuery.end());
+    int nStart=n;
+    vector<Point2f> cont;
+    for (int i=0; i<nStart; i++)
+    {
+        cont.push_back(contoursQuery[i]);
+    }
+    return cont;
+}
+
+void CV_ShapeTest::listShapeNames( vector<string> &listHeaders)
+{
+    listHeaders.push_back("apple"); //ok
+    listHeaders.push_back("children"); // ok
+    listHeaders.push_back("device7"); // ok
+    listHeaders.push_back("Heart"); // ok
+    listHeaders.push_back("teddy"); // ok
+}
+
+float CV_ShapeTest::computeShapeDistance(vector <Point2f>& query1, vector <Point2f>& query2,
+                                         vector <Point2f>& query3, vector <Point2f>& testq)
+{
+    //waitKey(0);
+    Ptr <ShapeContextDistanceExtractor> mysc = createShapeContextDistanceExtractor(angularBins, radialBins, minRad, maxRad);
+    //Ptr <HistogramCostExtractor> cost = createNormHistogramCostExtractor(cv::DIST_L1);
+    Ptr <HistogramCostExtractor> cost = createChiHistogramCostExtractor(30,0.15);
+    //Ptr <HistogramCostExtractor> cost = createEMDHistogramCostExtractor();
+    //Ptr <HistogramCostExtractor> cost = createEMDL1HistogramCostExtractor();
+    mysc->setIterations(1);
+    mysc->setCostExtractor( cost );
+    //mysc->setTransformAlgorithm(createAffineTransformer(true));
+    mysc->setTransformAlgorithm( createThinPlateSplineShapeTransformer() );
+    //mysc->setImageAppearanceWeight(1.6);
+    //mysc->setImageAppearanceWeight(0.0);
+    //mysc->setImages(im1,imtest);
+    return ( std::min( mysc->computeDistance(query1, testq),
+                       std::min(mysc->computeDistance(query2, testq), mysc->computeDistance(query3, testq) )));
+}
+
+void CV_ShapeTest::mpegTest()
+{
+    string baseTestFolder="shape/mpeg_test/";
+    string path = cvtest::TS::ptr()->get_data_path() + baseTestFolder;
+    vector<string> namesHeaders;
+    listShapeNames(namesHeaders);
+
+    // distance matrix //
+    Mat distanceMat=Mat::zeros(NSN*namesHeaders.size(), NSN*namesHeaders.size(), CV_32F);
+
+    // query contours (normal v flipped, h flipped) and testing contour //
+    vector<Point2f> contoursQuery1, contoursQuery2, contoursQuery3, contoursTesting;
+
+    // reading query and computing its properties //
+    int counter=0;
+    const int loops=NSN*namesHeaders.size()*NSN*namesHeaders.size();
+    for (size_t n=0; n<namesHeaders.size(); n++)
+    {
+        for (int i=1; i<=NSN; i++)
+        {
+            // read current image //
+            stringstream thepathandname;
+            thepathandname<<path+namesHeaders[n]<<"-"<<i<<".png";
+            Mat currentQuery, flippedHQuery, flippedVQuery;
+            currentQuery=imread(thepathandname.str(), IMREAD_GRAYSCALE);
+            Mat currentQueryBuf=currentQuery.clone();
+            flip(currentQuery, flippedHQuery, 0);
+            flip(currentQuery, flippedVQuery, 1);
+            // compute border of the query and its flipped versions //
+            vector<Point2f> origContour;
+            contoursQuery1=convertContourType(currentQuery, NP);
+            origContour=contoursQuery1;
+            contoursQuery2=convertContourType(flippedHQuery, NP);
+            contoursQuery3=convertContourType(flippedVQuery, NP);
+
+            // compare with all the rest of the images: testing //
+            for (size_t nt=0; nt<namesHeaders.size(); nt++)
+            {
+                for (int it=1; it<=NSN; it++)
+                {
+                    // skip self-comparisson //
+                    counter++;
+                    if (nt==n && it==i)
+                    {
+                        distanceMat.at<float>(NSN*n+i-1,
+                                              NSN*nt+it-1)=0;
+                        continue;
+                    }
+                    // read testing image //
+                    stringstream thetestpathandname;
+                    thetestpathandname<<path+namesHeaders[nt]<<"-"<<it<<".png";
+                    Mat currentTest;
+                    currentTest=imread(thetestpathandname.str().c_str(), 0);
+                    // compute border of the testing //
+                    contoursTesting=convertContourType(currentTest, NP);
+
+                    // compute shape distance //
+                    std::cout<<std::endl<<"Progress: "<<counter<<"/"<<loops<<": "<<100*double(counter)/loops<<"% *******"<<std::endl;
+                    std::cout<<"Computing shape distance between "<<namesHeaders[n]<<i<<
+                               " and "<<namesHeaders[nt]<<it<<": ";
+                    distanceMat.at<float>(NSN*n+i-1, NSN*nt+it-1)=
+                            computeShapeDistance(contoursQuery1, contoursQuery2, contoursQuery3, contoursTesting);
+                    std::cout<<distanceMat.at<float>(NSN*n+i-1, NSN*nt+it-1)<<std::endl;
+                }
+            }
+        }
+    }
+    // save distance matrix //
+    FileStorage fs(cvtest::TS::ptr()->get_data_path() + baseTestFolder + "distanceMatrixMPEGTest.yml", FileStorage::WRITE);
+    fs << "distanceMat" << distanceMat;
+}
+
+const int FIRST_MANY=2*NSN;
+void CV_ShapeTest::displayMPEGResults()
+{
+    string baseTestFolder="shape/mpeg_test/";
+    Mat distanceMat;
+    FileStorage fs(cvtest::TS::ptr()->get_data_path() + baseTestFolder + "distanceMatrixMPEGTest.yml", FileStorage::READ);
+    vector<string> namesHeaders;
+    listShapeNames(namesHeaders);
+
+    // Read generated MAT //
+    fs["distanceMat"]>>distanceMat;
+
+    int corrects=0;
+    int divi=0;
+    for (int row=0; row<distanceMat.rows; row++)
+    {
+        if (row%NSN==0) //another group
+        {
+            divi+=NSN;
+        }
+        for (int col=divi-NSN; col<divi; col++)
+        {
+            int nsmall=0;
+            for (int i=0; i<distanceMat.cols; i++)
+            {
+                if (distanceMat.at<float>(row,col)>distanceMat.at<float>(row,i))
+                {
+                    nsmall++;
+                }
+            }
+            if (nsmall<=FIRST_MANY)
+            {
+                corrects++;
+            }
+        }
+    }
+    float porc = 100*float(corrects)/(NSN*distanceMat.rows);
+    std::cout<<"%="<<porc<<std::endl;
+    if (porc >= CURRENT_MAX_ACCUR)
+        ts->set_failed_test_info(cvtest::TS::OK);
+    else
+        ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
+    //done
+}
+
+void CV_ShapeTest::run( int /*start_from*/ )
+{
+    mpegTest();
+    displayMPEGResults();
+    ts->set_failed_test_info(cvtest::TS::OK);
+}
+
+TEST(Shape_SCD, regression) { CV_ShapeTest test; test.safe_run(); }

samples/cpp/CMakeLists.txt

@@ -5,7 +5,7 @@
 
 SET(OPENCV_CPP_SAMPLES_REQUIRED_DEPS opencv_core opencv_flann opencv_imgproc
     opencv_highgui opencv_ml opencv_video opencv_objdetect opencv_photo opencv_nonfree opencv_softcascade
-    opencv_features2d opencv_calib3d opencv_legacy opencv_contrib opencv_stitching opencv_videostab opencv_bioinspired)
+    opencv_features2d opencv_calib3d opencv_legacy opencv_contrib opencv_stitching opencv_videostab opencv_bioinspired opencv_shape)
 
 ocv_check_dependencies(${OPENCV_CPP_SAMPLES_REQUIRED_DEPS})
 

samples/cpp/shape_example.cpp

@@ -0,0 +1,111 @@
+/*
+ * shape_context.cpp -- Shape context demo for shape matching
+ */
+
+#include "opencv2/shape.hpp"
+#include "opencv2/highgui.hpp"
+#include "opencv2/imgproc.hpp"
+#include <opencv2/core/utility.hpp>
+#include <iostream>
+#include <string>
+
+using namespace std;
+using namespace cv;
+
+static void help()
+{
+    printf("\n"
+            "This program demonstrates a method for shape comparisson based on Shape Context\n"
+            "You should run the program providing a number between 1 and 20 for selecting an image in the folder shape_sample.\n"
+            "Call\n"
+            "./shape_example [number between 1 and 20]\n\n");
+}
+
+static vector<Point> simpleContour( const Mat& currentQuery, int n=300 )
+{
+    vector<vector<Point> > _contoursQuery;
+    vector <Point> contoursQuery;
+    findContours(currentQuery, _contoursQuery, RETR_LIST, CHAIN_APPROX_NONE);
+    for (size_t border=0; border<_contoursQuery.size(); border++)
+    {
+        for (size_t p=0; p<_contoursQuery[border].size(); p++)
+        {
+            contoursQuery.push_back( _contoursQuery[border][p] );
+        }
+    }
+
+    // In case actual number of points is less than n
+    int dummy=0;
+    for (int add=contoursQuery.size()-1; add<n; add++)
+    {
+        contoursQuery.push_back(contoursQuery[dummy++]); //adding dummy values
+    }
+
+    // Uniformly sampling
+    random_shuffle(contoursQuery.begin(), contoursQuery.end());
+    vector<Point> cont;
+    for (int i=0; i<n; i++)
+    {
+        cont.push_back(contoursQuery[i]);
+    }
+    return cont;
+}
+
+int main(int argc, char** argv)
+{
+    help();
+    string path = "./shape_sample/";
+    int indexQuery = 1;
+    if( argc < 2 )
+    {
+        std::cout<<"Using first image as query."<<std::endl;
+    }
+    else
+    {
+        sscanf( argv[1], "%i", &indexQuery );
+    }
+    cv::Ptr <cv::ShapeContextDistanceExtractor> mysc = cv::createShapeContextDistanceExtractor();
+
+    Size sz2Sh(300,300);
+    stringstream queryName;
+    queryName<<path<<indexQuery<<".png";
+    Mat query=imread(queryName.str(), IMREAD_GRAYSCALE);
+    Mat queryToShow;
+    resize(query, queryToShow, sz2Sh);
+    imshow("QUERY", queryToShow);
+    moveWindow("TEST", 0,0);
+    vector<Point> contQuery = simpleContour(query);
+    int bestMatch = 0;
+    float bestDis=FLT_MAX;
+    for ( int ii=1; ii<=20; ii++ )
+    {
+        if (ii==indexQuery) continue;
+        waitKey(30);
+        stringstream iiname;
+        iiname<<path<<ii<<".png";
+        cout<<"name: "<<iiname.str()<<endl;
+        Mat iiIm=imread(iiname.str(), 0);
+        Mat iiToShow;
+        resize(iiIm, iiToShow, sz2Sh);
+        imshow("TEST", iiToShow);
+        moveWindow("TEST", sz2Sh.width+50,0);
+        vector<Point> contii = simpleContour(iiIm);
+        float dis = mysc->computeDistance( contQuery, contii );
+        if ( dis<bestDis )
+        {
+            bestMatch = ii;
+            bestDis = dis;
+        }
+        std::cout<<" distance between "<<queryName.str()<<" and "<<iiname.str()<<" is: "<<dis<<std::endl;
+    }
+    destroyWindow("TEST");
+    stringstream bestname;
+    bestname<<path<<bestMatch<<".png";
+    Mat iiIm=imread(bestname.str(), 0);
+    Mat bestToShow;
+    resize(iiIm, bestToShow, sz2Sh);
+    imshow("BEST MATCH", bestToShow);
+    moveWindow("BEST MATCH", sz2Sh.width+50,0);
+
+    return 0;
+}

samples/cpp/shape_transformation.cpp

@@ -0,0 +1,74 @@
+/*
+ * shape_context.cpp -- Shape context demo for shape matching
+ */
+
+#include "opencv2/shape.hpp"
+#include "opencv2/highgui.hpp"
+#include "opencv2/imgproc.hpp"
+#include "opencv2/features2d/features2d.hpp"
+#include "opencv2/nonfree/nonfree.hpp"
+#include <opencv2/core/utility.hpp>
+#include <iostream>
+#include <string>
+
+using namespace std;
+using namespace cv;
+
+static void help()
+{
+    printf("\nThis program demonstrates how to use common interface for shape transformers\n"
+           "Call\n"
+           "shape_transformation [image1] [image2]\n");
+}
+
+int main(int argc, char** argv)
+{
+    help();
+    Mat img1 = imread(argv[1], IMREAD_GRAYSCALE);
+    Mat img2 = imread(argv[2], IMREAD_GRAYSCALE);
+    if(img1.empty() || img2.empty() || argc<2)
+    {
+        printf("Can't read one of the images\n");
+        return -1;
+    }
+
+    // detecting keypoints
+    SurfFeatureDetector detector(5000);
+    vector<KeyPoint> keypoints1, keypoints2;
+    detector.detect(img1, keypoints1);
+    detector.detect(img2, keypoints2);
+
+    // computing descriptors
+    SurfDescriptorExtractor extractor;
+    Mat descriptors1, descriptors2;
+    extractor.compute(img1, keypoints1, descriptors1);
+    extractor.compute(img2, keypoints2, descriptors2);
+
+    // matching descriptors
+    BFMatcher matcher(NORM_L2);
+    vector<DMatch> matches;
+    matcher.match(descriptors1, descriptors2, matches);
+
+    // drawing the results
+    namedWindow("matches", 1);
+    Mat img_matches;
+    drawMatches(img1, keypoints1, img2, keypoints2, matches, img_matches);
+    imshow("matches", img_matches);
+
+    // extract points
+    vector<Point2f> pts1, pts2;
+    for (size_t ii=0; ii<keypoints1.size(); ii++)
+        pts1.push_back( keypoints1[ii].pt );
+    for (size_t ii=0; ii<keypoints2.size(); ii++)
+        pts2.push_back( keypoints2[ii].pt );
+
+    // Apply TPS
+    Ptr<ThinPlateSplineShapeTransformer> mytps = createThinPlateSplineShapeTransformer(25000); //TPS with a relaxed constraint
+    mytps->estimateTransformation(pts1, pts2, matches);
+    mytps->warpImage(img2, img2);
+
+    imshow("Tranformed", img2);
+    waitKey(0);
+
+    return 0;
+}