How to find the length of the upper and lower arcs from an ellipse image

Question

Here I am trying to find the upper arc and lower arc using the image vector (image contours), but it could not give the result of the extraction. Suggest any other method for finding the upper and lower arcs from images and their length.

Here is my code

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; }} } 

Above the code, I get the coordinates using coordinate points. I can find the length of the arc, but its inconsistency with the result of the extraction.

Here is the actual image:

Image1:

After thinning, I got:

Expected Result:

Answer 1

Since you cannot determine what the upper / lower arc is, I will assume that you cut the ellipse in half with a horizontal line through the midpoint of the ellipse. If this is not the case, you need to adapt it yourself ... Well, how to do it:

• binarization image

  When you provide jpg, the colors are distorted, so there is more to it than just black and white.

• reduce the border to 1 pixel

  Fill the inside with white, and then redraw all the white pixels that are not adjacent to any black pixels on some unused or black color. There are many other options how to achieve this ...

• find bounding box

  find all the pixels and remember the minimum, max x,y coordinates of all white pixels. Call them x0,y0,x1,y1 .

• ellipse computing center

  just find the midpoint of the bounding box

   cx=(x0+x1)/2 cy=(y0+y1)/2 

• count pixels for each elliptical arc

  have a counter for each arc and simply increase the counter of the upper arc for any white pixel with y<=cy and lower if y>=cy . If your coordinate system is different, then the conditions may be reversed.

• find ellipse parameters

  just find the white pixel closest to (cx,cy) , this will be the end point of the secondary semiaxis b , call it (bx,by) . Also find the farthest white pixel to (cx,cy) , which will become the main endpoint of the semi-axis (ax,ay) . The distance between them and the center will give you a,b , and their position, subtracted by the center, will give you vectors with the rotation of your ellipse. the angle can be obtained using atan2 or use basis vectors, like me. You can check orthogonality using a point product. For the nearest and closest point, there may be more than 2 points. in this case, you should find the middle of each group to increase accuracy.

• Integrated Ellipse Integration

  You need to first find the angle at which the points of the ellipse with y=cy , then integrate the ellipse between these two angles. The other half is the same as angles + PI integration. To determine which half is just to calculate the point in the middle between the range of angles and decide according to y>=cy ...

[Edit2] Here the C ++ code is updated, which I cheated for this:

  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; 

It uses my own image class with members:

• xs,ys image resolution
• p[y][x].dd pixel access as a 32-bit unsigned integer as color

• p[y][x].db[4] pixel access as an unsigned integer 4 * 8 bits as color channels

  You can see the picture::p member as a simple 2D array

   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... 

  where each cell can be accessed as a 32-bit pixel p[][].dd , since 0xAABBGGRR or 0xAARRGGBB now not sure what. You can also directly access channels using p [] []. Db [4] as 8-bit BYTE.

  The bmp element has a GDI bitmap, so bmp->Canvas-> get access to all GDI materials that are not important to you.

Result for the second image:

• The gray horizontal line is the boundary line of the arc through the center.
• Red, blue represent the halves of the arc (redrawn when counting)
• Green are semi-axial base vectors
• Aqua Dash Dash is an analytical ellipse overlay for matching comparisons.

As you can see, the fit is pretty close (+/- 1 pixel). The calculated arc length upper,lower quite close to the approximate average circumference of the circle (circle).

You must add a range check a0 to decide whether the start is the upper or lower half, because there is no guarantee which side of the main axis will be found. The integration of both halves is almost the same and saturated around the integration step of 0.001*M_PI around 307.3 pixels per arc length, which is equal to only 17 and 22 pixel differences from direct pixel counting, which is even better than what I expect to smooth ...

For more eccentric ellipses, the fit is not so good, but the results are still good enough: