I wanted to detect ellipse in an image. Since I was learning Mathematica at that time, I asked a question here and got a satisfactory result from the answer below, which used the RANSAC algorithm to detect ellipse.
我想检测图像中的椭圆。因为当时我正在学习Mathematica，所以我在这里问了一个问题，得到了一个满意的结果，下面的答案是RANSAC算法检测椭圆。
However, recently I need to port it to OpenCV, but there are some functions that only exist in Mathematica. One of the key function is the "GradientOrientationFilter" function.
然而，最近我需要将它移植到OpenCV，但是有些函数只存在于Mathematica中。其中一个关键的功能是“梯度方向过滤器”功能。
Since there are five parameters for a general ellipse, I need to sample five points to determine one. Howevere, the more sampling points indicates the lower chance to have a good guess, which leads to the lower success rate in ellipse detection. Therefore, the answer from Mathematica add another condition, that is the gradient of the image must be parallel to the gradient of the ellipse equation. Anyway, we'll only need three points to determine one ellipse using least square from the Mathematica approach. The result is quite good.
由于一般椭圆有5个参数，我需要采样5个点来确定一个。然而，采样点越多，猜测的概率越低，椭圆检测的成功率越低。因此，Mathematica给出的答案增加了另一个条件，即图像的梯度必须与椭圆方程的梯度平行。无论如何，我们只需要3个点就可以用最小二乘法确定一个椭圆。结果很好。
However, when I try to find the image gradient using Sobel or Scharr operator in OpenCV, it is not good enough, which always leads to the bad result.
但是，当我在OpenCV中尝试使用Sobel或Scharr算子来寻找图像梯度时，它是不够好的，这总是导致不好的结果。
How to calculate the gradient or the tangent of an image accurately? Thanks!
如何准确计算图像的梯度或正切?谢谢!
Result with gradient, three points 
结果与梯度，三点
Result without gradient, five points 
结果没有梯度，五点
----------updated----------
- - - - - - - - - - - - - - - - - - - - - -更新
I did some edge detect and median blur beforehand and draw the result on the edge image. My original test image is like this:
我事先做了一些边缘检测和中值模糊，并在边缘图像上绘制结果。我最初的测试图像是这样的:
In general, my final goal is to detect the ellipse in a scene or on an object. Something like this:
一般来说，我的最终目标是检测场景或对象中的椭圆。是这样的:
That's why I choose to use RANSAC to fit the ellipse from edge points.
这就是为什么我选择使用RANSAC从边缘点拟合椭圆。
3 个解决方案
#1
1
As for your final goal, you may try
至于你的最终目标，你可以试试
findContours and [fitEllipse] in OpenCV
OpenCV中的findContours和[fitEllipse]
The pseudo code will be
伪代码将是
1) some image process
1)一些形象的过程
2) find all contours
2)找到所有轮廓
3) fit each contours by fitEllipse
3)用fitEllipse对每个等高线进行拟合
here is part of code I use before
这是我以前使用的部分代码
[... image process ....you get a bwimage ]
vector<vector<Point> > contours;
findContours(bwimage, contours, CV_RETR_LIST, CV_CHAIN_APPROX_NONE);
for(size_t i = 0; i < contours.size(); i++)
{
size_t count = contours[i].size();
Mat pointsf;
Mat(contours[i]).convertTo(pointsf, CV_32F);
RotatedRect box = fitEllipse(pointsf);
/* You can put some limitation about size and aspect ratio here */
if( box.size.width > 20 &&
box.size.height > 20 &&
box.size.width < 80 &&
box.size.height < 80 )
{
if( MAX(box.size.width, box.size.height) > MIN(box.size.width, box.size.height)*30 )
continue;
//drawContours(SrcImage, contours, (int)i, Scalar::all(255), 1, 8);
ellipse(SrcImage, box, Scalar(0,0,255), 1, CV_AA);
ellipse(SrcImage, box.center, box.size*0.5f, box.angle, 0, 360, Scalar(200,255,255), 1, CV_AA);
}
}
imshow("result", SrcImage);
#2
0
You're using terms in an unusual way.
你使用术语的方式不同寻常。
Normally for images, the term "gradient" is interpreted as if the image is a mathematical function f(x,y). This gives us a (df/dx, df/dy) vector in each point.
通常对于图像，术语“梯度”被解释为如果图像是一个数学函数f(x,y)。这给了我们每个点的(df/dx, df/dy)向量。
Yet you're looking at the image as if it's a function y = f(x) and the gradient would be f(x)/dx.
你看着图像就好像它是一个函数y = f(x)梯度是f(x)/dx。
Now, if you look at your image, you'll see that the two interpretations are definitely related. Your ellipse is drawn as a set of contrasting pixels, and as a result there are two sharp gradients in the image - the inner and outer. These of course correspond to the two normal vectors, and therefore are in opposite directions.
现在，如果你看看你的图像，你会发现这两种解释是肯定相关的。你的椭圆被画成一组对比的像素，结果在图像中有两个清晰的渐变——内部和外部。这些当然与两个法向量对应，因此方向相反。
Also note that your image has pixels. The gradient is also pixelated. The way your ellipse is drawn, with a single pixel width means that your local gradient takes on only values that are a multiple of 45 degrees:
还要注意，您的图像有像素。渐变也被像素化。你绘制椭圆的方式，只有一个像素宽度，意味着你的局部梯度只能取45度的倍数:
▄▄ ▄▀ ▌ ▀▄
#3
0
If you focus on ellipse(no other shape), you can treat the value of the pixels of the ellipse as mass of the points.
如果您关注椭圆(没有其他形状)，您可以将椭圆的像素值视为点的质量。
Then you can calculate the moment of inertial Ixx, Iyy, Ixy to find out the angle, theta, which can rotate a general ellipse back to a canonical form (X-Xc)^2/a + (Y-Yc)^2/b = 1.
然后你可以计算惯性矩Ixx,Iyy,Ixy找出角,θ,一般椭圆旋转回规范形式(X-Xc)^ 2 b / a +(Y-Yc)^ 2 = 1。
Then you can find out Xc and Yc by the center of mass.
然后你可以通过质心求出Xc和Yc。
Then you can find out a and b by min X and min Y.
然后你可以用min X和min Y求出a和b。
--------------- update -----------
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -更新
This method can apply to filled ellipse too.
该方法也适用于填充椭圆。
More than one ellipse on a single image will fail unless you segment them first.
一个图像上的多个椭圆将失败，除非您首先对它们进行分割。
Let me explain more, I will use C to represent cos(theta) and S to represent sin(theta)
让我解释一下，我用C来表示cos和S表示sin
After rotation to canonical form, the new X is [eq0] X=xC-yS and Y is Y=xS+yC where x and y are original positions.
旋转到标准形式后，新的X为[eq0] X=xC-yS, Y为Y=xS+yC，其中X和Y为原始位置。
The rotation will give you min IYY.
旋转会得到最小值。
[eq1]
(eq1)
IYY= Sum(m*Y*Y) = Sum{m*(xS+yC)(xS+yC)} = Sum{ m(xxSS+yyCC+xySC) = Ixx*S^2 + Iyy*C^2 + Ixy*S*C
IYY =(m * * Y)之和=和{ m *(x + yC)(x + yC)} = {和m(xxs + yyCC + xySC)= Ixx * S ^ 2 + IYY * C ^ 2 + Ixy * * C
For min IYY, d(IYY)/d(theta) = 0 that is
对于最小IYY, d(IYY)/d() = 0
2IxxSC - 2IyySC + Ixy(CC-SS) = 0
2IxxSC - 2IyySC + Ixy(CC-SS) = 0
2(Ixx-Iyy)/Ixy = (SS-CC)/SC = S/C+C/S = Z+1/Z
2(Ixx-Iyy)/Ixy = (SS-CC)/SC = S/C+C/S = Z+1/Z
While programming, the LHS is just a number, let's said N
编程时，LHS只是一个数字，假设是N
Z^2 - NZ +1 =0
Z ^ 2 -新西兰+ 1 = 0
So there are two roots of Z hence theta, let's said Z1 and Z2, one will min the IYY and the other will max the IYY.
有两个根，Z1和Z2，一个是最小的IYY另一个是最大的IYY。
----------- pseudo code --------
- - - - - - - - - - - - - - - - - - - - - -伪代码
Compute Ixx, Iyy, Ixy for a hollow or filled ellipse.
计算空心椭圆或填充椭圆的Ixx、Iyy、Ixy。
Compute theta1=atan(Z1) and theta2=atan(Z2)
计算θ₁=(Z1)每股和θ=每股(Z2)
Put These two theta into eq1 find which is smaller. Then you get theta.
把这两项放到eq1中，发现哪个更小。然后你得到θ。
Go back to those non-zero pixels, transfer them to new X and Y by the theta you found.
回到那些非零像素点，把它们转换成新的X和Y。
Find center of mass Xc Yc and min X and min Y by sort().
通过排序()求出质量中心Xc Yc, min X, min Y。
-------------- by hand -----------
用手- - - - - - - - - - - - - - - - - - - - - - - - - - - -
If you need the original equation of the ellipse
如果你需要椭圆的原始方程。
Just put [eq0] into the canonical form
把[eq0]放入规范形式。
#1
1
As for your final goal, you may try
至于你的最终目标，你可以试试
findContours and [fitEllipse] in OpenCV
OpenCV中的findContours和[fitEllipse]
The pseudo code will be
伪代码将是
1) some image process
1)一些形象的过程
2) find all contours
2)找到所有轮廓
3) fit each contours by fitEllipse
3)用fitEllipse对每个等高线进行拟合
here is part of code I use before
这是我以前使用的部分代码
[... image process ....you get a bwimage ]
vector<vector<Point> > contours;
findContours(bwimage, contours, CV_RETR_LIST, CV_CHAIN_APPROX_NONE);
for(size_t i = 0; i < contours.size(); i++)
{
size_t count = contours[i].size();
Mat pointsf;
Mat(contours[i]).convertTo(pointsf, CV_32F);
RotatedRect box = fitEllipse(pointsf);
/* You can put some limitation about size and aspect ratio here */
if( box.size.width > 20 &&
box.size.height > 20 &&
box.size.width < 80 &&
box.size.height < 80 )
{
if( MAX(box.size.width, box.size.height) > MIN(box.size.width, box.size.height)*30 )
continue;
//drawContours(SrcImage, contours, (int)i, Scalar::all(255), 1, 8);
ellipse(SrcImage, box, Scalar(0,0,255), 1, CV_AA);
ellipse(SrcImage, box.center, box.size*0.5f, box.angle, 0, 360, Scalar(200,255,255), 1, CV_AA);
}
}
imshow("result", SrcImage);
#2
0
You're using terms in an unusual way.
你使用术语的方式不同寻常。
Normally for images, the term "gradient" is interpreted as if the image is a mathematical function f(x,y). This gives us a (df/dx, df/dy) vector in each point.
通常对于图像，术语“梯度”被解释为如果图像是一个数学函数f(x,y)。这给了我们每个点的(df/dx, df/dy)向量。
Yet you're looking at the image as if it's a function y = f(x) and the gradient would be f(x)/dx.
你看着图像就好像它是一个函数y = f(x)梯度是f(x)/dx。
Now, if you look at your image, you'll see that the two interpretations are definitely related. Your ellipse is drawn as a set of contrasting pixels, and as a result there are two sharp gradients in the image - the inner and outer. These of course correspond to the two normal vectors, and therefore are in opposite directions.
现在，如果你看看你的图像，你会发现这两种解释是肯定相关的。你的椭圆被画成一组对比的像素，结果在图像中有两个清晰的渐变——内部和外部。这些当然与两个法向量对应，因此方向相反。
Also note that your image has pixels. The gradient is also pixelated. The way your ellipse is drawn, with a single pixel width means that your local gradient takes on only values that are a multiple of 45 degrees:
还要注意，您的图像有像素。渐变也被像素化。你绘制椭圆的方式，只有一个像素宽度，意味着你的局部梯度只能取45度的倍数:
▄▄ ▄▀ ▌ ▀▄
#3
0
If you focus on ellipse(no other shape), you can treat the value of the pixels of the ellipse as mass of the points.
如果您关注椭圆(没有其他形状)，您可以将椭圆的像素值视为点的质量。
Then you can calculate the moment of inertial Ixx, Iyy, Ixy to find out the angle, theta, which can rotate a general ellipse back to a canonical form (X-Xc)^2/a + (Y-Yc)^2/b = 1.
然后你可以计算惯性矩Ixx,Iyy,Ixy找出角,θ,一般椭圆旋转回规范形式(X-Xc)^ 2 b / a +(Y-Yc)^ 2 = 1。
Then you can find out Xc and Yc by the center of mass.
然后你可以通过质心求出Xc和Yc。
Then you can find out a and b by min X and min Y.
然后你可以用min X和min Y求出a和b。
--------------- update -----------
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -更新
This method can apply to filled ellipse too.
该方法也适用于填充椭圆。
More than one ellipse on a single image will fail unless you segment them first.
一个图像上的多个椭圆将失败，除非您首先对它们进行分割。
Let me explain more, I will use C to represent cos(theta) and S to represent sin(theta)
让我解释一下，我用C来表示cos和S表示sin
After rotation to canonical form, the new X is [eq0] X=xC-yS and Y is Y=xS+yC where x and y are original positions.
旋转到标准形式后，新的X为[eq0] X=xC-yS, Y为Y=xS+yC，其中X和Y为原始位置。
The rotation will give you min IYY.
旋转会得到最小值。
[eq1]
(eq1)
IYY= Sum(m*Y*Y) = Sum{m*(xS+yC)(xS+yC)} = Sum{ m(xxSS+yyCC+xySC) = Ixx*S^2 + Iyy*C^2 + Ixy*S*C
IYY =(m * * Y)之和=和{ m *(x + yC)(x + yC)} = {和m(xxs + yyCC + xySC)= Ixx * S ^ 2 + IYY * C ^ 2 + Ixy * * C
For min IYY, d(IYY)/d(theta) = 0 that is
对于最小IYY, d(IYY)/d() = 0
2IxxSC - 2IyySC + Ixy(CC-SS) = 0
2IxxSC - 2IyySC + Ixy(CC-SS) = 0
2(Ixx-Iyy)/Ixy = (SS-CC)/SC = S/C+C/S = Z+1/Z
2(Ixx-Iyy)/Ixy = (SS-CC)/SC = S/C+C/S = Z+1/Z
While programming, the LHS is just a number, let's said N
编程时，LHS只是一个数字，假设是N
Z^2 - NZ +1 =0
Z ^ 2 -新西兰+ 1 = 0
So there are two roots of Z hence theta, let's said Z1 and Z2, one will min the IYY and the other will max the IYY.
有两个根，Z1和Z2，一个是最小的IYY另一个是最大的IYY。
----------- pseudo code --------
- - - - - - - - - - - - - - - - - - - - - -伪代码
Compute Ixx, Iyy, Ixy for a hollow or filled ellipse.
计算空心椭圆或填充椭圆的Ixx、Iyy、Ixy。
Compute theta1=atan(Z1) and theta2=atan(Z2)
计算θ₁=(Z1)每股和θ=每股(Z2)
Put These two theta into eq1 find which is smaller. Then you get theta.
把这两项放到eq1中，发现哪个更小。然后你得到θ。
Go back to those non-zero pixels, transfer them to new X and Y by the theta you found.
回到那些非零像素点，把它们转换成新的X和Y。
Find center of mass Xc Yc and min X and min Y by sort().
通过排序()求出质量中心Xc Yc, min X, min Y。
-------------- by hand -----------
用手- - - - - - - - - - - - - - - - - - - - - - - - - - - -
If you need the original equation of the ellipse
如果你需要椭圆的原始方程。
Just put [eq0] into the canonical form
把[eq0]放入规范形式。