1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkEllipsoidInteriorExteriorSpatialFunctionTest.cxx,v $
Language: C++
Date: $Date: 2005-08-09 18:01:44 $
Version: $Revision: 1.19 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkEllipsoidInteriorExteriorSpatialFunction.h"
#include "vnl/vnl_matrix.h"
int itkEllipsoidInteriorExteriorSpatialFunctionTest(int, char* [] )
{
std::cout << "itkEllipsoidInteriorExteriorSpatialFunction test start" << std::endl;
// Test will create an ellipsoid (3-dimensional)
const unsigned int dimension = 3;
// Ellipsoid spatial function typedef
typedef itk::EllipsoidInteriorExteriorSpatialFunction<3> TEllipsoidFunctionType;
// Point position typedef
typedef TEllipsoidFunctionType::InputType TEllipsoidFunctionVectorType;
// Create an ellipsoid spatial function for the source image
TEllipsoidFunctionType::Pointer spatialFunc = TEllipsoidFunctionType::New();
// Define and set the axes lengths for the ellipsoid
TEllipsoidFunctionVectorType axes;
axes[0] = 40;
axes[1] = 30;
axes[2] = 20;
spatialFunc->SetAxes(axes);
// Define function doitkEllipsoidInteriorExteriorSpatialFunctionTest, which encapsulates ellipsoid.
int xExtent = 50;
int yExtent = 50;
int zExtent = 50;
// Define and set the center of the ellipsoid in the center of
// the function doitkEllipsoidInteriorExteriorSpatialFunctionTest
TEllipsoidFunctionVectorType center;
center[0] = xExtent/2;
center[1] = yExtent/2;
center[2] = zExtent/2;
spatialFunc->SetCenter(center);
// Define the orientations of the ellipsoid axes
// (0,1,0) corresponds to the axes of length axes[0]
// (1,0,0) corresponds to the axes of length axes[1]
// (0,0,1) corresponds to the axes of length axes[2]
double data[] = {0, 1, 0, 1, 0, 0, 0, 0, 1};
vnl_matrix<double> orientations (data, 3, 3);
// Set the orientations of the ellipsoids
spatialFunc->SetOrientations(orientations);
// Evaluate all points in the spatial function and count the number of
// pixels that are inside the sphere.
double testPosition[dimension]; // position of a pixel in the function doitkEllipsoidInteriorExteriorSpatialFunctionTest
bool functionValue; // Value of pixel at a given position
int interiorPixelCounter = 0; // Count pixels inside ellipsoid
for(int x = 0; x < xExtent; x++)
{
for(int y = 0; y < yExtent; y++)
{
for(int z =0; z < zExtent; z++)
{
testPosition[0] = x;
testPosition[1] = y;
testPosition[2] = z;
functionValue = spatialFunc->Evaluate(testPosition);
if(functionValue == 1)
interiorPixelCounter ++;
}
}
}
// Evaluate the center of the ellipsoid, which is inside the ellipsoid and
// should equal 1.
testPosition[0] = center[0];
testPosition[1] = center[1];
testPosition[2] = center[2];
functionValue = spatialFunc->Evaluate(testPosition);
// Volume of ellipsoid using V=(4/3)*pi*(a/2)*(b/2)*(c/2)
double volume = 4.18879013333*(axes[0]/2)*(axes[1]/2)*(axes[2]/2);
// Percent difference in volume measurement and calculation
double volumeError = (fabs(volume - interiorPixelCounter)/volume)*100;
std::cout << spatialFunc;
// 5% error was randomly chosen as a successful ellipsoid fill.
// This should actually be some function of the image/ellipsoid size.
if(volumeError <= 5 || functionValue == 1)
{
// With testing settings, results should yield:
// calculated ellipsoid volume = 12566.4 pixels
// measured ellipsoid volume = 12428 pixels
// volume error = 1.10907%
// function value = 1
std::cout << "calculated ellipsoid volume = " << volume << std::endl
<< "measured ellipsoid volume = " << interiorPixelCounter << std::endl
<< "volume error = " << volumeError << "%" << std::endl
<< "function value = " << functionValue << std::endl
<< "center location = (" << spatialFunc->GetCenter()[0] << ", " << spatialFunc->GetCenter()[0]
<< ", " << spatialFunc->GetCenter()[2] << ")" << std::endl
<< "major axis length = " << spatialFunc->GetAxes()[0] << " minor axis 1 length = "
<< spatialFunc->GetAxes()[1] << " minor axis 2 length = " << spatialFunc->GetAxes()[2] << std::endl
<< "itkEllipsoidSpatialFunction ended successfully!" << std::endl;
return EXIT_SUCCESS;
}
//Default behavior is to fail
std::cerr << "calculated ellipsoid volume = " << volume << std::endl
<< "measured ellipsoid volume = " << interiorPixelCounter << std::endl
<< "volume error = " << volumeError << "%" << std::endl
<< "function value = " << functionValue << std::endl
<< "center location = (" << spatialFunc->GetCenter()[0] << ", " << spatialFunc->GetCenter()[0]
<< ", " << spatialFunc->GetCenter()[2] << ")" << std::endl
<< "major axis length = " << spatialFunc->GetAxes()[0] << " minor axis 1 length = "
<< spatialFunc->GetAxes()[1] << " minor axis 2 length = " << spatialFunc->GetAxes()[2] << std::endl
<< "itkEllipsoidSpatialFunction failed :(" << std::endl;
return EXIT_FAILURE;
}
|
