DGtal 1.3.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 ( auto 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 typedef std::vector<Point> PointRange;
101
102 // Generate a random polytope in the specified domain
103 Point lo = Point::zero;
104 Point hi = Point::diagonal( width );
105 DConvexity dconv( lo, hi );
106 std::vector< Point > X;
107 int nb = Space::dimension + rand() % 7;
108 makeRandomRange( X, nb, width );
109 auto C = dconv.StarCvxH( X );
110 auto E = dconv.envelope( X );
111 auto Y = C.extremaOfCells();
112 bool ok1 = dconv.isFullyConvex( E );
113 if ( ! ok1 )
114 trace.warning() << "FC*(X) is not fully convex !" << std::endl;
115 bool ok2 = dconv.isFullyConvex( Y );
116 if ( ! ok2 )
117 {
118 trace.warning() << "Extr(Star(CvxH(X))) is not fully convex !" << std::endl;
119 for ( auto p : Y ) std::cout << " " << p;
120 trace.warning() << std::endl;
121 }
122 bool ok3 = std::includes( Y.cbegin(), Y.cend(), E.cbegin(), E.cend() );
123 trace.info() << "#X=" << X.size()
124 << " #FC*(X)=" << E.size() << ( ok1 ? "/FC" : "/ERROR" )
125 << " #Extr(Star(CvxH(X)))=" << Y.size()
126 << ( ok2 ? "/FC" : "/ERROR" )
127 << ( ok3 ? " FC*(X) subset Extr(Star(CvxH(X)))" : " Inclusion ERROR" )
128 << std::endl;
129 return ok1 && ok2 && ok3;
130}
131
132// @param width the width of the domain
133template <typename Space>
134bool
135checkCvxHPlusHypercubeFullConvexity( int width )
136{
137 typedef typename Space::Integer Integer;
138 typedef DGtal::KhalimskySpaceND< Space::dimension, Integer > KSpace;
139 typedef DGtal::DigitalConvexity< KSpace > DConvexity;
140 typedef typename KSpace::Point Point;
141 typedef std::vector<Point> PointRange;
142
143 // Generate a random polytope in the specified domain
144 Point lo = Point::zero;
145 Point hi = Point::diagonal( width );
146 DConvexity dconv( lo, hi );
147 std::vector< Point > X, XpH, Y;
148 int nb = Space::dimension + rand() % 7;
149 makeRandomRange( X, nb, width );
150 std::sort( X.begin(), X.end() );
151 XpH = X;
152 for ( Dimension k = 0; k < Space::dimension; k++ )
153 XpH = dconv.U( k, XpH );
154 auto P = dconv.makePolytope( XpH );
155 P.getPoints( Y );
156 auto E = dconv.envelope( X );
157 bool ok1 = dconv.isFullyConvex( E );
158 if ( ! ok1 )
159 trace.warning() << "FC*(X) is not fully convex !" << std::endl;
160 bool ok2 = dconv.isFullyConvex( Y );
161 if ( ! ok2 )
162 {
163 trace.warning() << "CvxH(X+H) cap Z^d is not fully convex !" << std::endl;
164 for ( auto p : Y ) std::cout << " " << p;
165 trace.warning() << std::endl;
166 }
167 bool ok3 = std::includes( Y.cbegin(), Y.cend(), E.cbegin(), E.cend() );
168 trace.info() << "#X=" << X.size()
169 << " #CvxH(X+H) cap Z^d=" << Y.size()
170 << ( ok2 ? "/FC" : "/ERROR" )
171 << " #FC*(X)=" << E.size() << ( ok1 ? "/FC" : "/ERROR" )
172 << ( ok3 ? " FC*(X) subset CvxH(X+H) cap Z^d"
173 : " FC*(X) not subset CvxH(X+H) cap Z^d" )
174 << std::endl;
175 return ok1 && ok2;
176}
177
178// @param width the width of the domain
179template <typename Space>
180bool
181checkProjectionFullConvexity( int width )
182{
183 typedef typename Space::Integer Integer;
184 typedef DGtal::KhalimskySpaceND< Space::dimension, Integer > KSpace;
185 typedef DGtal::DigitalConvexity< KSpace > DConvexity;
186 typedef typename KSpace::Point Point;
187 typedef std::vector<Point> PointRange;
188 typedef DGtal::KhalimskySpaceND< Space::dimension-1, Integer > ProjKSpace;
189 typedef DGtal::DigitalConvexity< ProjKSpace > ProjDConvexity;
190 typedef typename ProjKSpace::Point ProjPoint;
191 typedef std::vector<ProjPoint> ProjPointRange;
192
193 // Generate a random polytope in the specified domain
194 Point lo = Point::zero;
195 Point hi = Point::diagonal( width );
196 DConvexity dconv( lo, hi );
197 std::vector< Point > X;
198 int n = Space::dimension + rand() % 7;
199 makeRandomRange( X, n, width );
200 auto E = dconv.envelope( X );
201 unsigned int nb = 0;
202 unsigned int nb_ok = 0;
203 for ( Dimension a = 0; a < Space::dimension; a++ )
204 {
205 ProjPoint plo, phi;
206 project( plo, lo, a );
207 project( phi, hi, a );
208 ProjDConvexity pdconv( plo, phi );
209 std::vector< ProjPoint > PE;
210 projectRange( PE, E, a );
211 bool ok = pdconv.isFullyConvex( PE );
212 if ( !ok )
213 trace.warning() << "Projection is not fully convex !" << std::endl;
214 nb_ok += ok ? 1 : 0;
215 nb += 1;
216 std::cout << "#E=" << E.size() << " #Proj(E)=" << PE.size()
217 << (ok ? "/FC" : "/ERROR" ) << std::endl;
218 }
219 return nb_ok == nb;
220}
221
222template <typename Point>
223void displayPointRange2D( const std::vector< Point >& X )
224{
225 if ( X.empty() ) return;
226 Point lo = X[ 0 ];
227 Point hi = X[ 0 ];
228 for ( auto&& p : X ) { lo = lo.inf( p ); hi = hi.sup( p ); }
229 auto w = hi[ 0 ] - lo[ 0 ] + 1;
230 auto h = hi[ 1 ] - lo[ 1 ] + 1;
231 for ( int y = 0; y < h; y++ )
232 {
233 for ( int x = 0; x < w; x++ )
234 {
235 Point q( lo[ 0 ] + x, lo[ 1 ] + y );
236 std::cout << ( std::binary_search( X.cbegin(), X.cend(), q ) ? "*" : "." );
237 }
238 std::cout << std::endl;
239 }
240}
241// @param width the width of the domain
242template <typename Space>
243bool
244checkFullConvexityCharacterization( int width )
245{
246 typedef typename Space::Integer Integer;
247 typedef DGtal::KhalimskySpaceND< Space::dimension, Integer > KSpace;
248 typedef DGtal::DigitalConvexity< KSpace > DConvexity;
249 typedef typename KSpace::Point Point;
250 typedef std::vector<Point> PointRange;
251 typedef DGtal::KhalimskySpaceND< Space::dimension-1, Integer > ProjKSpace;
252 typedef DGtal::DigitalConvexity< ProjKSpace > ProjDConvexity;
253 typedef typename ProjKSpace::Point ProjPoint;
254 typedef std::vector<ProjPoint> ProjPointRange;
255
256 // Generate a random polytope in the specified domain
257 Point lo = Point::zero;
258 Point hi = Point::diagonal( width );
259 DConvexity dconv( lo, hi );
260 std::vector< Point > X, Y;
261 int n = Space::dimension + rand() % 17;
262 makeRandomRange( X, n, width );
263 auto P = dconv.makePolytope( X );
264 P.getPoints( Y );
265 const bool cvx = dconv.is0Convex( Y );
266 const bool fc = dconv.isFullyConvex( Y );
267 bool proj_fc = true;
268 std::cout << "#X=" << Y.size()
269 << " " << ( cvx ? "X C" : "X NC" )
270 << "/" << ( fc ? "X FC" : "X NFC" );
271 for ( Dimension a = 0; a < Space::dimension; a++ )
272 {
273 ProjPoint plo, phi;
274 project( plo, lo, a );
275 project( phi, hi, a );
276 ProjDConvexity pdconv( plo, phi );
277 std::vector< ProjPoint > PE;
278 projectRange( PE, Y, a );
279 if ( Space::dimension == 3 )
280 {
281 std::cout << std::endl;
282 displayPointRange2D( PE );
283 }
284 bool ok = pdconv.isFullyConvex( PE );
285 bool ok0 = pdconv.is0Convex( PE );
286 std::cout << "/" << a << ( ok0 ? ( ok ? "FC" : "NFC" ) : "NC" );
287 proj_fc = proj_fc && ok;
288 if ( fc && !ok )
289 trace.warning() << "Projection is not fully convex !" << std::endl;
290 }
291 if ( fc != proj_fc )
292 trace.warning() << "X is " << ( fc ? "FCvx" : "not FCvx" )
293 << "proj(X) " << ( proj_fc ? "FCvx" : "not FCvx" )
294 << std::endl;
295 else std::cout << ( fc ? " => FC ok" : " => Not FC ok" ) << std::endl;
296 return fc == proj_fc;
297}
298
299int main( int argc, char* argv[] )
300{
301 int NB_TEST1 = 5;
302 int NB_TEST2 = 5;
303 int NB_TEST3 = 5;
304 int NB_TEST4 = 25;
305 {
306 trace.beginBlock( "Check SkelStarCvxH(X) full convexity 2D" );
307 typedef DGtal::SpaceND< 2, int > Space;
308 unsigned int nb = 0;
309 unsigned int nb_ok = 0;
310 for ( int i = 0; i < NB_TEST1; i++ )
311 {
312 nb_ok += checkSkelStarCvxHFullConvexity< Space >( 100 ) ? 1 : 0;
313 nb += 1;
314 }
315 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
316 trace.endBlock();
317 }
318 {
319 trace.beginBlock( "Check SkelStarCvxH(X) full convexity 3D" );
320 typedef DGtal::SpaceND< 3, int > Space;
321 unsigned int nb = 0;
322 unsigned int nb_ok = 0;
323 for ( int i = 0; i < NB_TEST1; i++ )
324 {
325 nb_ok += checkSkelStarCvxHFullConvexity< Space >( 30 ) ? 1 : 0;
326 nb += 1;
327 }
328 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
329 trace.endBlock();
330 }
331 {
332 trace.beginBlock( "Check SkelStarCvxH(X) full convexity 4D" );
333 typedef DGtal::SpaceND< 4, int > Space;
334 unsigned int nb = 0;
335 unsigned int nb_ok = 0;
336 for ( int i = 0; i < NB_TEST1; i++ )
337 {
338 nb_ok += checkSkelStarCvxHFullConvexity< Space >( 10 ) ? 1 : 0;
339 nb += 1;
340 }
341 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
342 trace.endBlock();
343 }
344 {
345 trace.beginBlock( "Check Projection full convexity 3D" );
346 typedef DGtal::SpaceND< 3, int > Space;
347 unsigned int nb = 0;
348 unsigned int nb_ok = 0;
349 for ( int i = 0; i < NB_TEST2; i++ )
350 {
351 nb_ok += checkProjectionFullConvexity< Space >( 30 ) ? 1 : 0;
352 nb += 1;
353 }
354 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
355 trace.endBlock();
356 }
357 {
358 trace.beginBlock( "Check Projection full convexity 4D" );
359 typedef DGtal::SpaceND< 4, int > Space;
360 unsigned int nb = 0;
361 unsigned int nb_ok = 0;
362 for ( int i = 0; i < NB_TEST2; i++ )
363 {
364 nb_ok += checkProjectionFullConvexity< Space >( 10 ) ? 1 : 0;
365 nb += 1;
366 }
367 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
368 trace.endBlock();
369 }
370 {
371 trace.beginBlock( "Check CvxH plus Hypercube full convexity 2D" );
372 typedef DGtal::SpaceND< 2, int > Space;
373 unsigned int nb = 0;
374 unsigned int nb_ok = 0;
375 for ( int i = 0; i < NB_TEST3; i++ )
376 {
377 nb_ok += checkCvxHPlusHypercubeFullConvexity< Space >( 100 ) ? 1 : 0;
378 nb += 1;
379 }
380 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
381 trace.endBlock();
382 }
383 {
384 trace.beginBlock( "Check CvxH plus Hypercube full convexity 3D" );
385 typedef DGtal::SpaceND< 3, int > Space;
386 unsigned int nb = 0;
387 unsigned int nb_ok = 0;
388 for ( int i = 0; i < NB_TEST3; i++ )
389 {
390 nb_ok += checkCvxHPlusHypercubeFullConvexity< Space >( 30 ) ? 1 : 0;
391 nb += 1;
392 }
393 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
394 trace.endBlock();
395 }
396 {
397 trace.beginBlock( "Check CvxH plus Hypercube full convexity 4D" );
398 typedef DGtal::SpaceND< 4, int > Space;
399 unsigned int nb = 0;
400 unsigned int nb_ok = 0;
401 for ( int i = 0; i < NB_TEST3; i++ )
402 {
403 nb_ok += checkCvxHPlusHypercubeFullConvexity< Space >( 10 ) ? 1 : 0;
404 nb += 1;
405 }
406 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
407 trace.endBlock();
408 }
409 {
410 trace.beginBlock( "Check full convexity characterization 3D" );
411 typedef DGtal::SpaceND< 3, int > Space;
412 unsigned int nb = 0;
413 unsigned int nb_ok = 0;
414 for ( int i = 0; i < NB_TEST4; i++ )
415 {
416 nb_ok += checkFullConvexityCharacterization< Space >( 10 ) ? 1 : 0;
417 nb += 1;
418 }
419 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
420 trace.endBlock();
421 }
422 {
423 trace.beginBlock( "Check full convexity characterization 4D" );
424 typedef DGtal::SpaceND< 4, int > Space;
425 unsigned int nb = 0;
426 unsigned int nb_ok = 0;
427 for ( int i = 0; i < NB_TEST4; i++ )
428 {
429 nb_ok += checkFullConvexityCharacterization< Space >( 10 ) ? 1 : 0;
430 nb += 1;
431 }
432 trace.info() << nb_ok << "/" << nb << " OK tests" << std::endl;
433 trace.endBlock();
434 }
435 return 0;
436}
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:181
checkFullConvexityCharacterization
bool checkFullConvexityCharacterization(int width)
Definition: checkFullConvexityTheorems.cpp:244
project
void project(ProjectedPoint &pp, const Point &p, Dimension a)
Definition: checkFullConvexityTheorems.cpp:71
checkCvxHPlusHypercubeFullConvexity
bool checkCvxHPlusHypercubeFullConvexity(int width)
Definition: checkFullConvexityTheorems.cpp:135
displayPointRange2D
void displayPointRange2D(const std::vector< Point > &X)
Definition: checkFullConvexityTheorems.cpp:223
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::SpaceND::Integer
TInteger Integer
Arithmetic ring induced by (+,-,*) and Integer numbers.
Definition: SpaceND.h:102
DGtal::SpaceND::dimension
static const Dimension dimension
static constants to store the dimension.
Definition: SpaceND.h:132
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()
PointRange
std::vector< Point > PointRange
Definition: digitalPolyhedronBuilder3D.cpp:84
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:137
DGtal::trace
Trace trace
Definition: Common.h:154
std
STL namespace.
main
int main()
Definition: testBits.cpp:56
Point
MyPointD Point
Definition: testClone2.cpp:383
Generated on Sun Nov 27 2022 15:15:32 for DGtal by doxygen 1.9.5