Skip to content

Commit a44a2ba

Browse files
committed
Merge pull request #1216 from LaurentBerger:FourierDescriptors
2 parents 1a61e8d + 8b733e0 commit a44a2ba

File tree

5 files changed

+668
-3
lines changed

5 files changed

+668
-3
lines changed

modules/ximgproc/doc/ximgproc.bib

Lines changed: 24 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -76,9 +76,9 @@ @incollection{PFF2004
7676
publisher={Springer}
7777
}
7878

79-
@article{deriche1987using,
79+
@article{deriche1987,
8080
title={Using Canny's criteria to derive a recursively implemented optimal edge detector},
81-
author={Deriche, Rachid},
81+
author={Deriche Rachid},
8282
journal={International journal of computer vision},
8383
volume={1},
8484
number={2},
@@ -259,3 +259,25 @@ @inproceedings{Khurshid2009
259259
year={2009},
260260
organization={International Society for Optics and Photonics}
261261
}
262+
263+
@article{BergerRaghunathan1998,
264+
title={Coalescence in 2 dimensions: experiments on thin copolymer films and numerical simulations},
265+
author={Berger, L and Raghunathan, V A and Launay, C and Ausserré, D and Gallot, Y},
266+
journal={The European Physical Journal B - Condensed Matter and Complex Systems},
267+
volume={2},
268+
number={1},
269+
pages={93-99},
270+
year={1998},
271+
publisher={Springer}
272+
}
273+
274+
@article{PersoonFu1977,
275+
title={Shape Discrimination Using Fourier Descriptors},
276+
author={E Persoon and King-Sun Fu},
277+
journal={IEEE Transactions on Pattern Analysis and Machine Intelligence},
278+
volume={7},
279+
number={3},
280+
pages={170-179},
281+
year={1977},
282+
publisher={IEEE Computer Society}
283+
}

modules/ximgproc/include/opencv2/ximgproc.hpp

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,8 @@
5252
#include "ximgproc/fast_line_detector.hpp"
5353
#include "ximgproc/deriche_filter.hpp"
5454
#include "ximgproc/peilin.hpp"
55+
#include "ximgproc/fourier_descriptors.hpp"
56+
5557

5658
/** @defgroup ximgproc Extended Image Processing
5759
@{
@@ -67,7 +69,9 @@ i.e. algorithms which somehow takes into account pixel affinities in natural ima
6769
@defgroup ximgproc_segmentation Image segmentation
6870
6971
@defgroup ximgproc_fast_line_detector Fast line detector
70-
@}
72+
73+
@defgroup ximgproc_fourier Fourier descriptors
74+
@}
7175
*/
7276

7377
namespace cv
Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
// This file is part of OpenCV project.
2+
// It is subject to the license terms in the LICENSE file found in the top-level directory
3+
// of this distribution and at http://opencv.org/license.html.
4+
5+
#ifndef __OPENCV_FOURIERDESCRIPTORS_HPP__
6+
#define __OPENCV_FOURIERDESCRIPTORS_HPP__
7+
8+
#include <opencv2/core.hpp>
9+
10+
namespace cv {
11+
namespace ximgproc {
12+
13+
//! @addtogroup ximgproc_fourier
14+
//! @{
15+
16+
/** @brief Class for ContourFitting algorithms.
17+
ContourFitting match two contours \f$ z_a \f$ and \f$ z_b \f$ minimizing distance
18+
\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
19+
*/
20+
class CV_EXPORTS_W ContourFitting : public Algorithm
21+
{
22+
int ctrSize;
23+
int fdSize;
24+
std::vector<std::complex<double> > b;
25+
std::vector<std::complex<double> > a;
26+
std::vector<double> frequence;
27+
std::vector<double> rho, psi;
28+
void frequencyInit();
29+
void fAlpha(double x, double &fn, double &df);
30+
double distance(std::complex<double> r, double alpha);
31+
double newtonRaphson(double x1, double x2);
32+
public:
33+
/** @brief Fit two closed curves using fourier descriptors. More details in @cite PersoonFu1977 and @cite BergerRaghunathan1998
34+
35+
* @param ctr number of Fourier descriptors equal to number of contour points after resampling.
36+
* @param fd Contour defining second shape (Target).
37+
*/
38+
ContourFitting(int ctr=1024,int fd=16):ctrSize(ctr),fdSize(fd){};
39+
/** @brief Fit two closed curves using fourier descriptors. More details in @cite PersoonFu1977 and @cite BergerRaghunathan1998
40+
41+
@param src Contour defining first shape.
42+
@param dst Contour defining second shape (Target).
43+
@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
44+
@param dist distance between src and dst after matching.
45+
@param fdContour false then src and dst are contours and true src and dst are fourier descriptors.
46+
*/
47+
CV_WRAP void estimateTransformation(InputArray src, InputArray dst, OutputArray alphaPhiST, double *dist = 0, bool fdContour = false);
48+
/** @brief Fit two closed curves using fourier descriptors. More details in @cite PersoonFu1977 and @cite BergerRaghunathan1998
49+
50+
@param src Contour defining first shape.
51+
@param dst Contour defining second shape (Target).
52+
@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
53+
@param dist distance between src and dst after matching.
54+
@param fdContour false then src and dst are contours and true src and dst are fourier descriptors.
55+
*/
56+
CV_WRAP void estimateTransformation(InputArray src, InputArray dst, OutputArray alphaPhiST, double &dist , bool fdContour = false);
57+
/** @brief set number of Fourier descriptors used in estimateTransformation
58+
59+
@param n number of Fourier descriptors equal to number of contour points after resampling.
60+
*/
61+
CV_WRAP void setCtrSize(int n);
62+
/** @brief set number of Fourier descriptors when estimateTransformation used vector<Point>
63+
64+
@param n number of fourier descriptors used for optimal curve matching.
65+
*/
66+
CV_WRAP void setFDSize(int n);
67+
/**
68+
@returns number of fourier descriptors
69+
*/
70+
CV_WRAP int getCtrSize() { return ctrSize; };
71+
/**
72+
@returns number of fourier descriptors used for optimal curve matching
73+
*/
74+
CV_WRAP int getFDSize() { return fdSize; };
75+
};
76+
/**
77+
* @brief Fourier descriptors for planed closed curves
78+
*
79+
* For more details about this implementation, please see @cite PersoonFu1977
80+
*
81+
* @param src contour type vector<Point> , vector<Point2f> or vector<Point2d>
82+
* @param dst Mat of type CV_64FC2 and nbElt rows A VERIFIER
83+
* @param nbElt number of rows in dst or getOptimalDFTSize rows if nbElt=-1
84+
* @param nbFD number of FD return in dst dst = [FD(1...nbFD/2) FD(nbFD/2-nbElt+1...:nbElt)]
85+
*
86+
*/
87+
CV_EXPORTS_W void fourierDescriptor(InputArray src, OutputArray dst, int nbElt=-1,int nbFD=-1);
88+
/**
89+
* @brief transform a contour
90+
*
91+
* @param src contour or Fourier Descriptors if fd is true
92+
* @param t transform Mat given by estimateTransformation
93+
* @param dst Mat of type CV_64FC2 and nbElt rows
94+
* @param fdContour true src are Fourier Descriptors. fdContour false src is a contour
95+
*
96+
*/
97+
CV_EXPORTS_W void transform(InputArray src, InputArray t,OutputArray dst, bool fdContour=true);
98+
/**
99+
* @brief Contour sampling .
100+
*
101+
* @param src contour type vector<Point> , vector<Point2f> or vector<Point2d>
102+
* @param out Mat of type CV_64FC2 and nbElt rows
103+
* @param nbElt number of points in out contour
104+
*
105+
*/
106+
CV_EXPORTS_W void contourSampling(InputArray src, OutputArray out, int nbElt);
107+
108+
/**
109+
* @brief create
110+
111+
* @param ctr number of Fourier descriptors equal to number of contour points after resampling.
112+
* @param fd Contour defining second shape (Target).
113+
*/
114+
CV_EXPORTS_W Ptr<ContourFitting> create(int ctr = 1024, int fd = 16);
115+
116+
//! @} ximgproc_fourier
117+
}
118+
}
119+
#endif
Lines changed: 186 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
#include <opencv2/core.hpp>
2+
#include <opencv2/core/utility.hpp>
3+
#include <opencv2/highgui.hpp>
4+
#include <opencv2/imgproc.hpp>
5+
#include <opencv2/ximgproc.hpp>
6+
#include <iostream>
7+
8+
using namespace cv;
9+
using namespace std;
10+
11+
struct ThParameters {
12+
int levelNoise;
13+
int angle;
14+
int scale10;
15+
int origin;
16+
int xg;
17+
int yg;
18+
bool update;
19+
} ;
20+
21+
static vector<Point> NoisyPolygon(vector<Point> pRef, double n);
22+
static void UpdateShape(int , void *r);
23+
static void AddSlider(String sliderName, String windowName, int minSlider, int maxSlider, int valDefault, int *valSlider, void(*f)(int, void *), void *r);
24+
25+
int main(void)
26+
{
27+
vector<Point> ctrRef;
28+
vector<Point> ctrRotate, ctrNoisy, ctrNoisyRotate, ctrNoisyRotateShift;
29+
// build a shape with 5 vertex
30+
ctrRef.push_back(Point(250,250)); ctrRef.push_back(Point(400, 250));
31+
ctrRef.push_back(Point(400, 300)); ctrRef.push_back(Point(250, 300));ctrRef.push_back(Point(180, 270));
32+
Point cg(0,0);
33+
for (int i=0;i<static_cast<int>(ctrRef.size());i++)
34+
cg+=ctrRef[i];
35+
cg.x /= static_cast<int>(ctrRef.size());
36+
cg.y /= static_cast<int>(ctrRef.size());
37+
ThParameters p;
38+
p.levelNoise=6;
39+
p.angle=45;
40+
p.scale10=5;
41+
p.origin=10;
42+
p.xg=150;
43+
p.yg=150;
44+
p.update=true;
45+
namedWindow("FD Curve matching");
46+
// A rotation with center at (150,150) of angle 45 degrees and a scaling of 5/10
47+
AddSlider("Noise", "FD Curve matching", 0, 20, p.levelNoise, &p.levelNoise, UpdateShape, &p);
48+
AddSlider("Angle", "FD Curve matching", 0, 359, p.angle, &p.angle, UpdateShape, &p);
49+
AddSlider("Scale", "FD Curve matching", 5, 100, p.scale10, &p.scale10, UpdateShape, &p);
50+
AddSlider("Origin%%", "FD Curve matching", 0, 100, p.origin, &p.origin, UpdateShape, &p);
51+
AddSlider("Xg", "FD Curve matching", 150, 450, p.xg, &p.xg, UpdateShape, &p);
52+
AddSlider("Yg", "FD Curve matching", 150, 450, p.yg, &p.yg, UpdateShape, &p);
53+
int code=0;
54+
double dist;
55+
vector<vector<Point> > c;
56+
Mat img;
57+
cout << "******************** PRESS g TO MATCH CURVES *************\n";
58+
do
59+
{
60+
code = waitKey(30);
61+
if (p.update)
62+
{
63+
Mat r = getRotationMatrix2D(Point(p.xg, p.yg), p.angle, 10.0/ p.scale10);
64+
ctrNoisy= NoisyPolygon(ctrRef,static_cast<double>(p.levelNoise));
65+
transform(ctrNoisy, ctrNoisyRotate,r);
66+
ctrNoisyRotateShift.clear();
67+
for (int i=0;i<static_cast<int>(ctrNoisy.size());i++)
68+
ctrNoisyRotateShift.push_back(ctrNoisyRotate[(i+(p.origin*ctrNoisy.size())/100)% ctrNoisy.size()]);
69+
// To draw contour using drawcontours
70+
c.clear();
71+
c.push_back(ctrRef);
72+
c.push_back(ctrNoisyRotateShift);
73+
p.update = false;
74+
Rect rglobal;
75+
for (int i = 0; i < static_cast<int>(c.size()); i++)
76+
{
77+
rglobal = boundingRect(c[i]) | rglobal;
78+
}
79+
rglobal.width += 10;
80+
rglobal.height += 10;
81+
img = Mat::zeros(2 * rglobal.height, 2 * rglobal.width, CV_8UC(3));
82+
drawContours(img, c, 0, Scalar(255,0,0));
83+
drawContours(img, c, 1, Scalar(0, 255, 0));
84+
circle(img, c[0][0], 5, Scalar(255, 0, 0));
85+
circle(img, c[1][0], 5, Scalar(0, 255, 0));
86+
imshow("FD Curve matching", img);
87+
}
88+
if (code == 'd')
89+
{
90+
destroyWindow("FD Curve matching");
91+
namedWindow("FD Curve matching");
92+
// A rotation with center at (150,150) of angle 45 degrees and a scaling of 5/10
93+
AddSlider("Noise", "FD Curve matching", 0, 20, p.levelNoise, &p.levelNoise, UpdateShape, &p);
94+
AddSlider("Angle", "FD Curve matching", 0, 359, p.angle, &p.angle, UpdateShape, &p);
95+
AddSlider("Scale", "FD Curve matching", 5, 100, p.scale10, &p.scale10, UpdateShape, &p);
96+
AddSlider("Origin%%", "FD Curve matching", 0, 100, p.origin, &p.origin, UpdateShape, &p);
97+
AddSlider("Xg", "FD Curve matching", 150, 450, p.xg, &p.xg, UpdateShape, &p);
98+
AddSlider("Yg", "FD Curve matching", 150, 450, p.yg, &p.yg, UpdateShape, &p);
99+
100+
}
101+
if (code == 'g')
102+
{
103+
ximgproc::ContourFitting fit;
104+
vector<Point2f> ctrRef2d, ctrRot2d;
105+
// sampling contour we want 256 points
106+
ximgproc::contourSampling(ctrRef, ctrRef2d, 256); // use a mat
107+
ximgproc::contourSampling(ctrNoisyRotateShift, ctrRot2d, 256); // use a vector og point
108+
fit.setFDSize(16);
109+
Mat t;
110+
fit.estimateTransformation(ctrRot2d, ctrRef2d, t, &dist, false);
111+
cout << "Transform *********\n "<<"Origin = "<< t.at<double>(0,0)*ctrNoisy.size() <<" expected "<< (p.origin*ctrNoisy.size()) / 100 <<" ("<< ctrNoisy.size()<<")\n";
112+
cout << "Angle = " << t.at<double>(0, 1) * 180 / M_PI << " expected " << p.angle <<"\n";
113+
cout << "Scale = " << t.at<double>(0, 2) << " expected " << p.scale10 / 10.0 << "\n";
114+
Mat dst;
115+
ximgproc::transform(ctrRot2d, t, dst, false);
116+
c.push_back(dst);
117+
drawContours(img, c, 2, Scalar(0,255,255));
118+
circle(img, c[2][0], 5, Scalar(0, 255, 255));
119+
imshow("FD Curve matching", img);
120+
}
121+
}
122+
while (code!=27);
123+
124+
return 0;
125+
}
126+
127+
vector<Point> NoisyPolygon(vector<Point> pRef, double n)
128+
{
129+
RNG rng;
130+
vector<Point> c;
131+
vector<Point> p = pRef;
132+
vector<vector<Point> > contour;
133+
for (int i = 0; i<static_cast<int>(p.size()); i++)
134+
p[i] += Point(Point2d(n*rng.uniform((double)-1, (double)1), n*rng.uniform((double)-1, (double)1)));
135+
if (n==0)
136+
return p;
137+
c.push_back(p[0]);
138+
int minX = p[0].x, maxX = p[0].x, minY = p[0].y, maxY = p[0].y;
139+
for (int i = 0; i <static_cast<int>(p.size()); i++)
140+
{
141+
int next = i + 1;
142+
if (next == static_cast<int>(p.size()))
143+
next = 0;
144+
Point2d u = p[next] - p[i];
145+
int d = static_cast<int>(norm(u));
146+
double a = atan2(u.y, u.x);
147+
int step = 1;
148+
if (n != 0)
149+
step = static_cast<int>(d / n);
150+
for (int j = 1; j<d; j += max(step, 1))
151+
{
152+
Point pNew;
153+
do
154+
{
155+
156+
Point2d pAct = (u*j) / static_cast<double>(d);
157+
double r = n*rng.uniform((double)0, (double)1);
158+
double theta = a + rng.uniform(0., 2 * CV_PI);
159+
pNew = Point(Point2d(r*cos(theta) + pAct.x + p[i].x, r*sin(theta) + pAct.y + p[i].y));
160+
} while (pNew.x<0 || pNew.y<0);
161+
if (pNew.x<minX)
162+
minX = pNew.x;
163+
if (pNew.x>maxX)
164+
maxX = pNew.x;
165+
if (pNew.y<minY)
166+
minY = pNew.y;
167+
if (pNew.y>maxY)
168+
maxY = pNew.y;
169+
c.push_back(pNew);
170+
}
171+
}
172+
return c;
173+
}
174+
175+
void UpdateShape(int , void *r)
176+
{
177+
((ThParameters *)r)->update = true;
178+
}
179+
180+
void AddSlider(String sliderName, String windowName, int minSlider, int maxSlider, int valDefault, int *valSlider, void(*f)(int, void *), void *r)
181+
{
182+
createTrackbar(sliderName, windowName, valSlider, 1, f, r);
183+
setTrackbarMin(sliderName, windowName, minSlider);
184+
setTrackbarMax(sliderName, windowName, maxSlider);
185+
setTrackbarPos(sliderName, windowName, valDefault);
186+
}

0 commit comments

Comments
 (0)