DGtal  1.5.beta
QuickHull.h
1 
17 #pragma once
18 
31 #if defined(QuickHull_RECURSES)
32 #error Recursive header files inclusion detected in QuickHull.h
33 #else // defined(QuickHull_RECURSES)
35 #define QuickHull_RECURSES
36 
37 #if !defined QuickHull_h
39 #define QuickHull_h
40 
42 // Inclusions
43 #include <iostream>
44 #include <string>
45 #include <vector>
46 #include <queue>
47 #include <set>
48 #include "DGtal/base/Common.h"
49 #include "DGtal/base/Clock.h"
50 #include "DGtal/geometry/tools/QuickHullKernels.h"
51 
52 namespace DGtal
53 {
55  // template class QuickHull
56 
138  template < typename TKernel >
139  struct QuickHull
140  {
141  typedef TKernel Kernel;
142  typedef typename Kernel::CoordinatePoint Point;
143  typedef typename Kernel::CoordinateVector Vector;
144  typedef typename Kernel::CoordinateScalar Scalar;
145  typedef typename Kernel::InternalScalar InternalScalar;
146  typedef std::size_t Index;
147  typedef std::size_t Size;
148  BOOST_STATIC_ASSERT(( Point::dimension == Vector::dimension ));
149  typedef std::vector< Index > IndexRange;
150  typedef typename Kernel::HalfSpace HalfSpace;
151  typedef typename Kernel::CombinatorialPlaneSimplex CombinatorialPlaneSimplex;
152  static const Size dimension = Point::dimension;
153 
155  enum { UNASSIGNED = (Index) -1 };
156 
160  struct Facet {
166 
167  Facet() = default;
168  Facet( const Facet& ) = default;
169  Facet( Facet&& ) = default;
170  Facet& operator=( Facet&& ) = default;
171  Facet& operator=( const Facet& ) = default;
172  Facet( const HalfSpace& aH, Index b )
173  : H( aH ), below( b ) {}
174 
175  void clear()
176  {
177  H = HalfSpace();
178  neighbors.clear();
179  outside_set.clear();
180  on_set.clear();
181  below = UNASSIGNED;
182  }
183  void addPointOn( Index p )
184  {
185  const auto it = std::find( on_set.cbegin(), on_set.cend(), p );
186  if ( it == on_set.cend() ) on_set.push_back( p );
187  }
188  void display( std::ostream& out ) const
189  {
190  const auto N = H.internalNormal();
191  out << "[Facet iN=(" << N[0];
192  for ( Dimension i = 1; i < N.dimension; i++ ) out << "," << N[ i ];
193  out << ") c=" << H.internalIntercept() << " b=" << below << " n={";
194  for ( auto&& n : neighbors ) out << " " << n;
195  out << " } #out=" << outside_set.size();
196  out << " on={";
197  for ( auto&& n : on_set ) out << " " << n;
198  out << " }]" << std::endl;
199  }
200 
201  void addNeighbor( Index n )
202  {
203  const auto it = std::find( neighbors.cbegin(), neighbors.cend(), n );
204  if ( it == neighbors.cend() ) neighbors.push_back( n );
205  }
206  void subNeighbor( Index n )
207  {
208  auto it = std::find( neighbors.begin(), neighbors.end(), n );
209  if ( it != neighbors.end() ) {
210  std::swap( *it, neighbors.back() );
211  neighbors.pop_back();
212  }
213  }
214  void swap( Facet& other )
215  {
216  if ( this != &other ) {
217  std::swap( H, other.H );
218  neighbors.swap ( other.neighbors );
219  outside_set.swap( other.outside_set );
220  on_set.swap ( other.on_set );
221  std::swap( below, other.below );
222  }
223  }
225  {
226  Size M;
227  M += neighbors.capacity() * sizeof( Index );
228  M += outside_set.capacity() * sizeof( Index );
229  M += on_set.capacity() * sizeof( Index );
230  return M;
231  }
232  };
233 
236  typedef std::pair< Index, Index > Ridge;
237 
239  enum class Status {
240  Uninitialized = 0,
241  InputInitialized = 1,
242  SimplexCompleted = 2,
243  FacetsCompleted = 3,
244  VerticesCompleted = 4,
245  AllCompleted = 5,
246  NotFullDimensional= 10,
247  InvalidRidge = 11,
248  InvalidConvexHull = 12
249  };
250 
251  // ----------------------- standard services --------------------------
252  public:
255 
259  QuickHull( const Kernel& K = Kernel(), int dbg = 0 )
260  : kernel( K ), debug_level( dbg ), myStatus( Status::Uninitialized )
261  {}
262 
266  Status status() const
267  { return myStatus; }
268 
270  void clear()
271  {
273  points.clear();
274  processed_points.clear();
275  input2comp.clear();
276  comp2input.clear();
277  assignment.clear();
278  facets.clear();
279  deleted_facets.clear();
280  p2v.clear();
281  v2p.clear();
282  timings.clear();
283  }
284 
285 
287  Size memory() const
288  {
289  // int debug_level;
290  Size M = sizeof( kernel ) + sizeof( int );
291  // std::vector< Point > points;
292  M += sizeof( std::vector< Point > )
293  + points.capacity() * sizeof( Point );
294  M += sizeof( std::vector< Index > )
295  + processed_points.capacity() * sizeof( Index );
296  M += sizeof( std::vector< Index > )
297  + input2comp.capacity() * sizeof( Index );
298  M += sizeof( std::vector< Index > )
299  + comp2input.capacity() * sizeof( Index );
300  // std::vector< Index > assignment;
301  M += sizeof( std::vector< Index > )
302  + assignment.capacity() * sizeof( Index );
303  // std::vector< Facet > facets;
304  M += sizeof( std::vector< Facet > )
305  + facets.capacity() * sizeof( Facet );
306  for ( const auto& f : facets ) M += f.variableMemory();
307  // std::set< Index > deleted_facets;
308  M += sizeof( std::set< Index > )
309  + deleted_facets.size() * ( sizeof( Index ) + 2*sizeof(Index*) );
310  // IndexRange p2v;
311  M += sizeof( std::vector< Index > )
312  + p2v.capacity() * sizeof( Index );
313  // IndexRange v2p;
314  M += sizeof( std::vector< Index > )
315  + v2p.capacity() * sizeof( Index );
316  // std::vector< double > timings;
317  M += sizeof( std::vector< double > )
318  + timings.capacity() * sizeof( double );
319  return M;
320  }
321 
323  Size nbPoints() const
324  { return points.size(); }
325 
328  Size nbFacets() const
329  {
330  ASSERT( status() >= Status::FacetsCompleted
331  && status() <= Status::AllCompleted );
332  return facets.size();
333  }
334 
337  Size nbVertices() const
338  {
339  ASSERT( status() >= Status::VerticesCompleted
340  && status() <= Status::AllCompleted );
341  return v2p.size();
342  }
343 
347  {
348  ASSERT( status() >= Status::VerticesCompleted
349  && status() <= Status::AllCompleted );
350  return nb_finite_facets;
351  }
352 
356  {
357  ASSERT( status() >= Status::VerticesCompleted
358  && status() <= Status::AllCompleted );
359  return nb_infinite_facets;
360  }
361 
363  // -------------------------- Convex hull services ----------------------------
364  public:
367 
381  template < typename InputPoint >
382  bool setInput( const std::vector< InputPoint >& input_points,
383  bool remove_duplicates = true )
384  {
385  Clock tic;
386  tic.startClock();
387  clear();
388  timings.clear();
389  kernel.makeInput( points, input2comp, comp2input,
390  input_points, remove_duplicates );
391  timings.push_back( tic.stopClock() );
392  if ( points.size() <= dimension ) {
394  return false;
395  }
397  return true;
398  }
399 
410  bool setInitialSimplex( const IndexRange& full_splx )
411  {
412  if ( status() != Status::InputInitialized ) return false;
413  if ( full_splx.size() != dimension + 1 )
414  {
415  trace.error() << "[QuickHull::setInitialSimplex]"
416  << " not a full dimensional simplex" << std::endl;
418  return false;
419  }
421  for ( Index j = 0; j < dimension; ++j )
422  splx[ j ] = full_splx[ j ];
423  const auto H = kernel.compute( points, splx, full_splx.back() );
424  const auto volume = kernel.volume( H, points[ full_splx.back() ] );
425  if ( volume > 0 )
426  return computeSimplexConfiguration( full_splx );
428  return false;
429  }
430 
432  // -------------------------- Convex hull services ----------------------------
433  public:
436 
461  {
462  if ( target < Status::InputInitialized || target > Status::AllCompleted )
463  return false;
464  Clock tic;
465  if ( status() == Status::InputInitialized )
466  { // Initialization
467  tic.startClock();
468  bool ok1 = computeInitialSimplex();
469  timings.push_back( tic.stopClock() );
470  if ( ! ok1 ) return false;
471  if ( status() == target ) return true;
472  }
473  if ( status() == Status::SimplexCompleted )
474  { // Computes facets
475  tic.startClock();
476  bool ok2 = computeFacets();
477  timings.push_back( tic.stopClock() );
478  if ( ! ok2 ) return false;
479  if ( status() == target ) return true;
480  }
481  if ( status() == Status::FacetsCompleted )
482  { // Computes vertices
483  tic.startClock();
484  bool ok3 = computeVertices();
485  timings.push_back( tic.stopClock() );
486  if ( ! ok3 ) return false;
487  if ( status() == target ) return true;
488  }
489  if ( target == Status::AllCompleted
491  { // for now, Status::VerticesCompleted and
492  // Status::AllCompleted are the same.
494  return true;
495  }
496  return false;
497  }
498 
505  {
506  const auto full_simplex = pickInitialSimplex();
507  if ( full_simplex.empty() ) {
509  return false;
510  }
511  return computeSimplexConfiguration( full_simplex );
512  }
513 
524  {
525  if ( status() != Status::SimplexCompleted ) return false;
526  std::queue< Index > Q;
527  for ( Index fi = 0; fi < facets.size(); ++fi )
528  Q.push( fi );
529  Index n = 0;
530  while ( processFacet( Q ) ) {
531  if ( debug_level >= 1 )
532  trace.info() << "---- Iteration " << n++ << " #Q=" << Q.size() << std::endl;
533  }
534  cleanFacets();
535  if ( debug_level >= 2 ) {
536  trace.info() << ".... #facets=" << facets.size()
537  << " #deleted=" << deleted_facets.size() << std::endl;
538  }
540  return true;
541  }
542 
554  {
555  static const int MAX_NB_VPF = 10 * dimension;
556  if ( status() != Status::FacetsCompleted ) return false;
557 
558  // Renumber infinite facets in case of Delaunay triangulation computation.
560 
561  // Builds the maps v2p: vertex -> point, and p2v : point -> vertex.
562  facet_counter = IndexRange( MAX_NB_VPF, 0 );
563  v2p.clear();
564  p2v.resize( points.size() );
565  std::vector< IndexRange > p2f( points.size() );
566  for ( Index f = 0; f < facets.size(); ++f ) {
567  for ( auto&& p : facets[ f ].on_set ) p2f[ p ].push_back( f );
568  }
569 
570  // vertices belong to at least d facets
571  Index v = 0;
572  for ( Index p = 0; p < points.size(); ++p ) {
573  const auto nbf = p2f[ p ].size();
574  facet_counter[ std::min( (int) nbf, MAX_NB_VPF-1 ) ] += 1;
575  if ( nbf >= dimension ) {
576  v2p.push_back( p );
577  p2v[ p ] = v++;
578  }
579  else p2v[ p ] = UNASSIGNED;
580  }
581 
582  // Display debug informations
583  if ( debug_level >= 1 ) {
584  trace.info() << "#vertices=" << v2p.size() << " #facets=" << facets.size()
585  << std::endl;
586  trace.info() << "#inc_facets/point= ";
587  for ( auto n : facet_counter ) trace.info() << n << " ";
588  trace.info() << std::endl;
589  }
591  return true;
592  }
593 
594 
596  // -------------------------- Output services ----------------------------
597  public:
600 
611  template < typename OutputPoint >
612  bool getVertexPositions( std::vector< OutputPoint >& vertex_positions )
613  {
614  vertex_positions.clear();
615  if ( ! ( status() >= Status::VerticesCompleted
616  && status() <= Status::AllCompleted ) ) return false;
617  vertex_positions.resize( v2p.size() );
618  for ( Index i = 0; i < v2p.size(); i++ ) {
619  kernel.convertPointTo( points[ v2p[ i ] ], vertex_positions[ i ] );
620  }
621  return true;
622  }
623 
632  bool getVertex2Point( IndexRange& vertex_to_point )
633  {
634  vertex_to_point.clear();
635  if ( ! ( status() >= Status::VerticesCompleted
636  && status() <= Status::AllCompleted ) ) return false;
637  vertex_to_point = v2p;
638  return true;
639  }
640 
651  bool getPoint2Vertex( IndexRange& point_to_vertex )
652  {
653  point_to_vertex.clear();
654  if ( ! ( status() >= Status::VerticesCompleted
655  && status() <= Status::AllCompleted ) ) return false;
656  point_to_vertex = p2v;
657  return true;
658  }
659 
666  bool getFacetVertices( std::vector< IndexRange >& facet_vertices ) const
667  {
668  facet_vertices.clear();
669  if ( ! ( status() >= Status::VerticesCompleted
670  && status() <= Status::AllCompleted ) ) return false;
671  facet_vertices.reserve( nbFacets() );
672  for ( Index f = 0; f < nbFacets(); ++f ) {
673  IndexRange ofacet = orientedFacetPoints( f );
674  for ( auto& v : ofacet ) v = p2v[ v ];
675  facet_vertices.push_back( ofacet );
676  }
677  return true;
678  }
679 
690  bool getFacetHalfSpaces( std::vector< HalfSpace >& facet_halfspaces )
691  {
692  facet_halfspaces.clear();
693  if ( ! ( status() >= Status::FacetsCompleted
694  && status() <= Status::AllCompleted ) ) return false;
695  facet_halfspaces.reserve( nbFacets() );
696  for ( Index f = 0; f < nbFacets(); ++f ) {
697  facet_halfspaces.push_back( facets[ f ].H );
698  }
699  return true;
700  }
701 
702 
704  // -------------------------- Check hull services ----------------------------
705  public:
708 
715  bool check()
716  {
717  bool ok = true;
719  trace.warning() << "[Quickhull::check] invalid status="
720  << (int)status() << std::endl;
721  ok = false;
722  }
723  if ( processed_points.size() != points.size() ) {
724  trace.warning() << "[Quickhull::check] not all points processed: "
725  << processed_points.size() << " / " << points.size()
726  << std::endl;
727  ok = false;
728  }
729  if ( ! checkHull() ) {
730  trace.warning() << "[Quickhull::check] Check hull is invalid. "
731  << std::endl;
732  ok = false;
733  }
734  if ( ! checkFacets() ) {
735  trace.warning() << "[Quickhull::check] Check facets is invalid. "
736  << std::endl;
737  ok = false;
738  }
739  return ok;
740  }
741 
743  bool checkFacets()
744  {
745  Size nb = 0;
746  Size nbok = 0;
747  for ( Index f = 0; f < facets.size(); ++f )
748  if ( deleted_facets.count( f ) == 0 ) {
749  bool ok = checkFacet( f );
750  nbok += ok ? 1 : 0;
751  nb += 1;
752  }
753  return nb == nbok;
754  }
755 
762  bool checkHull()
763  {
764  Index nbok = 0;
765  // for ( Index v = 0; v < points.size(); ++v ) {
766  for ( auto v : processed_points ) {
767  bool ok = true;
768  for ( Index f = 0; f < facets.size(); ++f )
769  if ( deleted_facets.count( f ) == 0 ) {
770  if ( above( facets[ f ], points[ v ] ) ) {
771  ok = false;
772  trace.error() << "- bad vertex " << v << " " << points[ v ]
773  << " dist="
774  << ( kernel.height( facets[ f ].H, points[ v ] ) )
775  << std::endl;
776  break;
777  }
778  }
779  nbok += ok ? 1 : 0;
780  }
781  if ( debug_level >= 2 ) {
782  trace.info() << nbok << "/"
783  << processed_points.size() << " vertices inside convex hull."
784  << std::endl;
785  }
786  if ( nbok != processed_points.size() ) myStatus = Status::InvalidConvexHull;
787  return nbok == processed_points.size();
788  }
789 
791  // ------------------------ public datas --------------------------
792  public:
795 
796  public:
798  mutable Kernel kernel;
802  std::vector< Point > points;
812  std::vector< Index > assignment;
814  std::vector< Facet > facets;
816  std::set< Index > deleted_facets;
826  std::vector< double > timings;
828  std::vector< Size > facet_counter;
829 
831  // ------------------------ protected datas --------------------------
832  protected:
835 
840 
842  // --------------------- protected services --------------------------
843  protected:
846 
850  InternalScalar height( const Facet& F, const Point& p ) const
851  { return kernel.height( F.H, p ); }
852 
856  bool above( const Facet& F, const Point& p ) const
857  { return kernel.above( F.H, p ); }
858 
862  bool aboveOrOn( const Facet& F, const Point& p ) const
863  { return kernel.aboveOrOn( F.H, p ); }
864 
868  bool on( const Facet& F, const Point& p ) const
869  { return kernel.on( F.H, p ); }
870 
873  void cleanFacets()
874  {
875  if ( deleted_facets.empty() ) return;
876  IndexRange renumbering( facets.size() );
877  Index i = 0;
878  Index j = 0;
879  for ( auto& l : renumbering ) {
880  if ( ! deleted_facets.count( j ) ) l = i++;
881  else l = UNASSIGNED;
882  j++;
883  }
884  const Index nf = facets.size() - deleted_facets.size();
885  deleted_facets.clear();
886  for ( Index f = 0; f < facets.size(); f++ )
887  if ( ( renumbering[ f ] != UNASSIGNED ) && ( f != renumbering[ f ] ) )
888  facets[ renumbering[ f ] ] = facets[ f ];
889  facets.resize( nf );
890  for ( auto& F : facets ) {
891  for ( auto& N : F.neighbors ) {
892  if ( renumbering[ N ] == UNASSIGNED )
893  trace.error() << "Invalid deleted neighboring facet." << std::endl;
894  else N = renumbering[ N ];
895  }
896  }
897  }
898 
902  {
903  if ( ! kernel.hasInfiniteFacets() ) return;
904  IndexRange renumbering( facets.size() );
905  Index i = 0;
906  Index j = 0;
907  for ( auto& l : renumbering ) {
908  if ( ! kernel.isHalfSpaceFacetInfinite( facets[ j ].H ) ) l = i++;
909  else l = UNASSIGNED;
910  j++;
911  }
912  const Index nf = i;
913  for ( Index f = 0; f < facets.size(); f++ )
914  if ( ( renumbering[ f ] != UNASSIGNED ) && ( f != renumbering[ f ] ) )
915  facets[ renumbering[ f ] ] = facets[ f ];
916  facets.resize( nf );
917  for ( auto& F : facets ) {
918  for ( auto& N : F.neighbors ) {
919  N = renumbering[ N ];
920  }
921  }
922  }
923 
927  {
928  nb_finite_facets = facets.size();
929  nb_infinite_facets = 0;
930  if ( ! kernel.hasInfiniteFacets() ) return;
931  IndexRange renumbering( facets.size() );
932  Index i = 0;
933  Index k = facets.size();
934  Index j = 0;
935  for ( auto& l : renumbering ) {
936  if ( ! kernel.isHalfSpaceFacetInfinite( facets[ j ].H ) ) l = i++;
937  else l = --k;
938  j++;
939  }
940  if ( i != k )
941  trace.error() << "[Quickhull::renumberInfiniteFacets]"
942  << " Error renumbering infinite facets "
943  << " up finite=" << i << " low infinite=" << k << std::endl;
944  std::vector< Facet > new_facets( facets.size() );
945  for ( Index f = 0; f < facets.size(); f++ )
946  new_facets[ renumbering[ f ] ].swap( facets[ f ] );
947  facets.swap( new_facets );
948  for ( auto& F : facets ) {
949  for ( auto& N : F.neighbors ) {
950  N = renumbering[ N ];
951  }
952  }
953  // Assign correct number of facets.
954  nb_finite_facets = i;
956  }
957 
958 
969  bool processFacet( std::queue< Index >& Q )
970  {
971  // If Q empty, we are done
972  if ( Q.empty() ) return false;
973  Index F = Q.front();
974  Q.pop();
975  // If F is already deleted, proceed to next in queue.
976  if ( deleted_facets.count( F ) ) return true;
977  // Take car of current facet.
978  const Facet& facet = facets[ F ];
979  if ( debug_level >= 3 ) {
980  trace.info() << "---------------------------------------------" << std::endl;
981  trace.info() << "---- ACTIVE FACETS---------------------------" << std::endl;
982  bool ok = true;
983  for ( Index i = 0; i < facets.size(); i++ )
984  if ( ! deleted_facets.count( i ) ) {
985  trace.info() << "- facet " << i << " ";
986  facets[ i ].display( trace.info() );
987  ok = ok && checkFacet( i );
988  }
989  if ( ! ok ) { // should not happen.
991  return false;
992  }
993  }
994  if ( debug_level >= 2 ) {
995  trace.info() << "---------------------------------------------" << std::endl;
996  trace.info() << "Processing facet " << F << " ";
997  facet.display( trace.info() );
998  }
999  if ( facet.outside_set.empty() ) return true;
1000  // Selects furthest vertex
1001  Index furthest_v = facet.outside_set[ 0 ];
1002  auto furthest_h = height( facet, points[ furthest_v ] );
1003  for ( Index v = 1; v < facet.outside_set.size(); v++ ) {
1004  auto h = height( facet, points[ facet.outside_set[ v ] ] );
1005  if ( h > furthest_h ) {
1006  furthest_h = h;
1007  furthest_v = facet.outside_set[ v ];
1008  }
1009  }
1010  const Point& p = points[ furthest_v ];
1011  // Extracts Visible facets V and Horizon Ridges H
1012  std::vector< Index > V; // visible facets
1013  std::set< Index > M; // marked facets (are in E or were in E)
1014  std::queue< Index > E; // queue to extract visible facets
1015  std::vector< Ridge > H; // visible facets
1016  E.push ( F );
1017  M.insert( F );
1018  while ( ! E.empty() ) {
1019  Index G = E.front(); E.pop();
1020  V.push_back( G );
1021  for ( auto& N : facets[ G ].neighbors ) {
1022  if ( aboveOrOn( facets[ N ], p ) ) {
1023  if ( M.count( N ) ) continue;
1024  E.push( N );
1025  } else {
1026  H.push_back( { G, N } );
1027  }
1028  M.insert( N );
1029  }
1030  } // while ( ! E.empty() )
1031  if ( debug_level >= 1 ) {
1032  trace.info() << "#Visible=" << V.size() << " #Horizon=" << H.size()
1033  << " furthest_v=" << furthest_v << std::endl;
1034  }
1035  // Create new facets
1036  IndexRange new_facets;
1037  // For each ridge R in H
1038  for ( Index i = 0; i < H.size(); i++ )
1039  {
1040  // Create a new facet from R and p
1041  IndexRange ridge = pointsOnRidge( H[ i ] );
1042  if ( debug_level >= 3 ) {
1043  trace.info() << "Ridge (" << H[i].first << "," << H[i].second << ") = {";
1044  for ( auto&& r : ridge ) trace.info() << " " << r;
1045  trace.info() << " } furthest_v=" << furthest_v << std::endl;
1046  }
1047  IndexRange base( 1 + ridge.size() );
1048  Index j = 0;
1049  base[ j++ ] = furthest_v;
1050  for ( auto&& v : ridge ) base[ j++ ] = v;
1051  if ( j < dimension ) {
1052  trace.error() << "Bad ridge between " << std::endl
1053  << "- facet " << H[i].first << " ";
1054  facets[ H[i].first ].display( trace.error() );
1055  trace.error() << "- facet " << H[i].second << " ";
1056  facets[ H[i].second ].display( trace.error() );
1057  }
1058  Index nf = newFacet();
1059  new_facets.push_back( nf );
1060  facets[ nf ] = makeFacet( base, facets[ H[i].first ].below );
1061  facets[ nf ].on_set = IndexRange { base.cbegin(), base.cend() };
1062  std::sort( facets[ nf ].on_set.begin(), facets[ nf ].on_set.end() );
1063  makeNeighbors( nf, H[ i ].second );
1064  if ( debug_level >= 3 ) {
1065  trace.info() << "* New facet " << nf << " ";
1066  facets[ nf ].display( trace.info() );
1067  }
1068  // Checks that the facet is not parallel to another in the Horizon
1069  for ( Index k = 0; k < new_facets.size() - 1; k++ )
1070  if ( areFacetsParallel( new_facets[ k ], nf ) ) {
1071  if ( debug_level >= 1 ) {
1072  trace.info() << "Facets " << new_facets[ k ] << " and " << nf
1073  << " are parallel => merge." << std::endl;
1074  }
1075  mergeParallelFacets( new_facets[ k ], nf );
1076  new_facets.pop_back();
1077  deleteFacet( nf );
1078  if ( debug_level >= 3 ) {
1079  facets[ new_facets[ k ] ].display( trace.info() );
1080  }
1081  }
1082  }
1083  // For each new facet
1084  for ( Index i = 0; i < new_facets.size(); i++ )
1085  { // link the new facet to its neighbors
1086  for ( Index j = i + 1; j < new_facets.size(); j++ )
1087  {
1088  const Index nfi = new_facets[ i ];
1089  const Index nfj = new_facets[ j ];
1090  if ( areFacetsNeighbor( nfi, nfj ) )
1091  makeNeighbors( nfi, nfj );
1092  }
1093  }
1094  // Extracts all outside points from visible facets V
1095  IndexRange outside_pts;
1096  for ( auto&& vf : V ) {
1097  for ( auto&& v : facets[ vf ].outside_set ) {
1098  if ( v != furthest_v ) {
1099  outside_pts.push_back( v );
1100  assignment[ v ] = UNASSIGNED;
1101  }
1102  }
1103  }
1104  // For each new facet F'
1105  for ( Index i = 0; i < new_facets.size(); i++ ) {
1106  Facet& Fp = facets[ new_facets[ i ] ];
1107  Index max_j = outside_pts.size();
1108  for ( Index j = 0; j < max_j; ) {
1109  const Index v = outside_pts[ j ];
1110  if ( above( Fp, points[ v ] ) ) {
1111  Fp.outside_set.push_back( v );
1112  assignment[ v ] = new_facets[ i ];
1113  outside_pts[ j ] = outside_pts.back();
1114  outside_pts.pop_back();
1115  max_j--;
1116  } else j++;
1117  }
1118  if ( debug_level >= 3 ) {
1119  trace.info() << "- New facet " << new_facets[ i ] << " ";
1120  Fp.display( trace.info() );
1121  }
1122  }
1123  // Update processed points
1124  processed_points.push_back( furthest_v );
1125  for ( auto v : outside_pts ) processed_points.push_back( v );
1126 
1127  // Delete the facets in V
1128  for ( auto&& v : V ) {
1129  if ( debug_level >= 2 ) {
1130  trace.info() << "Delete facet " << v << " ";
1131  facets[ v ].display( trace.info() );
1132  }
1133  deleteFacet( v );
1134  }
1135 
1136  // Add new facets to queue
1137  for ( Index i = 0; i < new_facets.size(); i++ )
1138  Q.push( new_facets[ i ] );
1139  if ( debug_level >= 1 ) {
1140  trace.info() << "#facets=" << facets.size()
1141  << " #deleted=" << deleted_facets.size() << std::endl;
1142  }
1143 
1144  // Checks that everything is ok.
1145  if ( debug_level >= 1 ) {
1146  trace.info() << "[CHECK INVARIANT] " << processed_points.size()
1147  << " / " << points.size() << " points processed." << std::endl;
1148  bool okh = checkHull();
1149  if ( ! okh )
1150  trace.error() << "[computeFacet] Invalid convex hull" << std::endl;
1151  bool okf = checkFacets();
1152  if ( ! okf )
1153  trace.error() << "[computeFacet] Invalid facets" << std::endl;
1154  if ( ! ( okh && okf ) ) myStatus = Status::InvalidConvexHull;
1155  }
1156 
1157  return status() == Status::SimplexCompleted;
1158  }
1159 
1161  bool checkFacet( Index f ) const
1162  {
1163  const Facet& F = facets[ f ];
1164  bool ok = F.on_set.size() >= dimension;
1165  for ( auto v : F.on_set )
1166  if ( ! on( F, points[ v ] ) ) {
1167  trace.error() << "[QuickHull::checkFacet( " << f
1168  << ") Invalid 'on' vertex " << v << std::endl;
1169  ok = false;
1170  }
1171  if ( F.neighbors.size() < dimension ) {
1172  trace.error() << "[QuickHull::checkFacet( " << f
1173  << ") Not enough neighbors " << F.neighbors.size() << std::endl;
1174  ok = false;
1175  }
1176  for ( auto nf : F.neighbors )
1177  if ( ! areFacetsNeighbor( f, nf ) ) {
1178  trace.error() << "[QuickHull::checkFacet( " << f
1179  << ") Invalid neighbor " << nf << std::endl;
1180  ok = false;
1181  }
1182  if ( aboveOrOn( F, points[ F.below ] ) ) {
1183  trace.error() << "[QuickHull::checkFacet( " << f
1184  << ") Bad below point " << F.below << std::endl;
1185  ok = false;
1186  }
1187  for ( auto ov : F.outside_set )
1188  if ( ! above( F, points[ ov ] ) ) {
1189  trace.error() << "[QuickHull::checkFacet( " << f
1190  << ") Bad outside point " << ov << std::endl;
1191  ok = false;
1192  }
1193  for ( auto && v : F.on_set ) {
1194  Size n = 0;
1195  for ( auto&& N : facets[ f ].neighbors )
1196  if ( on( facets[ N ], points[ v ] ) ) n += 1;
1197  if ( n < dimension-1 ) {
1198  trace.error() << "[QuickHull::checkFacet( " << f << ") 'on' point " << v
1199  << " is a vertex of " << n << " facets "
1200  << "(should be >= " << dimension-1 << ")" << std::endl;
1201  ok = false;
1202  }
1203  }
1204  return ok;
1205  }
1206 
1209  {
1210  // SLightly faster to postpone deletion of intermediate facets.
1211  const Index f = facets.size();
1212  facets.push_back( Facet() );
1213  return f;
1214  }
1215 
1218  void deleteFacet( Index f )
1219  {
1220  for ( auto n : facets[ f ].neighbors )
1221  facets[ n ].subNeighbor( f );
1222  deleted_facets.insert( f );
1223  facets[ f ].clear();
1224  }
1225 
1229  void makeNeighbors( const Index if1, const Index if2 )
1230  {
1231  facets[ if1 ].addNeighbor( if2 );
1232  facets[ if2 ].addNeighbor( if1 );
1233  }
1234 
1238  void unmakeNeighbors( const Index if1, const Index if2 )
1239  {
1240  facets[ if1 ].subNeighbor( if2 );
1241  facets[ if2 ].subNeighbor( if1 );
1242  }
1243 
1250  Index mergeParallelFacets( const Index if1, const Index if2 )
1251  {
1252  Facet& f1 = facets[ if1 ];
1253  Facet& f2 = facets[ if2 ];
1254  std::copy( f2.outside_set.cbegin(), f2.outside_set.cend(),
1255  std::back_inserter( f1.outside_set ) );
1256  IndexRange merge_idx;
1257  std::set_union( f1.on_set.cbegin(), f1.on_set.cend(),
1258  f2.on_set.cbegin(), f2.on_set.cend(),
1259  std::back_inserter( merge_idx ) );
1260  f1.on_set.swap( merge_idx );
1261  for ( auto && nf2 : f2.neighbors ) {
1262  if ( nf2 == if1 ) continue;
1263  facets[ nf2 ].subNeighbor( if2 );
1264  makeNeighbors( if1, nf2 );
1265  }
1266  return if2;
1267  }
1268 
1273  bool areFacetsParallel( const Index if1, const Index if2 ) const
1274  {
1275  ASSERT( if1 != if2 );
1276  const Facet& f1 = facets[ if1 ];
1277  const Facet& f2 = facets[ if2 ];
1278  if ( kernel.equal( f1.H, f2.H ) ) return true;
1279  // Need to check if one N is a multiple of the other.
1280  for ( auto&& v : f1.on_set )
1281  if ( ! on( f2, points[ v ] ) ) return false;
1282  return true;
1283  }
1284 
1292  bool areFacetsNeighbor( const Index if1, const Index if2 ) const
1293  {
1294  const Facet& f1 = facets[ if1 ];
1295  const Facet& f2 = facets[ if2 ];
1296  Index nb = 0;
1297  for ( Index i1 = 0, i2 = 0; i1 < f1.on_set.size() && i2 < f2.on_set.size(); )
1298  {
1299  if ( f1.on_set[ i1 ] == f2.on_set[ i2 ] ) { nb++; i1++; i2++; }
1300  else if ( f1.on_set[ i1 ] < f2.on_set[ i2 ] ) i1++;
1301  else i2++;
1302  }
1303  return nb >= ( dimension - 1 );
1304  }
1305 
1317  Facet makeFacet( const IndexRange& base, Index below ) const
1318  {
1319  CombinatorialPlaneSimplex simplex;
1320  for ( Size i = 0; i < dimension; i++ ) simplex[ i ] = base[ i ];
1321  auto plane = kernel.compute( points, simplex, below );
1322  return Facet( plane, below );
1323  }
1324 
1327  IndexRange pointsOnRidge( const Ridge& R ) const
1328  {
1329  // Slightly faster to look at points instead at true geometry.
1330  IndexRange result;
1331  const Facet& f1 = facets[ R.first ];
1332  const Facet& f2 = facets[ R.second ];
1333  std::set_intersection( f1.on_set.cbegin(), f1.on_set.cend(),
1334  f2.on_set.cbegin(), f2.on_set.cend(),
1335  std::back_inserter( result ) );
1336  return result;
1337  }
1338 
1349  {
1350  const Facet& F = facets[ f ];
1351  IndexRange result = F.on_set;
1352  // Sort a facet such that its points, taken in order, have
1353  // always the same orientation of the facet. More precisely,
1354  // facets span a `dimension-1` vector space. There are thus
1355  // dimension-2 fixed points, and the last ones (at least two)
1356  // may be reordered.
1358  for ( Dimension k = 0; k < dimension-2; k++ )
1359  splx[ k ] = result[ k ];
1360  std::sort( result.begin()+dimension-2, result.end(),
1361  [&] ( Index i, Index j )
1362  {
1363  splx[ dimension-2 ] = i;
1364  splx[ dimension-1 ] = j;
1365  const auto H = kernel.compute( points, splx );
1366  const auto orient = kernel.dot( F.H, H );
1367  return orient > 0;
1368  } );
1369  return result;
1370  }
1371 
1372 
1378  {
1379  auto & on_set = facets[ f ].on_set;
1380  for ( Index i = 0; i < on_set.size(); )
1381  {
1382  Index v = on_set[ i ];
1383  Size n = 0;
1384  for ( auto&& N : facets[ f ].neighbors )
1385  if ( on( facets[ N ], points[ v ] ) ) n += 1;
1386  if ( n >= dimension-1 ) i++;
1387  else {
1388  on_set[ i ] = on_set.back();
1389  on_set.pop_back();
1390  }
1391  }
1392  std::sort( on_set.begin(), on_set.end() );
1393  }
1394 
1398  {
1399  const Size nb = points.size();
1400  if ( nb < dimension + 1 ) return IndexRange();
1401  IndexRange best = pickIntegers( dimension + 1, nb );
1403  for ( Index j = 0; j < dimension; ++j ) splx[ j ] = best[ j ];
1404  const auto first_H = kernel.compute( points, splx, best.back() );
1405  auto best_volume = kernel.volume ( first_H, points[ best.back() ] );
1406  const Size nbtries = std::min( (Size) 10, 1 + nb / 10 );
1407  const Size max_nbtries = std::max( (Size) 10, 2 * nb );
1408  for ( Size i = 0; i < max_nbtries; i++ )
1409  {
1410  IndexRange tmp = pickIntegers( dimension + 1, nb );
1411  for ( Index j = 0; j < dimension; ++j ) splx[ j ] = tmp[ j ];
1412  const auto tmp_H = kernel.compute( points, splx, tmp.back() );
1413  const auto tmp_volume = kernel.volume ( tmp_H, points[ tmp.back() ] );
1414  if ( best_volume < tmp_volume ) {
1415  if ( debug_level >= 1 ) {
1416  trace.info() << "(" << i << ")"
1417  << " new_volume = " << tmp_volume
1418  << " > " << best_volume << std::endl;
1419  }
1420  best = tmp;
1421  best_volume = tmp_volume;
1422  }
1423  if ( i >= nbtries && best_volume > 0 ) return best;
1424  }
1425  return IndexRange();
1426  }
1427 
1431  static IndexRange pickIntegers( const Size d, const Size n )
1432  {
1433  IndexRange result( d );
1434  bool distinct = false;
1435  while ( ! distinct )
1436  {
1437  distinct = true;
1438  for ( Index i = 0; i < d; i++ ) result[ i ] = rand() % n;
1439  std::sort( result.begin(), result.end() );
1440  for ( Index i = 1; distinct && i < d; i++ )
1441  distinct = result[ i-1 ] != result[ i ];
1442  }
1443  return result;
1444  }
1445 
1454  bool computeSimplexConfiguration( const IndexRange& full_simplex )
1455  {
1456  assignment = std::vector< Index >( points.size(), UNASSIGNED );
1457  facets.resize( dimension + 1 );
1458  deleted_facets.clear();
1459  for ( Index j = 0; j < full_simplex.size(); ++j )
1460  {
1461  IndexRange lsimplex( dimension );
1462  IndexRange isimplex( dimension );
1463  Index s = 0;
1464  for ( Index i = 0; i <= dimension; i++ )
1465  if ( i != j ) {
1466  lsimplex[ s ] = i;
1467  isimplex[ s ] = full_simplex[ i ];
1468  s++;
1469  }
1470  facets[ j ] = makeFacet( isimplex, full_simplex[ j ] );
1471  facets[ j ].neighbors = lsimplex;
1472  for ( auto&& v : isimplex ) facets[ j ].on_set.push_back( v );
1473  std::sort( facets[ j ].on_set.begin(), facets[ j ].on_set.end() );
1474  }
1475  // List of unassigned vertices
1476  for ( Index fi = 0; fi < facets.size(); ++fi ) {
1477  Facet& f = facets[ fi ];
1478  for ( Index v = 0; v < points.size(); v++ )
1479  {
1480  if ( assignment[ v ] == UNASSIGNED && above( f, points[ v ] ) ) {
1481  f.outside_set.push_back( v );
1482  assignment[ v ] = fi;
1483  }
1484  }
1485  }
1486  for ( Index v = 0; v < points.size(); v++ )
1487  if ( assignment[ v ] == UNASSIGNED )
1488  processed_points.push_back( v );
1489 
1490  // Display some information
1491  if ( debug_level >= 2 ) {
1492  for ( auto&& f : facets ) f.display( trace.info() );
1493  }
1495  if ( debug_level >= 1 ) {
1496  trace.info() << "[CHECK INVARIANT] " << processed_points.size()
1497  << " / " << points.size() << " points processed." << std::endl;
1498  bool okh = checkHull();
1499  if ( ! okh )
1500  trace.error() << "[computeInitialSimplex] Invalid convex hull" << std::endl;
1501  bool okf = checkFacets();
1502  if ( ! okf )
1503  trace.error() << "[computeInitialSimplex] Invalid facets" << std::endl;
1504  if ( ! ( okh && okf ) ) myStatus = Status::InvalidConvexHull;
1505  }
1506  return status() == Status::SimplexCompleted;
1507  }
1508 
1510  // ----------------------- Interface --------------------------------------
1511  public:
1514 
1519  void selfDisplay ( std::ostream & out ) const
1520  {
1521  out << "[QuickHull " << dimension << "D"
1522  // << " status=" << status()
1523  << " #P=" << nbPoints();
1525  out << " #F=" << nbFacets();
1527  out << " #V=" << nbVertices();
1528  out << "]";
1530  && nbFacets() == 24 ) {
1531  for ( auto f : facets ) f.display( out );
1532  for ( auto v : v2p ) out << points[ v2p[ v ] ] << std::endl;
1533  }
1534  }
1535 
1540  bool isValid() const
1541  {
1542  return status() >= Status::Uninitialized
1543  && status() <= Status::AllCompleted;
1544  }
1546 
1547  };
1548 
1556  template < typename TKernel >
1557  std::ostream&
1558  operator<< ( std::ostream & out,
1559  const QuickHull< TKernel > & object )
1560  {
1561  object.selfDisplay( out );
1562  return out;
1563  }
1564 
1565 } // namespace DGtal
1566 
1568 // Includes inline functions.
1569 // //
1571 
1572 #endif // !defined QuickHull_h
1573 
1574 #undef QuickHull_RECURSES
1575 #endif // else defined(QuickHull_RECURSES)
1576 
void tic()
Starts timer.
std::ostream & error()
std::ostream & info()
std::ostream & warning()
DGtal is the top-level namespace which contains all DGtal functions and types.
void swap(UnorderedSetByBlock< Key, TSplitter, Hash, KeyEqual, UnorderedMapAllocator > &s1, UnorderedSetByBlock< Key, TSplitter, Hash, KeyEqual, UnorderedMapAllocator > &s2) noexcept
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
void subNeighbor(Index n)
Definition: QuickHull.h:206
void addPointOn(Index p)
Definition: QuickHull.h:183
Facet(const HalfSpace &aH, Index b)
Definition: QuickHull.h:172
void display(std::ostream &out) const
Definition: QuickHull.h:188
Facet(Facet &&)=default
void addNeighbor(Index n)
Definition: QuickHull.h:201
Facet & operator=(Facet &&)=default
Index below
index of point that is below this facet
Definition: QuickHull.h:165
Facet(const Facet &)=default
HalfSpace H
the facet geometry
Definition: QuickHull.h:161
Facet & operator=(const Facet &)=default
IndexRange outside_set
outside set, i.e. points above this facet
Definition: QuickHull.h:163
void swap(Facet &other)
Definition: QuickHull.h:214
IndexRange neighbors
neighbor facets
Definition: QuickHull.h:162
IndexRange on_set
on set, i.e. points on this facet, sorted
Definition: QuickHull.h:164
Size variableMemory() const
Definition: QuickHull.h:224
Aim: Implements the quickhull algorithm by Barber et al. , a famous arbitrary dimensional convex hull...
Definition: QuickHull.h:140
bool computeInitialSimplex()
Definition: QuickHull.h:504
bool computeFacets()
Definition: QuickHull.h:523
void unmakeNeighbors(const Index if1, const Index if2)
Definition: QuickHull.h:1238
bool setInput(const std::vector< InputPoint > &input_points, bool remove_duplicates=true)
Definition: QuickHull.h:382
Kernel::CoordinateVector Vector
Definition: QuickHull.h:143
bool areFacetsNeighbor(const Index if1, const Index if2) const
Definition: QuickHull.h:1292
QuickHull(const Kernel &K=Kernel(), int dbg=0)
Definition: QuickHull.h:259
Kernel kernel
Kernel that is duplicated for computing facet geometries.
Definition: QuickHull.h:798
int debug_level
debug_level from 0:no to 2
Definition: QuickHull.h:800
bool setInitialSimplex(const IndexRange &full_splx)
Definition: QuickHull.h:410
Kernel::InternalScalar InternalScalar
Definition: QuickHull.h:145
std::vector< Index > IndexRange
Definition: QuickHull.h:149
Size nb_infinite_facets
Number of infinite facets (!= 0 only for specific kernels)
Definition: QuickHull.h:824
bool on(const Facet &F, const Point &p) const
Definition: QuickHull.h:868
void filterVerticesOnFacet(const Index f)
Definition: QuickHull.h:1377
IndexRange processed_points
Points already processed (and within the convex hull).
Definition: QuickHull.h:810
Index newFacet()
Definition: QuickHull.h:1208
bool checkFacet(Index f) const
Definition: QuickHull.h:1161
BOOST_STATIC_ASSERT((Point::dimension==Vector::dimension))
Facet makeFacet(const IndexRange &base, Index below) const
Definition: QuickHull.h:1317
Size nb_finite_facets
Number of finite facets.
Definition: QuickHull.h:822
IndexRange pointsOnRidge(const Ridge &R) const
Definition: QuickHull.h:1327
IndexRange v2p
vertex index -> point index
Definition: QuickHull.h:820
Kernel::CombinatorialPlaneSimplex CombinatorialPlaneSimplex
Definition: QuickHull.h:151
bool computeSimplexConfiguration(const IndexRange &full_simplex)
Definition: QuickHull.h:1454
void cleanInfiniteFacets()
Definition: QuickHull.h:901
std::vector< Facet > facets
the current set of facets.
Definition: QuickHull.h:814
bool getFacetVertices(std::vector< IndexRange > &facet_vertices) const
Definition: QuickHull.h:666
std::size_t Index
Definition: QuickHull.h:146
TKernel Kernel
Definition: QuickHull.h:141
std::pair< Index, Index > Ridge
Definition: QuickHull.h:236
bool checkFacets()
Definition: QuickHull.h:743
Kernel::CoordinatePoint Point
Definition: QuickHull.h:142
std::vector< Point > points
the set of points, indexed as in the array.
Definition: QuickHull.h:802
IndexRange input2comp
Definition: QuickHull.h:805
std::set< Index > deleted_facets
set of deleted facets
Definition: QuickHull.h:816
void clear()
Clears the object.
Definition: QuickHull.h:270
Size nbFiniteFacets() const
Definition: QuickHull.h:346
bool areFacetsParallel(const Index if1, const Index if2) const
Definition: QuickHull.h:1273
bool computeVertices()
Definition: QuickHull.h:553
bool getVertex2Point(IndexRange &vertex_to_point)
Definition: QuickHull.h:632
IndexRange pickInitialSimplex() const
Definition: QuickHull.h:1397
bool isValid() const
Definition: QuickHull.h:1540
IndexRange orientedFacetPoints(Index f) const
Definition: QuickHull.h:1348
std::vector< Size > facet_counter
Counts the number of facets with a given number of vertices.
Definition: QuickHull.h:828
Status status() const
Definition: QuickHull.h:266
bool processFacet(std::queue< Index > &Q)
Definition: QuickHull.h:969
static IndexRange pickIntegers(const Size d, const Size n)
Definition: QuickHull.h:1431
Size nbVertices() const
Definition: QuickHull.h:337
std::size_t Size
Definition: QuickHull.h:147
void deleteFacet(Index f)
Definition: QuickHull.h:1218
bool getVertexPositions(std::vector< OutputPoint > &vertex_positions)
Definition: QuickHull.h:612
bool checkHull()
Definition: QuickHull.h:762
std::vector< Index > assignment
assignment of points to facets
Definition: QuickHull.h:812
static const Size dimension
Definition: QuickHull.h:152
Status
Represents the status of a QuickHull object.
Definition: QuickHull.h:239
@ InvalidRidge
Error: some ridge is not consistent (probably integer overflow).
@ InvalidConvexHull
Error: the convex hull does not contain all its points (probably integer overflow).
@ SimplexCompleted
An initial full-dimensional simplex has been found. QuickHull core algorithm can start.
@ FacetsCompleted
All facets of the convex hull are identified.
@ NotFullDimensional
Error: the initial simplex is not full dimensional.
@ VerticesCompleted
All vertices of the convex hull are determined.
@ InputInitialized
A range of input points has been given to QuickHull.
@ AllCompleted
Same as VerticesCompleted.
@ Uninitialized
QuickHull is empty and has just been instantiated.
Size nbInfiniteFacets() const
Definition: QuickHull.h:355
InternalScalar height(const Facet &F, const Point &p) const
Definition: QuickHull.h:850
bool getFacetHalfSpaces(std::vector< HalfSpace > &facet_halfspaces)
Definition: QuickHull.h:690
Index mergeParallelFacets(const Index if1, const Index if2)
Definition: QuickHull.h:1250
void makeNeighbors(const Index if1, const Index if2)
Definition: QuickHull.h:1229
void selfDisplay(std::ostream &out) const
Definition: QuickHull.h:1519
bool aboveOrOn(const Facet &F, const Point &p) const
Definition: QuickHull.h:862
IndexRange p2v
point index -> vertex index (or UNASSIGNED)
Definition: QuickHull.h:818
Size memory() const
Definition: QuickHull.h:287
Size nbFacets() const
Definition: QuickHull.h:328
bool getPoint2Vertex(IndexRange &point_to_vertex)
Definition: QuickHull.h:651
void renumberInfiniteFacets()
Definition: QuickHull.h:926
bool computeConvexHull(Status target=Status::VerticesCompleted)
Definition: QuickHull.h:460
Kernel::HalfSpace HalfSpace
Definition: QuickHull.h:150
Kernel::CoordinateScalar Scalar
Definition: QuickHull.h:144
IndexRange comp2input
Definition: QuickHull.h:808
std::vector< double > timings
Timings of the different phases: 0: init, 1: facets, 2: vertices.
Definition: QuickHull.h:826
bool above(const Facet &F, const Point &p) const
Definition: QuickHull.h:856
Size nbPoints() const
Definition: QuickHull.h:323
void cleanFacets()
Definition: QuickHull.h:873
int max(int a, int b)
KSpace K