1

Here i try to find the upper arc and lower arc using image vector(contours of images) But It could n't gave Extract result. Suggest any other method to find upper and lower arc from images and their length.

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

By Above code I get Coordinates, by using Coordinates points I can find the length of arc but its mismatch with extract result.

Here is actual image:

Image1:

enter image description here

After thinning i got:

enter image description here

Expected Output:

enter image description here

Bala Kumaran
  • 131
  • 3
  • 14
  • 1
    add specification what you mean by upper/lower arc (you mean you cut the ellipse in halves by horizontal line going through ellipse middle point?) also add the expected result and the incorrect output of yours – Spektre Dec 21 '15 at 09:46
  • @Spektre Spektre I can't divides into two halves because upper arc is slightly bigger than lower arc due to that I need to find values of two arc. And then Images could n't come in extract position . Its looks like ellipse but not a perfect ellipse, so i unable to use any formula regards to ellipse properties . Here i need to find its ellipse or not. Is there any other way?. thanks for suggestion. – Bala Kumaran Dec 21 '15 at 12:02
  • well what exactly are the arc then how are they defined where are the boundary (periaxises perhaps)? If you just want to know if it is ellipse or not then find (avg) mid point. then the farest and closest points will give you peri-axises and then just compute the avg distance of your points from that analytical ellipse. the bigger the distance the far is this off an ellipse shape. Also if the periaxises are not perpendicular then it is skewed ... so also not an ellipse – Spektre Dec 21 '15 at 13:14
  • btw your ellipse looks like normaly rendered circle with diameter `D=212` pixels without any deviations so where is the problem? – Spektre Dec 21 '15 at 13:20
  • @Spektre yes, It's gave Diameter `D=(210-212)`. check with second image(image2) its not a perfect ellipse.but it so gave similar value with image1 ellipse. How to differentiate those images(not with size). – Bala Kumaran Dec 22 '15 at 04:29
  • Your last image is pretty confusing you really have to declare the definition of upper and lower arc you want to measure because it is clearly not obvious.(may be the border is slope +/-1? but that is just speculation on my side) – Spektre Dec 23 '15 at 08:51
  • @Spectre I finds the upper and lower arc using contours and `x=320` (midpoint of image not ellipse),and i tried to find contours for y values in midpoint, by using contours midpoint value (320,154 and 320,355) i could finds the upper arc and lower arc by adding y values with semi major/2. but upper arc and lower arc are same due to camera optics. – Bala Kumaran Dec 23 '15 at 09:23
  • that is really unclear. ... mid point of image is `(320,240)` ... `x=320` is vertical line splitting image to left and right sides not up/down. Also your ellipse is not centered in mid of image instead it is centered around `(317,254)`. How did you came up with the horizontal lines? where they come from ? What does optics does with this? – Spektre Dec 23 '15 at 10:58
  • @spectre Check my actual image(image1). – Bala Kumaran Dec 24 '15 at 05:05

1 Answers1

2

As you are unable to define what exactly is upper/lower arc then I will assume you cut the ellipse in halves by horizontal line going through the ellipse's middle point. If that is not the case then you have to adapt this on your own... Ok now how to do it:

  1. binarize image

    As you provide JPG the colors are distorted so there is more then just black and white

  2. thin the border to 1 pixel

    Fill the inside with white and then recolor all white pixels not neighboring any black pixels to some unused or black color. There are many other variation how to achieve this...

  3. find the bounding box

    search all pixels and remember min,max x,y coordinates of all white pixels. Let call them x0,y0,x1,y1.

  4. compute center of ellipse

    simply find middle point of bounding box

    cx=(x0+x1)/2
    cy=(y0+y1)/2
    
  5. count the pixels for each elliptic arc

    have counter for each arc and simply increment upper arc counter for any white pixel that have y<=cy and lower if y>=cy. If your coordinate system is different then the conditions can be reverse.

  6. find ellipse parameters

    simply find white pixel closest to (cx,cy) this will be endpoint of minor semi-axis b let call it (bx,by). Also find the most far white pixel to (cx,cy) that will be the major semi axis endpoint (ax,ay). The distances between them and center will give you a,b and their position substracted by center will give you vectors with rotation of your ellipse. the angle can be obtained by atan2 or use basis vectors as I do. You can test ortogonality by dot product. There can be more then 2 points for closest and farest point. in that case you should find the middle of each group to enhance precision.

  7. Integrate fitted ellipse

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

