From f6fc39ce8ff602302b6f2be1890bf1fd96ea5868 Mon Sep 17 00:00:00 2001
From: Juan Manuel Perez <juanmanpr@gmail.com>
Date: Wed, 25 Sep 2013 23:25:10 +0200
Subject: [PATCH] Putting definitions of SCD and SCDMatcher separated from
 sc_dis.cpp file

---
 modules/shape/src/sc_dis.cpp  | 989 ++++++++++++++++++++----------------------
 modules/shape/src/scd_def.hpp | 128 ++++++
 2 files changed, 589 insertions(+), 528 deletions(-)
 create mode 100644 modules/shape/src/scd_def.hpp

diff --git a/modules/shape/src/sc_dis.cpp b/modules/shape/src/sc_dis.cpp
index f19080b..fc5ed42 100644
--- a/modules/shape/src/sc_dis.cpp
+++ b/modules/shape/src/sc_dis.cpp
@@ -46,532 +46,7 @@
  */
 #include "precomp.hpp"
 #include "opencv2/core.hpp"
-/*
- * 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)
-    {
-        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)++;
-                }
-            }
-        }
-    }
-
-    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
-    {
-        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 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 buildNormalizedDistanceMatrix(cv::Mat& contour,
-                          cv::Mat& disMatrix, const std::vector<int> &queryInliers,
-                          const float _meanDistance=-1)
-    {
-        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 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);
-                      //angleMatrix.at<float>(i,j) = 1+floor( angleMatrix.at<float>(i,j)*nAngularBins/(2*CV_PI) );
-                  }
-              }
-          }
-      }
-};
-
-/*
- * 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)
-    {
-        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);
-    }
-
-    // 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
-    {
-        comparer->buildCostMatrix(descriptors1, descriptors2, costMatrix);
-    }
-
-    void hungarian(cv::Mat& costMatrix, std::vector<cv::DMatch>& outMatches, std::vector<int> &inliers1,
-                   std::vector<int> &inliers2, int sizeScd1=0, int sizeScd2=0)
-    {
-        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;
-        }
-    }
-
-};
-
-/*
- *
- */
+#include "scd_def.hpp"
 
 namespace cv
 {
@@ -734,9 +209,9 @@ float ShapeContextDistanceExtractorImpl::computeDistance(InputArray contour1, In
     }
 
     // Initializing Extractor, Descriptor structures and Matcher //
-    SCD set1SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, false);
+    SCD set1SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, rotationInvariant);
     Mat set1SCD;
-    SCD set2SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, false);
+    SCD set2SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, rotationInvariant);
     Mat set2SCD;
     SCDMatcher matcher;
     std::vector<DMatch> matches;
@@ -846,3 +321,461 @@ Ptr <ShapeContextDistanceExtractor> createShapeContextDistanceExtractor(int nAng
 }
 
 } // cv
+
+//! 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;
+    }
+}
diff --git a/modules/shape/src/scd_def.hpp b/modules/shape/src/scd_def.hpp
new file mode 100644
index 0000000..c7476ea
--- /dev/null
+++ b/modules/shape/src/scd_def.hpp
@@ -0,0 +1,128 @@
+/*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>
+
+/*
+ * 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);
+
+};
-- 
2.7.4