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