plugins/Experiments/SkeletonHausdorffDistance.cxx

   1 #include <Experiments/SkeletonHausdorffDistance.h>
   2 #include <cpInstances/DataObjects/Image.h>
   3 #include <cpInstances/DataObjects/Mesh.h>
   4 #include <cpInstances/DataObjects/Skeleton.h>
   5 #include <cpExtensions/DataStructures/Skeleton.h>
   6
   7 #include <vtkImageData.h>
   8 #include <vtkPolyData.h>
   9
  10 // -------------------------------------------------------------------------
  11 fpaPluginsExperiments::SkeletonHausdorffDistance::
  12 SkeletonHausdorffDistance( )
  13   : Superclass( )
  14 {
  15   this->_ConfigureInput< cpInstances::DataObjects::Image >( "DistanceMap", true, false );
  16   this->_ConfigureInput< cpInstances::DataObjects::Mesh >( "Seeds", true, false );
  17   this->_ConfigureInput< cpInstances::DataObjects::Skeleton >( "Skeleton1", true, false );
  18   this->_ConfigureInput< cpInstances::DataObjects::Skeleton >( "Skeleton2", true, false );
  19 }
  20
  21 // -------------------------------------------------------------------------
  22 fpaPluginsExperiments::SkeletonHausdorffDistance::
  23 ~SkeletonHausdorffDistance( )
  24 {
  25 }
  26
  27 // -------------------------------------------------------------------------
  28 void fpaPluginsExperiments::SkeletonHausdorffDistance::
  29 _GenerateData( )
  30 {
  31   typedef cpExtensions::DataStructures::Skeleton< 3 > _TSkeleton;
  32
  33   auto dmap = this->GetInputData< vtkImageData >( "DistanceMap" );
  34   auto seeds = this->GetInputData< vtkPolyData >( "Seeds" );
  35   auto sk1 = this->GetInputData< _TSkeleton >( "Skeleton1" );
  36   auto sk2 = this->GetInputData< _TSkeleton >( "Skeleton2" );
  37
  38   double buf[ 3 ], pcoords[ 3 ];
  39   seeds->GetPoint( 0, buf );
  40   int ijk[ 3 ];
  41   dmap->ComputeStructuredCoordinates( buf, ijk, pcoords );
  42   double radius =
  43     dmap->GetScalarComponentAsDouble( ijk[ 0 ], ijk[ 1 ], ijk[ 2 ], 0 );
  44
  45   double d1 = this->_Distance( sk1, sk2, NULL, 0 );
  46   double d2 = this->_Distance( sk1, sk2, buf, radius );
  47
  48   std::cout
  49     << std::endl
  50     << "-------------------------------" << std::endl
  51     << "D1    : " << d1 << std::endl
  52     << "D2    : " << d2 << std::endl
  53     << "Radius: " << radius << std::endl
  54     << "Seed  : " << buf[ 0 ] << " " << buf[ 1 ] << " " << buf[ 2 ] << std::endl
  55     << "-------------------------------" << std::endl;
  56 }
  57
  58 // -------------------------------------------------------------------------
  59 template< class _TSkeleton >
  60 double fpaPluginsExperiments::SkeletonHausdorffDistance::
  61 _Distance( _TSkeleton* sk1, _TSkeleton* sk2, double* center, double radius )
  62 {
  63   auto lst1 = this->_PointList( sk1, center, radius );
  64   auto lst2 = this->_PointList( sk2, center, radius );
  65   double dist = -std::numeric_limits< double >::max( );
  66   typename _TSkeleton::TPath::TPoint point;
  67   for( auto p1 : lst1 )
  68   {
  69     double ldist = std::numeric_limits< double >::max( );
  70     typename _TSkeleton::TPath::TPoint lpoint;
  71     for( auto p2 : lst2 )
  72     {
  73       double d = p1.EuclideanDistanceTo( p2 );
  74       if( d < ldist )
  75       {
  76         ldist = d;
  77         lpoint = p2;
  78
  79       } // fi
  80
  81     } // rof
  82     if( ldist > dist )
  83     {
  84       dist = ldist;
  85       point = lpoint;
  86
  87     } // fi
  88
  89   } // rof
  90   return( dist );
  91 }
  92
  93 // -------------------------------------------------------------------------
  94 template< class _TSkeleton >
  95 std::vector< typename _TSkeleton::TPath::TPoint >
  96 fpaPluginsExperiments::SkeletonHausdorffDistance::
  97 _PointList( _TSkeleton* sk, double* center, double radius )
  98 {
  99   typename _TSkeleton::TPath::TPoint p_center;
 100   if( center != NULL )
 101     for( unsigned int d = 0; d < _TSkeleton::Dimension; ++d )
 102       p_center[ d ] = center[ d ];
 103
 104   std::vector< typename _TSkeleton::TPath::TPoint > lst;
 105   auto mIt = sk->BeginEdgesRows( );
 106   for( ; mIt != sk->EndEdgesRows( ); ++mIt )
 107   {
 108     auto rIt = mIt->second.begin( );
 109     for( ; rIt != mIt->second.end( ); ++rIt )
 110     {
 111       auto eIt = rIt->second.begin( );
 112       for( ; eIt != rIt->second.end( ); ++eIt )
 113       {
 114         auto path = *eIt;
 115         for( unsigned int i = 0; i < path->GetSize( ); ++i )
 116         {
 117           auto p = path->GetPoint( i );
 118           if( center != NULL )
 119           {
 120             if( p_center.EuclideanDistanceTo( p ) > radius )
 121               lst.push_back( p );
 122           }
 123           else
 124             lst.push_back( p );
 125
 126         } // rof
 127
 128       } // rof
 129
 130     } // rof
 131
 132   } // rof
 133   return( lst );
 134 }
 135
 136 // eof - $RCSfile$