[Insight-users] Problems with fixedReader->Update();
motes motes
mort.motes at gmail.com
Tue Jul 21 18:15:19 EDT 2009
I tried adding:
std::cout << "COMPONENT: " << this->GetComponentType() << std::endl;
before the switch as you recommended but nothing gets printed. Should I
rebuild ITK after adding that line?
Below is the code which is a slightly modified version of the
BSplineWarping1.cxx example. I call it with:
std::string configFile = "BSplineDisplacements1.txt";
std::string fixed = "c:/test/lena.png";
std::string moving = "c:/test/lena.png";
std::string deformed = "c:/test/lena_deformed.png";
std::string deformed_field = "c:/test/lena_deformed_field.png";
BSplineWarping2D warp2D;
warp2D.setup();
warp2D.deform(configFile, fixed, moving, deformed, deformed_field);
here is the code (the fieldWriter is in the bottom):
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImage.h"
#include "itkResampleImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkBSplineDeformableTransform.h"
#include <fstream>
#include "itkCommand.h"
class CommandProgressUpdate : public itk::Command
{
public:
typedef CommandProgressUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandProgressUpdate() {};
public:
void Execute(itk::Object *caller, const itk::EventObject & event)
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object * object, const itk::EventObject & event)
{
const itk::ProcessObject * filter =
dynamic_cast< const itk::ProcessObject * >( object );
if( ! itk::ProgressEvent().CheckEvent( &event ) )
{
return;
}
std::cout << filter->GetProgress() << std::endl;
}
};
class BSplineWarping2D {
public:
// Types
static const unsigned int ImageDimension = 2;
static const unsigned int SplineOrder = 3;
static const unsigned int SpaceDimension = ImageDimension;
typedef unsigned
char
PixelType;
//typedef unsigned
short
PixelType;
typedef itk::Image<PixelType,
ImageDimension> FixedImageType;
typedef itk::Image<PixelType,
ImageDimension> MovingImageType;
typedef
itk::ImageFileReader<FixedImageType>
FixedReaderType;
typedef
itk::ImageFileReader<MovingImageType>
MovingReaderType;
typedef
itk::ImageFileWriter<MovingImageType>
MovingWriterType;
typedef itk::ResampleImageFilter<MovingImageType,
FixedImageType> FilterType;
typedef itk::LinearInterpolateImageFunction< MovingImageType,
double> InterpolatorType;
typedef
double
CoordinateRepType;
typedef itk::BSplineDeformableTransform<CoordinateRepType, SpaceDimension,
SplineOrder> TransformType;
typedef
TransformType::RegionType
RegionType;
typedef itk::Point<float,
ImageDimension> PointType;
typedef itk::Vector<float,
ImageDimension> VectorType;
typedef itk::Image<VectorType,
ImageDimension>
DeformationFieldType;
typedef
TransformType::SpacingType
SpacingType;
typedef
TransformType::ParametersType
ParametersType;
typedef
TransformType::OriginType
OriginType;
typedef
itk::ImageRegionIterator<DeformationFieldType>
FieldIterator;
typedef
itk::ImageFileWriter<DeformationFieldType>
FieldWriterType;
// Instances
FixedReaderType::Pointer fixedReader;
MovingReaderType::Pointer movingReader;
MovingWriterType::Pointer movingWriter;
FilterType::Pointer resampler;
InterpolatorType::Pointer interpolator;
TransformType::Pointer bsplineTransform;
DeformationFieldType::Pointer field;
FieldWriterType::Pointer fieldWriter;
BSplineWarping2D(){}
void setup(){
createComponents();
}
void createComponents(){
fixedReader = FixedReaderType::New();
movingReader = MovingReaderType::New();
movingWriter = MovingWriterType::New();
resampler = FilterType::New();
interpolator = InterpolatorType::New();
bsplineTransform = TransformType::New();
field = DeformationFieldType::New();
fieldWriter = FieldWriterType::New();
}
int deform(std::string in_configFile,
std::string in_fixed,
std::string in_moving,
std::string in_deformed,
std::string in_field)
{
fixedReader->SetFileName(in_fixed);
try
{
fixedReader->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
movingReader->SetFileName(in_moving);
movingWriter->SetFileName(in_deformed);
FixedImageType::ConstPointer fixedImage = fixedReader->GetOutput();
resampler->SetInterpolator(interpolator);
FixedImageType::SpacingType fixedSpacing =
fixedImage->GetSpacing();
FixedImageType::PointType fixedOrigin = fixedImage->GetOrigin();
FixedImageType::DirectionType fixedDirection =
fixedImage->GetDirection();
resampler->SetOutputSpacing(fixedSpacing);
resampler->SetOutputOrigin( fixedOrigin );
resampler->SetOutputDirection( fixedDirection );
FixedImageType::RegionType fixedRegion =
fixedImage->GetBufferedRegion();
FixedImageType::SizeType fixedSize = fixedRegion.GetSize();
resampler->SetSize( fixedSize );
resampler->SetOutputStartIndex( fixedRegion.GetIndex() );
resampler->SetInput( movingReader->GetOutput() );
movingWriter->SetInput( resampler->GetOutput() );
// We instantiate now the type of the \code{BSplineDeformableTransform}
using
// as template parameters the type for coordinates representation, the
// dimension of the space, and the order of the B-spline.
//
// \index{BSplineDeformableTransform!New}
// \index{BSplineDeformableTransform!Instantiation}
//
// Since we are using a B-spline of order 3, the coverage of the
BSpling grid
// should exceed by one the spatial extent of the image on the lower
region of
// image indices, and by two grid points on the upper region of image
indices.
// We choose here to use a $8 \times 8$ B-spline grid, from which only
a $5
// \times 5$ sub-grid will be covering the input image. If we use an
input
// image of size $500 \times 500$ pixels, and pixel spacing $2.0 \times
2.0$
// then we need the $5 \times 5$ B-spline grid to cover a physical
extent of $1000
// \times 1000$ mm. This can be achieved by setting the pixel spacing
of the
// B-spline grid to $250.0 \times 250.0$ mm. The origin of the B-spline
grid
// must be set at one grid position away from the origin of the desired
output
// image. All this is done with the following lines of code.
//
// \index{BSplineDeformableTransform}
RegionType bsplineRegion;
RegionType::SizeType size;
const unsigned int numberOfGridNodesOutsideTheImageSupport = 3;
const unsigned int numberOfGridNodesInsideTheImageSupport = 5;
const unsigned int numberOfGridNodes =
numberOfGridNodesInsideTheImageSupport +
numberOfGridNodesOutsideTheImageSupport;
const unsigned int numberOfGridCells =
numberOfGridNodesInsideTheImageSupport - 1;
size.Fill(numberOfGridNodes);
bsplineRegion.SetSize(size);
SpacingType spacing;
OriginType origin;
for( unsigned int i=0; i< SpaceDimension; i++ )
{
spacing[i] = fixedSpacing[i] * (fixedSize[i] - 1) / numberOfGridCells;
}
origin = fixedOrigin - fixedDirection * spacing;
bsplineTransform->SetGridSpacing( spacing );
bsplineTransform->SetGridOrigin( origin );
bsplineTransform->SetGridRegion( bsplineRegion );
bsplineTransform->SetGridDirection( fixedDirection );
const unsigned int numberOfParameters =
bsplineTransform->GetNumberOfParameters();
const unsigned int numberOfNodes = numberOfParameters / SpaceDimension;
ParametersType parameters( numberOfParameters );
//
// The B-spline grid should now be fed with coeficients at each node.
Since
// this is a two dimensional grid, each node should receive two
coefficients.
// Each coefficient pair is representing a displacement vector at this
node.
// The coefficients can be passed to the B-spline in the form of an
array where
// the first set of elements are the first component of the
displacements for
// all the nodes, and the second set of elemets is formed by the second
// component of the displacements for all the nodes.
//
// In this example we read such displacements from a file, but for
convinience
// we have written this file using the pairs of $(x,y)$ displacement
for every
// node. The elements read from the file should therefore be
reorganized when
// assigned to the elements of the array. We do this by storing all the
odd
// elements from the file in the first block of the array, and all the
even
// elements from the file in the second block of the array. Finally the
array
// is passed to the B-spline transform using the
\code{SetParameters()}.
//
// Read nodes from the config file into the parameter array.
std::ifstream infile;
std::wstring widestr =
std::wstring(in_configFile.begin(),in_configFile.end());
const wchar_t* widecstr = widestr.c_str();
infile.open(widecstr);
for( unsigned int n=0; n < numberOfNodes; n++ )
{
infile >> parameters[n];
infile >> parameters[n+numberOfNodes];
}
infile.close();
// Finally the array is passed to the B-spline transform using the
// \code{SetParameters()}.
bsplineTransform->SetParameters(parameters);
CommandProgressUpdate::Pointer observer = CommandProgressUpdate::New();
resampler->AddObserver(itk::ProgressEvent(), observer);
// At this point we are ready to use the transform as part of the
resample
// filter. We trigger the execution of the pipeline by invoking
// \code{Update()} on the last filter of the pipeline, in this case
writer.
resampler->SetTransform( bsplineTransform );
try
{
movingWriter->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Store deformation field.
field->SetRegions(fixedRegion);
field->SetOrigin(fixedOrigin);
field->SetSpacing(fixedSpacing);
field->SetDirection(fixedDirection);
field->Allocate();
FieldIterator fi(field, fixedRegion);
fi.GoToBegin();
TransformType::InputPointType fixedPoint;
TransformType::OutputPointType movingPoint;
DeformationFieldType::IndexType index;
VectorType displacement;
while(!fi.IsAtEnd())
{
index = fi.GetIndex();
field->TransformIndexToPhysicalPoint( index, fixedPoint );
movingPoint = bsplineTransform->TransformPoint( fixedPoint );
displacement = movingPoint - fixedPoint;
fi.Set( displacement );
++fi;
}
if( !in_field.empty() )
{
fieldWriter->SetInput( field );
fieldWriter->SetFileName(in_field);
try
{
fieldWriter->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
}
return EXIT_SUCCESS;
}
};
2009/7/22 Ramón Casero Cañas <ramon.casero at comlab.ox.ac.uk>
> motes motes wrote:
>
>> Ok the error has been related to this line:
>>
>> fieldWriter->Update();
>>
>> all the time, but I still cannot understand why. Have made sure to specify
>> :
>>
>> typedef unsigned char PixelType;
>>
>>
>> Any help appreciated!
>>
>
>
> Hi, I cannot solve the problem, but maybe some indications will help.
>
> Your code is failing because the writer is calling PNGImageIO::WriteSlice,
> and it gets to the following switch statement in itkPNGImageIO.cxx:
>
> <CODE>
> switch (this->GetComponentType())
> {
> case UCHAR:
> bitDepth = 8;
> break;
>
> case USHORT:
> bitDepth = 16;
> break;
>
> default:
> {
> // [cropped]
> ::itk::ExceptionObject excp(__FILE__, __LINE__, "PNG supports unsigned
> char and unsigned short", ITK_LOCATION);
> throw excp;
> }
> }
> </CODE>
>
>
> The component type is defined in Code/IO/itkImageIOBase.h
>
> <CODE>
> /** Enums used to manipulate the component type. The component type
> * refers to the actual storage class associated with either a
> * SCALAR pixel type or elements of a compound pixel.
> */
> typedef enum {UNKNOWNCOMPONENTTYPE,UCHAR,CHAR,USHORT,SHORT,UINT,INT,
> ULONG,LONG, FLOAT,DOUBLE} IOComponentType;
> </CODE>
>
>
> Can you add a line before the switch statement in itkPNGImageIO.cxx
>
> std::cout << "COMPONENT: " << this->GetComponentType() << std::endl;
>
> to see what component type it thinks it is?
>
> Besides, can you submit your code?
>
> Cheers,
>
>
> R.
>
> --
> Ramón Casero Cañas, DPhil
>
> Computational Biology, Computing Laboratory
> University of Oxford
> Wolfson Building, Parks Rd
> Oxford OX1 3QD
>
> tlf +44 (0) 1865 610807
> web http://web.comlab.ox.ac.uk/people/Ramon.CaseroCanas
> photos http://www.flickr.com/photos/rcasero/
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090722/8b5fec92/attachment-0001.htm>
More information about the Insight-users
mailing list