diff --git a/ChangeLog.md b/ChangeLog.md index dbfd6ff6b2..1c9f56aa4e 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,6 @@ # DGtal 1.3 (beta) + ## New features / critical changes - *Docker File* - A Dockerfile is added to create a Docker image to have a base to start development diff --git a/src/DGtal/geometry/doc/images/Alcapone3D.png b/src/DGtal/geometry/doc/images/Alcapone3D.png new file mode 100644 index 0000000000..dfb6204f6e Binary files /dev/null and b/src/DGtal/geometry/doc/images/Alcapone3D.png differ diff --git a/src/DGtal/geometry/doc/images/Alcapone3D_reduced.png b/src/DGtal/geometry/doc/images/Alcapone3D_reduced.png new file mode 100644 index 0000000000..378f33a222 Binary files /dev/null and b/src/DGtal/geometry/doc/images/Alcapone3D_reduced.png differ diff --git a/src/DGtal/geometry/doc/images/flower_5.png b/src/DGtal/geometry/doc/images/flower_5.png new file mode 100644 index 0000000000..cc9c7c1131 Binary files /dev/null and b/src/DGtal/geometry/doc/images/flower_5.png differ diff --git a/src/DGtal/geometry/doc/images/flower_5_reduced.png b/src/DGtal/geometry/doc/images/flower_5_reduced.png new file mode 100644 index 0000000000..d1b8376ce3 Binary files /dev/null and b/src/DGtal/geometry/doc/images/flower_5_reduced.png differ diff --git a/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.1.png b/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.1.png new file mode 100644 index 0000000000..ea0df51bac Binary files /dev/null and b/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.1.png differ diff --git a/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.3.png b/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.3.png new file mode 100644 index 0000000000..a187d7c867 Binary files /dev/null and b/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.3.png differ diff --git a/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.7.png b/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.7.png new file mode 100644 index 0000000000..0900dc7648 Binary files /dev/null and b/src/DGtal/geometry/doc/images/scale_Alcapone3D_1.7.png differ diff --git a/src/DGtal/geometry/doc/images/scale_Alcapone3D_4.png b/src/DGtal/geometry/doc/images/scale_Alcapone3D_4.png new file mode 100644 index 0000000000..ebeb003e3a Binary files /dev/null and b/src/DGtal/geometry/doc/images/scale_Alcapone3D_4.png differ diff --git a/src/DGtal/geometry/doc/moduleMedialAxis.dox b/src/DGtal/geometry/doc/moduleMedialAxis.dox new file mode 100644 index 0000000000..e9f9c29e1c --- /dev/null +++ b/src/DGtal/geometry/doc/moduleMedialAxis.dox @@ -0,0 +1,513 @@ +/** + * @file moduleMedialAxis.dox + * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * + * @date 15/06/2021 + * + * Documentation file for feature medialAxis files + * + * This file is part of the DGtal library. + */ + +/* + * Useful to avoid writing DGtal:: in front of every class. + * Do not forget to add an entry in src/DGtal/base/Config.h.in ! + */ +namespace DGtal { +//---------------------------------------- +/*! +@page moduleMedialAxis Discrete Medial Axis and its filtration + +[TOC] + +@author Robin Lamy, David Coeurjolly, Isabelle Sivignon + +This part of the manual describes the DGtal Medial Axis module. We +detail how to extract a discrete medial axis @f$ MA(X)@f$ (set of +Euclidean balls) of a given binary object @f$ X @f$. There are many +definitions of such strucutre, both in continuous and discrete +settings, please refer to @cite TODO and @cite TODO for an overview. + + +In our discrete setting, the discrete medial axis will be defined such that @f[ +X = \cup_{b\in MA(X)} b\cap \mathbb{Z}^d\,.@f] Furthermore, all balls +in @f$ MA(X)@f$ are @e maximal @e in the sense that for all @f$ b\in +MA(X)@f$ there is no @f$ b'\in MA(X) @f$ such that @f$ b\subset b@f$. + + + + +\section Medial Axis Introduction + +The Medial Axis of a digital binary shape is a subset +of points from the original shape. Coming from a PowerMap, we only keep +the subset of points which associated balls is not fully contained into +another subset of balls. +(Hence there is no loss of informations when time of reconstruction comes) +The Medial Axis as been used in many situations such as shape analysis, shape matching, +shape-based interpolation, motion planning, image registration, or +differential measurement estimation. + +The Reduced Medial Axis algorithm is available in all dimensions + +\section Reduced Medial Axis + +For now, we get the Reduced Medial Axis from a +static function of the ReducedMedialAxis struct. + +In DGtal, the reduced medial axis is parametrized by the following objects : + +Template : +- A PowerMap Type + +Arguments : +- A powermap + +@code + typedef PowerMap< ... > PMap; + typedef ReducedMedialAxis RMA; + + PMap power( ... ); // A power map + + RMA::Type rdma = RMA::getReducedMedialAxisFromPowerMap(power); +@endcode + +This provides us the following 2D images : + +@image html flower_5.png +@image latex flower_5.png + +@image html flower_5_reduced.png +@image latex flower_5_reduced.png + +And 3D images : + +@image html Alcapone3D.png +@image latex Alcapone3D.png + +@image html Alcapone3D_reduced.png +@image latex Alcapone3D_reduced.png + +\section Reduced Medial Axis Filter + +The Reduced Medial Axis is really sensitive to the shape border, +and noise might make it a lot more furnished that it could be. +Hence, we might want to filter it in order to trade a little bit +of details on the border for simplicity. +Of course, the reconstruction will not be exactly the same, but +with a well-choosen parameter we can keep the original shape and +simplify a lot our set of balls. + +\section Scale Medial Axis + +In order to compute the Scale Medial Axis, we start from a +Distance Transformation (refers to @ref distanceTransformation) +of the considered shape. The goal is to multiply all balls radius +by a coefficient alpha, then compute the Reduced Medial Axis on +this set of balls, and finally divide all the Reduced Medial Axis +balls radius by alpha. It helps to erase the shape border imperfections. + +In DGtal, the Scale Medial Axis is parametrized by the following objects : + +Template : +- A Distance Transformation Type +- A Power Separable Metric Type + +Arguments : +- A distance transformation +- coefficient alpha (double) +- A power separable metric + +@code + typedef DistanceTransformation< ... > DTL2; + typedef ScaledMedialAxis SMA; + + DTL2 dt( ... ) // A distance transformation + double alpha = 1.3 + + SMA::Type sma = SMA::getScaleAxisFromWeightedMap(dt, alpha, Z2i::l2PowerMetric); + +@endcode + +Examples of Scale Medial Axis with different alpha : + +Original Medial Axis (alpha = 1) : + +@image html Alcapone3D_reduced.png +@image latex Alcapone3D_reduced.png + +alpha = 1.1 : + +@image html scale_Alcapone3D_1.1.png +@image latex scale_Alcapone3D_1.1.png + +alpha = 1.3 : + +@image html scale_Alcapone3D_1.3.png +@image latex scale_Alcapone3D_1.3.png + +alpha = 1.7 : + +@image html scale_Alcapone3D_1.7.png +@image latex scale_Alcapone3D_1.7.png + +alpha = 4 : + +@image html scale_Alcapone3D_4.png +@image latex scale_Alcapone3D_4.png + +@warning : According to the shape considered, a relatively too big +coefficient alpha might completely destroy the shape, especially +if there are thin and close details that you want to keep. + +\section Lambda Medial Axis + +In order to filter what can be consider secondary ball, you can +use the Lambda Medial Axis Filter, which is based on the idea +that the further its contact points with outside are, the more +important the ball becomes. Hence, we filter the balls computed by +the reduced medial axis by the size of the ball enclosing its contact +points. + +Here are example of the shapes computed : + + + + + + + + + + + + +In DGtal, we have chosen to implement such volumetric tools such that +the underlying metric could be specified independently. As an +example, we illustrate the distance maps from a single source point +for various metrics in 2D using the generic DistanceTransformation method: + +@image html DTmetrics-small.png "Distance maps for metrics L_1, L_2, L_4, L_80 and Chamfer masks M_{3−4} and M_{5−7−11}." +@image latex DTmetrics-small.png "Distance maps for metrics L_1, L_2, L_4, L_80 and Chamfer masks M_{3−4} and M_{5−7−11}." + +For a complete discussion of metric concepts in DGtal, please refer to +@ref moduleMetrics. + + +\section voronoiSect Digital Voronoi Diagram Computation + +The generic distance transformation is based on a prior Voronoi map +construction. Indeed, if we compute the Voronoi diagram of background +points, the distance transformation at an object point is exactly its +distance to the site associated with the Voronoi cell it belongs to. + +The core of the implementation is based on a separable approach: For +example, in dimension 2, partial digital Voronoi maps of dimension one are +first computed in each row independently. Then such partial Voronoi +maps are updated using independent processes along the columns, +leading to a valid Voronoi map of dimension 2. In an algorithmic point +of view, the 1D processes used for both columns and rows are +the same. In higher dimensions, the other dimensions are processed +similarly. + +@note We say digital Voronoi map instead of Voronoi diagram since the +output of the result is the intersection between the Voronoi diagram +of exterior points with @f$ \mathbb{Z}^d @f$. Furthermore, along +Voronoi faces (@e e.i. when two sites are equidistant), only one sites +is considered when intersection with @f$ \mathbb{Z}^d @f$. + + +In DGtal, the class VoronoiMap implements such digital Vornoi map +extraction. Such class is parametrized by the following types: +- a type representing the underlying digital space (model of CSpace); +- a type representing the object @f$ X @f$ as a point predicate (model + of concepts::CPointPredicate) ; +- a type representing the underlying metric (model of +CSeparableMetric, see below) +- and an optional image container to store the resulting Voronoi map +(by default, the type is +ImageContainerBySTLVector,typename +TSpace::Vector>). + +The VoronoiMap constructor is parametrized by +- an instance of Domain (the Domain type associated with the image +container); +- an instance of the PointPredicate ; +- and an instance of the separable metric. + +The VoronoiMap will be computed on the specified and will use the +point predicate to decide if a point of such domain is in the object +or note. + +@warning The point predicate must be valid for each point in the +specified domain. Basically, the domain can a sub-domain of the point +predicate definition domain. + + +Once the VoronoiMap object is created, the voronoi map is computed and +the class itself is a model of CConstImage. In other words, you can +access to the VoronoiMap value at a point @a p and iterate of values +using image ranges (see @ref moduleImages). For example + +@code + VoronoiMap<....> myVoronoiMap( .... ); //object construction + VoronoiMap<....>::Point p(12,34); + + VoronoiMap<....>::Value vectorToClosestSiteAtP = myVoronoiMap( p ); + for(VoronoiMap<....>::Domain::ConstIterator it = myVoronoiMap.domain().begin() , itend = myVoronoiMap.domain().end(); + it != itend ; ++it) + //do something on myVoronoiMap(it) +@endcode + +Since we are constructing a VoronoiMap, the value type of the map is a +vector (pointing to the closest site) of type Space::Vector (if Space +was the underlying digital space type used when specifying VoronoiMap +template parameters). + +Let us illustrate the construction in dimension 2 (see +voronoimap2D.cpp). Other examples can be found in distancetransform2D.cpp and distancetransform3D.cpp. + + First of all, we need couple of includes: + +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-header + +We will discuss later about the metric definition but let us consider +a classical Euclidean @f$ l_2 @f$ metric: + +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-Metric + +We now consider an object in a [0,0]x[16,16] domain with three +background points. To construct such point predicate, we first define +a set containing the three points, then we consider the point +predicate defined on this set (which returns true at a point if the +point is inside the set) and we consider the negation of such predicate +in order to return true for object points. Here you have the +construction: + +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-SmallImage + +and the resulting set: + +@image html voronoimap-inputset.png "Input set." +@image latex voronoimap-inputset.png "Input set." + + +The voronoi map is simply given by: + +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-Voro + +At each point of the object, we thus have a vector to the closest +background point. We can display this information as follows: + +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-trace + +To obtain: + +@image html voronoimap-voro.png "Voronoi map." +@image latex voronoimap-voro.png "Voronoi map." + +Changing the board output, we can see the Voronoi cells accordingly: + +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-traceCell + +To get: + +@image html voronoimap-cells.png "Voronoi map cells." +@image latex voronoimap-cells.png "Voronoi map cells." + +We could easily change the metric (here to the @f$ l_8 @f$) and get a +new Voronoi map: +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-l8Metric + +@image html voronoimap-vorol8.png "Voronoi map for the l_8 metric." +@image latex voronoimap-vorol8.png "Voronoi map for the l_8 metric." + +@section DTsec Distance Transformation + +As discussed earlier, the distance transformation is given by +computing distances once the Voronoi map is obtained. In DGtal, the +class DistanceTransformation simply adapts the VoronoiMap class in +order to override output image getters to return the distance for the +given metric to the closest site instead of the vector. + +As a consequence, the DistanceTransformation class simply inherits from the +VoronoiMap class and overrides methods required by +the concepts::CConstImage concept. Note that the DistanceTransformation::Value +type is double. If you want to get the underlying vector instead of +the distance to perform exact computations, you can use the +DistanceTransformation::getVoronoiVector method. + +In the following example, we consider the previous small image and use +a colormap to display distance values for the @f$ l_2 @f$ mertic: + +@snippet geometry/volumes/distance/voronoimap2D.cpp Voro2D-DT + +@image html voronoimap-dt.png "Distance transformation for the l_2 metric." +@image latex voronoimap-dt.png "Distance transformation for the l_2 metric." + + + +@section RDTSec Digital Power Map and Reverse Distance Transformation + +Similarly to Voronoi diagram and digital Voronoi maps, digital Power +maps are defined as the intersection between the integer grid and a +power diagram. Given a set of weighed points, power diagram can be +seen as Voronoi diagram where the metric is modified with additive +weights. For example, considering the @f$ l_2@f$ metric, the power +distance between a point @f$p@f$ and a weighted point @f$(q,w)@f$ is +defined by +@f[ pow(p,q) = \| p - q\|_2^2 - w @f] + +Hence, similarly to Voronoi diagram, the power diagram is a +decomposition of the space ino cells from weighed sites where each +cell (maybe empty) is associated with a site and each point in the +cell has got minimal power distance to the cell site (compared to its +power distance to all other sites) @cite Aurenhammer1987 . + + +Separable algorithms similar to VoronoiMap/DistanceTransformation can +be designed to compute respectively PowerMap and +ReverseDistanceTransformation. The only difference is that the input of +PowerMap is a weighted set of points instead of a point predicate. + +@note for @f$l_p@f$ metrics, the power distance is defined by @f$ +pow(p,q) = \| p - q\|_p^p - w @f$. Hence, both the distance and the weight value type +capacity must be able to represent d-sums of power p numbers (if d is +the dimension of the space). + + +Hence such class is parametrized by the following types: +- a type representing the mapping between points and weights (WeightImage, model of concepts::CImage); +- a type representing the underlying power metric (model of concepts::CPowerSeparableMetric, see below) +- and an optional image container to store the resulting Power map (by default, the type is ImageContainerBySTLVector). + +The PowerMap constructor is parametrized by +- an instance of Domain (the Domain type associated with the image container); +- an instance of PowerMap::WeightImage; +- and an instance of the power separable metric. + +Similarly to DistanceTransformation, ReverseDistanceTransformation +remaps the PowerMap vectors to map the power metric to the closest +weighted site. + +As a consequence, for the Euclidean @f$ l_2 @f$ metric, if we consider +a set of balls @f$ B_i(p_i,r_i) @f$ and if we create an WeightImage +whose domain contains points @f$\{ p_i \}@f$ and with +values @f$ r_i^2@f$, negative (strictly) values of the +ReverseDistanceTransformation will correspond to digital points +belonging to the union @f$ \bigcup \{B_i\}@f$ (see @cite dcoeurjo_pami_RDMA). + +ReverseDistanceTransformation can thus be used to reconstructed a +binary shape from a given Medial Axis or any set of balls. Another +consequence is that given a binary shape, the pipeline +@f[ Shape \rightarrow DT \rightarrow ReverseDT \rightarrow \text{ strictly negative values }@f] +for the same metric/power metric, returns the input binary shape. + + + +@note Power separable metrics are formalized in concepts::CPowerMetric and +concepts::CPowerSeparableMetric concepts whose main model is +ExactPredicateLpPowerSeparableMetric, see @ref moduleMetrics + +@section toricVol Volumetric Analysis on Toric domains + +In some specific applications, toric domains and volumetric analysis +of shapes in toric domains are crucial. Thanks to the separability +property of VoronoiMap, PowerMap (and their associated wrappers +DistanceTransformation and ReverseDistanceTransformation), one can +easily consider volumetric transformation in arbitrary dimension on +toric domains @cite Coeurjo2008. Note that changing the periodicity +property of the domain has no impact on the computational cost. + +More precisely, VoronoiMap and DistanceTransformation classes have +constructors allowing to specify periodicity information. In dimension +@e d, such periodicity information is an array with @e d boolean flags +where the i-th value is \c true if the i-th dimension of the space is +periodic, \c false otherwise. Hence, computation can be performed +either on a full toric domain, or on a domain with toricity property +along only one axis. Similar extensions to toric domains have been +implemented for the PowerMap and ReverseDistanceTransformation +classes. + +As illustrated in the example toricdomainvolumetric.cpp, given an +input: + +@image html toric-inputShape-red.png "Input domain and sites." + +We consider the following distance transformation objects: + +@snippet geometry/volumes/distance/toricdomainvolumetric.cpp DTComputeToric + +The following results illustrate both distance transformation and +Voronoi maps. For the VoronoiMap results, points may be attached to +sites exterior to the initial domain. In fact such sites correspond to +toric replicas of existing sites within the domain. + +
+ + + + + + + + + + +
+Classical domain + +Toric domain +
+@image html toric-example-DT-L2-red.png "DT values." + +@image html toric-example-DT-L2-toric-red.png "Toric DT values." +
+@image html toric-example-Voro-L2-red.png "Voronoi map as vectors." + +@image html toric-example-Voro-L2-toric-red.png "Voronoi map as vectors." +
+
+ +Using VoronoiMap::projectPoint(Point) const, site's coordinates can be projected +into the initial domain, even for VoronoiMap calculated on toric domains: + +@image html toric-example-Voro-L2-toric-projected.png "Voronoi map on toric domain with projected sites." + +With partial periodicity specification (along the first or second dimension only): + +@snippet geometry/volumes/distance/toricdomainvolumetric.cpp DTComputePartialToric + +we obtain the following VoronoiMap: + +
+ + + + + + + + + + +
+Periodic domain along the 1th dimension. + +Periodic domain along the 2th dimension. +
+@image html toric-example-DT-L2-toricX-red.png "DT values." + +@image html toric-example-DT-L2-toricY-red.png "DT values." +
+@image html toric-example-Voro-L2-toricX-red.png "Voronoi map as vectors." + +@image html toric-example-Voro-L2-toricY-red.png "Voronoi map as vectors." +
+
+ +*/ + + +} diff --git a/src/DGtal/geometry/doc/packageGeometry.dox b/src/DGtal/geometry/doc/packageGeometry.dox index 526bd2d339..803a229bbd 100644 --- a/src/DGtal/geometry/doc/packageGeometry.dox +++ b/src/DGtal/geometry/doc/packageGeometry.dox @@ -77,6 +77,7 @@ of arbitrary dimension, by the means of separable and incremental distance trans - \ref moduleVolumetric (David Coeurjolly, Roland Denis, Robin Lamy (VoronoiComplete)) - \ref moduleFMM (Tristan Roussillon) - \ref moduleDigitalConvexity (Jacques-Olivier Lachaud) + - \ref moduleMedialAxis (Robin Lamy, David Coeurjolly, Isabelle Sivignon) - Geometry Processing - \ref moduleRegularization (David Coeurjolly, Jacques-Olivier Lachaud, Pierre Gueth) diff --git a/src/DGtal/geometry/tools/ImplicitBallTools.h b/src/DGtal/geometry/tools/ImplicitBallTools.h new file mode 100644 index 0000000000..bc79a07882 --- /dev/null +++ b/src/DGtal/geometry/tools/ImplicitBallTools.h @@ -0,0 +1,207 @@ +#include +#include +#include "DGtal/base/Common.h" +#include "DGtal/helpers/StdDefs.h" +#include "DGtal/kernel/NumberTraits.h" +#include "DGtal/kernel/domains/HyperRectDomain.h" +#include "DGtal/images/ImageContainerBySTLMap.h" +#include "DGtal/kernel/sets/CDigitalSet.h" +#include "DGtal/shapes/implicit/ImplicitBall.h" + +#ifdef WITH_EIGEN +#include +#endif + +namespace DGtal +{ + namespace functions + { + /** + * Returns the smallest nD ball going through two points. + * + * @param aPoint one point + * @param otherPoint other point + * + * @return the smallest ImplicitBall that goes + * through the 2 points + */ + template + ImplicitBall > ballFrom2Points(const TPoint & aPoint, + const TPoint & otherPoint) + { + typename SpaceND::RealPoint center(aPoint); + center += otherPoint; + double radius = (aPoint - otherPoint).norm()/2; + return ImplicitBall >(0.5*center, radius); + } + + + /** + * Return the smallest enclosing ball from 4 points in 3D + */ + template + ImplicitBall > + getBallFrom4Point3D(const TPoint & A, const TPoint & B, const TPoint & C, const TPoint & D) + { + BOOST_ASSERT_MSG(TPoint::dimension == 3, "Point dimension must be 3"); + Eigen::Matrix3f M; + Eigen::Vector3f b; + int b1, b2, b3; + + b1 = -(pow(B[0], 2) - pow(A[0], 2) + pow(B[1], 2) - pow(A[1], 2) + pow(B[2], 2) - pow(A[2], 2)); + b2 = -(pow(C[0], 2) - pow(A[0], 2) + pow(C[1], 2) - pow(A[1], 2) + pow(C[2], 2) - pow(A[2], 2)); + b3 = -(pow(D[0], 2) - pow(A[0], 2) + pow(D[1], 2) - pow(A[1], 2) + pow(D[2], 2) - pow(A[2], 2)); + b << b1, b2, b3; + + M << 2*(A[0] - B[0]) , 2*(A[1] - B[1]) , 2*(A[2] - B[2]) , + 2*(A[0] - C[0]) , 2*(A[1] - C[1]) , 2*(A[2] - C[2]) , + 2*(A[0] - D[0]) , 2*(A[1] - D[1]) , 2*(A[2] - D[2]); + + Eigen::Vector3f centerVector = M.colPivHouseholderQr().solve(b); + typename SpaceND::RealPoint center(centerVector(0), centerVector(1), centerVector(2)); + float x, y, z; + x = center[0]; + y = center[1]; + z = center[2]; + float radius = sqrt(pow(x-A[0], 2) + pow(y-A[1], 2) + pow(z-A[2], 2)); + ImplicitBall > ball(center, radius); + return ball; + } + + + /** + * @param p1 one point of dimension 2 + * @param p2 second point of dimension 2 + * @param p3 third point of dimension 2 + * + * @return the smallest ImplicitBall of dimension 2 that goes + * through the 3 points + */ + template + ImplicitBall > ballFrom3Points(const TPoint & p1, + const TPoint & p2, + const TPoint & p3) + { + BOOST_ASSERT(TPoint::dimension == 2 || TPoint::dimension == 3); + + if (TPoint::dimension == 3) + { + return getBallFrom4Point3D(p1, p2, p3, p3); + } + /** + * Formule d'après : Centre et rayon d’un cerclepassant par trois points donnés(Phm 2006/02/05) + * https://cral-perso.univ-lyon1.fr/labo/fc/Ateliers_archives/ateliers_2005-06/cercle_3pts.pdf + */ + double x1, x2, x3, y1, y2, y3; + // Cas d'un côté parallèle à l'axe des abscisses pour éviter la division par 0, ici [P1, P2] + if (!(p1[1] - p2[1])) { + x1 = p1[0]; + y1 = p1[1]; + x3 = p2[0]; + y3 = p2[1]; + x2 = p3[0]; + y2 = p3[1]; + } else if (!(p2[1] - p3[1])) { // Ici [P2, P3] + x1 = p2[0]; + y1 = p2[1]; + x3 = p3[0]; + y3 = p3[1]; + x2 = p1[0]; + y2 = p1[1]; + } else { + x1 = p1[0]; + y1 = p1[1]; + x2 = p2[0]; + y2 = p2[1]; + x3 = p3[0]; + y3 = p3[1]; + } + //Fin de gestion des cas particuliers + + double X1 = pow(x1, 2); + double X2 = pow(x2, 2); + double X3 = pow(x3, 2); + double Y1 = pow(y1, 2); + double Y2 = pow(y2, 2); + double Y3 = pow(y3, 2); + + double centerX = -((X3-X2+Y3-Y2)/(2*(y3-y2)) - (X2-X1+Y2-Y1)/(2*(y2-y1))) / ((x2-x1)/(y2-y1) - (x3-x2)/(y3-y2)); + double centerY = -((x2-x1)/(y2-y1)) * centerX + (X2-X1+Y2-Y1)/(2*(y2-y1)); + double radius = sqrt(pow(centerX - x1, 2) + pow(centerY - y1, 2)); + typename SpaceND::RealPoint center(centerX, centerY); + return ImplicitBall >(center, radius); + } + + + + /** + * @param R a set (DigitalSet, SmallSet etc) of size 3 or less + * + * @return the smallest ball with R on its border + */ + template + ImplicitBall trivialCircle(TSet & R) { + BOOST_ASSERT_MSG(R.size() < 5, "Set size must be 4 or less to find the trivial ball"); + BOOST_CONCEPT_ASSERT(( concepts::CDigitalSet )); + typedef typename TSet::Space::RealPoint RealPoint; + typedef ImplicitBall ImplicitBall; + double x, y, r; + + switch (R.size()) { + case 1 : + return ImplicitBall(*(R.begin()), 0); + case 2 : + return ballFrom2Points(*(R.begin()), *(R.begin()+1)); + case 3 : + return ballFrom3Points(*(R.begin()), *(R.begin()+1), *(R.begin()+2)); + case 4 : + return getBallFrom4Point3D(*(R.begin()), *(R.begin()+1), *(R.begin()+2), *(R.begin()+3)); + default : + return ImplicitBall(RealPoint(0.5, 0.5), 0); + } + } + + /** + * @param P the points we want to put in the + * smallest enclosing ball. The type is necessary to use + * random access, to generalize later... + * @param R the points that are on the border of + * the ball, initially empty. + * + * @return the smallest enclosing ImplicitBall of Point dimension + */ + template + ImplicitBall smallestEuclideanEnclosingBall(std::vector & P, TSet & R) + { + //This methods has trivialCircle code only in 2d. + //In higher dimensions, use ........ + typedef typename TSet::Point Point; + typedef typename TSet::Space::RealPoint RealPoint; + typedef ImplicitBall ImplicitBall; + + BOOST_STATIC_ASSERT_MSG( Point::dimension == 2, "Point dimension must be 2"); + BOOST_CONCEPT_ASSERT(( concepts::CDigitalSet )); + + double x, y, r; + if (P.empty() || R.size() == Point::dimension) { + ImplicitBall ball(trivialCircle(R)); + return ball; + } + int random = std::rand()%P.size(); + Point p = P[random]; + P.erase(P.begin()+random); + ImplicitBall D = smallestEuclideanEnclosingBall(P, R); + RealPoint center(D.center()); + r = D.radius(); + if ( (p-center).squaredNorm() < r*r ) { + P.insert(P.begin(), p); + return D; + } + R.insert(p); + ImplicitBall ball(smallestEuclideanEnclosingBall(P, R)); + R.erase(p); + P.insert(P.begin(), p); + return ball; + } + } +} diff --git a/src/DGtal/geometry/volumes/distance/LambdaMedialAxis.h b/src/DGtal/geometry/volumes/distance/LambdaMedialAxis.h new file mode 100644 index 0000000000..f18143ec8d --- /dev/null +++ b/src/DGtal/geometry/volumes/distance/LambdaMedialAxis.h @@ -0,0 +1,282 @@ +/** + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + **/ + +#pragma once + +/** + * @file LambdaMedialAxis.h + * @brief Linear in time distance transformation + * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * + * @date 2012/12/08 + * + * Header file for module + * + * This file is part of the DGtal library. + * + * @see testReducedMedialAxis.cpp, testReducedMedialAxisND.cpp, testReverseDT.cpp + */ + +#if defined(LambdaMedialAxis_RECURSES) +#error Recursive header files inclusion detected in ReducedMedialAxis.h +#else // defined(ReducedMedialAxis_RECURSES) +/** Prevents recursive inclusion of headers. */ +#define LambdaMedialAxis_RECURSES + +#if !defined LambdaMedialAxis_h +/** Prevents repeated inclusion of headers. */ +#define LambdaMedialAxis_h + +////////////////////////////////////////////////////////////////////////////// +// Inclusions +#include +#include +#include "DGtal/base/Common.h" +#include "DGtal/helpers/StdDefs.h" +#include "DGtal/kernel/NumberTraits.h" +#include "DGtal/geometry/volumes/distance/PowerMap.h" +#include "DGtal/kernel/domains/HyperRectDomain.h" +#include "DGtal/images/ImageContainerBySTLMap.h" +#include "DGtal/images/Image.h" +#include "DGtal/geometry/volumes/distance/DistanceTransformation.h" +#include "DGtal/geometry/volumes/distance/VoronoiMap.h" +#include "DGtal/kernel/sets/CDigitalSet.h" +#include "DGtal/geometry/tools/ImplicitBallTools.h" +////////////////////////////////////////////////////////////////////////////// + +namespace DGtal +{ + + ///////////////////////////////////////////////////////////////////////////// + // template class LambdaMedialAxis + /** + * Description of template class 'LambdaMedialAxis'

+ * \brief Aim: Implementation of the separable lambda medial axis + * extraction. + * + * This utility struct extract lambda medial axis balls from a + * Set. Basically, each (weighted) site of the PowerMap defines + * a digital maximal ball if its digital power cell restricted to + * the input shape is not empty @cite dcoeurjo_pami_RDMA . + * + * Optimal Separable Algorithms to Compute the Reverse + * Euclidean Distance Transformation and Discrete Medial Axis in + * Arbitrary Dimension, D. Coeurjolly and A. Montanvert, IEEE + * Transactions on Pattern Analysis and Machine Intelligence, + * 29(3):437-448, 2007. + * + * The output is an image associating ball radii (weight of the + * power map site) to maximal ball centers. Most methods output a + * lightweight proxy to an image container (of type ImageContainer, + * see below). + * + * @note Following ReverseDistanceTransformation, the input shape is + * defined as points with negative power distance. + * + * @tparam TSet any set type (to be generalized as PointPredicate) + * @tparam TDistanceTransformation any DistanceTransformation type + * that has SmallSet as its PointPredicate + * @tparam TPowerMap any specialized PowerMap type + * @tparam TSmallObject a SmallObject type in dimension 2 or 3 + * @tparam TTopology a topology type that match the SmallObject + * @tparam TImageContainer any model of CImage to store the medial axis + * points (default: ImageContainerBySTLVector). + * + * @see testLambdaMedialAxis.cpp //TODO + */ + + + template , double> > + struct LambdaMedialAxis + { + typedef TSpace Space; + + typedef DGtal::int32_t Integer; + + typedef typename Space::Point Point; + + typedef typename Space::RealPoint RealPoint; + + typedef HyperRectDomain Domain; + + typedef ExactPredicateLpSeparableMetric L2Metric; + + typedef ExactPredicateLpPowerSeparableMetric L2PowerMetric; + + typedef typename DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet; + + typedef typename DigitalSetSelector < Domain, SMALL_DS + HIGH_ITER_DS >::Type SmallSet; + + typedef Image Type; + + typedef DistanceTransformation DT; + + typedef PowerMap PMap; + + typedef VoronoiMap VMap; + + typedef ReverseDistanceTransformation RDT; + + typedef std::pair Ball; + + /** + * Detail on algorithm + * + * @param set a set of points + * @param lambda the minimal radius to accept a ball in the lambda medial axis + * + * @return a lightweight proxy to the ImageContainer specified in + * template arguments. (soon to be a vector > instead) + */ + static + std::vector computeLambdaAxisFromSet(DigitalSet & set, double lambda) { + // Initialize metrics and topology + L2Metric l2Metric; + L2PowerMetric l2PowerMetric; + + std::vector medialAxis; + + TImageContainer *computedMA = new TImageContainer( set.domain() ); + TImageContainer *squaredDT = new TImageContainer( set.domain() ); + + // Building the VoronoiMap, DistanceTransformation, + // SquaredDistanceTransformation and PowerMap for the next steps + VMap vmap(set.domain(), set, l2Metric); + DT dt(set.domain(), set, l2Metric); + for (typename DT::Domain::ConstIterator it = dt.domain().begin(); it != dt.domain().end(); it++) { + squaredDT->setValue((*it), pow(dt(*it), 2)); + } + PMap aPowerMap(squaredDT->domain(), squaredDT, l2PowerMetric); + + for (typename PMap::Domain::ConstIterator it = aPowerMap.domain().begin(), + itend = aPowerMap.domain().end(); it != itend; ++it) + { + const auto v = aPowerMap( *it ); + const auto pv = aPowerMap.projectPoint( v ); + + if ( aPowerMap.metricPtr()->powerDistance( *it, v, aPowerMap.weightImagePtr()->operator()( pv ) ) + < NumberTraits::ZERO ) { + + //Lambda Medial Axis filter + Point voronoiPoint(vmap(*it)); + double voronoiDistance = (voronoiPoint - *it).norm(); //Squared Distance to its Voronoi map point + std::vector otherVoronoiPoints; + DigitalSet box(set.domain()); + Point center = *it; + subBox(center, 2, set, box); // point neighborhood + for (auto neighbor : box) { + Point neighborVoronoiPoint(vmap(neighbor)); // Voronoi point of the neighbor + double distanceToTest = (neighborVoronoiPoint - *it).norm(); // Distance between the point, and the neighbor voronoi point + if (distanceToTest == voronoiDistance) { + otherVoronoiPoints.push_back(neighborVoronoiPoint); + } + } + SmallSet R(set.domain()); + ImplicitBall enclosing(smallestEuclideanEnclosingBall(otherVoronoiPoints, R)); + if (enclosing.radius() > lambda) { + Ball ball(v, aPowerMap.weightImagePtr()->operator()( pv ) ); + if (std::find(medialAxis.begin(), medialAxis.end(), ball) == medialAxis.end()) { + medialAxis.push_back(ball); + } + } + // End of filter + + } + } + + return medialAxis; + } + + static std::vector computeLambdaAxisFromVector(std::vector & ballSet, Domain & domain, double lambda) { + + L2Metric l2Metric; + L2PowerMetric l2PowerMetric; + + //Reconstruction of the shape + TImageContainer *weightedMap = new TImageContainer( domain ); + for (Ball ball : ballSet) { + weightedMap->setValue(ball.first, ball.second); + } + + RDT rdt(domain, *weightedMap, l2PowerMetric); + DigitalSet set(domain); + for (typename RDT::Domain::ConstIterator it = rdt.domain().begin(); it != rdt.domain().end(); it++) { + if (rdt(*it) < 0) { + set.insert(*it); + } + } + + std::vector medialAxis; + + medialAxis = computeLambdaAxisFromSet(set, lambda); + + return medialAxis; + } + + + /** + * @param center + * @param size the box goes from -size to +size in all dimensions + * @param set the set we want a subBox of + * + * @return a subset of type DigitalSet that contains the box + * + */ + static void subBox(Point & center, + Integer size, + DigitalSet & set, + DigitalSet & box, + Integer currentDimension = 0, + std::vector indexes = std::vector()) { + Integer dimension = Point::dimension; + if (indexes.size() == 0) { + for (Integer i = 0; i < dimension; i++) { + indexes.push_back(-size); + } + } + Point point; + for (Integer i = 0; i < dimension; i++) { + point[i] = indexes[i]; + } + + while(indexes[currentDimension] <= size) { + if (currentDimension == dimension - 1) { + if (set(center + point)) { + box.insert(center + point); + } + } else { + subBox(center, size, set, box, currentDimension+1, indexes); + } + indexes[currentDimension] += 1; + point[currentDimension] += 1; + } + + return; + } + + }; // end of class LambdaMedialAxis + +} // namespace DGtal + +// // +/////////////////////////////////////////////////////////////////////////////// + +#endif // !defined LambdaMedialAxis_h + +#undef LambdaMedialAxis_RECURSES +#endif // else defined(LambdaMedialAxisdesign pa_RECURSES) diff --git a/src/DGtal/geometry/volumes/distance/ReducedMedialAxis.h b/src/DGtal/geometry/volumes/distance/ReducedMedialAxis.h index 6e152f2ebb..db5b2cb258 100644 --- a/src/DGtal/geometry/volumes/distance/ReducedMedialAxis.h +++ b/src/DGtal/geometry/volumes/distance/ReducedMedialAxis.h @@ -97,7 +97,8 @@ namespace DGtal struct ReducedMedialAxis { //MA Container - typedef Image Type; + typedef std::pair Ball; + typedef typename TPowerMap::Point Point; /** * Extract reduced medial axis from a power map. @@ -109,23 +110,33 @@ namespace DGtal * template arguments. */ static - Type getReducedMedialAxisFromPowerMap(const TPowerMap &aPowerMap) + std::vector + getReducedMedialAxisFromPowerMap(const TPowerMap &aPowerMap) { - TImageContainer *computedMA = new TImageContainer( aPowerMap.domain() ); + std::vector medialAxis; + medialAxis.push_back(Ball(*(aPowerMap.domain().begin()), 0)); //To have a non empty vector for the loop for (typename TPowerMap::Domain::ConstIterator it = aPowerMap.domain().begin(), itend = aPowerMap.domain().end(); it != itend; ++it) { + const auto v = aPowerMap( *it ); const auto pv = aPowerMap.projectPoint( v ); if ( aPowerMap.metricPtr()->powerDistance( *it, v, aPowerMap.weightImagePtr()->operator()( pv ) ) - < NumberTraits::ZERO ) - computedMA->setValue( v, aPowerMap.weightImagePtr()->operator()( pv ) ); + < NumberTraits::ZERO ) { + + Ball ball(v, aPowerMap.weightImagePtr()->operator()( pv ) ); + if (std::find(medialAxis.begin(), medialAxis.end(), ball) == medialAxis.end()) { + medialAxis.push_back(ball); + } + } } - return Type( computedMA ); + medialAxis.erase(medialAxis.begin()); + return medialAxis; } + }; // end of class ReducedMedialAxis diff --git a/src/DGtal/geometry/volumes/distance/ScaleMedialAxis.h b/src/DGtal/geometry/volumes/distance/ScaleMedialAxis.h new file mode 100644 index 0000000000..fdc4f2b187 --- /dev/null +++ b/src/DGtal/geometry/volumes/distance/ScaleMedialAxis.h @@ -0,0 +1,213 @@ +/** + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + **/ + +#pragma once + +/** + * @file ScaledMedialAxis.h + * @brief Linear in time distance transformation + * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * + * @date 2012/12/08 + * + * Header file for module + * + * This file is part of the DGtal library. + * + * @see testReducedMedialAxis.cpp, testReducedMedialAxisND.cpp, testReverseDT.cpp + */ + +#if defined(ScaledMedialAxis_RECURSES) +#error Recursive header files inclusion detected in ScaledMedialAxis.h +#else // defined(ScaledMedialAxis_RECURSES) +/** Prevents recursive inclusion of headers. */ +#define ScaledMedialAxis_RECURSES + +#if !defined ScaledMedialAxis_h +/** Prevents repeated inclusion of headers. */ +#define ScaledMedialAxis_h + +////////////////////////////////////////////////////////////////////////////// +// Inclusions +#include +#include +#include +#include "DGtal/base/Common.h" +#include "DGtal/kernel/NumberTraits.h" +#include "DGtal/geometry/volumes/distance/CPowerSeparableMetric.h" +#include "DGtal/geometry/volumes/distance/PowerMap.h" +#include "DGtal/geometry/volumes/distance/DistanceTransformation.h" +#include "DGtal/geometry/volumes/distance/ReverseDistanceTransformation.h" +#include "DGtal/geometry/volumes/distance/ReducedMedialAxis.h" +#include "DGtal/images/DefaultConstImageRange.h" +#include "DGtal/kernel/domains/HyperRectDomain.h" +#include "DGtal/images/ImageContainerBySTLMap.h" +#include "DGtal/images/CImage.h" +#include "DGtal/images/Image.h" +#include "DGtal/images/ImageSelector.h" +////////////////////////////////////////////////////////////////////////////// + +namespace DGtal +{ + + ///////////////////////////////////////////////////////////////////////////// + // template class ReducedMedialAxis + /** + * Description of template class 'ReducedMedialAxis'

