Merge pull request #1216 from LaurentBerger:FourierDescriptors

modules/ximgproc/doc/ximgproc.bib

@@ -76,9 +76,9 @@
   publisher={Springer}
 }
 
-@article{deriche1987using,
+@article{deriche1987,
   title={Using Canny's criteria to derive a recursively implemented optimal edge detector},
-  author={Deriche, Rachid},
+  author={Deriche Rachid},
   journal={International journal of computer vision},
   volume={1},
   number={2},
@@ -259,3 +259,25 @@
   year={2009},
   organization={International Society for Optics and Photonics}
 }
+
+@article{BergerRaghunathan1998,
+  title={Coalescence in 2 dimensions: experiments on thin copolymer films and numerical simulations},
+  author={Berger, L and Raghunathan, V A and Launay, C and  Ausserré, D and Gallot, Y},
+  journal={The European Physical Journal B - Condensed Matter and Complex Systems},
+  volume={2},
+  number={1},
+  pages={93-99},
+  year={1998},
+  publisher={Springer}
+}
+
+@article{PersoonFu1977,
+  title={Shape Discrimination Using Fourier Descriptors},
+  author={E  Persoon and King-Sun  Fu},
+  journal={IEEE Transactions on Pattern Analysis and Machine Intelligence},
+  volume={7},
+  number={3},
+  pages={170-179},
+  year={1977},
+  publisher={IEEE Computer Society}
+}

modules/ximgproc/include/opencv2/ximgproc.hpp

@@ -52,6 +52,8 @@
 #include "ximgproc/fast_line_detector.hpp"
 #include "ximgproc/deriche_filter.hpp"
 #include "ximgproc/peilin.hpp"
+#include "ximgproc/fourier_descriptors.hpp"
+
 
 /** @defgroup ximgproc Extended Image Processing
   @{
@@ -67,7 +69,9 @@ i.e. algorithms which somehow takes into account pixel affinities in natural ima
     @defgroup ximgproc_segmentation Image segmentation
 
     @defgroup ximgproc_fast_line_detector Fast line detector
-  @}
+
+    @defgroup ximgproc_fourier Fourier descriptors
+    @}
 */
 
 namespace cv

modules/ximgproc/include/opencv2/ximgproc/fourier_descriptors.hpp

+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#ifndef __OPENCV_FOURIERDESCRIPTORS_HPP__
+#define __OPENCV_FOURIERDESCRIPTORS_HPP__
+
+#include <opencv2/core.hpp>
+
+namespace cv {
+namespace ximgproc {
+
+    //! @addtogroup ximgproc_fourier
+    //! @{
+
+    /** @brief Class for ContourFitting algorithms.
+    ContourFitting match two contours \f$ z_a \f$ and \f$ z_b \f$ minimizing distance
+    \f[ d(z_a,z_b)=\sum (a_n - s  b_n e^{j(n \alpha +\phi )})^2 \f] where \f$ a_n \f$ and \f$ b_n \f$ are Fourier descriptors of \f$ z_a \f$ and \f$ z_b \f$ and s is a scaling factor and  \f$ \phi \f$ is angle rotation and \f$ \alpha \f$ is starting point factor adjustement
+    */
+    class CV_EXPORTS_W ContourFitting : public Algorithm
+    {
+        int ctrSize;
+        int fdSize;
+        std::vector<std::complex<double> > b;
+        std::vector<std::complex<double> > a;
+        std::vector<double> frequence;
+        std::vector<double> rho, psi;
+        void frequencyInit();
+        void fAlpha(double x, double &fn, double &df);
+        double distance(std::complex<double> r, double alpha);
+        double  newtonRaphson(double x1, double x2);
+    public:
+        /** @brief Fit two closed curves using fourier descriptors. More details in @cite PersoonFu1977 and @cite BergerRaghunathan1998
+
+        * @param ctr number of Fourier descriptors equal to number of contour points after resampling.
+        * @param fd Contour defining second shape (Target).
+        */
+        ContourFitting(int ctr=1024,int fd=16):ctrSize(ctr),fdSize(fd){};
+        /** @brief Fit two closed curves using fourier descriptors. More details in @cite PersoonFu1977 and @cite BergerRaghunathan1998
+
+        @param src Contour defining first shape.
+        @param dst Contour defining second shape (Target).
+        @param alphaPhiST : \f$ \alpha \f$=alphaPhiST(0,0), \f$ \phi \f$=alphaPhiST(0,1) (in radian), s=alphaPhiST(0,2), Tx=alphaPhiST(0,3), Ty=alphaPhiST(0,4) rotation center
+        @param dist distance between src and dst after matching.
+        @param fdContour false then src and dst are contours and true src and dst are fourier descriptors.
+        */
+        CV_WRAP void estimateTransformation(InputArray src, InputArray dst, OutputArray alphaPhiST, double *dist = 0, bool fdContour = false);
+        /** @brief Fit two closed curves using fourier descriptors. More details in @cite PersoonFu1977 and @cite BergerRaghunathan1998
+
+        @param src Contour defining first shape.
+        @param dst Contour defining second shape (Target).
+        @param alphaPhiST : \f$ \alpha \f$=alphaPhiST(0,0), \f$ \phi \f$=alphaPhiST(0,1) (in radian), s=alphaPhiST(0,2), Tx=alphaPhiST(0,3), Ty=alphaPhiST(0,4) rotation center
+        @param dist distance between src and dst after matching.
+        @param fdContour false then src and dst are contours and true src and dst are fourier descriptors.
+        */
+        CV_WRAP void estimateTransformation(InputArray src, InputArray dst, OutputArray alphaPhiST, double &dist , bool fdContour = false);
+        /** @brief set number of Fourier descriptors used in estimateTransformation
+
+        @param n number of Fourier descriptors equal to number of contour points after resampling.
+        */
+        CV_WRAP void setCtrSize(int n);
+        /** @brief set number of Fourier descriptors when estimateTransformation used vector<Point>
+
+        @param n number of fourier descriptors used for optimal curve matching.
+        */
+        CV_WRAP void setFDSize(int n);
+        /**
+        @returns number of fourier descriptors
+        */
+        CV_WRAP int getCtrSize() { return ctrSize; };
+        /**
+        @returns number of fourier descriptors used for optimal curve matching
+        */
+        CV_WRAP int getFDSize() { return fdSize; };
+    };
+    /**
+    * @brief   Fourier descriptors for planed closed curves
+    *
+    * For more details about this implementation, please see @cite PersoonFu1977
+    *
+    * @param   src   contour type vector<Point> , vector<Point2f>  or vector<Point2d>
+    * @param   dst   Mat of type CV_64FC2 and nbElt rows A VERIFIER
+    * @param   nbElt number of rows in dst or getOptimalDFTSize rows if nbElt=-1
+    * @param   nbFD number of FD return in dst dst = [FD(1...nbFD/2) FD(nbFD/2-nbElt+1...:nbElt)]
+    *
+    */
+    CV_EXPORTS_W void fourierDescriptor(InputArray src, OutputArray dst, int nbElt=-1,int nbFD=-1);
+    /**
+    * @brief   transform a contour
+    *
+    * @param   src   contour or Fourier Descriptors if fd is true
+    * @param   t   transform Mat given by estimateTransformation
+    * @param   dst   Mat of type CV_64FC2 and nbElt rows
+    * @param   fdContour true src are Fourier Descriptors. fdContour false src is a contour
+    *
+    */
+    CV_EXPORTS_W void transform(InputArray src, InputArray t,OutputArray dst, bool fdContour=true);
+    /**
+    * @brief   Contour sampling .
+    *
+    * @param   src   contour type vector<Point> , vector<Point2f>  or vector<Point2d>
+    * @param   out   Mat of type CV_64FC2 and nbElt rows
+    * @param   nbElt number of points in out contour
+    *
+    */
+    CV_EXPORTS_W void contourSampling(InputArray src, OutputArray out, int nbElt);
+
+    /**
+    * @brief create
+
+    * @param ctr number of Fourier descriptors equal to number of contour points after resampling.
+    * @param fd Contour defining second shape (Target).
+    */
+    CV_EXPORTS_W Ptr<ContourFitting> create(int ctr = 1024, int fd = 16);
+
+    //! @} ximgproc_fourier
+}
+}
+#endif

modules/ximgproc/samples/fourier_descriptors_demo.cpp

+#include <opencv2/core.hpp>
+#include <opencv2/core/utility.hpp>
+#include <opencv2/highgui.hpp>
+#include <opencv2/imgproc.hpp>
+#include <opencv2/ximgproc.hpp>
+#include <iostream>
+
+using namespace cv;
+using namespace std;
+
+struct ThParameters {
+    int levelNoise;
+    int angle;
+    int scale10;
+    int origin;
+    int xg;
+    int yg;
+    bool update;
+} ;
+
+static vector<Point> NoisyPolygon(vector<Point> pRef, double n);
+static void UpdateShape(int , void *r);
+static void AddSlider(String sliderName, String windowName, int minSlider, int maxSlider, int valDefault, int *valSlider, void(*f)(int, void *), void *r);
+
+int main(void)
+{
+    vector<Point> ctrRef;
+    vector<Point> ctrRotate, ctrNoisy, ctrNoisyRotate, ctrNoisyRotateShift;
+    // build a shape with 5 vertex
+    ctrRef.push_back(Point(250,250)); ctrRef.push_back(Point(400, 250));
+    ctrRef.push_back(Point(400, 300)); ctrRef.push_back(Point(250, 300));ctrRef.push_back(Point(180, 270));
+    Point cg(0,0);
+    for (int i=0;i<static_cast<int>(ctrRef.size());i++)
+        cg+=ctrRef[i];
+    cg.x /= static_cast<int>(ctrRef.size());
+    cg.y /= static_cast<int>(ctrRef.size());
+    ThParameters p;
+    p.levelNoise=6;
+    p.angle=45;
+    p.scale10=5;
+    p.origin=10;
+    p.xg=150;
+    p.yg=150;
+    p.update=true;
+    namedWindow("FD Curve matching");
+    // A rotation with center at (150,150) of angle 45 degrees and a scaling of 5/10
+    AddSlider("Noise", "FD Curve matching", 0, 20, p.levelNoise, &p.levelNoise, UpdateShape, &p);
+    AddSlider("Angle", "FD Curve matching", 0, 359, p.angle, &p.angle, UpdateShape, &p);
+    AddSlider("Scale", "FD Curve matching", 5, 100, p.scale10, &p.scale10, UpdateShape, &p);
+    AddSlider("Origin%%", "FD Curve matching", 0, 100, p.origin, &p.origin, UpdateShape, &p);
+    AddSlider("Xg", "FD Curve matching", 150, 450, p.xg, &p.xg, UpdateShape, &p);
+    AddSlider("Yg", "FD Curve matching", 150, 450, p.yg, &p.yg, UpdateShape, &p);
+    int code=0;
+    double dist;
+    vector<vector<Point> > c;
+    Mat img;
+    cout << "******************** PRESS g TO MATCH CURVES *************\n";
+    do
+    {
+        code = waitKey(30);
+        if (p.update)
+        {
+            Mat r = getRotationMatrix2D(Point(p.xg, p.yg), p.angle, 10.0/ p.scale10);
+            ctrNoisy= NoisyPolygon(ctrRef,static_cast<double>(p.levelNoise));
+            transform(ctrNoisy, ctrNoisyRotate,r);
+            ctrNoisyRotateShift.clear();
+            for (int i=0;i<static_cast<int>(ctrNoisy.size());i++)
+                ctrNoisyRotateShift.push_back(ctrNoisyRotate[(i+(p.origin*ctrNoisy.size())/100)% ctrNoisy.size()]);
+            // To draw contour using drawcontours
+            c.clear();
+            c.push_back(ctrRef);
+            c.push_back(ctrNoisyRotateShift);
+            p.update = false;
+            Rect rglobal;
+            for (int i = 0; i < static_cast<int>(c.size()); i++)
+            {
+                rglobal = boundingRect(c[i]) | rglobal;
+            }
+            rglobal.width += 10;
+            rglobal.height += 10;
+            img = Mat::zeros(2 * rglobal.height, 2 * rglobal.width, CV_8UC(3));
+            drawContours(img, c, 0, Scalar(255,0,0));
+            drawContours(img, c, 1, Scalar(0, 255, 0));
+            circle(img, c[0][0], 5, Scalar(255, 0, 0));
+            circle(img, c[1][0], 5, Scalar(0, 255, 0));
+            imshow("FD Curve matching", img);
+        }
+        if (code == 'd')
+        {
+            destroyWindow("FD Curve matching");
+            namedWindow("FD Curve matching");
+            // A rotation with center at (150,150) of angle 45 degrees and a scaling of 5/10
+            AddSlider("Noise", "FD Curve matching", 0, 20, p.levelNoise, &p.levelNoise, UpdateShape, &p);
+            AddSlider("Angle", "FD Curve matching", 0, 359, p.angle, &p.angle, UpdateShape, &p);
+            AddSlider("Scale", "FD Curve matching", 5, 100, p.scale10, &p.scale10, UpdateShape, &p);
+            AddSlider("Origin%%", "FD Curve matching", 0, 100, p.origin, &p.origin, UpdateShape, &p);
+            AddSlider("Xg", "FD Curve matching", 150, 450, p.xg, &p.xg, UpdateShape, &p);
+            AddSlider("Yg", "FD Curve matching", 150, 450, p.yg, &p.yg, UpdateShape, &p);
+
+        }
+        if (code == 'g')
+        {
+            ximgproc::ContourFitting fit;
+            vector<Point2f> ctrRef2d, ctrRot2d;
+            // sampling contour we want 256 points
+            ximgproc::contourSampling(ctrRef, ctrRef2d, 256); // use a mat
+            ximgproc::contourSampling(ctrNoisyRotateShift, ctrRot2d, 256); // use a vector og point
+            fit.setFDSize(16);
+            Mat t;
+            fit.estimateTransformation(ctrRot2d, ctrRef2d, t, &dist, false);
+            cout << "Transform *********\n "<<"Origin = "<< t.at<double>(0,0)*ctrNoisy.size() <<" expected "<< (p.origin*ctrNoisy.size()) / 100 <<" ("<< ctrNoisy.size()<<")\n";
+            cout << "Angle = " << t.at<double>(0, 1) * 180 / M_PI << " expected " << p.angle  <<"\n";
+            cout << "Scale = " << t.at<double>(0, 2) << " expected " << p.scale10 / 10.0 << "\n";
+            Mat dst;
+            ximgproc::transform(ctrRot2d, t, dst, false);
+            c.push_back(dst);
+            drawContours(img, c, 2, Scalar(0,255,255));
+            circle(img, c[2][0], 5, Scalar(0, 255, 255));
+            imshow("FD Curve matching", img);
+        }
+    }
+    while (code!=27);
+
+    return 0;
+}
+
+vector<Point> NoisyPolygon(vector<Point> pRef, double n)
+{
+    RNG rng;
+    vector<Point> c;
+    vector<Point> p = pRef;
+    vector<vector<Point> > contour;
+    for (int i = 0; i<static_cast<int>(p.size()); i++)
+        p[i] += Point(Point2d(n*rng.uniform((double)-1, (double)1), n*rng.uniform((double)-1, (double)1)));
+    if (n==0)
+        return p;
+    c.push_back(p[0]);
+    int minX = p[0].x, maxX = p[0].x, minY = p[0].y, maxY = p[0].y;
+    for (int i = 0; i <static_cast<int>(p.size()); i++)
+    {
+        int next = i + 1;
+        if (next == static_cast<int>(p.size()))
+            next = 0;
+        Point2d u = p[next] - p[i];
+        int d = static_cast<int>(norm(u));
+        double a = atan2(u.y, u.x);
+        int step = 1;
+        if (n != 0)
+            step = static_cast<int>(d / n);
+        for (int j = 1; j<d; j += max(step, 1))
+        {
+            Point pNew;
+            do
+            {
+
+                Point2d pAct = (u*j) / static_cast<double>(d);
+                double r = n*rng.uniform((double)0, (double)1);
+                double theta = a + rng.uniform(0., 2 * CV_PI);
+                pNew = Point(Point2d(r*cos(theta) + pAct.x + p[i].x, r*sin(theta) + pAct.y + p[i].y));
+            } while (pNew.x<0 || pNew.y<0);
+            if (pNew.x<minX)
+                minX = pNew.x;
+            if (pNew.x>maxX)
+                maxX = pNew.x;
+            if (pNew.y<minY)
+                minY = pNew.y;
+            if (pNew.y>maxY)
+                maxY = pNew.y;
+            c.push_back(pNew);
+        }
+    }
+    return c;
+}
+
+void UpdateShape(int , void *r)
+{
+    ((ThParameters *)r)->update = true;
+}
+
+void AddSlider(String sliderName, String windowName, int minSlider, int maxSlider, int valDefault, int *valSlider, void(*f)(int, void *), void *r)
+{
+    createTrackbar(sliderName, windowName, valSlider, 1, f, r);
+    setTrackbarMin(sliderName, windowName, minSlider);
+    setTrackbarMax(sliderName, windowName, maxSlider);
+    setTrackbarPos(sliderName, windowName, valDefault);
+}

modules/ximgproc/src/fourier_descriptors.cpp

+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#include "precomp.hpp"
+#include <math.h>
+#include <vector>
+#include <iostream>
+
+/*
+If you use this code please cite this @cite BergerRaghunathan1998
+Coalescence in 2 dimensions: experiments on thin copolymer films and numerical simulations  The European Physical Journal B - Condensed Matter and Complex Systems  (Volume:2 ,  Issue: 1 )  1998
+*/
+
+namespace cv {
+namespace ximgproc {
+
+
+void ContourFitting::setCtrSize(int n)
+{
+    CV_Assert(n>0);
+    ctrSize = n;
+}
+
+void ContourFitting::setFDSize(int n)
+{
+    CV_Assert(n>0);
+    fdSize=n;
+}
+
+void ContourFitting::frequencyInit()
+{
+    int nbElt = ctrSize;
+    frequence.resize(ctrSize);
+
+    for (int i = 0; i <= nbElt / 2; i++)
+        frequence[i] = 2 * M_PI*(float)i / nbElt;
+    for (int i = nbElt / 2 + 1; i<nbElt; i++)
+        frequence[i] = 2 * M_PI*(float)(i - nbElt) / nbElt;
+}
+
+void ContourFitting::fAlpha(double x, double &fn, double &df)
+{
+    int 	nbElt = static_cast<int>(rho.size());
+    double	s1 = 0, s2 = 0, s3 = 0, s4 = 0;
+    double	ds1 = 0, ds2 = 0, ds3 = 0, ds4 = 0;
+
+    for (int n = 1; n <= fdSize; n++)
+    {
+        s1 += rho[n] * sin(psi[n] + frequence[n] * x) +
+            rho[nbElt - n] * sin(psi[nbElt - n] + frequence[nbElt - n] * x);
+        s2 += frequence[n] * rho[n] * cos(psi[n] + frequence[n] * x) +
+            frequence[nbElt - n] * rho[nbElt - n] * cos(psi[nbElt - n] + frequence[nbElt - n] * x);
+        s3 += rho[n] * cos(psi[n] + frequence[n] * x) +
+            rho[nbElt - n] * cos(psi[nbElt - n] + frequence[nbElt - n] * x);
+        s4 += frequence[n] * rho[n] * sin(psi[n] + frequence[n] * x) +
+            frequence[nbElt - n] * rho[nbElt - n] * sin(psi[nbElt - n] + frequence[nbElt - n] * x);
+        ds1 += frequence[n] * rho[n] * cos(psi[n] + frequence[n] * x) +
+            frequence[nbElt - n] * rho[nbElt - n] * cos(psi[nbElt - n] + frequence[nbElt - n] * x);
+        ds2 += -frequence[n] * frequence[n] * rho[n] * sin(psi[n] + frequence[n] * x) -
+            frequence[nbElt - n] * frequence[nbElt - n] * rho[nbElt - n] * sin(psi[nbElt - n] + frequence[nbElt - n] * x);
+        ds3 += -frequence[n] * rho[n] * sin(psi[n] + frequence[n] * x) -
+            frequence[nbElt - n] * rho[nbElt - n] * sin(psi[nbElt - n] + frequence[nbElt - n] * x);
+        ds4 += frequence[n] * frequence[n] * rho[n] * cos(psi[n] + frequence[n] * x) +
+            frequence[nbElt - n] * frequence[nbElt - n] * rho[nbElt - n] * cos(psi[nbElt - n] + frequence[nbElt - n] * x);
+    }
+    fn = s1 * s2 - s3 *s4;
+    df = ds1 * s2 + s1 * ds2 - ds3 * s4 - s3 * ds4;
+}
+
+double ContourFitting::distance(std::complex<double> r, double alpha)
+{
+    std::complex<double>		j(0, 1);
+    double 		d = 0;
+
+    for (int i = 1; i <= fdSize; i++)
+        d += abs(a[i] - b[i] * r*exp(j*alpha*frequence[i])) + abs(a[a.size() - i] - b[a.size() - i] * r*exp(j*alpha*frequence[a.size() - i]));
+    return d/fdSize/2;
+};
+
+double  ContourFitting::newtonRaphson(double x1, double x2)
+{
+    double f1,df1;
+    fAlpha(x1,f1,df1);
+    if (f1 < 0)
+    {
+        x1=x2;
+        fAlpha(x1, f1, df1);
+    }
+    CV_Assert(f1>=0);
+    if (f1==0)
+        return x1;
+    for (int i = 0; i < 5; i++)
+    {
+        x1 = x1 -f1/df1;
+        fAlpha(x1, f1, df1);
+        if (f1 == 0)
+            return x1;
+    }
+    return x1;
+}
+
+void ContourFitting::estimateTransformation(InputArray _src, InputArray _ref, OutputArray _alphaPhiST, double &distFin, bool fdContour)
+{
+    estimateTransformation(_src, _ref, _alphaPhiST, &distFin, fdContour);
+}
+void ContourFitting::estimateTransformation(InputArray _src, InputArray _ref, OutputArray _alphaPhiST,double *distFin, bool fdContour)
+{
+    if (!fdContour)
+        CV_Assert( _src.kind() == _InputArray::STD_VECTOR && _ref.kind() == _InputArray::STD_VECTOR);
+    else
+        CV_Assert(fdContour && _src.kind() == _InputArray::MAT && _ref.kind() == _InputArray::MAT);
+    CV_Assert(_src.channels() == 2 && _ref.channels() == 2);
+    Mat fdCtr1,fdCtr2;
+    if (!fdContour)
+    {
+        Mat newCtr1,newCtr2;
+        contourSampling(_src, newCtr1, ctrSize);
+        contourSampling(_ref, newCtr2, ctrSize);
+        fourierDescriptor(newCtr1, fdCtr1);
+        fourierDescriptor(newCtr2, fdCtr2);
+    }
+    else
+    {
+        fdCtr1=_src.getMat();
+        fdCtr2= _ref.getMat();
+        CV_Assert(fdCtr1.rows == fdCtr2.rows);
+    }
+    CV_Assert( fdSize<= ctrSize / 2 - 1);
+    if (fdCtr1.type() != CV_64FC2)
+        fdCtr1.convertTo(fdCtr1, CV_64F);
+    if (fdCtr2.type() != CV_64FC2)
+        fdCtr2.convertTo(fdCtr2, CV_64F);
+    rho.resize(ctrSize);
+    psi.resize(ctrSize);
+    b.resize(ctrSize);
+    a.resize(ctrSize);
+    frequencyInit();
+    double  alphaMin, phiMin, sMin;
+    long		    n, nbElt = ctrSize;
+    double 		    s1, s2, sign1, sign2, df, x1 = nbElt, x2 = nbElt, dx;
+    double		    dist, distMin = 10000, alpha, s, phi;
+    std::complex<double> 	j(0, 1), zz;
+
+    for (n = 0; n<nbElt; n++)
+    {
+        b[n] = std::complex<double>(fdCtr1.at<Vec2d>(n,0)[0], fdCtr1.at<Vec2d>(n, 0)[1]);
+        a[n] = std::complex<double>(fdCtr2.at<Vec2d>(n, 0)[0], fdCtr2.at<Vec2d>(n, 0)[1]);
+        zz = conj(a[n])*b[n];
+        rho[n] = abs(zz);
+        psi[n] = arg(zz);
+    }
+
+    x1 = nbElt, x2 = nbElt;
+    sMin = 1;
+    alphaMin = 0;
+    phiMin = arg(a[1] / b[1]);
+    do
+    {
+        x2 = x1;
+        fAlpha(x2, sign2, df);
+        dx = 1;
+        x1 = x2;
+        do
+        {
+            x2 = x1;
+            x1 -= dx;
+            fAlpha(x1, sign1, df);
+        }
+        while ((sign1*sign2>0) && (x1>-nbElt));
+        if (sign1*sign2<0)
+        {
+            alpha=newtonRaphson(x1,x2);
+            s1 = 0;
+            s2 = 0;
+            for (n = 1; n<nbElt; n++)
+            {
+                s1 += rho[n] * sin(psi[n] + frequence[n] * alpha);
+                s2 += rho[n] * cos(psi[n] + frequence[n] * alpha);
+            }
+            phi = -atan2(s1, s2);
+            s1 = 0;
+            s2 = 0;
+            for (n = 1; n < nbElt; n++)
+            {
+                s1 += rho[n] * cos(psi[n] + frequence[n] * alpha + phi);
+                s2 += abs(b[n] * conj(b[n]));
+            }
+            s = s1 / s2;
+            zz = s*exp(j * phi);
+            if (s>0)
+                dist = distance(zz, alpha);
+            else
+                dist = 10000;
+            if (dist<distMin)
+            {
+                distMin = dist;
+                alphaMin = alpha;
+                phiMin = phi;
+                sMin = s;
+            }
+        }
+    }
+    while ((x1>-nbElt));
+    Mat x=(Mat_<double>(1,5)<<alphaMin/ nbElt,phiMin,sMin, fdCtr2.at<Vec2d>(0, 0)[0]- fdCtr1.at<Vec2d>(0, 0)[0], fdCtr2.at<Vec2d>(0, 0)[1]- fdCtr1.at<Vec2d>(0, 0)[1]);
+    if (distFin)
+        *distFin= distMin;
+    x.copyTo(_alphaPhiST);
+}
+
+void fourierDescriptor(InputArray _src, OutputArray _dst, int nbElt, int nbFD)
+{
+    CV_Assert(_src.kind() == _InputArray::MAT || _src.kind() == _InputArray::STD_VECTOR);
+    CV_Assert(_src.empty() || (_src.channels() == 2 && (_src.depth() == CV_32S || _src.depth() == CV_32F || _src.depth() == CV_64F)));
+    Mat  z = _src.getMat();
+    CV_Assert(z.rows == 1 || z.cols == 1);
+    if (nbElt==-1)
+        nbElt = getOptimalDFTSize(max(z.rows, z.cols));
+    CV_Assert((nbFD >= 1 && nbFD <=nbElt/2) || nbFD==-1);
+    Mat  Z;
+    if (z.rows*z.cols!=nbElt)
+        contourSampling(_src, z,nbElt);
+    dft(z, Z, DFT_SCALE | DFT_REAL_OUTPUT);
+    if (nbFD == -1)
+    {
+        Z.copyTo(_dst);
+    }
+    else
+    {
+        int n1 = nbFD / 2, n2 = nbElt - n1;
+        Mat d(nbFD, 1, Z.type());
+        Z.rowRange(Range(1, n1+1)).copyTo(d.rowRange(Range(0, n1)));
+        if (n2>0)
+            Z.rowRange(Range(n2, Z.rows)).copyTo(d.rowRange(Range(n1, nbFD)));
+        d.copyTo(_dst);
+    }
+}
+
+void contourSampling(InputArray _src, OutputArray _out, int nbElt)
+{
+    CV_Assert(_src.kind() == _InputArray::STD_VECTOR || _src.kind() == _InputArray::MAT);
+    CV_Assert(_src.empty() || (_src.channels() == 2 && (_src.depth() == CV_32S || _src.depth() == CV_32F || _src.depth() == CV_64F)));
+    CV_Assert(nbElt>0);
+    Mat ctr;
+    _src.getMat().convertTo(ctr,CV_32F);
+    if (ctr.rows*ctr.cols == 0)
+    {
+        _out.release();
+        return;
+    }
+    CV_Assert(ctr.rows==1 || ctr.cols==1);
+    double		l1 = 0, l2, p, d, s;
+ //   AutoBuffer<Point2d> _buf(nbElt);
+    Mat r;
+    if (ctr.rows==1)
+        ctr=ctr.t();
+    int j = 0;
+    int 		nb = ctr.rows;
+    p = arcLength(_src, true);
+    l2 = norm(ctr.row(j) - ctr.row(j + 1)) / p;
+    for (int i = 0; i<nbElt; i++)
+    {
+        s = (float)i / (float)nbElt;
+        while (s >= l2)
+        {
+            j++;
+            l1 = l2;
+            d = norm(ctr.row(j % nb) - ctr.row((j + 1) % nb));
+            l2 = l1 + d / p;
+        }
+        if ((s >= l1) && (s < l2))
+        {
+            Mat d1=ctr.row((j + 1) % nb);
+            Mat d0=ctr.row(j % nb);
+            Mat d10 = d1 - d0;
+            Mat pn = d0 + d10 * (s - l1) / (l2 - l1);
+            r.push_back(pn);
+ //           _buf[i]=Point2d(pn.at<Point2f>(0,0));
+        }
+    }
+    r.copyTo(_out);
+}
+
+void transform(InputArray _src, InputArray _t,OutputArray _dst,  bool fdContour)
+{
+    if (!fdContour)
+        CV_Assert(_src.kind() == _InputArray::STD_VECTOR);
+    else
+        CV_Assert( _src.kind() == _InputArray::MAT );
+    CV_Assert(_src.channels() == 2);
+    CV_Assert(_t.kind() == _InputArray::MAT);
+    Mat t=_t.getMat();
+    CV_Assert(t.rows == 1 && t.cols==5 && t.depth()==CV_64F);
+    Mat Z;
+    if (!fdContour)
+    {
+        Mat ctr1 = _src.getMat();
+        if (ctr1.rows==1)
+            ctr1 = ctr1.t();
+        Mat newCtr1;
+        int M = getOptimalDFTSize(ctr1.rows);
+        contourSampling(ctr1, newCtr1, M);
+        fourierDescriptor(newCtr1, Z);
+    }
+    else
+        Z = _src.getMat();
+    if (Z.type()!=CV_64FC2)
+        Z.convertTo(Z,CV_64F);
+    std::complex<double> expitheta = t.at<double>(0,2)  * std::complex<double>(cos(t.at<double>(0, 1)), sin(t.at<double>(0, 1)));
+    for (int j = 1; j<Z.rows; j++)
+    {
+        std::complex<double> zr(Z.at<Vec2d>(j, 0)[0], Z.at<Vec2d>(j,0)[1]);
+        if (j<=Z.rows/2)
+            zr = zr*expitheta*exp(t.at<double>(0, 0) * 2 * (M_PI*j) *  std::complex<double>(0, 1));
+        else
+            zr = zr*expitheta*exp(t.at<double>(0, 0)* 2 * (M_PI*(j - Z.rows))    * std::complex<double>(0, 1));
+        Z.at<Vec2d>(j, 0) = Vec2d(zr.real(),zr.imag());
+    }
+    Z.at<Vec2d>(0, 0) += Vec2d(t.at<double>(0, 3), t.at<double>(0, 4));
+    std::vector<Point2d> z;
+    dft(Z, z, DFT_INVERSE);
+    std::vector<Point> c;
+    for (int j = 0; j<static_cast<int>(z.size()); j++)
+        c.push_back(Point2d(z[j].x, z[j].y));
+    Mat(z).copyTo(_dst);
+}
+
+cv::Ptr<ContourFitting> create(int ctr, int fd)
+{
+    return makePtr<ContourFitting>(ctr, fd);
+}
+
+}
+}