[Insight-users] Problem with output image in custom filter
Peter Cech
pcech at vision.ee.ethz.ch
Wed Dec 8 07:39:43 EST 2004
Hi,
I tried to write an image filter that performs labels per-slice
connected components (independent labeling in each slice of 3D image).
My approach is to loop through slices and process each by 2D
ConnectedComponetsImageFilter.
The problem occurs when trying to copy labeled slice into output image.
Application segfaults because of null pointer access. After adding index
of output iterator into debugging output, application terminates with
"Floating exception" which seems to be caused by integer division by
zero (zero comes from ImageBase::m_OffsetTable). It seems like output
image was not allocated/initialized prior to calling GenerateData().
I compared my filter to other filters in ITK, but did not notice they do
anything more to initialize output image. Could somebody have a look at
the attached sources, please?
I'm using CVS version of ITK as of November 15, 2004.
Regards,
Peter Cech
-------------- next part --------------
CMAKE_MINIMUM_REQUIRED(VERSION 1.8)
FIND_PACKAGE(ITK)
IF(ITK_FOUND)
INCLUDE(${ITK_USE_FILE})
ELSE(ITK_FOUND)
MESSAGE(FATAL_ERROR
"Cannot build SpineSeg without ITK. Please set ITK_DIR.")
ENDIF(ITK_FOUND)
ADD_EXECUTABLE(PerSliceConnectedComponentImageFilterTest
PerSliceConnectedComponentImageFilterTest.cpp)
TARGET_LINK_LIBRARIES(PerSliceConnectedComponentImageFilterTest
ITKCommon)
-------------- next part --------------
#ifndef PerSliceConnectedComponentImageFilter_h
#define PerSliceConnectedComponentImageFilter_h
#include "PerSliceImageFilter.h"
template <class TInputImage, class TOutputImage>
class PerSliceConnectedComponentImageFilter
: public PerSliceImageFilter<TInputImage, TOutputImage>
{
public:
/** Standard typedefs. */
typedef PerSliceConnectedComponentImageFilter<TInputImage, TOutputImage> Self;
typedef PerSliceImageFilter<TInputImage, TOutputImage> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Runtime information support. */
itkTypeMacro(PerSliceConnectedComponentImageFilter, PerSliceImageFilter);
/** Image related typedefs. */
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::InputImagePointer InputImagePointer;
typedef typename Superclass::InputImageConstPointer InputImageConstPointer;
typedef typename Superclass::InputImageRegionType InputImageRegionType;
typedef typename Superclass::InputImagePixelType InputImagePixelType;
typedef typename Superclass::OutputImageType OutputImageType;
// typedef typename Superclass::OutputImagePointer OutputImagePointer;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
// typedef typename Superclass::OutputImagePixelType OutputImagePixelType;
/** Slice number type. */
typedef typename Superclass::SliceNumberType SliceNumberType;
/** Slice types. */
typedef typename Superclass::InputSliceImageType InputSliceImageType;
typedef typename Superclass::InputSliceImagePointer InputSliceImagePointer;
typedef typename Superclass::OutputSliceImageType OutputSliceImageType;
/** Pixel types. */
typedef typename Superclass::InputPixelType InputPixelType;
typedef typename Superclass::OutputPixelType OutputPixelType;
protected:
PerSliceConnectedComponentImageFilter();
~PerSliceConnectedComponentImageFilter() {}
/** Actual slice processing. Called once for each slice of image. */
virtual void ProcessSlice(SliceNumberType slice_number);
private:
PerSliceConnectedComponentImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
};
#ifndef ITK_MANUAL_INSTANTIATION
#include "PerSliceConnectedComponentImageFilter.txx"
#endif
#endif // PerSliceConnectedComponentImageFilter_h
-------------- next part --------------
#include "PerSliceConnectedComponentImageFilter.h"
#include <itkConnectedComponentImageFilter.h>
template <class TInputImage, class TOutputImage>
PerSliceConnectedComponentImageFilter<TInputImage, TOutputImage>
::PerSliceConnectedComponentImageFilter()
: Superclass()
{
}
template <class TInputImage, class TOutputImage>
void
PerSliceConnectedComponentImageFilter<TInputImage, TOutputImage>
::ProcessSlice(SliceNumberType slice_number)
{
typedef itk::ConnectedComponentImageFilter<InputSliceImageType, OutputSliceImageType> ConnectedFilterType;
InputSliceImagePointer input_slice = GetInputSlice(slice_number);
typename ConnectedFilterType::Pointer connected_filter = ConnectedFilterType::New();
connected_filter->SetInput(input_slice);
connected_filter->SetFullyConnected(true);
connected_filter->Update();
std::cout << connected_filter << std::endl;
std::cout << connected_filter->GetOutput() << std::endl;
SetOutputSlice(connected_filter->GetOutput(), slice_number);
}
-------------- next part --------------
#include "PerSliceConnectedComponentImageFilter.h"
#include <itkImageRegionIteratorWithIndex.h>
#include <itkRandomImageSource.h>
int PerSliceConnectedComponentImageFilterTest(int, char **)
{
// Define the dimension of the images
const unsigned int ImageDimension = 3;
// Size of test image
const unsigned int size = 5;
// slicing direction to test
const unsigned int slicing_direction = 1;
// Declare the types of the images
typedef itk::Image<unsigned char, ImageDimension> InputImageType;
typedef itk::Image<unsigned int, ImageDimension> OutputImageType;
typedef InputImageType::PixelType InputPixelType;
typedef OutputImageType::PixelType OutputPixelType;
// Declare iterator type
typedef itk::ImageRegionIteratorWithIndex<
InputImageType> InputIteratorType;
typedef itk::ImageRegionIteratorWithIndex<
OutputImageType> OutputIteratorType;
// Use a random image source as input
typedef itk::RandomImageSource<InputImageType> SourceType;
SourceType::Pointer source = SourceType::New();
unsigned long sizeArray[ImageDimension] = { size, size, size };
source->SetMin( itk::NumericTraits<InputPixelType>::Zero );
source->SetMax( itk::NumericTraits<InputPixelType>::One );
source->SetSize( sizeArray );
// Declare the type for the binary threshold filter
typedef PerSliceConnectedComponentImageFilter< InputImageType,
OutputImageType > FilterType;
// Create a filter
FilterType::Pointer filter = FilterType::New();
// Setup filter
filter->SetSlicingDirection(slicing_direction);
filter->Print( std::cout );
// exercise Get methods
std::cout << "SlicingDirection: " << filter->GetSlicingDirection() << std::endl;
// Connect the input images
InputImageType::Pointer input_image = source->GetOutput();
filter->SetInput(input_image);
// Get the Smart Pointer to the Filter Output
OutputImageType::Pointer outputImage = filter->GetOutput();
// Execute the filter
try
{
filter->Update();
}
catch(...)
{
std::cerr << "Caught an unexpected exception. " << std::endl;
std::cerr << "Test failed. " << std::endl;
return EXIT_FAILURE;
}
/* TODO: somehow to test for correct result
// Create an iterator for going through the image output
InputIteratorType it( source->GetOutput(), source->GetOutput()->GetRequestedRegion() );
OutputIteratorType ot(outputImage, outputImage->GetRequestedRegion());
// Check the content of the result image
std::cout << "Verification of the output " << std::endl;
ot.GoToBegin();
it.GoToBegin();
while( !ot.IsAtEnd() )
{
const InputPixelType input = it.Get();
const OutputPixelType output = ot.Get();
std::cout << (double) input << " " << (double) output << std::endl;
// compute mean
double mean = 0;
InputIteratorType::IndexType index = it.GetIndex();
for (signed int i = -(signed int)radius; i <= (signed int)radius; ++i)
for (signed int j = -(signed int)radius; j <= (signed int)radius; ++j)
{
// boundary condition
InputIteratorType::IndexType sample_index = index;
if (sample_index[0]+i < 0)
sample_index[0] = 0;
else if (sample_index[0]+i >= size)
sample_index[0] = size-1;
else
sample_index[0] += i;
if (sample_index[1]+j < 0)
sample_index[1] = 0;
else if (sample_index[1]+j >= size)
sample_index[1] = size-1;
else
sample_index[1] += j;
mean += input_image->GetPixel(sample_index);
}
mean /= diameter*diameter;
// compute standard deviation
double sigma = 0;
for (signed int i = -(signed int)radius; i <= (signed int)radius; ++i)
for (signed int j = -(signed int)radius; j <= (signed int)radius; ++j)
{
// boundary condition
InputIteratorType::IndexType sample_index = index;
if (sample_index[0]+i < 0)
sample_index[0] = 0;
else if (sample_index[0]+i >= size)
sample_index[0] = size-1;
else
sample_index[0] += i;
if (sample_index[1]+j < 0)
sample_index[1] = 0;
else if (sample_index[1]+j >= size)
sample_index[1] = size-1;
else
sample_index[1] += j;
InputPixelType pixel = input_image->GetPixel(sample_index);
sigma += (pixel-mean)*(pixel-mean);
}
sigma /= diameter*diameter;
using namespace std;
sigma = sqrt(sigma);
// expected threshloding result for this pixel
OutputPixelType expected_output_pixel;
if (input >= mean+K*sigma)
expected_output_pixel = inside;
else
expected_output_pixel = outside;
if(ot.Get() != expected_output_pixel)
{
std::cerr << "Error in NiblackThresholdImageFilterTest " << std::endl;
std::cout << "OutsideValue: " << filter->GetOutsideValue() << std::endl;
std::cout << "InsideValue: " << filter->GetInsideValue() << std::endl;
std::cout << "K: " << filter->GetK() << std::endl;
std::cout << "NeighborhoodRadius: " << filter->GetNeighborhoodRadius() << std::endl;
std::cout << "Expected output: " << expected_output_pixel << std::endl;
std::cout << std::endl;
std::cout << " mean: " << mean << std::endl;
std::cout << " sigma: " << sigma << std::endl;
std::cout << " input = " << (double) input << std::endl;
std::cout << " output = " << (double) output;
std::cout << std::endl;
return EXIT_FAILURE;
}
++ot;
++it;
}
*/
std::cout << "Test passsed. " << std::endl;
return EXIT_SUCCESS;
}
int main(int argc, char ** argv)
{
try {
return PerSliceConnectedComponentImageFilterTest(argc, argv);
} catch (std::exception & e) {
printf("Exception caught: %s\n", e.what());
return 99;
} catch (...) {
printf("Exception caught!");
return 99;
}
}
-------------- next part --------------
#ifndef PerSliceImageFilter_h
#define PerSliceImageFilter_h
#include <itkImageToImageFilter.h>
template <class TInputImage, class TOutputImage>
class PerSliceImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>
{
public:
/** Standard typedefs. */
typedef PerSliceImageFilter<TInputImage, TOutputImage> Self;
typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Runtime information support. */
itkTypeMacro(PerSliceImageFilter, ImageToImageFilter);
/** Image related typedefs. */
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::InputImagePointer InputImagePointer;
typedef typename Superclass::InputImageConstPointer InputImageConstPointer;
typedef typename Superclass::InputImageRegionType InputImageRegionType;
typedef typename Superclass::InputImagePixelType InputImagePixelType;
typedef typename Superclass::OutputImageType OutputImageType;
// typedef typename Superclass::OutputImagePointer OutputImagePointer;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
// typedef typename Superclass::OutputImagePixelType OutputImagePixelType;
/** ImageDimension constants. */
itkStaticConstMacro(InputImageDimension, unsigned int,
TInputImage::ImageDimension);
itkStaticConstMacro(OutputImageDimension, unsigned int,
TOutputImage::ImageDimension);
/** The input and output image types must have the same dimension. */
itkConceptMacro(SameDimension,
(itk::Concept::SameDimension<itkGetStaticConstMacro(InputImageDimension), itkGetStaticConstMacro(OutputImageDimension)>));
// TODO: Input/Output image dimension has to be at least 1
/** Pixel types. */
typedef typename TInputImage::PixelType InputPixelType;
typedef typename TOutputImage::PixelType OutputPixelType;
/** Slice number type. */
typedef typename InputImageType::SizeValueType SliceNumberType;
/** Slice types. */
typedef itk::Image<InputPixelType, InputImageDimension-1> InputSliceImageType;
typedef typename InputSliceImageType::Pointer InputSliceImagePointer;
typedef itk::Image<OutputPixelType, TOutputImage::ImageDimension-1> OutputSliceImageType;
/** Set/get slicing direction. Default slicing direction is
* InputImageDimension-1. */
void SetSlicingDirection(unsigned int slicing_direction)
{
itkDebugMacro("setting SlicingDirection to " << slicing_direction);
if (slicing_direction >= InputImageDimension)
itkExceptionMacro( << "SlicingDirection too high (" << slicing_direction <<", max. is " << InputImageDimension << ")");
if (this->m_SlicingDirection != slicing_direction)
{
this->m_SlicingDirection = slicing_direction;
this->Modified();
}
}
itkGetMacro(SlicingDirection, unsigned int);
/** So far, only sigle threaded. Need to properly partition input (in slicing
* direction) to run more slices concurrently. */
void GenerateData();
protected:
PerSliceImageFilter();
~PerSliceImageFilter() {}
virtual void PrintSelf(std::ostream& os, itk::Indent indent);
/** Get region of slice in input/output image. */
InputImageRegionType GetInputSliceRegion(SliceNumberType slice_number);
OutputImageRegionType GetOutputSliceRegion(SliceNumberType slice_number);
/** Get slice of input image. */
InputSliceImagePointer GetInputSlice(SliceNumberType slice_number);
/** Set slice in output image. */
void SetOutputSlice(const OutputSliceImageType * slice, SliceNumberType slice_number);
/** Slice processing. Called once for each slice of image. */
virtual void ProcessSlice(SliceNumberType slice_number) = 0;
private:
PerSliceImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
unsigned int m_SlicingDirection;
};
#ifndef ITK_MANUAL_INSTANTIATION
#include "PerSliceImageFilter.txx"
#endif
#endif // PerSliceImageFilter_h
-------------- next part --------------
#include "PerSliceImageFilter.h"
#include <itkExtractImageFilter.h>
#include <itkImageRegionIterator.h>
template <class TInputImage, class TOutputImage>
PerSliceImageFilter<TInputImage, TOutputImage>
::PerSliceImageFilter()
: Superclass(), m_SlicingDirection(InputImageDimension-1)
{
}
template <class TInputImage, class TOutputImage>
void
PerSliceImageFilter<TInputImage, TOutputImage>
::GenerateData()
{
//TODO precondition: GetInput() != 0
// Workaround: ComputeOffsetTable on output image was never called, force call it
this->GetOutput()->SetBufferedRegion(this->GetOutput()->GetLargestPossibleRegion());
this->GetOutput()->Initialize();
SliceNumberType slice_count = this->GetInput()->GetLargestPossibleRegion().GetSize(this->m_SlicingDirection);
for (SliceNumberType slice = 0; slice < slice_count; ++slice)
this->ProcessSlice(slice);
}
template <class TInputImage, class TOutputImage>
typename PerSliceImageFilter<TInputImage, TOutputImage>::InputImageRegionType
PerSliceImageFilter<TInputImage, TOutputImage>
::GetInputSliceRegion(SliceNumberType slice_number)
{
//TODO precondition: GetInput() != 0
//TODO precondition: slice_number < GetInput()->GetLargestPossibleRegion().GetSize(this->m_SlicingDirection)
InputImageRegionType region(this->GetInput()->GetLargestPossibleRegion());
region.SetSize(this->m_SlicingDirection, 1);
region.SetIndex(this->m_SlicingDirection, slice_number);
return region;
}
template <class TInputImage, class TOutputImage>
typename PerSliceImageFilter<TInputImage, TOutputImage>::OutputImageRegionType
PerSliceImageFilter<TInputImage, TOutputImage>
::GetOutputSliceRegion(SliceNumberType slice_number)
{
//TODO precondition: GetOutput() != 0
//TODO precondition: slice_number < GetOutput()->GetLargestPossibleRegion().GetSize(this->m_SlicingDirection)
OutputImageRegionType region(this->GetOutput()->GetLargestPossibleRegion());
region.SetSize(this->m_SlicingDirection, 1);
region.SetIndex(this->m_SlicingDirection, slice_number);
return region;
}
template <class TInputImage, class TOutputImage>
typename PerSliceImageFilter<TInputImage, TOutputImage>::InputSliceImagePointer
PerSliceImageFilter<TInputImage, TOutputImage>
::GetInputSlice(SliceNumberType slice_number)
{
typedef itk::ExtractImageFilter<InputImageType, InputSliceImageType> ExtractorType;
typename ExtractorType::Pointer extractor = ExtractorType::New();
extractor->SetInput(this->GetInput());
InputImageRegionType region(this->GetInputSliceRegion(slice_number));
region.SetSize(this->m_SlicingDirection, 0); // collapse in slicing direction
extractor->SetExtractionRegion(region);
extractor->Update();
return extractor->GetOutput();
}
template <class TInputImage, class TOutputImage>
void
PerSliceImageFilter<TInputImage, TOutputImage>
::SetOutputSlice(const OutputSliceImageType * slice, SliceNumberType slice_number)
{
//TODO precondition: slice != 0
//TODO precondition: GetOutput() != 0
//TODO precondition: slice_number < GetOutput()->GetLargestPossibleRegion().GetSize(this->m_SlicingDirection)
//TODO precondition: slice size == output slice region size
printf("Slice: %p\n", slice);
printf("Output: %p\n", this->GetOutput());
std::cout << "Output region: " << this->GetOutputSliceRegion(slice_number).GetIndex() << std::endl;
std::cout << "Output region: " << this->GetOutputSliceRegion(slice_number).GetSize() << std::endl;
std::cout << "Output maximum: " << this->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl;
typedef itk::ImageRegionIterator<const OutputSliceImageType> SliceIterator;
typedef itk::ImageRegionIterator<OutputImageType> ImageIterator;
ImageIterator oit(this->GetOutput(), this->GetOutputSliceRegion(slice_number));
SliceIterator sit(slice, slice->GetLargestPossibleRegion());
for (oit.GoToBegin(), sit.GoToBegin(); !oit.IsAtEnd(); ++oit, ++sit) {
std::cout << "InIt: " << sit.GetIndex() << std::endl;
std::cout << "OutIt: " << oit.GetIndex()[0] << "," << oit.GetIndex()[1] << std::endl;
std::cout << "oit image: " << oit.GetImage() << std::endl;
std::cout << "oit region: " << oit.GetRegion().GetIndex() << std::endl;
std::cout << "oit region: " << oit.GetRegion().GetSize() << std::endl;
oit.Set(sit.Get());
}
}
template <class TInputImage, class TOutputImage>
void
PerSliceImageFilter<TInputImage, TOutputImage>
::PrintSelf(std::ostream& os, itk::Indent indent)
{
// print inherited stuff
Superclass::PrintSelf(os, indent);
// print own stuff
os << indent << "SlicingDirection: "
<< static_cast<typename itk::NumericTraits<unsigned int>::PrintType>(this->m_SlicingDirection) << std::endl;
}
More information about the Insight-users
mailing list