[ITK-users] Linear case of itk::EllipsoidInteriorExteriorSpatialFunction
Daoued23
mlt.khawla at gmail.com
Sun Jul 5 17:13:05 EDT 2015
Hello everyone,
I am currently working on a project in which I am using the class
"itk::EllipsoidInteriorExteriorSpatialFunction".
I am trying to get ellipsoids with different orientation and different
sizes.
In the volumic case, it works of course.
In the Planar case, it works perfectly, I just had to set the length of one
of the axes to one voxel, and I get my ellipse just fine.
But in the linear case, which means I set the length of two axes to one
voxel and the third one to a value>1, it works only in the vertical and the
horizontal directions, and I get as a result a vertical or horizontal line.
But with an arbitrary direction, it doesn't work and I don't get anything,
can you please explain to me why? and what should i do to make it work?
I used the examle of public wiki, and i made some adjustments. I just pasted
a part of the code that deals with the ellipsoid creation,I set an
arbitrary direction and the length of the axes as [1,1,7].
Thanks
{//*********************ELLIPSOID DRAWING****************
// Image size and spacing parameters
unsigned long xExtent = 9;
unsigned long yExtent = 9;
unsigned long zExtent = 9;
unsigned long structuringElementSize[] = { xExtent, yExtent, zExtent };
double structuringElementSpacing[] = { 1.0, 1.0, 1.0 };
double structuringElementOrigin[] = { 0, 0, 0 };
// Calculate image volume
unsigned long imageVolume = xExtent * yExtent * zExtent;
// Image typedef
typedef itk::Image< unsigned char, Dimension> ImageType;
// Creates the structuringElement (but doesn't set the size or allocate
memory)
ImageType::Pointer structuringElement = ImageType::New();
structuringElement->SetOrigin(structuringElementOrigin);
structuringElement->SetSpacing(structuringElementSpacing);
std::cout << "New physical structuringElement created\n";
//-----The following block allocates the structuringElement-----
// Create a size object native to the structuringElement type
ImageType::SizeType structuringElementSizeObject;
// Set the size object to the array defined earlier
structuringElementSizeObject.SetSize(structuringElementSize);
// Create a region object native to the structuringElement type
ImageType::RegionType largestPossibleRegion;
// Resize the region
largestPossibleRegion.SetSize(structuringElementSizeObject);
// Set the largest legal region size (i.e. the size of the whole
structuringElement) to what we just defined
structuringElement->SetLargestPossibleRegion(largestPossibleRegion);
// Set the buffered region
structuringElement->SetBufferedRegion(largestPossibleRegion);
// Set the requested region
structuringElement->SetRequestedRegion(largestPossibleRegion);
// Now allocate memory for the structuringElement
structuringElement->Allocate();
std::cout << "New physical structuringElement allocated\n";
// Initialize the image to hold all 0
itk::ImageRegionIterator<ImageType> it
=itk::ImageRegionIterator<ImageType>(structuringElement,
largestPossibleRegion);
unsigned long numImagePixels = 0;
unsigned char exteriorPixelValue = 0;
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
it.Set(exteriorPixelValue);
++numImagePixels;
}
//-----Create ellipsoid in structuringElement-----------------
// Ellipsoid spatial function typedef
typedef itk::EllipsoidInteriorExteriorSpatialFunction<Dimension>
EllipsoidFunctionType;
// Create an ellipsoid spatial function for the source image
EllipsoidFunctionType::Pointer spatialFunc = EllipsoidFunctionType::New();
// Define and set the axes lengths for the ellipsoid
EllipsoidFunctionVectorType axes;
axes[0] =1;
axes[1] =1;
axes[2] = 7;
cout << "Axes: " << axes << endl;
spatialFunc->SetAxes(axes);
// Define and set the center of the ellipsoid in physical space
EllipsoidFunctionVectorType center;
center[0] = 4.5;
center[1] = 4.5;
center[2] = 4.5;
spatialFunc->SetCenter(center);
std::cout << center << std::endl;
// Define the orientations of the ellipsoid axes, vectors must be normalized
// (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 lenght axes[2]
//double data[] = {0, 1, 0, 1, 0, 0, 0, 0, 1};
double data[] ={0.7,-0.7,0,0.408,0.408,-0.816, 0.57735,0.57735,0.57735}; //
I just took an exemple
vnl_matrix<double> orientations(data, 3, 3);
// Set the orientations of the ellipsoids
spatialFunc->SetOrientations(orientations);
//cout<<"spatialfunction"<<spatialFunc<<endl;
typedef ImageType::IndexType IndexType;
typedef IndexType::IndexValueType IndexValueType;
IndexType seedPos;
seedPos[0] = static_cast<IndexValueType>(center[0]);
seedPos[1] = static_cast<IndexValueType>(center[1]);
seedPos[2] = static_cast<IndexValueType>(center[2]);
itk::FloodFilledSpatialFunctionConditionalIterator<ImageType,
EllipsoidFunctionType>
sfi = itk::FloodFilledSpatialFunctionConditionalIterator<ImageType,
EllipsoidFunctionType>(structuringElement, spatialFunc, seedPos);
// Iterate through the entire image and set interior pixels to 255
int numInteriorPixels1 = 0;
unsigned char interiorPixelValue = 255;
for (; !sfi.IsAtEnd(); ++sfi)
{
sfi.Set(interiorPixelValue);
++numInteriorPixels1;
}
cout << "numInteriorPixels1 : " << numInteriorPixels1 << endl;
ImageType::PixelType apixel;
int numExteriorPixels = 0; // Number of pixels not filled by spatial
function
int numInteriorPixels2 = 0; // Number of pixels filled by spatial function
int numErrorPixels = 0; // Number of pixels not set by spatial function
ImageType::IndexValueType indexarray[3] = { 0, 0, 0 };
// Iterate through source image and get pixel values and count pixels
// iterated through, not filled by spatial function, filled by spatial
// function, and not set by the spatial function.
for (unsigned long x = 0; x < xExtent; x++)
{
for (unsigned long y = 0; y < yExtent; y++)
{
for (unsigned long z = 0; z < zExtent; z++)
{
indexarray[0] = x;
indexarray[1] = y;
indexarray[2] = z;
ImageType::IndexType index;
index.SetIndex(indexarray);
apixel =
structuringElement->GetPixel(index);
if (apixel == exteriorPixelValue)
{
++numExteriorPixels;
}
else
{
if (apixel == interiorPixelValue)
{
++numInteriorPixels2;
}
else
{
if (apixel !=
interiorPixelValue || apixel != exteriorPixelValue)
{
++numErrorPixels;
}
}
}
}
}
}
cout << "numInteriorPixels2 : " << numInteriorPixels2 << endl;
std::cout << "Number of interior, exterior and error pixels: "
<< numInteriorPixels2 << ", "
<< numExteriorPixels << ", "
<< numErrorPixels << std::endl;
// Check to see that number of pixels within ellipsoid are equal
// for different iteration loops.
int numInteriorPixels;
if(numInteriorPixels1 == numInteriorPixels2)
{
std::cerr << "numInteriorPixels1 != numInteriorPixels2" <<
std::endl;
return EXIT_FAILURE;
}
else
{
numInteriorPixels = numInteriorPixels1;
}
// Volume of ellipsoid using V=(4/3)*pi*(a/2)*(b/2)*(c/2)
double ellipsoidVolume = 4.18879013333*(axes[0] / 2)*(axes[1] / 2)*(axes[2]
/ 2);
// Percent difference in volume measurement and calculation.
double volumeError = (fabs(ellipsoidVolume - numInteriorPixels2) /
ellipsoidVolume) * 100;
// Test the center of the ellipsoid which should be within the sphere
// and return 1.
double testPosition[Dimension];
bool functionValue;
testPosition[0] = center[0];
testPosition[1] = center[1];
testPosition[2] = center[2];
functionValue = spatialFunc->Evaluate(testPosition);
// 20% error was randomly chosen as a successful ellipsoid fill.
// This should actually be some function of the image/ellipsoid size.
if (volumeError > 20 || functionValue == 0)
{
std::cerr << std::endl << "calculated ellipsoid volume = " <<
ellipsoidVolume << std::endl
<< "measured ellipsoid volume = " <<
numInteriorPixels << std::endl
<< "volume error = " << volumeError << "%" <<
std::endl
<< "function value = " << functionValue << std::endl
<< "itkEllipsoidInteriorExteriorSpatialFunction
failed :(" << std::endl;
return EXIT_FAILURE;
}
else if (numImagePixels != (imageVolume))
{
// Make sure that the number of pixels iterated through from source
image
// is equal to the pre-defined image size.
std::cerr << "Number of pixels iterated through in
structuringElement = "
<< numImagePixels << std::endl;
return EXIT_FAILURE;
}
else
{
std::cout << std::endl << "calculated ellipsoid volume = " <<
ellipsoidVolume << std::endl
<< "measured ellipsoid volume = " <<
numInteriorPixels << std::endl
<< "volume error = " << volumeError << "%" <<
std::endl
<< "function value = " << functionValue << std::endl
<< "itkEllipsoidInteriorExteriorSpatialFunction
ended succesfully!" <<
std::endl;
}
Thank you very much for your help.
Best Regards
Daoud
--
View this message in context: http://itk-users.7.n7.nabble.com/Linear-case-of-itk-EllipsoidInteriorExteriorSpatialFunction-tp35872.html
Sent from the ITK - Users mailing list archive at Nabble.com.
More information about the Insight-users
mailing list