如何从椭圆图像中找到上下弧的长度答案

【问题标题】：How to find length of upper and lower arc from ellipse image如何从椭圆图像中找到上下弧的长度
【发布时间】：2016-03-27 04:13:46
【问题描述】：

在这里，我尝试使用图像矢量（图像的轮廓）找到上弧和下弧，但它无法给出提取结果。建议任何其他方法从图像及其长度中找到上弧和下弧。

这是我的代码

    Mat image =cv::imread("thinning/20d.jpg");
    int i=0,j=0,k=0,x=320;
    for(int y = 0; y < image.rows; y++)
    {
    if(image.at<Vec3b>(Point(x, y))[0] >= 250 && image.at<Vec3b>(Point(x, y))[1] >= 250 && image.at<Vec3b>(Point(x, y))[2] >= 250){
              qDebug()<<x<<y;
              x1[i]=x;
              y1[i]=y;
              i=i+1;
    }
    }
    for(i=0;i<=1;i++){
      qDebug()<<x1[i]<<y1[i];
    }
    qDebug()<<"UPPER ARC";
    for(int x = 0; x < image.cols; x++)
    {
      for(int y = 0; y <= (y1[0]+20); y++)
      {
          if(image.at<Vec3b>(Point(x, y))[0] >= 240 && image.at<Vec3b>(Point(x, y))[1] >= 240 && image.at<Vec3b>(Point(x, y))[2] >= 240){
              x2[j]=x;
              y2[j]=y;
              j=j+1;
            qDebug()<<x<<y;
          }}
    }   
    qDebug()<<"Lower ARC";
    for(int x = 0; x < image.cols; x++)
    {
      for(int y = (y1[1]-20); y <= image.rows; y++)
      {
          if(image.at<Vec3b>(Point(x, y))[0] >= 240 && image.at<Vec3b>(Point(x, y))[1] >= 240 && image.at<Vec3b>(Point(x, y))[2] >= 240){
              x3[k]=x;
              y3[k]=y;
              k=k+1;
             qDebug()<<x<<y;
          }}
   }

通过上面的代码我得到坐标，通过使用坐标点我可以找到弧的长度，但它与提取结果不匹配。

这是实际图像：

图片1：

瘦身后得到：

预期输出：

【问题讨论】：

添加规范你所说的上/下弧是什么意思（你的意思是你用穿过椭圆中点的水平线将椭圆切成两半？）还添加预期结果和你的错误输出跨度>
@Spektre Spektre 我不能分成两半，因为上弧略大于下弧，因为我需要找到两个弧的值。然后图像无法进入提取位置。它看起来像椭圆但不是一个完美的椭圆，所以我无法使用任何关于椭圆属性的公式。在这里我需要找到它的椭圆。有没有别的办法。。谢谢你的建议。
那么弧到底是什么，那么它们如何定义边界在哪里（也许是轴心）？如果您只想知道它是否为椭圆，则找到（平均）中点。然后最远和最近的点将为您提供周轴，然后仅计算您的点与该分析椭圆的平均距离。距离越大，这个椭圆形越远。此外，如果 periaxises 不垂直，那么它是倾斜的......所以也不是椭圆
顺便说一句，您的椭圆看起来像正常渲染的圆，直径为D=212 像素，没有任何偏差，那么问题出在哪里？
@Spektre 是的，它给出了直径D=(210-212)。检查第二个图像（image2）它不是一个完美的椭圆。但它与 image1 椭圆给出了相似的值。如何区分这些图像（不是大小）。

标签： c++ qt opencv image-processing ellipse

【解决方案1】：

由于您无法准确定义什么是上/下弧，因此我假设您通过穿过椭圆中点的水平线将椭圆切成两半。如果不是这种情况，那么你必须自己调整它......现在好，该怎么做：

二值化图像

当您提供 JPG 时，颜色会失真，因此不仅仅是黑白
将边框缩小到 1 个像素

用白色填充内部，然后将不与任何黑色像素相邻的所有白色像素重新着色为一些未使用的颜色或黑色。还有许多其他变体如何实现这一点...
找到边界框

搜索所有像素并记住所有白色像素的最小、最大x,y 坐标。让他们叫他们x0,y0,x1,y1。
计算椭圆的中心

只需找到边界框的中点
```
cx=(x0+x1)/2
cy=(y0+y1)/2
```
计算每个椭圆弧的像素数

