Putting definitions of SCD and SCDMatcher separated from sc_dis.cpp file

modules/shape/src/sc_dis.cpp

@@ -46,165 +46,398 @@
  */
 #include "precomp.hpp"
 #include "opencv2/core.hpp"
-/*
- * ShapeContextDescriptor class
- */
-class SCD
+#include "scd_def.hpp"
+
+namespace cv
+{
+class ShapeContextDistanceExtractorImpl : public ShapeContextDistanceExtractor
 {
 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)
+    /* Constructors */
+    ShapeContextDistanceExtractorImpl(int _nAngularBins, int _nRadialBins, float _innerRadius, float _outerRadius, int _iterations,
+                                      const Ptr<HistogramCostExtractor> &_comparer, const Ptr<ShapeTransformer> &_transformer)
     {
-        setAngularBins(_nAngularBins);
-        setRadialBins(_nRadialBins);
-        setInnerRadius(_innerRadius);
-        setOuterRadius(_outerRadius);
-        setRotationInvariant(_rotationInvariant);
+        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";
     }
 
-    void extractSCD(cv::Mat& contour, cv::Mat& descriptors,
-                    const std::vector<int>& queryInliers=std::vector<int>(),
-                    const float _meanDistance=-1)
+    /* Destructor */
+    ~ShapeContextDistanceExtractorImpl()
     {
-        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);
+    virtual AlgorithmInfo* info() const { return 0; }
 
-        // Now, build the descriptor matrix (each row is a point) //
-        descriptors = cv::Mat::zeros(contourMat.cols, descriptorSize(), CV_32F);
+    //! the main operator
+    virtual float computeDistance(InputArray contour1, InputArray contour2);
 
-        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)
+    //! 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)
     {
-                    if (queryInliers[ptidx]==0 || queryInliers[cmp]==0) continue; //avoid outliers
+        Mat image1_=_image1.getMat(), image2_=_image2.getMat();
+        CV_Assert((image1_.depth()==0) && (image2_.depth()==0));
+        image1=image1_;
+        image2=image2_;
     }
 
-                int angidx=-1, radidx=-1;
-                for (int i=0; i<nRadialBins; i++)
-                {
-                    if (disMatrix.at<float>(ptidx, cmp)<logspaces[i])
+    virtual void getImages(OutputArray _image1, OutputArray _image2) const
     {
-                        radidx=i;
-                        break;
-                    }
+        CV_Assert((!image1.empty()) && (!image2.empty()));
+        _image1.create(image1.size(), image1.type());
+        _image2.create(image2.size(), image2.type());
+        _image1.getMat()=image1;
+        _image2.getMat()=image2;
     }
-                for (int i=0; i<nAngularBins; i++)
-                {
-                    if (angleMatrix.at<float>(ptidx, cmp)<angspaces[i])
+
+    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
     {
-                        angidx=i;
-                        break;
-                    }
+        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;
     }
-                if (angidx!=-1 && radidx!=-1)
+
+    virtual void read(const FileNode& fn)
     {
-                    int idx = angidx+radidx*nAngularBins;
-                    descriptors.at<float>(ptidx, idx)++;
-                }
-            }
-        }
+        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"];
     }
 
-    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;
+    float innerRadius;
+    float outerRadius;
     bool rotationInvariant;
-    float meanDistance;
+    int costFlag;
+    int iterations;
+    Ptr<ShapeTransformer> transformer;
+    Ptr<HistogramCostExtractor> comparer;
+    Mat image1;
+    Mat image2;
+    float bendingEnergyWeight;
+    float imageAppearanceWeight;
+    float shapeContextWeight;
+    float sigma;
 
 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;
+    String name_;
+};
 
-        for (int i=0; i<nRadialBins; i++)
-        {
-            double val = std::pow(10,logmin+accdelta);
-            vecSpaces.push_back(val);
-            accdelta += delta;
-        }
-    }
+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);
 
-    void angularSpaces(std::vector<double>& vecSpaces) const
-    {
-        double delta=2*CV_PI/nAngularBins;
-        double val=0;
+    if (set2.type() != CV_32F)
+        sset2.convertTo(set2, CV_32F);
+    else
+        sset1.copyTo(set2);
 
-        for (int i=0; i<nAngularBins; i++)
+    CV_Assert((set1.channels()==2) && (set1.cols>0));
+    CV_Assert((set2.channels()==2) && (set2.cols>0));
+    if (imageAppearanceWeight!=0)
     {
-            val += delta;
-            vecSpaces.push_back(val);
-        }
+        CV_Assert((!image1.empty()) && (!image2.empty()));
     }
 
