[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&lt;&lt;endl;
typedef   ImageType::IndexType       IndexType;
typedef   IndexType::IndexValueType   IndexValueType;

IndexType seedPos;
seedPos[0] = static_cast&lt;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