Commit f941d542 authored by Kasper Marstal's avatar Kasper Marstal Committed by GitHub

Merge pull request #369 from SuperElastix/SELX-351-Fix-ComposeDisplacementFieldComponent

ENH: Add support for displacement fields of different sizes and in di…
parents dbd019cd 87d26a8f
......@@ -344,13 +344,16 @@ def compose_displacement_fields(superelastix, disp_field_file_name_0, disp_field
output_image_file_name = output_image_base_name + '_inverse_consistency' + "." + output_image_ext
try:
stdout = subprocess.check_output([superelastix,
'--conf', os.path.join(os.path.dirname(os.path.realpath(__file__)), 'compose_displacement_fields.json'),
'--in', 'WarpingDisplacementField=%s' % disp_field_file_name_0,
'DisplacementField=%s' % disp_field_file_name_1,
'--out', 'DisplacementFieldSink=%s' % output_image_file_name,
'--loglevel', 'trace',
'--logfile', os.path.splitext(output_image_file_name)[0] + '.log'])
cmd = [superelastix,
'--conf', os.path.join(os.path.dirname(os.path.realpath(__file__)), 'compose_displacement_fields.json'),
'--in', 'WarpingDisplacementField=%s' % disp_field_file_name_0,
'DisplacementField=%s' % disp_field_file_name_1,
'--out', 'DisplacementFieldSink=%s' % output_image_file_name,
'--loglevel', 'trace',
'--logfile', os.path.splitext(output_image_file_name)[0] + '.log']
logging.debug(' '.join(cmd))
stdout = subprocess.check_output(cmd)
except Exception as e:
logging.error('Failed to compose blueprints %s and %s' % (disp_field_file_name_0, disp_field_file_name_1))
raise e
......
#ifndef selxComposeDisplacementFieldsImageFilter_h
#define selxComposeDisplacementFieldsImageFilter_h
#include "itkComposeDisplacementFieldsImageFilter.h"
#include "itkVectorInterpolateImageFunction.h"
namespace selx
{
/**
* \class ComposeDisplacementFieldsImageFilter
*
* \brief Compose two displacement fields.
*
* \author Nick Tustison
* \author Brian Avants
*
* \ingroup ITKDisplacementField
*
*/
template <typename TInputImage, typename TOutputImage = TInputImage>
class ComposeDisplacementFieldsImageFilter
: public itk::ComposeDisplacementFieldsImageFilter<TInputImage, TOutputImage>
{
public:
typedef ComposeDisplacementFieldsImageFilter Self;
typedef itk::ComposeDisplacementFieldsImageFilter<TInputImage, TOutputImage>
Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro( Self );
/** Extract dimension from input image. */
using Superclass::ImageDimension;
using typename Superclass::InputFieldType;
using typename Superclass::OutputFieldType;
/** Image typedef support. */
using typename Superclass::PixelType;
using typename Superclass::VectorType;
using typename Superclass::RegionType;
using typename Superclass::IndexType;
using typename Superclass::PointType;
using typename Superclass::SpacingType;
using typename Superclass::OriginType;
using typename Superclass::SizeType;
using typename Superclass::DirectionType;
/** Other typedef */
using typename Superclass::RealType;
using typename Superclass::InterpolatorType;
/** Get the interpolator. */
using Superclass::GetInterpolator;
/** Set the deformation field */
using Superclass::SetDisplacementField;
/**
* Get the deformation field.
*/
using Superclass::GetDisplacementField;
/** Set the warping field */
using Superclass::SetWarpingField;
/**
* Get the warping field.
*/
using Superclass::GetWarpingField;
/* Set the interpolator. */
using Superclass::SetInterpolator;
/* Override default behaviour (do not check that images have the same
* world coordinate systems). This is why we make this new class. The
* superclass should have done this from the beginning */
void VerifyInputInformation() override {};
/* This function implements the default behaviour of
* itk::ComposeDisplacementFieldsImageFilter::GenerateInputRequestedRegion()
* EXCEPT we let the warping field keep its own requested region (i.e. do not copy
* requested region from output to warping field) */
void GenerateInputRequestedRegion() override {
typedef itk::ImageBase< ImageDimension > ImageBaseType;
ImageBaseType* displacementField = dynamic_cast< ImageBaseType * >( const_cast< InputFieldType * >( this->GetInput( 0 ) ) );
if( displacementField )
{
typename ImageBaseType::RegionType inputRegion;
this->CallCopyOutputRegionToInputRegion( inputRegion, this->GetOutput()->GetRequestedRegion() );
displacementField->SetRequestedRegion(inputRegion);
}
}
/* Same as the base class except we crop the region to ensure it is contained in both output and
* warping field. Actually this is due to a bug in the Superclass. The output gets its domain from
* the displacement field so the superclass _should_ loop over the displacement field and the output
* (since these two domains are the same) NOT loop can be over the warping field and the output.
*/
void ThreadedGenerateData( const RegionType & region, itk::ThreadIdType itkNotUsed( threadId ) ) override
{
typename OutputFieldType::Pointer output = this->GetOutput();
typename InputFieldType::ConstPointer warpingField = this->GetWarpingField();
auto croppedRegion = const_cast< RegionType & >( region );
croppedRegion.Crop( warpingField->GetRequestedRegion() );
itk::ImageRegionConstIteratorWithIndex<InputFieldType> warpingFieldIterator( warpingField, croppedRegion );
itk::ImageRegionIterator<OutputFieldType> outputIterator( output, croppedRegion );
PointType pointIn1;
PointType pointIn2;
PointType pointIn3;
typename OutputFieldType::PixelType outDisplacement;
for( warpingFieldIterator.GoToBegin(), outputIterator.GoToBegin(); !warpingFieldIterator.IsAtEnd();
++warpingFieldIterator, ++outputIterator )
{
warpingField->TransformIndexToPhysicalPoint( warpingFieldIterator.GetIndex(), pointIn1 );
VectorType warpVector = warpingFieldIterator.Get();
for( unsigned int d = 0; d < ImageDimension; d++ )
{
pointIn2[d] = pointIn1[d] + warpVector[d];
}
typename InterpolatorType::OutputType displacement( 0.0 );
if( this->GetInterpolator()->IsInsideBuffer( pointIn2 ) )
{
displacement = this->GetInterpolator()->Evaluate( pointIn2 );
}
for( unsigned int d = 0; d < ImageDimension; d++ )
{
pointIn3[d] = pointIn2[d] + displacement[d];
}
outDisplacement = pointIn3 - pointIn1;
outputIterator.Set( outDisplacement );
}
}
protected:
/** Constructor */
ComposeDisplacementFieldsImageFilter() {};
/** Deconstructor */
~ComposeDisplacementFieldsImageFilter() {};
private:
ITK_DISALLOW_COPY_AND_ASSIGN(ComposeDisplacementFieldsImageFilter);
};
} // end namespace itk
#endif
......@@ -23,7 +23,7 @@
#include "selxSuperElastixComponent.h"
#include "selxItkObjectInterfaces.h"
#include "itkComposeDisplacementFieldsImageFilter.h"
#include "selxComposeDisplacementFieldsImageFilter.h"
namespace selx {
......@@ -55,7 +55,7 @@ public:
using WarpingDisplacementFieldType = typename WarpingDisplacementFieldInterfaceType::ItkDisplacementFieldType;
using WarpingDisplacementFieldPointer = typename WarpingDisplacementFieldType::Pointer;
using ComposeDisplacementFieldsImageFilterType = typename itk::ComposeDisplacementFieldsImageFilter< DisplacementFieldType, DisplacementFieldType >;
using ComposeDisplacementFieldsImageFilterType = ComposeDisplacementFieldsImageFilter< DisplacementFieldType, DisplacementFieldType >;
using ComposeDisplacementFieldsImageFilterPointer = typename ComposeDisplacementFieldsImageFilterType::Pointer;
// Accept interfaces
......
......@@ -55,9 +55,14 @@ void
ItkDisplacementFieldComposerComponent< Dimensionality, TPixel, CoordRepType >
::BeforeUpdate()
{
this->m_DisplacementField->UpdateOutputInformation();
this->m_DisplacementField->Update();
this->m_WarpingDisplacementField->Update();
// TODO: By default, the output region is copied to all inputs. Find out where it happens,
// and disable this behaviour
this->m_ComposeDisplacementFieldsImageFilter->SetDisplacementField( this->m_DisplacementField );
this->m_WarpingDisplacementField->UpdateOutputInformation();
this->m_ComposeDisplacementFieldsImageFilter->SetWarpingField( this->m_WarpingDisplacementField );
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment