DGtal 2.1.0
Loading...
Searching...
No Matches
checkFullConvexityTheorems.cpp
Go to the documentation of this file.
1
44#include <iostream>
45#include <queue>
46#include "DGtal/base/Common.h"
47#include "DGtal/helpers/StdDefs.h"
48#include "DGtal/geometry/volumes/DigitalConvexity.h"
49
51
52using namespace std;
53using namespace DGtal;
54
55template <typename Point>
56void makeRandom( Point& p, int width )
57{
58 for ( Dimension i = 0; i < Point::dimension; i++ )
59 p[ i ] = rand() % width;
60}
61
62template <typename Point>
63void makeRandomRange( std::vector< Point >& X, int nb, int width )
64{
65 X.resize( nb );
66 for ( int i = 0; i < nb; i++ )
67 makeRandom( X[ i ], width );
68}
69
70template <typename ProjectedPoint, typename Point >
71void project( ProjectedPoint& pp, const Point& p, Dimension a )
72{
73 Dimension j = 0;
74 for ( Dimension i = 0; i < Point::dimension; i++ )
75 if ( i != a ) pp[ j++ ] = p[ i ];
76}
77
78template <typename ProjectedPoint, typename Point >
79void projectRange( std::vector< ProjectedPoint >& pp,
80 const std::vector< Point > & p, Dimension a )
81{
82 pp.resize( p.size() );
83 for ( size_t i = 0; i < p.size(); i++ )
84 project( pp[ i ], p[ i ], a );
85 std::sort( pp.begin(), pp.end() );
86 auto last = std::unique( pp.begin(), pp.end() );
87 pp.erase( last, pp.end() );
88}
89
90
91// @param width the width of the domain
92template <typename Space>
93bool
94checkSkelStarCvxHFullConvexity( int width )
95{
96 typedef typename Space::Integer Integer;
97 typedef DGtal::KhalimskySpaceND< Space::dimension, Integer > KSpace;
98 typedef DGtal::DigitalConvexity< KSpace > DConvexity;
99 typedef typename KSpace::Point Point;
100
101 // Generate a random polytope in the specified domain
102 Point lo = Point::zero;
103 Point hi = Point::diagonal( width );
104 DConvexity dconv( lo, hi );
105 std::vector< Point > X;
106 int nb = Space::dimension + rand() % 7;
107 makeRandomRange( X, nb, width );
108 auto C = dconv.StarCvxH( X );
109 auto E = dconv.envelope( X );
110 auto Y = C.extremaOfCells();
111 bool ok1 = dconv.isFullyConvex( E );
112 if ( ! ok1 )
113 trace.warning() << "FC*(X) is not fully convex !" << std::endl;
114 bool ok2 = dconv.isFullyConvex( Y );
115 if ( ! ok2 )
116 {
117 trace.warning() << "Extr(Star(CvxH(X))) is not fully convex !" << std::endl;
118 for ( auto p : Y ) std::cout << " " << p;
119 trace.warning() << std::endl;
120 }
121 bool ok3 = std::includes( Y.cbegin(), Y.cend(), E.cbegin(), E.cend() );
122 trace.info() << "#X=" << X.size()
123 << " #FC*(X)=" << E.size() << ( ok1 ? "/FC" : "/ERROR" )
124 << " #Extr(Star(CvxH(X)))=" << Y.size()
125 << ( ok2 ? "/FC" : "/ERROR" )
126 << ( ok3 ? " FC*(X) subset Extr(Star(CvxH(X)))" : " Inclusion ERROR" )
127 << std::endl;
128 return ok1 && ok2 && ok3;
129}
130
131// @param width the width of the domain
132template <typename Space>
133bool
134checkCvxHPlusHypercubeFullConvexity( int width )
135{
136 typedef typename Space::Integer Integer;
137 typedef DGtal::KhalimskySpaceND< Space::dimension, Integer > KSpace;
138 typedef DGtal::DigitalConvexity< KSpace > DConvexity;
139 typedef typename KSpace::Point Point;
140 // typedef std::vector<Point> PointRange;
141
142 // Generate a random polytope in the specified domain
143 Point lo = Point::zero;
144 Point hi = Point::diagonal( width );
145 DConvexity dconv( lo, hi );
146 std::vector< Point > X, XpH, Y;
147 int nb = Space::dimension + rand() % 7;
148 makeRandomRange( X, nb, width );
149 std::sort( X.begin(), X.end() );
150 XpH = X;
151 for ( Dimension k = 0; k < Space::dimension; k++ )
152 XpH = dconv.U( k, XpH );
153 auto P = dconv.makePolytope( XpH );
154 P.getPoints( Y );
155 auto E = dconv.envelope( X );
156 bool ok1 = dconv.isFullyConvex( E );
157 if ( ! ok1 )
158 trace.warning() << "FC*(X) is not fully convex !" << std::endl;
159 bool ok2 = dconv.isFullyConvex( Y );
160 if ( ! ok2 )
161 {
162 trace.warning() << "CvxH(X+H) cap Z^d is not fully convex !" << std::endl;
163 for ( auto p : Y ) std::cout << " " << p;
164 trace.warning() << std::endl;
165 }
166 bool ok3 = std::includes( Y.cbegin(), Y.cend(), E.cbegin(), E.cend() );
167 trace.info() << "#X=" << X.size()
168 << " #CvxH(X+H) cap Z^d=" << Y.size()
169 << ( ok2 ? "/FC" : "/ERROR" )
170 << " #FC*(X)=" << E.size() << ( ok1 ? "/FC" : "/ERROR" )
171 << ( ok3 ? " FC*(X) subset CvxH(X+H) cap Z^d"
172 : " FC*(X) not subset CvxH(X+H) cap Z^d" )
173 << std::endl;
174 return ok1 && ok2;
175}
176
177// @param width the width of the domain
178template <typename Space>
179bool
180checkProjectionFullConvexity( int width )
181{
182 typedef typename Space::Integer Integer;
183 typedef DGtal::KhalimskySpaceND< Space::dimension, Integer > KSpace;
184 typedef DGtal::DigitalConvexity< KSpace > DConvexity;
185 typedef typename KSpace::Point Point;
186 // typedef std::vector<Point> PointRange;
187 typedef DGtal::KhalimskySpaceND< Space::dimension-1, Integer > ProjKSpace;
188 typedef DGtal::DigitalConvexity< ProjKSpace > ProjDConvexity;
189 typedef typename ProjKSpace::Point ProjPoint;
190 // typedef std::vector<ProjPoint> ProjPointRange;
191
192 // Generate a random polytope in the specified domain
193 Point lo = Point::zero;
194 Point hi = Point::diagonal( width );
195 DConvexity dconv( lo, hi );
196 std::vector< Point > X;
197 int n = Space::dimension + rand() % 7;
198 makeRandomRange( X, n, width );
199 auto E = dconv.envelope( X );
200 unsigned int nb = 0;
201 unsigned int nb_ok = 0;
202 for ( Dimension a = 0; a < Space::dimension; a++ )
203 {
204 ProjPoint plo, phi;
205 project( plo, lo, a );
206 project( phi, hi, a );
207 ProjDConvexity pdconv( plo, phi );
208 std::vector< ProjPoint > PE;
209 projectRange( PE, E, a );
210 bool ok = pdconv.isFullyConvex( PE );
211 if ( !ok )
212 trace.warning() << "Projection is not fully convex !" << std::endl;
213 nb_ok += ok ? 1 : 0;
214 nb += 1;
215 std::cout << "#E=" << E.size() << " #Proj(E)=" << PE.size()
216 << (ok ? "/FC" : "/ERROR" ) << std::endl;
217 }
218 return nb_ok == nb;
219}
220
221template <typename Point>
222void displayPointRange2D( const std::vector< Point >& X )
223{
224 if ( X.empty() ) return;
225 Point lo = X[ 0 ];
226 Point hi = X[ 0 ];
227 for ( auto&& p : X ) { lo = lo.inf( p ); hi = hi.sup( p ); }
228 auto w = hi[ 0 ] - lo[ 0 ] + 1;
229 auto h = hi[ 1 ] - lo[ 1 ] + 1;
230 for ( int y = 0; y < h; y++ )
231 {
232 for ( int x = 0; x < w; x++ )
233 {
234 Point q( lo[ 0 ] + x, lo[ 1 ] + y );
235 std::cout << ( std::binary_search( X.cbegin(), X.cend(), q ) ? "*" : "." );
236 }
237 std::cout << std::endl;
238 }
239}
240// @param width the width of the domain
241template <typename Space>
242bool
243checkFullConvexityCharacterization( int width )
244{
245 typedef typename Space::Integer Integer;
246 typedef DGtal::KhalimskySpaceND< Space::dimension, Integer > KSpace;
247 typedef DGtal::DigitalConvexity< KSpace > DConvexity;
248 typedef typename KSpace::Point Point;
249 // typedef std::vector<Point> PointRange;
250 typedef DGtal::KhalimskySpaceND< Space::dimension-1, Integer > ProjKSpace;
251 typedef DGtal::DigitalConvexity< ProjKSpace > ProjDConvexity;
252 typedef typename ProjKSpace::Point ProjPoint;
253 // typedef std::vector<ProjPoint> ProjPointRange;
254
255 // Generate a random polytope in the specified domain
256 Point lo = Point::zero;
257 Point hi = Point::diagonal( width );
258 DConvexity dconv( lo, hi );
259 std::vector< Point > X, Y;
260 int n = Space::dimension + rand() % 17;
261 makeRandomRange( X, n, width );
262 auto P = dconv.makePolytope( X );
263 P.getPoints( Y );
264 const bool cvx = dconv.is0Convex( Y );
265 const bool fc = dconv.isFullyConvex( Y );
266 bool proj_fc = true;
267 std::cout << "#X=" << Y.size()
268 << " " << ( cvx ? "X C" : "X NC" )
269 << "/" << ( fc ? "X FC" : "X NFC" );
270 for ( Dimension a = 0; a < Space::dimension; a++ )
271 {
272 ProjPoint plo, phi;
273 project( plo, lo, a );
274 project( phi, hi, a );
275 ProjDConvexity pdconv( plo, phi );
276 std::vector< ProjPoint > PE;
277 projectRange( PE, Y, a );
278 if ( Space::dimension == 3 )
279 {
280 std::cout << std::endl;
281 displayPointRange2D( PE );
282 }
283 bool ok = pdconv.isFullyConvex( PE );
284 bool ok0 = pdconv.is0Convex( PE );
285 std::cout << "/" << a << ( ok0 ? ( ok ? "FC" : "NFC" ) : "NC" );
286 proj_fc = proj_fc && ok;
287 if ( fc && !ok )
288 trace.warning() << "Projection is not fully convex !" << std::endl;
289 }
290 if ( fc != proj_fc )
291 trace.warning() << "X is " << ( fc ? "FCvx" : "not FCvx" )
292 << "proj(X) " << ( proj_fc ? "FCvx" : "not FCvx" )
293 << std::endl;
294 else std::cout << ( fc ? " => FC ok" : " => Not FC ok" ) << std::endl;
295 return fc == proj_fc;
296}
297
298int main( )
299{
300 int NB_TEST1 = 5;
301 int NB_TEST2 = 5;
302 int NB_TEST3 = 5;
303 int NB_TEST4 = 25;
304 {
305 trace.beginBlock( "Check SkelStarCvxH(X) full convexity 2D" );
306 typedef DGtal::SpaceND< 2, int > Space;
307 unsigned int nb = 0;
308 unsigned int nb_ok = 0;
309 for ( int i = 0; i < NB_TEST1; i++ )
310 {
311 nb_ok += checkSkelStarCvxHFullConvexity< Space >( 100 ) ? 1 : 0;
312 nb += 1;
313 }
314 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
315 trace.endBlock();
316 }
317 {
318 trace.beginBlock( "Check SkelStarCvxH(X) full convexity 3D" );
319 typedef DGtal::SpaceND< 3, int > Space;
320 unsigned int nb = 0;
321 unsigned int nb_ok = 0;
322 for ( int i = 0; i < NB_TEST1; i++ )
323 {
324 nb_ok += checkSkelStarCvxHFullConvexity< Space >( 30 ) ? 1 : 0;
325 nb += 1;
326 }
327 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
328 trace.endBlock();
329 }
330 {
331 trace.beginBlock( "Check SkelStarCvxH(X) full convexity 4D" );
332 typedef DGtal::SpaceND< 4, int > Space;
333 unsigned int nb = 0;
334 unsigned int nb_ok = 0;
335 for ( int i = 0; i < NB_TEST1; i++ )
336 {
337 nb_ok += checkSkelStarCvxHFullConvexity< Space >( 10 ) ? 1 : 0;
338 nb += 1;
339 }
340 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
341 trace.endBlock();
342 }
343 {
344 trace.beginBlock( "Check Projection full convexity 3D" );
345 typedef DGtal::SpaceND< 3, int > Space;
346 unsigned int nb = 0;
347 unsigned int nb_ok = 0;
348 for ( int i = 0; i < NB_TEST2; i++ )
349 {
350 nb_ok += checkProjectionFullConvexity< Space >( 30 ) ? 1 : 0;
351 nb += 1;
352 }
353 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
354 trace.endBlock();
355 }
356 {
357 trace.beginBlock( "Check Projection full convexity 4D" );
358 typedef DGtal::SpaceND< 4, int > Space;
359 unsigned int nb = 0;
360 unsigned int nb_ok = 0;
361 for ( int i = 0; i < NB_TEST2; i++ )
362 {
363 nb_ok += checkProjectionFullConvexity< Space >( 10 ) ? 1 : 0;
364 nb += 1;
365 }
366 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
367 trace.endBlock();
368 }
369 {
370 trace.beginBlock( "Check CvxH plus Hypercube full convexity 2D" );
371 typedef DGtal::SpaceND< 2, int > Space;
372 unsigned int nb = 0;
373 unsigned int nb_ok = 0;
374 for ( int i = 0; i < NB_TEST3; i++ )
375 {
376 nb_ok += checkCvxHPlusHypercubeFullConvexity< Space >( 100 ) ? 1 : 0;
377 nb += 1;
378 }
379 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
380 trace.endBlock();
381 }
382 {
383 trace.beginBlock( "Check CvxH plus Hypercube full convexity 3D" );
384 typedef DGtal::SpaceND< 3, int > Space;
385 unsigned int nb = 0;
386 unsigned int nb_ok = 0;
387 for ( int i = 0; i < NB_TEST3; i++ )
388 {
389 nb_ok += checkCvxHPlusHypercubeFullConvexity< Space >( 30 ) ? 1 : 0;
390 nb += 1;
391 }
392 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
393 trace.endBlock();
394 }
395 {
396 trace.beginBlock( "Check CvxH plus Hypercube full convexity 4D" );
397 typedef DGtal::SpaceND< 4, int > Space;
398 unsigned int nb = 0;
399 unsigned int nb_ok = 0;
400 for ( int i = 0; i < NB_TEST3; i++ )
401 {
402 nb_ok += checkCvxHPlusHypercubeFullConvexity< Space >( 10 ) ? 1 : 0;
403 nb += 1;
404 }
405 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
406 trace.endBlock();
407 }
408 {
409 trace.beginBlock( "Check full convexity characterization 3D" );
410 typedef DGtal::SpaceND< 3, int > Space;
411 unsigned int nb = 0;
412 unsigned int nb_ok = 0;
413 for ( int i = 0; i < NB_TEST4; i++ )
414 {
415 nb_ok += checkFullConvexityCharacterization< Space >( 10 ) ? 1 : 0;
416 nb += 1;
417 }
418 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
419 trace.endBlock();
420 }
421 {
422 trace.beginBlock( "Check full convexity characterization 4D" );
423 typedef DGtal::SpaceND< 4, int > Space;
424 unsigned int nb = 0;
425 unsigned int nb_ok = 0;
426 for ( int i = 0; i < NB_TEST4; i++ )
427 {
428 nb_ok += checkFullConvexityCharacterization< Space >( 10 ) ? 1 : 0;
429 nb += 1;
430 }
431 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
432 trace.endBlock();
433 }
434 return 0;
435}
makeRandomRange
void makeRandomRange(std::vector< Point > &X, int nb, int width)
Definition checkFullConvexityTheorems.cpp:63
projectRange
void projectRange(std::vector< ProjectedPoint > &pp, const std::vector< Point > &p, Dimension a)
Definition checkFullConvexityTheorems.cpp:79
checkSkelStarCvxHFullConvexity
bool checkSkelStarCvxHFullConvexity(int width)
Definition checkFullConvexityTheorems.cpp:94
makeRandom
void makeRandom(Point &p, int width)
Definition checkFullConvexityTheorems.cpp:56
checkProjectionFullConvexity
bool checkProjectionFullConvexity(int width)
Definition checkFullConvexityTheorems.cpp:180
checkFullConvexityCharacterization
bool checkFullConvexityCharacterization(int width)
Definition checkFullConvexityTheorems.cpp:243
project
void project(ProjectedPoint &pp, const Point &p, Dimension a)
Definition checkFullConvexityTheorems.cpp:71
checkCvxHPlusHypercubeFullConvexity
bool checkCvxHPlusHypercubeFullConvexity(int width)
Definition checkFullConvexityTheorems.cpp:134
displayPointRange2D
void displayPointRange2D(const std::vector< Point > &X)
Definition checkFullConvexityTheorems.cpp:222
DGtal::BoundedLatticePolytope::getPoints
void getPoints(std::vector< Point > &pts) const
DGtal::DigitalConvexity::U
PointRange U(Dimension i, const PointRange &X) const
DGtal::DigitalConvexity::isFullyConvex
bool isFullyConvex(const PointRange &X, bool convex0=false) const
DGtal::DigitalConvexity::is0Convex
bool is0Convex(const PointRange &X) const
DGtal::DigitalConvexity::envelope
PointRange envelope(const PointRange &Z, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
DGtal::DigitalConvexity::StarCvxH
LatticeSet StarCvxH(const PointRange &X, Dimension axis=dimension) const
DGtal::DigitalConvexity::makePolytope
LatticePolytope makePolytope(const PointRange &X, bool make_minkowski_summable=false) const
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition KhalimskySpaceND.h:394
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::warning
std::ostream & warning()
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::endBlock
double endBlock()
phi
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector phi(const Face f)
Definition dgtalCalculus-bunny.cpp:69
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition ClosedIntegerHalfPlane.h:49
DGtal::Dimension
DGtal::uint32_t Dimension
Definition Common.h:119
DGtal::trace
Trace trace
std
STL namespace.
Point
MyPointD Point
Definition testClone2.cpp:381