[ITK] [ITK-users] itk::EllipsoidInteriorExteriorSpatialFunction

Daoued23 mlt.khawla at gmail.com
Fri Jul 3 08:55:25 EDT 2015


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,10].

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] = 10;
				
			
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;
}














siavashk wrote
> It would help if you pasted your code. I think you might have a small bug
> that you overlooked.
> 
> 
> 
> --
> View this message in context:
> http://itk-insight-users.2283740.n2.nabble.com/itk-EllipsoidInteriorExteriorSpatialFunction-tp7587522p7587523.html
> Sent from the ITK Insight Users mailing list archive at Nabble.com.
> _____________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
> 
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
> 
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users





--
View this message in context: http://itk-users.7.n7.nabble.com/Re-ITK-users-itk-EllipsoidInteriorExteriorSpatialFunction-tp35851p35864.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users


More information about the Community mailing list