为每个弧设置计数器，并简单地为任何具有y<=cy 的白色像素增加弧计数器上限，如果y>=cy 则为更低。如果您的坐标系不同，则条件可以相反。
查找椭圆参数

只需找到最接近(cx,cy) 的白色像素，这将是小半轴的端点b，我们称之为(bx,by)。还要找到距离(cx,cy) 最远的白色像素，这将是长半轴端点(ax,ay)。它们与中心之间的距离将为您提供a,b，并且它们的位置减去中心将为您提供椭圆旋转的向量。角度可以通过 atan2 获得，也可以像我一样使用基向量。您可以通过点积测试正交性。最近点和最远点可以有超过 2 个点。在这种情况下，您应该找到每个组的中间以提高精度。
积分拟合椭圆

您需要首先找到椭圆点与y=cy 的角度，然后在这两个角度之间整合椭圆。另一半也一样，只是整合angles + PI。确定哪一半只是计算角度范围中间的点并根据y>=cy决定...

[Edit2] 这里更新了我为此破解的 C++ 代码：

    picture pic0,pic1,pic2;
        // pic0 - source
        // pic1 - output
    float a,b,a0,a1,da,xx0,xx1,yy0,yy1,ll0,ll1;
    int x,y,i,threshold=127,x0,y0,x1,y1,cx,cy,ax,ay,bx,by,aa,bb,dd,l0,l1;
    pic1=pic0;
    // bbox,center,recolor (white,black)
    x0=pic1.xs; x1=0;
    y0=pic1.ys; y1=0;
    for (y=0;y<pic1.ys;y++)
     for (x=0;x<pic1.xs;x++)
      if (pic1.p[y][x].db[0]>=threshold)
        {
        if (x0>x) x0=x;
        if (y0>y) y0=y;
        if (x1<x) x1=x;
        if (y1<y) y1=y;
        pic1.p[y][x].dd=0x00FFFFFF;
        } else pic1.p[y][x].dd=0x00000000;
    cx=(x0+x1)/2; cy=(y0+y1)/2;
    // fill inside (gray) left single pixel width border (thining)
    for (y=y0;y<=y1;y++)
        {
        for (x=x0;x<=x1;x++) if (pic1.p[y][x].dd)
            {
            for (i=x1;i>=x;i--) if (pic1.p[y][i].dd)
                {
                for (x++;x<i;x++) pic1.p[y][x].dd=0x00202020;
                break;
                }
            break;
            }
        }
    for (x=x0;x<=x1;x++)
        {
        for (y=y0;y<=y1;y++) if (pic1.p[y][x].dd) { pic1.p[y][x].dd=0x00FFFFFF; break; }
        for (y=y1;y>=y0;y--) if (pic1.p[y][x].dd) { pic1.p[y][x].dd=0x00FFFFFF; break; }
        }
    // find min,max radius (periaxes)
    bb=pic1.xs+pic1.ys; bb*=bb; aa=0;
    ax=cx; ay=cy; bx=cx; by=cy;
    for (y=y0;y<=y1;y++)
     for (x=x0;x<=x1;x++)
      if (pic1.p[y][x].dd==0x00FFFFFF)
        {
        dd=((x-cx)*(x-cx))+((y-cy)*(y-cy));
        if (aa<dd) { ax=x; ay=y; aa=dd; }
        if (bb>dd) { bx=x; by=y; bb=dd; }
        }
    aa=sqrt(aa); ax-=cx; ay-=cy;
    bb=sqrt(bb); bx-=cx; by-=cy;
    //a=float((ax*bx)+(ay*by))/float(aa*bb);    // if (fabs(a)>zero_threshold) not perpendicular semiaxes

    // separate/count upper,lower arc by horizontal line
    l0=0; l1=0;
    for (y=y0;y<=y1;y++)
     for (x=x0;x<=x1;x++)
      if (pic1.p[y][x].dd==0x00FFFFFF)
        {
        if (y>=cy) { l0++; pic1.p[y][x].dd=0x000000FF; } // red
        if (y<=cy) { l1++; pic1.p[y][x].dd=0x00FF0000; } // blue
        }
    // here is just VCL/GDI info layer output so you can ignore it...

    // arc separator axis
    pic1.bmp->Canvas->Pen->Color=0x00808080;
    pic1.bmp->Canvas->MoveTo(x0,cy);
    pic1.bmp->Canvas->LineTo(x1,cy);

    // draw analytical ellipse to compare
    pic1.bmp->Canvas->Pen->Color=0x0000FF00;
    pic1.bmp->Canvas->MoveTo(cx,cy);
    pic1.bmp->Canvas->LineTo(cx+ax,cy+ay);
    pic1.bmp->Canvas->MoveTo(cx,cy);
    pic1.bmp->Canvas->LineTo(cx+bx,cy+by);
    pic1.bmp->Canvas->Pen->Color=0x00FFFF00;
    da=0.01*M_PI;   // dash step [rad]
    a0=0.0;         // start
    a1=2.0*M_PI;    // end
    for (i=1,a=a0;i;)
        {
        a+=da; if (a>=a1) { a=a1; i=0; }
        x=cx+(ax*cos(a))+(bx*sin(a));
        y=cy+(ay*cos(a))+(by*sin(a));
        pic1.bmp->Canvas->MoveTo(x,y);
        a+=da; if (a>=a1) { a=a1; i=0; }
        x=cx+(ax*cos(a))+(bx*sin(a));
        y=cy+(ay*cos(a))+(by*sin(a));
        pic1.bmp->Canvas->LineTo(x,y);
        }

    // integrate the arclengths from fitted ellipse
    da=0.001*M_PI;      // integration step [rad] (accuracy)
    // find start-end angles
    ll0=M_PI; ll1=M_PI;
    for (i=1,a=0.0;i;)
        {
        a+=da; if (a>=2.0*M_PI) { a=0.0; i=0; }
        xx1=(ax*cos(a))+(bx*sin(a));
        yy1=(ay*cos(a))+(by*sin(a));
        b=atan2(yy1,xx1);
        xx0=fabs(b-0.0); if (xx0>M_PI) xx0=2.0*M_PI-xx0;
        xx1=fabs(b-M_PI);if (xx1>M_PI) xx1=2.0*M_PI-xx1;
        if (ll0>xx0) { ll0=xx0; a0=a; }
        if (ll1>xx1) { ll1=xx1; a1=a; }
        }
    // [upper half]
    ll0=0.0;
    xx0=cx+(ax*cos(a0))+(bx*sin(a0));
    yy0=cy+(ay*cos(a0))+(by*sin(a0));
    for (i=1,a=a0;i;)
        {
        a+=da; if (a>=a1) { a=a1; i=0; }
        xx1=cx+(ax*cos(a))+(bx*sin(a));
        yy1=cy+(ay*cos(a))+(by*sin(a));
        // sum arc-line sizes
        xx0-=xx1; xx0*=xx0;
        yy0-=yy1; yy0*=yy0;
        ll0+=sqrt(xx0+yy0);
//      pic1.p[int(yy1)][int(xx1)].dd=0x0000FF00; // recolor for visualy check the right arc selection
        xx0=xx1; yy0=yy1;
        }
    // lower half
    a0+=M_PI; a1+=M_PI; ll1=0.0;
    xx0=cx+(ax*cos(a0))+(bx*sin(a0));
    yy0=cy+(ay*cos(a0))+(by*sin(a0));
    for (i=1,a=a0;i;)
        {
        a+=da; if (a>=a1) { a=a1; i=0; }
        xx1=cx+(ax*cos(a))+(bx*sin(a));
        yy1=cy+(ay*cos(a))+(by*sin(a));
        // sum arc-line sizes
        xx0-=xx1; xx0*=xx0;
        yy0-=yy1; yy0*=yy0;
        ll1+=sqrt(xx0+yy0);
//      pic1.p[int(yy1)][int(xx1)].dd=0x00FF00FF; // recolor for visualy check the right arc selection
        xx0=xx1; yy0=yy1;
        }
    // handle if the upper/lower parts are swapped
    a=a0+0.5*(a1-a0);
    if ((ay*cos(a))+(by*sin(a))<0.0) { a=ll0; ll0=ll1; ll1=a; }
    // info texts
    pic1.bmp->Canvas->Font->Color=0x00FFFF00;
    pic1.bmp->Canvas->Brush->Style=bsClear;
    x=5; y=5; i=16; y-=i;
    pic1.bmp->Canvas->TextOutA(x,y+=i,AnsiString().sprintf("center = (%i,%i) px",cx,cy));
    pic1.bmp->Canvas->TextOutA(x,y+=i,AnsiString().sprintf("a = %i px",aa));
    pic1.bmp->Canvas->TextOutA(x,y+=i,AnsiString().sprintf("b = %i px",bb));
    pic1.bmp->Canvas->TextOutA(x,y+=i,AnsiString().sprintf("upper = %i px",l0));
    pic1.bmp->Canvas->TextOutA(x,y+=i,AnsiString().sprintf("lower = %i px",l1));
    pic1.bmp->Canvas->TextOutA(x,y+=i,AnsiString().sprintf("upper`= %.3lf px",ll0));
    pic1.bmp->Canvas->TextOutA(x,y+=i,AnsiString().sprintf("lower`= %.3lf px",ll1));
    pic1.bmp->Canvas->Brush->Style=bsSolid;

