Get best of the week

Using the OpenCV library to recognize elliptical arcs in 2D sections of 3D point clouds

In connection with the wider spread of affordable laser scanners (lidars) capable of receiving 3D point clouds ( 3dOT ) and the wider application of this technology in various fields (from mechanical engineering to security, from the oil industry to architecture), interest in processing algorithms has revived clouds of points.

One of the popular applications of 3d in industry is the creation of design documentation for only constructed, old or converted equipment, which usually consists of pipelines and other structures of cylindrical geometry.

To detect geometric primitives in 3dOT , specialized 3D libraries, for example, Microsoft PCL , are usually used. The approach with the use of ready-made libraries, along with advantages, has disadvantages. For example, it is difficult to incorporate them into existing Kadov processing schemes, which usually have a 2D dimension.

Let's consider how it would be possible to process 3dOT , for example, a pumping station, starting with 2D sections and using the whole arsenal of 2D processing, which is available in reliable and optimized image processing libraries, for example OpenCV .


Figure 1. 3D OT model of a pumping station

The main element of the sections obtained by scanning various pipe structures are elliptical arcs .


Figure 2. Horizontal cross section of a 3D model of a pumping station at an average level.

For this article, we restrict our consideration to one key algorithm that allows us to detect arbitrary elliptical arcs - this is an iterative algorithm for the growth of arc segments and region binding ( region growth and edge linking ).

Growth algorithms are the most obvious and easily verifiable, albeit time-consuming, in comparison with statistical algorithms, which are better suited to the case when the scene contains loosely coupled, distant objects that belong to one ellipse. These algorithms will be discussed in future articles.


For now, for simplicity, we omit the procedure for obtaining a section from the source 3dOT file , preprocessing a section, clustering it to isolate geometric primitives, as well as subsequent binding, rectification, and other photogrammetric operations needed to obtain model parameters. We will not discuss the parameterization of heuristic search algorithms in the same way. Let us describe all the basic operations from which the algorithm is built.

We assume that we need to detect (recognize, classify) an elliptical arc (that is, calculate the parameters of the ellipse, as well as the initial and final angle of the elliptical arc) in this image, cut out from the horizontal section of the point cloud.


Figure 3. One of the elliptical arcs of the cross section of the 3D model (after smoothing)

In order to minimize the work with the raster blindly, we will carry out all operations with the raster through outlining .

The OpenCV findContours procedure finds on the raster mat all external (without internal shapes) contours in the form of a vector of integer point vectors (in raster coordinates):
 Mat mat(size);
 vector<vector<Point>> contours;
 findContours(mat, contours, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE);

This is our key operation, which in some simple cases completely solves the task . But since degenerate cases are not always found, let us consider in more detail the processing technology through contouring.

The reverse operation, generating a raster according to the existing external circuit using the OpenCV function, also looks simple:
 drawContours(mat, contours, -1, Scalar(255), -1);

It is also very often used to mask contours, draw, or to calculate area.

Thus, at the initial stage, we have a set of patches (pieces of a certain curve) that need to be connected into an elliptical arc, eliminating parts of other components of the structure (for example, fasteners) or optical noise caused by shadowing during scanning and other reasons.

Let's create a discriminant function that will return the type of contour (ellipse, linear segment, hatching or something else), as well as the end points of the contour and its rotated outline rectangle:
 contourTypeSearch(
   const vector<Point> &contour, Vec4i &def, RotatedRect &rc);

The ratio of the length and width of the rectangle helps to quickly discriminate contours close to linear segments , as well as small noise contours .

The rotated rectangle in OpenCV has a complex coordinate system. If it is not the angle itself that is required, but its trigonometric functions , it is all the more less obvious from the context. If the absolute value of the angle is used , it must be taken into account that the angle is counted from the horizontal to the first edge of the rectangle counterclockwise and has a negative value .

The endpoints of the elliptical contours are found using our procedure, which receives the mat rasterwith a discriminated outline extracted from the original image by masking and returns the maximum defect :
 contourConvFeature(mat, &def, … );

The main code for this function is to call two OpenCV procedures :

 vector<int> *hull = new vector<int>();
 convexHull(contour, *hull);
 vector<Vec4i> *defs = new vector<Vec4i>();
 convexityDefects(contour, *hull, *defs);

The first procedure finds a convex polygon for the contour under study, the second - calculates all the convexity defects .

We take only the largest defect in terms of convexity, considering that it determines the end points of the contour. This may not be the case if the external or internal boundaries of the contour have features . In order to smooth them , we apply additional smoothing to the contour under study (and not to the entire image so as not to “blur” the isthmuses between the contours and not violate the original topology).


Figure 4. Calculation of the bulge defect.

Option (a) mistakenly defines the red endpoint. Option (b)correctly defines endpoints. Option (c) redefines the endpoints on the original shape.

Since in the technology we adopted, the circuit is regenerated each time , we have to re-search for the points of correspondence (or rather, their indices) by the exhaustive search procedure :
 nearestContourPtIdx(const vector<Point> &contour, const Point& pt);