-    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);
+    // 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;
 
-        for (int i=0; i<contourMat.cols; i++)
+    // 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;
+    for (int ii=0; ii<iterations; ii++)
     {
-          for (int j=0; j<contourMat.cols; j++)
+        // 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=std::pow(set1SCE.getMeanDistance(),2);
+
+        // 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)
         {
-              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)
+            // Have to accumulate the transformation along all the iterations
+            if (ii==0)
             {
-                  if (queryInliers.size()>0)
+                if ( !transDown.empty() )
                 {
-                      mask.at<char>(i,j)=char(queryInliers[j] && queryInliers[i]);
+                    image2.copyTo(warpedImage);
                 }
                 else
                 {
-                      mask.at<char>(i,j)=1;
+                    image1.copyTo(warpedImage);
                 }
             }
+            transformer->warpImage(warpedImage, warpedImage);
         }
     }
 
-        if (_meanDistance<0)
+    Mat gaussWindow, diffIm;
+    if (imageAppearanceWeight!=0)
     {
-          meanDistance=(float)mean(disMatrix, mask)[0];
+        // compute appearance cost
+        if ( !transDown.empty() )
+        {
+            resize(warpedImage, warpedImage, image1.size());
+            Mat temp=(warpedImage-image1);
+            multiply(temp, temp, diffIm);
         }
         else
         {
-          meanDistance=_meanDistance;
-        }
-        disMatrix/=meanDistance+FLT_EPSILON;
+            resize(warpedImage, warpedImage, image2.size());
+            Mat temp=(warpedImage-image2);
+            multiply(temp, temp, diffIm);
+        }
+        gaussWindow = Mat::zeros(warpedImage.rows, warpedImage.cols, CV_32F);
+        for (int pt=0; pt<sset1.cols; pt++)
+        {
+            for (int ii=0; ii<diffIm.rows; ii++)
+            {
+                for (int jj=0; jj<diffIm.cols; jj++)
+                {
+                    float xx = sset1.at<Point2f>(0,pt).x;
+                    float yy = sset1.at<Point2f>(0,pt).y;
+                    float val = float(std::exp( -float( (xx-jj)*(xx-jj) + (yy-ii)*(yy-ii) )/(2*sigma*sigma) ) / (sigma*sigma*2*CV_PI));
+                    gaussWindow.at<float>(ii,jj) += val;
+                }
+            }
+        }
+
+        Mat appIm(diffIm.rows, diffIm.cols, CV_32F);
+        for (int ii=0; ii<diffIm.rows; ii++)
+        {
+            for (int 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) );
+}
+
+} // 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;
+              }
+          }
+      }
     }
 
-    void buildAngleMatrix(cv::Mat& contour,
-                              cv::Mat& angleMatrix) const
+    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 //
@@ -240,28 +473,15 @@ protected:
                     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
+//! 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)
 {
-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 //
@@ -270,24 +490,17 @@ public:
 
     // 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
-    {
+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 hungarian(cv::Mat& costMatrix, std::vector<cv::DMatch>& outMatches, std::vector<int> &inliers1,
-                   std::vector<int> &inliers2, int sizeScd1=0, int sizeScd2=0)
-    {
+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);
@@ -565,284 +778,4 @@ protected:
         else
             inliers2[kc]=0;
     }
-    }
-
-};
-
-/*
- *
- */
-
-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"];
-    }
-
-private:
-    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;
-
-protected:
-    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, false);
-    Mat set1SCD;
-    SCD set2SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, false);
-    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;
-    for (int 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=std::pow(set1SCE.getMeanDistance(),2);
-
-        // 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 (int pt=0; pt<sset1.cols; pt++)
-        {
-            for (int ii=0; ii<diffIm.rows; ii++)
-            {
-                for (int jj=0; jj<diffIm.cols; jj++)
-                {
-                    float xx = sset1.at<Point2f>(0,pt).x;
-                    float yy = sset1.at<Point2f>(0,pt).y;
-                    float val = float(std::exp( -float( (xx-jj)*(xx-jj) + (yy-ii)*(yy-ii) )/(2*sigma*sigma) ) / (sigma*sigma*2*CV_PI));
-                    gaussWindow.at<float>(ii,jj) += val;
-                }
-            }
-        }
-
-        Mat appIm(diffIm.rows, diffIm.cols, CV_32F);
-        for (int ii=0; ii<diffIm.rows; ii++)
-        {
-            for (int 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) );
-}
-
-} // cv

modules/shape/src/scd_def.hpp

+/*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);
+
+};