[Edit2] Here updated C++ code I busted 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 use my own picture class with members:

  • xs,ys resolution of image
  • p[y][x].dd pixel access as 32bit unsigned integer as color
  • p[y][x].db[4] pixel access as 4*8bit unsigned integer as color channels

    You can look at picture::p member as simple 2D array of

    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 32 bit pixel p[][].dd as 0xAABBGGRR or 0xAARRGGBB not sure now which. Also you can access the channels directly with p[][].db[4] as 8bit BYTEs.

    The bmp member is GDI bitmap so bmp->Canvas-> access all the GDI stuff which is not important for you.

Here result for your second image:

example

  • Gray horizontal line is the arc boundary line going through center
  • Red,Blue are the arc halves (recolored during counting)
  • Green are the semi-axes basis vectors
  • Aqua dash-dash is analytical ellipse overlay to compare the fit.

As you can see the fit is pretty close (+/-1 pixel). The counted arc-lengths upper,lower are pretty close to approximated average circle half perimeter(circumference).

You should add a0 range check to decide if the start is upper or lower half because there is no quarantee which side of major axis this will find. The integration of both halves are almost the same and saturated around integration step 0.001*M_PI around 307.3 pixels per arc-length which is only 17 and 22 pixels difference from the direct pixel count which is even better then I anticipate due to aliasing ...

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

example

eddie
  • 1,252
  • 3
  • 15
  • 20
Spektre
  • 49,595
  • 11
  • 110
  • 380
  • @Spectre, yes man I got similar result but that's a problem. I get upper arc and lower arc of 1/3th from yours. here know of image is exactly at x=320 and y is variant,images moves over y range. And I use x value to detect arc. And semi major and semi minor axis are unable to calculate due to their wrong images,if compute means its lay on error. Thanks for yours great efforts. – Bala Kumaran Dec 22 '15 at 09:33
  • @BalaKumaran your image is JPG so if you are comparing directly to `0x00FFFFFF` or `0x00000000` then you ignore all points with dstorted color due to JPG compression hence the error ... try to binarize by treshold so set all pixels which have `r+g+b>=3*128` as white and the rest as black and that should do the trick unless you got some other bug somewhere ...In my example I treshold just single channel >= 127 but the idea is the same ... – Spektre Dec 22 '15 at 10:29
  • I tried with grayscale images also,Please put some details about your pictures class. i tried compile your code it raise error in `p[y][x].db[]>0` . and `p[y][x].dd` – Bala Kumaran Dec 22 '15 at 10:35
  • @Spectre add code about picture class, I try your code it got error in image reading,P[][] and db[] declaration thanks in advance. – Bala Kumaran Dec 22 '15 at 10:55
  • @BalaKumaran Added the integration and also the meaning of `p[][]` from mine picture class. The class it self is quite big containing many unimportant things for this task sorry too lazy to cut it dow ... instead just change mine image access to yours. The only thing important is the direct pixel access so it should not be too hard... Also my nik is `Spektre` not `Spectre` The notification luckily works anyway as you comment mine answer instead of your question ... – Spektre Dec 22 '15 at 11:22
  • @BalaKumaran BTW the picture class is based on this [GDI Bitmap ScanLine](http://stackoverflow.com/a/21699076/2521214) – Spektre Dec 22 '15 at 11:28
  • @BalaKumaran reedited answer ... removed old code and added the new with hopefully correct integration now. btw for your almost circular ellipse the result was different about `0.001` pixel from the buggy integration... – Spektre Dec 22 '15 at 12:25