ITK  4.9.0
Insight Segmentation and Registration Toolkit
WikiExamples/Morphology/FlatStructuringElementRadiusIsParametric.cxx
// Helper function
template< class SEType>
bool ComputeAreaError(SEType k, unsigned int thickness = 0);
int main(int, char *[])
{
int scalarRadius = 5;
int scalarThickness = 2;
bool radiusIsParametric = true;
SE2Type::RadiusType r2;
r2.Fill( scalarRadius );
SE2Type k2;
std::cout << "2D ball of radius " << scalarRadius
<< " with radiusIsParametric mode off:" << std::endl;
k2 = SE2Type::Ball( r2 );
ComputeAreaError(k2);
// Test the radiusIsParametric mode.
std::cout << "2D ball of radius " << scalarRadius
<< " with radiusIsParametric mode on:" << std::endl;
k2 = SE2Type::Ball(r2, radiusIsParametric);
ComputeAreaError(k2);
std::cout << "2D annulus of radius " << scalarRadius
<< " and thickness " << scalarThickness
<< " with radiusIsParametric mode off:" << std::endl;
k2 = SE2Type::Annulus(r2,scalarThickness,false);
ComputeAreaError(k2,scalarThickness);
// Test the radiusIsParametric mode.
std::cout << "2D annulus of radius " << scalarRadius
<< " and thickness " << scalarThickness
<< " with radiusIsParametric mode on:" << std::endl;
k2 = SE2Type::Annulus(r2,scalarThickness,false,radiusIsParametric);
ComputeAreaError(k2,scalarThickness);
SE3Type::RadiusType r3;
r3.Fill( scalarRadius );
SE3Type k3;
std::cout << "3D ball of radius " << scalarRadius
<< " with radiusIsParametric mode off:" << std::endl;
k3 = SE3Type::Ball( r3 );
ComputeAreaError(k3);
// Test the radiusIsParametric mode.
std::cout << "3D ball of radius " << scalarRadius
<< " with radiusIsParametric mode on:" << std::endl;
k3 = SE3Type::Ball(r3, radiusIsParametric);
ComputeAreaError(k3);
std::cout << "3D annulus of radius " << scalarRadius
<< " and thickness " << scalarThickness
<< " with radiusIsParametric mode off:" << std::endl;
k3 = SE3Type::Annulus(r3,scalarThickness,false);
ComputeAreaError(k3,scalarThickness);
// Test the radiusIsParametric mode.
std::cout << "3D annulus of radius " << scalarRadius
<< " and thickness " << scalarThickness
<< " with radiusIsParametric mode on:" << std::endl;
k3 = SE3Type::Annulus(r3,scalarThickness,false,radiusIsParametric);
ComputeAreaError(k3,scalarThickness);
SE4Type::RadiusType r4;
r4.Fill( scalarRadius );
SE4Type k4;
std::cout << "4D ball of radius " << scalarRadius
<< " with radiusIsParametric mode off:" << std::endl;
k4 = SE4Type::Ball( r4 );
ComputeAreaError(k4);
// Test the radiusIsParametric mode.
std::cout << "4D ball of radius " << scalarRadius
<< " with radiusIsParametric mode on:" << std::endl;
k4 = SE4Type::Ball(r4, radiusIsParametric);
ComputeAreaError(k4);
std::cout << "4D annulus of radius " << scalarRadius
<< " and thickness " << scalarThickness
<< " with radiusIsParametric mode off:" << std::endl;
k4 = SE4Type::Annulus(r4,scalarThickness,false);
ComputeAreaError(k4,scalarThickness);
// Test the radiusIsParametric mode.
std::cout << "4D annulus of radius " << scalarRadius
<< " and thickness " << scalarThickness
<< " with radiusIsParametric mode on:" << std::endl;
k4 = SE4Type::Annulus(r4,scalarThickness,false,radiusIsParametric);
ComputeAreaError(k4,scalarThickness);
return EXIT_SUCCESS;
}
template< class SEType >
bool ComputeAreaError(SEType k, unsigned int thickness)
{
float expectedOuterForegroundArea = 1;
float expectedInnerForegroundArea;
if( thickness == 0 )
{
// Circle/Ellipse has no inner area to subract.
expectedInnerForegroundArea = 0;
}
else
{
// Annulus does have inner area to subract.
expectedInnerForegroundArea = 1;
}
if( SEType::NeighborhoodDimension == 2)
{
expectedOuterForegroundArea *= vnl_math::pi;
expectedInnerForegroundArea *= vnl_math::pi;
}
else if( SEType::NeighborhoodDimension == 3 )
{
expectedOuterForegroundArea *= 4.0/3.0 * vnl_math::pi;
expectedInnerForegroundArea *= 4.0/3.0 * vnl_math::pi;
}
else if ( SEType::NeighborhoodDimension == 4 )
{
expectedOuterForegroundArea *= 0.5 * vnl_math::pi * vnl_math::pi;
expectedInnerForegroundArea *= 0.5 * vnl_math::pi * vnl_math::pi;
}
else
{
return EXIT_FAILURE;
}
for( unsigned int i = 0; i < SEType::NeighborhoodDimension; i++ )
{
expectedOuterForegroundArea *= k.GetRadius()[i];
expectedInnerForegroundArea *= (k.GetRadius()[i] - thickness);
}
float expectedForegroundArea = expectedOuterForegroundArea - expectedInnerForegroundArea;
// Show the neighborhood if it is 2D.
typename SEType::Iterator SEIt;
if( SEType::NeighborhoodDimension == 2 )
{
for( SEIt = k.Begin(); SEIt != k.End(); ++SEIt )
{
std::cout << *SEIt << "\t";
if( (SEIt - k.Begin()+1) % k.GetSize()[0] == 0 )
{
std::cout << std::endl;
}
}
}
// Compute the area/volume.
float computedForegroundArea = 0;
for( SEIt = k.Begin(); SEIt != k.End(); ++SEIt )
{
if( *SEIt )
{
computedForegroundArea++;
}
}
std::cout << "Expected foreground area: " << expectedForegroundArea << std::endl;
std::cout << "Computed foreground area: " << computedForegroundArea << std::endl;
std::cout << "Foreground area error: "
<< 100 * vnl_math_abs(expectedForegroundArea-computedForegroundArea)/expectedForegroundArea
<< "%" << "\n\n";
return EXIT_FAILURE;
}