它使用我自己的图片类和成员：

xs,ys图片分辨率
p[y][x].dd 像素访问作为 32 位无符号整数作为颜色
p[y][x].db[4] 像素访问为 4*8bit 无符号整数作为颜色通道

您可以将picture::p 成员视为简单的二维数组
```
union color
    {
    DWORD dd; WORD dw[2]; byte db[4];
    int i; short int ii[2];
    color(){}; color(color& a){ *this=a; }; ~color(){}; color* operator = (const color *a) { dd=a->dd; return this; }; /*color* operator = (const color &a) { ...copy... return this; };*/
    };
int xs,ys;
color p[ys][xs];
Graphics::TBitmap *bmp; // VCL GDI Bitmap object you do not need this...
```
每个单元格可以作为 32 位像素p[][].dd 访问，0xAABBGGRR 或 0xAARRGGBB 现在不确定哪个。您也可以使用 p[][].db[4] 作为 8 位字节直接访问通道。

bmp 成员是 GDI 位图，因此bmp->Canvas-> 可以访问所有对您不重要的 GDI 内容。

第二张图片的结果如下：

灰色水平线是通过中心的圆弧边界线
红色、蓝色是弧的一半（计数时重新着色）
绿色是半轴基向量
Aqua dash-dash 是用于比较拟合的解析椭圆叠加。