+ * \brief Aim: Implementation of the separable medial axis + * extraction. + * + * This utility struct extract scale medial axis balls from a + * DistanceTransformation. Basically, each (weighted) site of the + * DistanceTransformation defines a digital ball, we make it + * alpha-bigger, extract the medial axis, and then make it + * alpha-smaller. + * + * Optimal Separable Algorithms to Compute the Reverse + * Euclidean Distance Transformation and Discrete Medial Axis in + * Arbitrary Dimension, D. Coeurjolly and A. Montanvert, IEEE + * Transactions on Pattern Analysis and Machine Intelligence, + * 29(3):437-448, 2007. + * + * The output is an image associating ball radii (weight of the + * power map site) to maximal ball centers. Most methods output a + * lightweight proxy to an image container (of type ImageContainer, + * see below). + * + * @note Following ReverseDistanceTransformation, the input shape is + * defined as points with negative power distance. + * + * @tparam TSpace any Space (usually SpaceND) + * @tparam TImageContainer any model of CImage to store the medial axis + * points (default: ImageContainerBySTLVector). + * + * @see testScaleMedialAxis.cpp + */ + template , double> > + struct ScaleMedialAxis + { + //MA Container + typedef TSpace Space; + + typedef typename Space::Point Point; + + typedef typename Space::RealPoint RealPoint; + + typedef HyperRectDomain Domain; + + typedef ExactPredicateLpSeparableMetric L2Metric; + + typedef ExactPredicateLpPowerSeparableMetric L2PowerMetric; + + typedef typename DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet; + + typedef DistanceTransformation DT; + + typedef ReverseDistanceTransformation RDT; + + typedef std::pair Ball; + + /** + * Extract the scaled medial axis from a distance transformation. + * This methods is in @f$ O(|distanceTransformation|)@f$. + * + * @param dt the input DistanceTransformation + * @param alpha the double ratio we use to expand and compress the RDMA + * + * @return a std::vector > that contains centers and + * radius of the Scale Medial Axis balls. + */ + + static std::vector computeScaleAxisFromDT(DT & dt, double alpha) + { + std::vector medialAxis; + L2PowerMetric l2PowerMetric; + + // SquaredDT construction, there might be a better way to skip this step + // All balls get alpha times bigger + TImageContainer *bigdt = new TImageContainer( dt.domain() ); + for (auto it : dt.domain()) { + bigdt->setValue(it, pow(dt(it), 2) * alpha); + } + + // Medial Axis Extraction + PowerMap power(bigdt->domain(), bigdt, l2PowerMetric); + std::vector BigRdma = ReducedMedialAxis, TImageContainer >::getReducedMedialAxisFromPowerMap(power); + + // alpha time reduction of all balls size + for (Ball ball : BigRdma) { + const auto v = ball.second; // Value of distance in BigRdma + const auto compressedValue = v/alpha; + Ball compressedBall(ball.first, compressedValue); + + if (std::find(medialAxis.begin(), medialAxis.end(), compressedBall) == medialAxis.end()) { + medialAxis.push_back(compressedBall); + } + } + return medialAxis; + } + + /** + * Extract the scaled medial axis from a distance transformation. + * This methods is in @f$ O(|distanceTransformation|)@f$. + * + * @param ballSet the input std::vector > set of balls + * @param alpha the double ratio we use to expand and compress the RDMA + * + * @return a lightweight proxy to the ImageContainer specified in + * template arguments. + */ + + static std::vector computeScaleAxisFromVector(std::vector & ballSet, Domain & domain, double alpha) { + + L2Metric l2Metric; + L2PowerMetric l2PowerMetric; + + //Reconstruction of the shape + TImageContainer *weightedMap = new TImageContainer( domain ); + for (Ball ball : ballSet) { + weightedMap->setValue(ball.first, ball.second); + } + + RDT rdt(domain, *weightedMap, l2PowerMetric); + DigitalSet set(domain); + for (typename RDT::Domain::ConstIterator it = rdt.domain().begin(); it != rdt.domain().end(); it++) { + if (rdt(*it) < 0) { + set.insert(*it); + } + } + DT dt(domain, set, l2Metric); + + std::vector medialAxis; + + medialAxis = computeScaleAxisFromDT(dt, alpha); + + return medialAxis; + } + + }; // end of class ScaledMedialAxis + +} // namespace DGtal + +// // +/////////////////////////////////////////////////////////////////////////////// + +#endif // !defined ScaledMedialAxis_h + +#undef ScalededMedialAxis_RECURSES +#endif // else defined(ScaledMedialAxisdesign pa_RECURSES) diff --git a/tests/geometry/tools/testImplicitBallTools.cpp b/tests/geometry/tools/testImplicitBallTools.cpp new file mode 100644 index 0000000000..e9893d0e9d --- /dev/null +++ b/tests/geometry/tools/testImplicitBallTools.cpp @@ -0,0 +1,98 @@ +/** + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + **/ + +/** + * @file + * @ingroup Tests + * @author + * + * + * @date 30/06/2021 + * + * This file is part of the DGtal library + */ + +/** + * Description of testPointVector-catch'

+ * Aim: simple test of \ref ImplicitBallTools functions with Catch unit test framework. + */ + +#include +#include +#include "DGtal/base/Common.h" +#include "DGtal/helpers/StdDefs.h" +#include "DGtal/kernel/NumberTraits.h" +#include "DGtal/kernel/domains/HyperRectDomain.h" +#include "DGtal/images/ImageContainerBySTLMap.h" +#include "DGtal/kernel/sets/CDigitalSet.h" +#include "DGtal/shapes/implicit/ImplicitBall.h" +#include "DGtal/geometry/tools/ImplicitBallTools.h" + +#include "DGtalCatch.h" + + +using namespace DGtal; + +TEST_CASE( "Testing ImplicitBallTools functions" ) { + + typedef ImplicitBall > ImplicitBall2D; + typedef ImplicitBall > ImplicitBall3D; + + typedef Z2i::Point Point2D; + typedef Z2i::RealPoint RealPoint2D; + typedef Z3i::Point Point3D; + typedef Z3i::RealPoint RealPoint3D; + + SECTION( "Ball From 2 Points Dimension 2 and 3" ) + { + ImplicitBall2D ball = functions::ballFrom2Points(Point2D(0, 0), Point2D(5, 5)); + CHECK(ball.center() == RealPoint2D(2.5, 2.5)); + + ImplicitBall2D ball2 = functions::ballFrom2Points(Point2D(4, 4), Point2D(8, 8)); + CHECK(ball2.center() == RealPoint2D(6, 6)); + + ImplicitBall2D ball3 = functions::ballFrom2Points(Point2D(2, 2), Point2D(2, 2)); + CHECK(ball3.center() == RealPoint2D(2, 2)); + + ImplicitBall3D ball3D = functions::ballFrom2Points(Point3D(0, 0, 0), Point3D(5, 5, 5)); + CHECK(ball3D.center() == RealPoint3D(2.5, 2.5, 2.5)); + + ImplicitBall3D ball4 = functions::ballFrom2Points(Point3D(4, 5, 1), Point3D(2, 3, 7)); + CHECK(ball4.center() == RealPoint3D(3, 4, 4)); + + ImplicitBall3D ball3D2 = functions::ballFrom2Points(Point3D(2, 2, 2), Point3D(2, 2, 2)); + CHECK(ball3D2.center() == RealPoint3D(2, 2, 2)); + } + + SECTION( "Ball From 3 Points Dimension 2 and 3" ) + { + ImplicitBall2D ball = functions::ballFrom3Points(Point2D(0, 0), Point2D(1, 0), Point2D(0, 1)); + CHECK(ball.center() == RealPoint2D(0.5, 0.5)); + + ImplicitBall3D ball3D3 = functions::ballFrom3Points(Point3D(0, 0, 0), Point3D(1, 0, 0), Point3D(0, 1, 0)); + CHECK(ball3D3.center() == RealPoint3D(0.5, 0.5, 0.5)); + } + + /*SECTION( "Ball From 4 Points") + { + ImplicitBall3D ball3D4 = functions::ballFrom4Points(Point3D(0, 0, 0), Point3D(1, 1, 1), Point3D(1, 0, 0), Point3D(0, 1, 0)); + CHECK(ball3D4.center() == RealPoint3D(0.5, 0.5, 0.5)); + }*/ + + SECTION( "Trivial Circle" ) { + + } +} diff --git a/tests/geometry/volumes/distance/testReducedMedialAxis.cpp b/tests/geometry/volumes/distance/testReducedMedialAxis.cpp index b755288acd..741382b82c 100644 --- a/tests/geometry/volumes/distance/testReducedMedialAxis.cpp +++ b/tests/geometry/volumes/distance/testReducedMedialAxis.cpp @@ -113,24 +113,11 @@ bool testReducedMedialAxis( std::array const& aPeriodicity = {{ false, trace.info()< >::Type rdma = ReducedMedialAxis< PowerMap >::getReducedMedialAxisFromPowerMap(power); - - //Reconstruction - for(unsigned int i=0; i<11; i++) - { - for(unsigned int j=0; j<11; j++) - { - Z2i::Point p(i,j); - if (rdma.domain().isInside(p) ) - trace.info()<< rdma(p); - else - trace.info()<< " - "; - - } - std::cerr< >::getReducedMedialAxisFromPowerMap(power); + for(auto &ball: rdma) + trace.info() <<"Got ball: "<< ball.first<<" r="<< ball.second< const& aPeriodicity = {{ false, << "true == true" << std::endl; trace.endBlock(); - bool isEqual = true; - for ( auto const & pt : domain ) - { - const Image::Value a = image.domain().isInside(pt) ? image(pt) : 0; - const Image::Value b = rdma.domain().isInside(pt) ? rdma(pt) : 0; - if ( a != b ) - { - isEqual = false; - break; - } - } - - trace.info() << "Equality ? " << isEqual << std::endl; - return nbok == nb; }