Generate Structuring Elements With Accurate Area

Synopsis

Generate structuring elements with accurate area.

Results

Output:

2D ball of radius 5 with radiusIsParametric mode off:
0       0       0       1       1       1       1       1       0       0       0
0       0       1       1       1       1       1       1       1       0       0
0       1       1       1       1       1       1       1       1       1       0
1       1       1       1       1       1       1       1       1       1       1
1       1       1       1       1       1       1       1       1       1       1
1       1       1       1       1       1       1       1       1       1       1
1       1       1       1   1   1       1       1       1       1       1
1       1       1       1       1       1       1       1       1       1       1
0       1       1       1       1       1       1       1       1       1       0
0       0       1       1   1   1       1       1       1       0       0
0       0       0       1       1       1       1       1       0       0       0
Expected foreground area: 78.5398
Computed foreground area: 97
Foreground area error: 23.5042%

2D ball of radius 5 with radiusIsParametric mode on:
0       0       0       0       0       1       0       0       0       0       0
0       0       1       1       1       1       1       1       1       0       0
0       1       1       1       1       1       1       1       1       1       0
0       1       1       1       1       1       1       1   1   1       0
0       1       1       1       1       1       1       1       1       1       0
1       1   1   1       1       1       1       1       1       1       1
0       1       1       1       1       1       1       1       1       1       0
0       1       1       1       1       1       1       1       1       1       0
0       1       1       1       1       1       1       1       1       1       0
0       0       1       1       1       1       1       1       1       0       0
0       0       0       0       0       1       0       0       0       0       0
Expected foreground area: 78.5398
Computed foreground area: 81
Foreground area error: 3.1324%

2D annulus of radius 5 and thickness 2 with radiusIsParametric mode off:
0       0       0       1       1       1       1       1       0       0       0
0       0       1       1       1       1       1       1       1       0       0
0   1   1       1   0   0       0       1       1       1       0
1       1       1       0       0       0       0       0       1       1       1
1       1       0       0       0       0       0       0       0       1       1
1       1       0       0       0       0       0       0       0       1       1
1       1       0       0   0   0       0       0       0       1       1
1       1       1       0       0       0       0       0       1       1       1
0       1       1       1       0       0       0       1       1       1       0
0       0       1       1       1       1       1       1       1       0       0
0   0   0       1       1       1       1       1       0       0       0
Expected foreground area: 50.2655
Computed foreground area: 60
Foreground area error: 19.3662%

2D annulus of radius 5 and thickness 2 with radiusIsParametric mode on:
0   0   0       0       0       1       0       0       0       0       0
0       0       1       1       1   1   1       1       1       0       0
0       1       1       1       1       0       1       1       1       1       0
0       1       1   0   0       0       0       0       1   1   0
0       1       1       0       0       0       0       0       1       1       0
1       1       0       0       0       0       0       0   0   1       1
0       1       1       0       0       0       0       0       1       1       0
0       1       1       0       0       0       0       0       1       1       0
0       1       1       1       1       0       1       1       1       1       0
0       0       1       1       1       1       1       1       1       0       0
0       0       0       0       0       1       0       0       0       0       0
Expected foreground area: 50.2655
Computed foreground area: 52
Foreground area error: 3.45071%

3D ball of radius 5 with radiusIsParametric mode off:
Expected foreground area: 523.599
Computed foreground area: 739
Foreground area error: 41.1386%

3D ball of radius 5 with radiusIsParametric mode on:
Expected foreground area: 523.599
Computed foreground area: 515
Foreground area error: 1.64224%

3D annulus of radius 5 and thickness 2 with radiusIsParametric mode off:
Expected foreground area: 410.501
Computed foreground area: 560
Foreground area error: 36.4185%

3D annulus of radius 5 and thickness 2 with radiusIsParametric mode on:
Expected foreground area: 410.501
Computed foreground area: 392
Foreground area error: 4.50703%

4D ball of radius 5 with radiusIsParametric mode off:
Expected foreground area: 3084.25
Computed foreground area: 4785
Foreground area error: 55.143%

4D ball of radius 5 with radiusIsParametric mode on:
Expected foreground area: 3084.25
Computed foreground area: 2929
Foreground area error: 5.03368%

4D annulus of radius 5 and thickness 2 with radiusIsParametric mode off:
Expected foreground area: 2684.53
Computed foreground area: 4024
Foreground area error: 49.8957%

4D annulus of radius 5 and thickness 2 with radiusIsParametric mode on:
Expected foreground area: 2684.53
Computed foreground area: 2504
Foreground area error: 6.72491%

Code

C++

#include "itkFlatStructuringElement.h"

// 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;

  using SE2Type = itk::FlatStructuringElement<2>;
  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);

  using SE3Type = itk::FlatStructuringElement<3>;
  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);

  using SE4Type = itk::FlatStructuringElement<4>;
  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 *= itk::Math::pi;
    expectedInnerForegroundArea *= itk::Math::pi;
  }
  else if (SEType::NeighborhoodDimension == 3)
  {
    expectedOuterForegroundArea *= 4.0 / 3.0 * itk::Math::pi;
    expectedInnerForegroundArea *= 4.0 / 3.0 * itk::Math::pi;
  }
  else if (SEType::NeighborhoodDimension == 4)
  {
    expectedOuterForegroundArea *= 0.5 * itk::Math::pi * itk::Math::pi;
    expectedInnerForegroundArea *= 0.5 * itk::Math::pi * itk::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 * std::abs(expectedForegroundArea - computedForegroundArea) / expectedForegroundArea << "%"
            << "\n\n";

  return EXIT_FAILURE;
}

Classes demonstrated