[Insight-users] discrete gaussian filter causes exception
Toon Huysmans
denhuys at hotmail.com
Thu May 27 06:17:48 EDT 2004
Hi,
I have got a strange error. I simply want to blur a 3D image. When I do
this using the recursive gaussian filter, there are no errors, but the
smoothing is only in one dimension.
So now I use the discrete gaussian filter, and I didn`t change my code
except for the setting of the parameter Variance. But the filter throws an
exception in itkUnaryFunctorImageFilter.txx at line 147.
extract from itkUnaryFunctorImageFilter.txx:
----------------------------------------------------------------------------
----------------------------------------------------------------------------
------------
/**
* ThreadedGenerateData Performs the pixel-wise addition
*/
template <class TInputImage, class TOutputImage, class TFunction >
void
UnaryFunctorImageFilter<TInputImage,TOutputImage,TFunction>::ThreadedGenerat
eData( const OutputImageRegionType &outputRegionForThread,int threadId)
{
InputImagePointer inputPtr = this->GetInput();
OutputImagePointer outputPtr = this->GetOutput(0);
// Define the portion of the input to walk for this thread, using
// the CallCopyOutputRegionToInputRegion method allows for the input
// and output images to be different dimensions
InputImageRegionType inputRegionForThread;
this->CallCopyOutputRegionToInputRegion(inputRegionForThread,
outputRegionForThread);
// Define the iterators
ImageRegionConstIterator<TInputImage> inputIt(inputPtr,
inputRegionForThread);
ImageRegionIterator<TOutputImage> outputIt(outputPtr,
outputRegionForThread);
ProgressReporter progress(this, threadId,
outputRegionForThread.GetNumberOfPixels());
inputIt.GoToBegin();
outputIt.GoToBegin();
while( !inputIt.IsAtEnd() )
{
outputIt.Set( m_Functor( inputIt.Get() ) );
(/******HERE******/)
++inputIt;
++outputIt;
progress.CompletedPixel(); // potential exception thrown here
}
}
----------------------------------------------------------------------------
----------------------------------------------------------------------------
------------
The Exception:
Unhandled exception at 0x0042a848 in THCochleaSegmentationLevelSet.exe:
0xC0000005: Access violation reading location 0x0f171000.
The callstack on exception is:
----------------------------------------------------------------------------
----------------------------------------------------------------------------
------------
THCochleaSegmentationLevelSet.exe!itk::UnaryFunctorImageFilter<itk::Image<fl
oat,3>,itk::Image<unsigned char,3>,itk::Functor::Cast<float,unsigned char>
>::ThreadedGenerateData(const itk::ImageRegion<3> &
outputRegionForThread={...}, int threadId=1) Line 147 + 0x4 C++
THCochleaSegmentationLevelSet.exe!itk::ImageSource<itk::Image<unsigned
char,3> >::ThreaderCallback(void * arg=0x04f283d8) Line 281 C++
msvcr71.dll!_threadstartex(void * ptd=0x04f62a78) Line 241 + 0x6 C
kernel32.dll!77e7d33b()
THCochleaSegmentationLevelSet.exe!itk::Object::~Object() Line 456 + 0x12
C++
----------------------------------------------------------------------------
----------------------------------------------------------------------------
------------
I have attached my code, the gaussian filter is named blur in my code.
It seems the code crashes when doing itkExporter->Update(). The pipeline on
that moment is:
ImageFileReader<ImageType>
CastImageFilter< ImageType, RealImageType >
DiscreteGaussianImageFilter<RealImageType,RealImageType>
BinaryThresholdImageFilter<RealImageType,RealImageType>
CastImageFilter< RealImageType, ImageType >
ImageToVTKImageFilter<ImageType>
vtkOrthogonalImagePlaneSet
I know this probably is a fuzzy post but the reason I post this, is that by
replacing a simple filter I get an Exception, so that it could be a bug in
ITK.
Thanks,
Toon.
-------------- next part --------------
//------------------------------------------------------------------------------
#ifndef THCochleaSegmentationBase_h
#define THCochleaSegmentationBase_h
//------------------------------------------------------------------------------
#include "vtkRenderer.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderWindow.h"
#include "vtkOrthogonalImagePlaneSet.h"
#include "vtkOutlineFilter.h"
#include "vtkActor.h"
#include "vtkConeSource.h"
#include "vtkPolyDataMapper.h"
#include "vtkDataSet.h"
#include "vtkImageData.h"
//------------------------------------------------------------------------------
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkExceptionObject.h"
#include "itkCastImageFilter.h"
#include "itkShrinkImageFilter.h"
#include "itkImageToVTKImageFilter.h"
#include "itkWatershedImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeImageFilter.h"
#include "itkBilateralImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
//------------------------------------------------------------------------------
#include <string>
using std::string;
//------------------------------------------------------------------------------
class THCochleaSegmentationLevelSetBase
{
private:
typedef itk::Image< unsigned char , 3 > ImageType;
typedef itk::Image< float , 3 > RealImageType;
typedef itk::ImageFileReader<ImageType> ReaderType;
typedef itk::ImageFileWriter<ImageType> WriterType;
typedef itk::WatershedImageFilter<RealImageType>
WatershedFilterType;
typedef itk::Image<unsigned long, 3> LabeledImageType;
typedef itk::CastImageFilter< ImageType, RealImageType >
CastFilterCToF;
typedef itk::CastImageFilter< RealImageType, ImageType >
CastFilterFToC;
typedef itk::CastImageFilter< LabeledImageType, ImageType >
CastFilterLToC;
typedef itk::ImageToVTKImageFilter<ImageType>
ExporterType;
typedef itk::CurvatureAnisotropicDiffusionImageFilter<RealImageType,RealImageType>
DiffusionFilterType;
typedef itk::BilateralImageFilter<RealImageType,RealImageType>
BilateralFilterType;
typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<RealImageType,RealImageType>
GradientMagnitudeFilterType;
typedef itk::SigmoidImageFilter<RealImageType,RealImageType>
SigmoidFilterType;
typedef itk::BinaryThresholdImageFilter<RealImageType,RealImageType>
BinaryThresholdFilterType;
typedef itk::GeodesicActiveContourLevelSetImageFilter<RealImageType,RealImageType>
GACFilterType;
typedef itk::DiscreteGaussianImageFilter<RealImageType,RealImageType>
GaussianFilterType;
public:
THCochleaSegmentationLevelSetBase(void);
~THCochleaSegmentationLevelSetBase(void);
void SetRenderer(vtkRenderer* ren);
void Initialize();
void ReadVolume(string filename);
//Thresholding
void SetLowerThreshold(float t);
void SetUpperThreshold(float t);
//smoothing
void SetSigmaInitial(float s);
//smoothing
void SetDiffusionTimeStep(float t);
void SetDiffusionIterations(int i);
void SetDiffusionConductance(float c);
void UpdateDiffusion();
//gradient
void SetGradientSigma(float s);
//sigmoid
void SetSigmoidAlpha(float a);
void SetSigmoidBeta(float b);
//GAC
void SetSigmaGAC(float s);
void SetPropagationScale(float s);
void SetAdvectionScale(float s);
void SetCurvatureScale(float s);
void SetMaximumIterations(int i);
void SetMaximumRMSError(float f);
void WriteVolume(string filename);
void view(int i);
protected:
DiffusionFilterType::Pointer diffusion;
ExporterType::Pointer itkExporter;
vtkOrthogonalImagePlaneSet *planeSet;
CastFilterCToF::Pointer castFilter1;
CastFilterLToC::Pointer castFilter2;
CastFilterFToC::Pointer castFilter3;
vtkRenderer* ren;
ReaderType::Pointer reader;
WriterType::Pointer writer;
BinaryThresholdFilterType::Pointer threshold;
BinaryThresholdFilterType::Pointer finalThreshold;
GradientMagnitudeFilterType::Pointer gradient;
SigmoidFilterType::Pointer sigmoid;
GACFilterType::Pointer gac;
GaussianFilterType::Pointer blur;
};
//------------------------------------------------------------------------------
#endif //THCochleaSegmentationBase_h
//------------------------------------------------------------------------------
-------------- next part --------------
//------------------------------------------------------------------------------
#include ".\thcochleasegmentationLevelSetbase.h"
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
THCochleaSegmentationLevelSetBase::THCochleaSegmentationLevelSetBase(void)
{
reader = ReaderType::New();
diffusion = DiffusionFilterType::New();
castFilter3 = CastFilterFToC::New();
castFilter2 = CastFilterLToC::New();
itkExporter = ExporterType::New();
planeSet = vtkOrthogonalImagePlaneSet::New();
writer = WriterType::New();
threshold = BinaryThresholdFilterType::New();
gradient = GradientMagnitudeFilterType::New();
sigmoid = SigmoidFilterType::New();
gac = GACFilterType::New();
castFilter1 = CastFilterCToF::New();
finalThreshold = BinaryThresholdFilterType::New();
blur = GaussianFilterType::New();
}
//------------------------------------------------------------------------------
THCochleaSegmentationLevelSetBase::~THCochleaSegmentationLevelSetBase(void)
{
reader->Delete();
diffusion->Delete();
castFilter3->Delete();
castFilter2->Delete();
itkExporter->Delete();
planeSet->Delete();
writer->Delete();
threshold->Delete();
gradient->Delete();
sigmoid->Delete();
gac->Delete();
castFilter1->Delete();
finalThreshold->Delete();
blur->Delete();
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetRenderer(vtkRenderer* ren)
{
this->ren = ren;
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::ReadVolume(string filename)
{
//read
reader->SetFileName( filename.c_str() );
reader->Update();
//cast from char to float
castFilter1->SetInput(reader->GetOutput());
//Edge Potential Image...
//denoise
diffusion->SetNumberOfIterations( 0 );
diffusion->SetConductanceParameter( 2.0 );
diffusion->SetTimeStep(0.125);
diffusion->SetInput(castFilter1->GetOutput());
//gradient
gradient->SetInput(diffusion->GetOutput());
gradient->SetSigma(1.0);
//sigmoid
sigmoid->SetInput(gradient->GetOutput());
sigmoid->SetOutputMaximum(255.0);
sigmoid->SetOutputMinimum(0.0);
sigmoid->SetBeta(5.0);
sigmoid->SetAlpha(-1.0);
//Initial Image...
//blur
blur->SetInput(castFilter1->GetOutput());
GaussianFilterType::ArrayType v;
v[0] = 1.0;
v[1] = 1.0;
v[2] = 1.0;
blur->SetVariance(v);
//threshold
threshold->SetInput(blur->GetOutput());
threshold->SetUpperThreshold(255.0);
threshold->SetLowerThreshold(0.0);
threshold->SetOutsideValue(1.0);
threshold->SetInsideValue(-1.0);
//Geodesic Active Contours...
gac->SetDerivativeSigma(1.0);
gac->SetPropagationScaling(1.0);
gac->SetAdvectionScaling(1.0);
gac->SetCurvatureScaling(1.0);
gac->SetNumberOfIterations(100);
gac->SetMaximumRMSError(0.010);
gac->SetInput(threshold->GetOutput());
gac->SetFeatureImage(sigmoid->GetOutput());
//get the negative part only
finalThreshold->SetInput(gac->GetOutput());
finalThreshold->SetLowerThreshold( -1000.0 );
finalThreshold->SetUpperThreshold( 0.0 );
finalThreshold->SetOutsideValue( 0 );
finalThreshold->SetInsideValue( 255 );
castFilter3->SetInput(finalThreshold->GetOutput());
itkExporter->SetInput(castFilter3->GetOutput());
planeSet->Register(ren);
planeSet->SetInteractor(ren->GetRenderWindow()->GetInteractor());
//outline
{
vtkOutlineFilter *outline = vtkOutlineFilter::New();
outline->SetInput((itkExporter->GetOutput()));
vtkPolyDataMapper *outlineMapper = vtkPolyDataMapper::New();
outlineMapper->SetInput(outline->GetOutput());
vtkActor *outlineActor = vtkActor::New();
outlineActor->SetMapper(outlineMapper);
ren->AddActor(outlineActor);
outline->Delete();
outlineMapper->Delete();
outlineActor->Delete();
}
view(0);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::WriteVolume(string filename)
{
//write...
writer->SetFileName( filename.c_str() );
//...current volume
writer->SetInput(castFilter3->GetOutput());
writer->Write();
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetLowerThreshold(float t)
{
threshold->SetLowerThreshold(t);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetUpperThreshold(float t)
{
threshold->SetUpperThreshold(t);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetDiffusionTimeStep(float t)
{
diffusion->SetTimeStep(t);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetDiffusionIterations(int i)
{
diffusion->SetNumberOfIterations( i );
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetDiffusionConductance(float c)
{
diffusion->SetConductanceParameter( c );
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::UpdateDiffusion()
{
diffusion->Update();
ren->GetRenderWindow()->GetInteractor()->Render();
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetGradientSigma(float s)
{
gradient->SetSigma(s);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetSigmoidAlpha(float a)
{
sigmoid->SetAlpha(a);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetSigmoidBeta(float b)
{
sigmoid->SetBeta(b);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetSigmaGAC(float s)
{
gac->SetDerivativeSigma(s);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetPropagationScale(float s)
{
gac->SetPropagationScaling(s);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetAdvectionScale(float s)
{
gac->SetAdvectionScaling(s);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetCurvatureScale(float s)
{
gac->SetCurvatureScaling(s);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetMaximumIterations(int i)
{
gac->SetNumberOfIterations(i);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetMaximumRMSError(float f)
{
gac->SetMaximumRMSError(f);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::SetSigmaInitial(float s)
{
GaussianFilterType::ArrayType v;
v[0] = s;
v[1] = s;
v[2] = s;
blur->SetVariance(v);
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::view(int i)
{
switch(i)
{
case 0:
{
//Original
cout << "Viewing Orignal Image" << endl;
this->castFilter3->SetInput(this->castFilter1->GetOutput());
break;
}
case 1:
{
cout << "Viewing Diffused Image" << endl;
this->castFilter3->SetInput(this->diffusion->GetOutput());
itkExporter->Update();
break;
}
case 2:
{
cout << "Viewing Gradient Image" << endl;
this->castFilter3->SetInput(this->gradient->GetOutput());
itkExporter->Update();
break;
}
case 3:
{
cout << "Viewing Sigmoid Image" << endl;
this->castFilter3->SetInput(this->sigmoid->GetOutput());
itkExporter->Update();
break;
}
case 4:
{
cout << "Viewing Initial Segmentation Image" << endl;
this->castFilter3->SetInput(this->threshold->GetOutput());
itkExporter->Update();
break;
}
case 5:
{
cout << "Viewing Final Segmentation Image" << endl;
this->castFilter3->SetInput(this->finalThreshold->GetOutput());
itkExporter->Update();
break;
}
case 6:
{
cout << "Viewing Final Segmentation Image 3D model" << endl;
cout << "TODO!" << endl;
break;
}
case 7:
{
cout << "Viewing blurred initial Image" << endl;
this->castFilter3->SetInput(this->blur->GetOutput());
itkExporter->Update();
break;
}
default:
{
cout << "Error: view number does not exist..." << endl;
}
}
planeSet->SetInput(itkExporter->GetOutput());
planeSet->SetLUTToRandom(false);
planeSet->On();
ren->GetRenderWindow()->GetInteractor()->Render();
}
//------------------------------------------------------------------------------
void
THCochleaSegmentationLevelSetBase::Initialize()
{
}
//------------------------------------------------------------------------------
More information about the Insight-users
mailing list