[ITK Community] [Insight-users] issues with inverting deformation fields
Tim Bhatnagar
tim.bhatnagar at gmail.com
Wed Dec 4 13:39:29 EST 2013
Hello all,
Sorry to make a repeat request for this issue.. I have made a project file
for Visual C++ Express 2010 using the attached cmakelists and
InvertingDefFields.cxx files, which are linked to the ITK libraries.
The cxx file (and its corresponding hxx file) was downloaded from the
github recently, but I am not able to build my solution without getting the
following error:
InvertingDefFields.cxx
2>c:\itkrap\insighttoolkit-4.2.0\modules\filtering\displacementfield\include\itkInverseDisplacementFieldImageFilter.hxx(353):
error C2143: syntax error : missing ';' before
'itk::InverseDisplacementFieldImageFilter<TInputImage,TOutputImage>::GetMTime'
2>c:\itkrap\insighttoolkit-4.2.0\modules\filtering\displacementfield\include\itkInverseDisplacementFieldImageFilter.hxx(353):
error C4430: missing type specifier - int assumed. Note: C++ does not
support default-int
I am, by no means, an apt c++ programmer, and I have tried everything I can
to address this error, but to no avail. Is there anything that anyone can
suggest? I'll try anything.. this is most likely the last applictiton I
need before my dissertation can be written up (hopefully!), and this
problem has been dragging on for a while.
I've attached the related files - any help is much, much appreciated. OR,
if you have access to a different set of files to acquire an the inverse of
a displacement field generated by Demon's Registration, I'd be happy to try
that too!
Thanks so much,
--
Tim Bhatnagar
PhD Candidate
Orthopaedic Injury Biomechanics Group
Department of Mechanical Engineering
University of British Columbia
Rm 5000 - 818 West 10th Ave.
Vancouver, BC
Canada
V5Z 1M9
Ph: (604) 675-8845
Fax: (604) 675-8820
Web: oibg.mech.ubc.ca
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20131204/656e657d/attachment.html>
-------------- next part --------------
PROJECT( InvertingDefFields )
FIND_PACKAGE ( ITK )
IF ( ITK_FOUND )
INCLUDE( ${ITK_USE_FILE} )
ENDIF( ITK_FOUND )
ADD_EXECUTABLE( InvertingDefFields InvertingDefFields.cxx )
TARGET_LINK_LIBRARIES ( InvertingDefFields ITKCommon ITKIO ITKRegistration)
-------------- next part --------------
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "itkInverseDisplacementFieldImageFilter.h"
#include "itkImageFileWriter.h"
#include "itkFilterWatcher.h"
int itkInverseDisplacementFieldImageFilterTest( int argc, char * argv[] )
{
if( argc < 2 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " outputImage" << std::endl;
return EXIT_FAILURE;
}
const unsigned int Dimension = 2;
typedef float VectorComponentType;
typedef itk::Vector< VectorComponentType, Dimension > VectorType;
typedef itk::Image< VectorType, Dimension > DisplacementFieldType;
typedef itk::InverseDisplacementFieldImageFilter<
DisplacementFieldType,
DisplacementFieldType
> FilterType;
FilterType::Pointer filter = FilterType::New();
FilterWatcher watcher(filter);
// Creating an input displacement field
DisplacementFieldType::Pointer field = DisplacementFieldType::New();
DisplacementFieldType::SpacingType spacing;
spacing.Fill( 1.0 );
DisplacementFieldType::PointType origin;
origin.Fill( 0.0 );
DisplacementFieldType::RegionType region;
DisplacementFieldType::SizeType size;
DisplacementFieldType::IndexType start;
size[0] = 128;
size[1] = 128;
start[0] = 0;
start[1] = 0;
region.SetSize( size );
region.SetIndex( start );
field->SetOrigin( origin );
field->SetSpacing( spacing );
field->SetRegions( region );
field->Allocate();
VectorType pixelValue;
itk::ImageRegionIteratorWithIndex< DisplacementFieldType > it( field, region );
// Fill the field with some vectors
it.GoToBegin();
while( !it.IsAtEnd() )
{
DisplacementFieldType::IndexType index = it.GetIndex();
pixelValue[0] = index[0] * 2.0 - index[0];
pixelValue[1] = index[1] * 2.0 - index[1];
it.Set( pixelValue );
++it;
}
// Since the tested transform is upsampling by a factor of two, the
// size of the inverse field should be twice the size of the input
// field. All other geomtry parameters are the same.
filter->SetOutputSpacing( spacing );
// keep the origin
filter->SetOutputOrigin( origin );
// set the size
DisplacementFieldType::SizeType invFieldSize;
invFieldSize[0] = size[0] * 2;
invFieldSize[1] = size[1] * 2;
filter->SetSize( invFieldSize );
filter->SetInput( field );
filter->SetSubsamplingFactor( 16 );
try
{
filter->UpdateLargestPossibleRegion();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
}
// Write an image for regression testing
typedef itk::ImageFileWriter< DisplacementFieldType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetInput (filter->GetOutput() );
writer->SetFileName( argv[1] );
try
{
writer->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception thrown by writer" << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Now, test for loop invariant (acts as filter validation)
// f^-1(f(p1) + p1 ) - f(p1) = 0
it.GoToBegin();
while( !it.IsAtEnd() )
{
DisplacementFieldType::PointType p1;
field->TransformIndexToPhysicalPoint(it.GetIndex(), p1);
DisplacementFieldType::PixelType fp1 = it.Get();
DisplacementFieldType::PointType p2;
p2[0] = p1[0] + fp1[0];
p2[1] = p1[1] + fp1[1];
DisplacementFieldType::IndexType id2;
filter->GetOutput()->TransformPhysicalPointToIndex(p2,id2);
DisplacementFieldType::PixelType fp2 = filter->GetOutput()->GetPixel(id2);
if(vcl_abs(fp2[0] + fp1[0]) > 0.001
|| vcl_abs(fp2[1] + fp1[1]) > 0.001)
{
std::cerr<<"Loop invariant not satisfied for index "<<it.GetIndex()<<" : f^-1(f(p1) + p1 ) + f(p1) = 0"<< std::endl;
std::cerr<<"f(p1) = "<<fp1<<std::endl;
std::cerr<<"f^-1(f(p1) + p1 ) = "<<fp2<<std::endl;
std::cerr<<"diff: "<<fp1[0]+fp2[0]<<", "<<fp1[1]+fp2[1]<<std::endl;
return EXIT_FAILURE;
}
++it;
}
return EXIT_SUCCESS;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: itkInverseDisplacementFieldImageFilter.h
Type: text/x-chdr
Size: 9272 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/community/attachments/20131204/656e657d/attachment-0002.h>
-------------- next part --------------
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef __itkInverseDisplacementFieldImageFilter_hxx
#define __itkInverseDisplacementFieldImageFilter_hxx
#include "itkInverseDisplacementFieldImageFilter.h"
#include "itkObjectFactory.h"
#include "itkProgressReporter.h"
#include "itkThinPlateSplineKernelTransform.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkVectorResampleImageFilter.h"
#include "itkImageIOBase.h"
#include "itkImageIOFactory.h"
#include "itkCommand.h"
#include "itkTimeStamp.h"
#include "itkProcessObject.h"
#include "itkImageAdaptor.h"
namespace itk
{
/**
* Initialize new instance
*/
template< class TInputImage, class TOutputImage >
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::InverseDisplacementFieldImageFilter()
{
m_OutputSpacing.Fill(1.0);
m_OutputOrigin.Fill(0.0);
for ( unsigned int i = 0; i < ImageDimension; i++ )
{
m_Size[i] = 0;
}
typedef ThinPlateSplineKernelTransform<
double,
itkGetStaticConstMacro(ImageDimension) > DefaultTransformType;
m_KernelTransform = DefaultTransformType::New();
m_SubsamplingFactor = 16;
}
/**
* Print out a description of self
*
* \todo Add details about this class
*/
template< class TInputImage, class TOutputImage >
void
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::PrintSelf(std::ostream & os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "Size: " << m_Size << std::endl;
os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl;
os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl;
os << indent << "KernelTransform: " << m_KernelTransform.GetPointer() << std::endl;
os << indent << "SubsamplingFactor: " << m_SubsamplingFactor << std::endl;
return;
}
/**
* Set the output image spacing.
*/
template< class TInputImage, class TOutputImage >
void
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::SetOutputSpacing(const double *spacing)
{
SpacingType s(spacing);
this->SetOutputSpacing(s);
}
/**
* Set the output image origin.
*/
template< class TInputImage, class TOutputImage >
void
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::SetOutputOrigin(const double *origin)
{
OriginPointType p(origin);
this->SetOutputOrigin(p);
}
/**
* Sub-sample the input displacement field and prepare the KernelBase
* BSpline
*/
template< class TInputImage, class TOutputImage >
void
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::PrepareKernelBaseSpline()
{
typedef typename KernelTransformType::PointsContainer LandmarkContainer;
typedef typename LandmarkContainer::Pointer LandmarkContainerPointer;
// Source contains points with physical coordinates of the
// destination displacement fields (the inverse field)
LandmarkContainerPointer source = LandmarkContainer::New();
// Target contains vectors (stored as points) indicating
// displacement in the inverse direction.
LandmarkContainerPointer target = LandmarkContainer::New();
typedef itk::VectorResampleImageFilter<
InputImageType,
InputImageType > ResamplerType;
typename ResamplerType::Pointer resampler = ResamplerType::New();
const InputImageType *inputImage = this->GetInput();
resampler->SetInput(inputImage);
resampler->SetOutputOrigin( inputImage->GetOrigin() );
typename InputImageType::SpacingType spacing = inputImage->GetSpacing();
typedef typename InputImageType::RegionType InputRegionType;
typedef typename InputImageType::SizeType InputSizeType;
typedef typename InputImageType::IndexType InputIndexType;
InputRegionType region;
region = inputImage->GetLargestPossibleRegion();
InputSizeType size = region.GetSize();
for ( unsigned int i = 0; i < ImageDimension; i++ )
{
size[i] = static_cast< SizeValueType >( size[i] / m_SubsamplingFactor );
spacing[i] *= m_SubsamplingFactor;
}
InputRegionType subsampledRegion;
subsampledRegion.SetSize(size);
subsampledRegion.SetIndex( region.GetIndex() );
resampler->SetSize(size);
resampler->SetOutputStartIndex( subsampledRegion.GetIndex() );
resampler->SetOutputSpacing(spacing);
resampler->Update();
// allocate a landmark pair for each
// pixel in the subsampled field
const SizeValueType numberOfLandmarks = subsampledRegion.GetNumberOfPixels();
source->Reserve(numberOfLandmarks);
target->Reserve(numberOfLandmarks);
const InputImageType *sampledInput = resampler->GetOutput();
typedef ImageRegionConstIteratorWithIndex< InputImageType > IteratorType;
unsigned int landmarkId = 0;
IteratorType ot(sampledInput, subsampledRegion);
ot.GoToBegin();
OutputPixelType value;
Point< double, ImageDimension > sourcePoint;
Point< double, ImageDimension > targetPoint;
while ( !ot.IsAtEnd() )
{
value = ot.Get();
// Here we try to evaluate the inverse transform, so points from
// input displacement field are actually the target points
sampledInput->TransformIndexToPhysicalPoint(ot.GetIndex(), targetPoint);
target->InsertElement(landmarkId, targetPoint);
for ( unsigned int i = 0; i < ImageDimension; i++ )
{
sourcePoint[i] = targetPoint[i] + value[i];
}
source->InsertElement(landmarkId, sourcePoint); // revert direction of
// displacement
++landmarkId;
++ot;
}
itkDebugMacro(<< "Number of Landmarks created = " << numberOfLandmarks);
m_KernelTransform->GetModifiableTargetLandmarks()->SetPoints(target);
m_KernelTransform->GetModifiableSourceLandmarks()->SetPoints(source);
itkDebugMacro(<< "Before ComputeWMatrix() ");
m_KernelTransform->ComputeWMatrix();
itkDebugMacro(<< "After ComputeWMatrix() ");
}
/**
* GenerateData
*/
template< class TInputImage, class TOutputImage >
void
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::GenerateData()
{
// First subsample the input displacement field in order to create
// the KernelBased spline.
this->PrepareKernelBaseSpline();
itkDebugMacro(<< "Actually executing");
// Get the output pointers
OutputImageType *outputPtr = this->GetOutput();
outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
outputPtr->Allocate();
// Create an iterator that will walk the output region for this thread.
typedef ImageRegionIteratorWithIndex<
TOutputImage > OutputIterator;
OutputImageRegionType region = outputPtr->GetRequestedRegion();
OutputIterator outIt(outputPtr, region);
// Define a few indices that will be used to translate from an input pixel
// to an output pixel
IndexType outputIndex; // Index to current output pixel
typedef typename KernelTransformType::InputPointType InputPointType;
typedef typename KernelTransformType::OutputPointType OutputPointType;
InputPointType outputPoint; // Coordinates of current output pixel
// Support for progress methods/callbacks
ProgressReporter progress(this, 0, region.GetNumberOfPixels(), 10);
outIt.GoToBegin();
// Walk the output region
while ( !outIt.IsAtEnd() )
{
// Determine the index of the current output pixel
outputIndex = outIt.GetIndex();
outputPtr->TransformIndexToPhysicalPoint(outputIndex, outputPoint);
// Compute corresponding inverse displacement vector
OutputPointType interpolation =
m_KernelTransform->TransformPoint(outputPoint);
OutputPixelType inverseDisplacement;
for ( unsigned int i = 0; i < ImageDimension; i++ )
{
inverseDisplacement[i] = interpolation[i] - outputPoint[i];
}
outIt.Set(inverseDisplacement); // set inverse displacement.
++outIt;
progress.CompletedPixel();
}
return;
}
/**
* Inform pipeline of necessary input image region
*
* Determining the actual input region is non-trivial, especially
* when we cannot assume anything about the transform being used.
* So we do the easy thing and request the entire input image.
*/
template< class TInputImage, class TOutputImage >
void
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::GenerateInputRequestedRegion()
{
// call the superclass's implementation of this method
Superclass::GenerateInputRequestedRegion();
if ( !this->GetInput() )
{
return;
}
// get pointers to the input and output
InputImagePointer inputPtr =
const_cast< InputImageType * >( this->GetInput() );
// Request the entire input image
InputImageRegionType inputRegion;
inputRegion = inputPtr->GetLargestPossibleRegion();
inputPtr->SetRequestedRegion(inputRegion);
return;
}
/**
* Inform pipeline of required output region
*/
template< class TInputImage, class TOutputImage >
void
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::GenerateOutputInformation()
{
// call the superclass' implementation of this method
Superclass::GenerateOutputInformation();
// get pointers to the input and output
OutputImagePointer outputPtr = this->GetOutput();
if ( !outputPtr )
{
return;
}
// Set the size of the output region
typename TOutputImage::RegionType outputLargestPossibleRegion;
outputLargestPossibleRegion.SetSize(m_Size);
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
// Set spacing and origin
outputPtr->SetSpacing(m_OutputSpacing);
outputPtr->SetOrigin(m_OutputOrigin);
return;
}
/**
* Verify if any of the components has been modified.
*/
template< class TInputImage, class TOutputImage >
ModifiedTimeType
InverseDisplacementFieldImageFilter< TInputImage, TOutputImage >
::GetMTime( void ) const
{
ModifiedTimeType latestTime = Object::GetMTime();
if ( m_KernelTransform )
{
if ( latestTime < m_KernelTransform->GetMTime() )
{
latestTime = m_KernelTransform->GetMTime();
}
}
return latestTime;
}
} // end namespace itk
#endif
-------------- next part --------------
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users
More information about the Community
mailing list