For cases when it is not possible to completely get rid of the features, an additional mode of arc separation was also implemented (work separately with the internal / external arc). This is important, for example, in cases where the external arc of the contour is in contact with other objects or is noisy . In this case, you can go to work with the internal arc. In this case, it is not necessary to process the external and internal arcs separately.

Also, according to the well-known formula for the ratio of the convexity of the arc , the radius of the circle is approximately estimated and too large ellipses are rejected:
R = bulge / 2 + SQR(hypot) / (8 * bulge);

Thus, for all the contours, their convexity defect metric was found (or they are classified as linear or small and removed from the procedure). At the last stage, additional parameters are added to the original metric, such as the rotated dimension parameter, etc., and the full set of metrics under study is ordered by size .
 typedef tuple<int , // 
   RotatedRect, //  
   Vec4i, //  
   int> // 
   RectDefMetric;


Algorithm for linking arc segments at end points


The growth algorithm is clear and obvious: we take the largest contour as a seed and try to grow it, that is, find and attach the nearest patches to its end points that satisfy the growth conditions. In the grown figure, we enter the desired elliptical arc . Mask and subtract the figure from the original set. We repeat the growth procedure until the initial set runs out .

The basic procedure of the growth algorithm looks like this:
 vector<Point> *patch =
    growingContours(contour, def, tmp, hull);

where contour is the contour under study, def is its convexity defect, hull is the convex polygon of the entire region, tmp is the auxiliary buffer matrix. At the output, we get a vector grown contour.

The procedure consists of a cycle of attempts to seed growth, ending either by exhausting the available patches for growth, or limited by the parameter of the maximum number of iterations .


Figure 5. Many patches for growth without seed The

main difficulty is to select the nearest patches to the end points of the contour, so that the figure grows only forward . For tangential directionwe take the averaged line belonging to the arc in the vicinity of the endpoint. In Figure 6 displays the candidates for connection to the seed at a certain iteration.


Figure 6. Seed surrounded by a plurality of growth candidate patches.

For each candidate patch, the following metric is calculated:
typedef tuple<
   double, //    2      2   
   bool, bool, //,  4   
   int, // 
   Vec4i> //  
   DistMetric;

Only patches that fall into the tangential cone are taken into account . Then the patch with the smallest distance is selected and, by imprinting the connecting section into the raster, connects to the corresponding end of the seed. For the other end of the seed, a patch matching the parameters is searched, and if found, is also connected to the seed. Then the seed is masked and subtracted from the many patches. The procedure is repeated from the very beginning.

At the end of the growth procedure, we got an elliptical arc , which remains to be verified.

To get started, using the standard OpenCVthe procedure that our patch receives (in the form of a path, we remember that the path and raster are interchangeable with us) and returns the rotated dimension , that is, a complete ellipse.
 RotatedRect box = fitEllipse(patch);

Then, we reject too large and too small ellipses, and then we apply our original procedure for comparing the areas of the resulting elliptical arc and the initial growth patch in raster form . This procedure includes some disguised tricks, so we will omit its description for now.

And finally, we find the remaining parameters of the detected ellipse - the start and end angles (we already know the half-axes from fitEllipse ).

To determine the starting and ending angles, we proceed as follows: we transform our complete ellipse, back to the polygon, and by direct enumeration we find its points closest to our ends. Their angular coordinates(in fact, indices) and will be the starting and ending angles of an elliptical arc. In code, it looks like this (a bit simplified):
 pair<int, int>
   ellipseAngles(const RotatedRect &box,
   vector<Point> &ell, const Point &ps, 
   const Point &pe, const Point &pm) 
 {
    vector<Point> ell0;
    ellipse2Poly(Point(box.center.x, box.center.y), 
      Size(box.size.width / 2, box.size.height / 2),
      box.angle, 0, 355, 5, ell0);
    int i0 = nearestContourPtIdx(ell0, ps);
    int i1 = nearestContourPtIdx(ell0, pe);
    cutSides(ell0, i0, i1, i2, ell, nullptr);
    return pair<int, int>(i0, i1);
}

Our cutSides procedure takes into account the elliptical arc traversal topology . In total, eight possible cases of bypassing the indices i0, i1, i2 should be considered . Do we go along the outer contour or the inner one, and which of the indices is greater, initial or final?

Easier to see code:
 void cutSides(
   const vector<Point> &ell0, int i0, int i1, int i2, 
   vector<Point> *ell_in, vector<Point> *ell_out)
 {
   if (i0 < i1) {
      if (i2 > i0 && i2 < i1) {
         if (ell_in) {...}
            if (ell_out) {...}
        } else {
            if (ell_in) {...}
            if (ell_out) {...}
        }}
    else {
        if (i2 > i1 && i2 < i0) {
            if (ell_in) {...}
            if (ell_out) {...}
        } else {
            if (ell_in) {...}
            if (ell_out) {...}
        }}}

Some results of detecting ellipses in complex cases are shown in Figure 7 .



In the following articles, statistical detection methods will be considered.

More articles:


All Articles