[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://www.itk.org/pipermail/insight-users/attachments/20131204/656e657d/attachment.htm>
-------------- 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: 9271 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20131204/656e657d/attachment.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


More information about the Insight-users mailing list