271 views
Practical: 2D Estimation == > Tristan Roussillon In this practical, we will see how one can get several estimates of geometrical quantities from the contour of a 2D digitized shape: curvatures values, but also unit normal vectors and, as a by-product, the contribution of every contour element to the total perimeter. Several approaches are implemented in [DGtal](https://dgtal.org) to get such estimates. In what follows, we will focus on a purely geometrical approach, based on maximal digital circular arcs (MDCA for short) [^1]. The objective is to write a program to check if the integral of MDCA-based curvature estimates along the contour of a digitized shape tends towards $2\pi$ as the grid step tends towards $0$. It should be, because the MDCA-based curvature estimation has been proven to be pointwise multigrid-convergent [^2]. # Getting Started In the [DGtal-DGMM](https://github.com/DGtal-team/DGtal-Tutorials-DGMM2022) repository, have a look to the `2D-estimation-template.cpp` file. It should compile as is. To run it, just type in: ``` ./2D-estimation-template ``` This program discretizes an ellipse and retrieves the boundary points. Note that the grid step is read from the command line and can be given as an input argument ($0.1$ is the default value), for example: ``` ./2D-estimation-template 0.05 ``` Note that in this template, we use the namespace `Z2i`. It gathers several useful types when working in the digital space $\mathbb{Z}^2$ or its associated cubical complex $\mathbb{F}^2$. :::success :::spoiler `Z2i` types to DGtal classes `Z2i` also simplifies the names for conveniency. For that reason, it is not always easy to find the corresponding class in the documentation. That's why we provide below this small dictionnary: - Space: SpaceND - K2: KhalimskySpaceND - Point, Vector, RealPoint, RealVector: PointVector - Curve: GridCurve ::: :::success :::spoiler Digital space and Khalimsky space In digital images, pixels are often identified by a couple of integer coordinates, i.e., an element of the digital space $\mathbb{Z}^2$. Any finite set of pixels may represent a region in an image. If one is interested in the interface between two regions, it is often convenient to represent not only the pixels, i.e., the unit squares, but also their edges and vertices. In DGtal, we therefore also consider the associated cubical complex $\mathbb{F}^2$, which models cells of different dimension: - 2-dimensional cells (closed unit square) = pixels - 1-dimensional cells (closed unit segment) = linels or grid edge - 0-dimensional cells (closed point) = pointels or grid point A Khalimsky space is just an efficient representation of the cubic complex $\mathbb{F}^2$ by couples of integers. It is heavily used is DGtal, e.g., to track the topological border of a region. See the page [Cellular grid space and topology](https://dgtal.org/doc/stable/moduleCellularTopology.html) for more details. ::: # Grid Curve A first thing to add is the creation of a grid curve: ``` c++ Curve gridcurve( kspace ); gridcurve.initFromVector( points ); ``` If your curve is not closed, it is best to return because something's going wrong. Check the documentation of the class [GridCurve](https://dgtal.org/doc/stable/classDGtal_1_1GridCurve.html) to find out how one can check whether a curve is closed or not. Then you can see what your curve looks like using [Board2D](https://dgtal.org/doc/stable/moduleBoard2D.html): ```c++ Board2D board; // draw the grid curve onto the board board << SetMode( "PointVector", "Grid" ) << gridcurve; // save the drawing board.saveSVG("gridcurve.svg"); ``` You can view the svg file in your favorite browser. You will see that your grid curve is displayed as a sequence of 1-dimensional cells (also called grid edges). A grid curve is actually internally represented in that way but can offer other views, as shown in the following. <img alt="grid curve" src="https://perso.liris.cnrs.fr/tristan.roussillon/DGtal-Tuto/gridcurve.svg" width="250"> # Estimation In this section, we will consider the 2-dimensional cells that are incident to every grid edges and that are conveniently represented as pairs of integer points. <img alt="incident points" src="https://perso.liris.cnrs.fr/tristan.roussillon/DGtal-Tuto/incident-points.svg" width="250"> That range can be obtained from the grid curve as follows: ```c++ using Range = Curve::IncidentPointsRange; Range range = gridcurve.getIncidentPointsRange(); ``` Note that, since the curve is closed, we will use circulators: ```c++ using Iterator = Range::ConstCirculator; ``` The MDCA approach is based on a routine that searches for the longest circle that separate two point sets, lying on either side of the grid curve. You can have a look at the function `getEstimates`. It is based on the class [StabbingCircleComputer](https://dgtal.org/doc/stable/classDGtal_1_1StabbingCircleComputer.html#details), but also on a functor, which determines which geometrical quantity must be returned from a separating circle. The goal of this section is to call `getEstimates` with the appropriate arguments. For curvature values, you must used: ```c++ using CFunctor = CurvatureFromDCAEstimator<StabbingCircleComputer<Iterator>,false>; CFunctor cf; ``` For normal values, you must used: ```c++ using NFunctor = NormalFromDCAEstimator<StabbingCircleComputer<Iterator> >; NFunctor nf; ``` :::success :::spoiler Need more help? ```c++ std::vector<double> curvatures; using CFunctor = CurvatureFromDCAEstimator<StabbingCircleComputer<Iterator>,false>; CFunctor cf; getEstimates(range.c(), range.c(), cf, std::back_inserter(curvatures), h); ``` You can have similarly the normal vectors as `RealVector` instances. ::: You could now print to the standard output your estimates and use your favorite plot tool to get the curvature graph. ![incident points](https://perso.liris.cnrs.fr/tristan.roussillon/DGtal-Tuto/curvature-graph.svg) # Integral of curvature In this section, you need first to evaluate the contribution of every grid edge to the total perimeter, which can be done by computing the dot product between the previously estimated normal vector and the unit vector orthogonal to the grid edge. You have now all the required material to compute the integral of curvature along the grid curve. You could compare your results with ours: | h | 0.01 | 0.005 | 0.001 | | -------------- | -------- | -------- | -------- | | curv. int. | 6.27058 | 6.28821 | 6.28 | We indeed experimentally observe a convergence towards $2\pi$ as the grid step tends towards $0$. # References [^1]: Roussillon, T., Lachaud, JO. (2011). Accurate Curvature Estimation along Digital Contours with Maximal Digital Circular Arcs. In: Aggarwal, J.K., Barneva, R.P., Brimkov, V.E., Koroutchev, K.N., Korutcheva, E.R. (eds) Combinatorial Image Analysis. IWCIA 2011. Lecture Notes in Computer Science, vol 6636. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-21073-0_7 [^2]: Schindele, A., Massopust, P. & Forster, B. Multigrid Convergence for the MDCA Curvature Estimator. J Math Imaging Vis 57, 423–438 (2017). https://doi.org/10.1007/s10851-016-0685-1.