如您所见，拟合非常接近（+/-1 像素）。计算的弧长upper,lower 非常接近近似的平均圆半周长（周长）。

您应该添加a0 范围检查来确定起点是上半部分还是下半部分，因为没有保证会找到主轴的哪一侧。两半的积分几乎相同，并且在积分步长0.001*M_PI 附近饱和，每个弧长只有17 和22 与直接像素计数的像素差异，甚至比我预期的要好到别名...

对于更偏心的椭圆，拟合效果不是很好，但结果仍然足够好：

【讨论】：

@Spectre，是的，我得到了类似的结果，但这是个问题。我从你那里得到了 1/3 的上弧和下弧。这里知道图像正好在 x=320 并且 y 是变体，图像在 y 范围内移动。我使用 x 值来检测弧。而半长轴和半短轴由于图像错误而无法计算，如果计算意味着它的放置错误。感谢您的努力。
@BalaKumaran 您的图像是 JPG，因此如果您直接与 0x00FFFFFF 或 0x00000000 进行比较，那么由于 JPG 压缩，您会忽略所有带有 dstorted 颜色的点，因此会出现错误...尝试通过二值化阈值因此将所有具有r+g+b>=3*128 的像素设置为白色，其余像素设置为黑色，除非您在某处遇到其他错误，否则应该可以解决问题...在我的示例中，我将单个通道的阈值设置为> = 127，但想法是相同的...
我也尝试了灰度图像，请提供有关您的图片类的一些详细信息。我尝试编译您的代码，它在 p[y][x].db[]>0 中引发错误。和p[y][x].dd
@Spectre 添加关于图片类的代码，我试试你的代码它在图像读取时出错，P[][] 和 db[] 声明提前谢谢。
@BalaKumaran 添加了我的图片类中p[][] 的集成和含义。它自己的课程很大，包含许多对于这项任务不重要的事情，抱歉懒得把它删掉......相反，只需将我的图像访问权限更改为您的。唯一重要的是直接像素访问，所以它不应该太难......而且我的 nik 是 Spektre 而不是 Spectre 幸运的是，当你评论我的答案而不是你的问题时，通知仍然有效......