Unverified Commit dcba82dc authored by N-Dekker's avatar N-Dekker Committed by GitHub
Browse files

Merge pull request #102 from SuperElastix/SELX-87-NiftyReg-Deformation-Field-Output

Selx 87 nifty reg deformation field output
parents ec470c9d 64c566dc
......@@ -50,6 +50,7 @@ set( ${MODULE}_SOURCE_FILES
set( ${MODULE}_TEST_SOURCE_FILES
${${MODULE}_SOURCE_DIR}/test/selxNiftyregComponentTest.cxx
${${MODULE}_SOURCE_DIR}/test/selxNiftiItkConversionsTest.cxx
)
set( ${MODULE}_LIBRARIES
......
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* 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 selxDisplacementFieldNiftiToItkImageSinkComponent_h
#define selxDisplacementFieldNiftiToItkImageSinkComponent_h
#include "selxSuperElastixComponent.h"
#include "selxInterfaces.h"
#include "selxNiftyregInterfaces.h"
#include "selxNiftiToItkImage.h"
#include "selxSinksAndSourcesInterfaces.h"
#include "selxItkObjectInterfaces.h"
#include <string.h>
#include "itkImageFileWriter.h"
#include "itkImportImageFilter.h"
#include "selxAnyFileWriter.h"
#include "selxFileWriterDecorator.h"
namespace selx
{
template< int Dimensionality, class TPixel >
class DisplacementFieldNiftiToItkImageSinkComponent :
public SuperElastixComponent<
Accepting< NiftyregDisplacementFieldImageInterface< TPixel >, itkImageDomainFixedInterface< Dimensionality >>,
Providing< SinkInterface, AfterRegistrationInterface >
>
{
public:
/** Standard ITK typedefs. */
typedef DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel > Self;
typedef SuperElastixComponent<
Accepting< NiftyregDisplacementFieldImageInterface< TPixel >, itkImageDomainFixedInterface< Dimensionality >>,
Providing< SinkInterface, AfterRegistrationInterface >
> Superclass;
typedef std::shared_ptr< Self > Pointer;
typedef std::shared_ptr< const Self > ConstPointer;
typedef itkImageDomainFixedInterface< Dimensionality > ItkImageDomainInterfaceType;
typedef std::shared_ptr< nifti_image > NiftiImagePointer;
typedef NiftyregDisplacementFieldImageInterface< TPixel > NiftiDisplacementFieldInterfaceType;
typedef itk::Image< itk::Vector< TPixel, Dimensionality >, Dimensionality > ItkDisplacementFieldType;
typedef typename ItkDisplacementFieldType::Pointer ItkDisplacementFieldPointer;
typedef itk::ImageFileWriter< ItkDisplacementFieldType > ItkDisplacementFieldWriterType;
typedef FileWriterDecorator< ItkDisplacementFieldWriterType > DecoratedWriterType;
typedef itk::ImportImageFilter< itk::Vector< TPixel, Dimensionality >, Dimensionality > ImportFilterType;
DisplacementFieldNiftiToItkImageSinkComponent( const std::string & name, LoggerImpl & logger );
virtual ~DisplacementFieldNiftiToItkImageSinkComponent();
// accepting interfaces
virtual int Set( typename NiftiDisplacementFieldInterfaceType::Pointer ) override;
virtual int Set( typename ItkImageDomainInterfaceType::Pointer ) override;
// providing interfaces
virtual void SetMiniPipelineOutput( itk::DataObject::Pointer ) override;
virtual itk::DataObject::Pointer GetMiniPipelineOutput( void ) override;
virtual AnyFileWriter::Pointer GetOutputFileWriter( void ) override;
virtual itk::DataObject::Pointer GetInitializedOutput( void ) override;
virtual void AfterRegistration() override;
virtual bool MeetsCriterion( const ComponentBase::CriterionType & criterion ) override;
static const char * GetDescription() { return "DisplacementFieldNiftiToItkImageSink Component"; }
private:
typename NiftiDisplacementFieldInterfaceType::Pointer m_DisplacementFieldInterface;
typename ItkImageDomainInterfaceType::Pointer m_ImageDomainInterface;
typename ImportFilterType::Pointer m_ImportFilter;
typename ItkDisplacementFieldType::Pointer m_MiniPipelineOutputImage;
typename ItkDisplacementFieldType::Pointer m_NetworkBuilderOutputImage;
//NiftiImagePointer m_MiniPipelineOutputImage;
protected:
// return the class name and the template arguments to uniquely identify this component.
static inline const std::map< std::string, std::string > TemplateProperties()
{
return { { keys::NameOfClass, "DisplacementFieldNiftiToItkImageSinkComponent" }, { keys::PixelType, PodString< TPixel >::Get() }, { keys::Dimensionality, std::to_string( Dimensionality ) } };
}
};
} //end namespace selx
#ifndef ITK_MANUAL_INSTANTIATION
#include "selxDisplacementFieldNiftiToItkImageSinkComponent.hxx"
#endif
#endif // #define selxDisplacementFieldNiftiToItkImageSinkComponent_h
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* 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 "selxDisplacementFieldNiftiToItkImageSinkComponent.h"
#include "selxCheckTemplateProperties.h"
namespace selx
{
template< int Dimensionality, class TPixel >
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::DisplacementFieldNiftiToItkImageSinkComponent( const std::string & name, LoggerImpl & logger ) :
Superclass( name, logger ), m_NetworkBuilderOutputImage( nullptr ), m_DisplacementFieldInterface( nullptr ), m_ImageDomainInterface( nullptr )
{
m_MiniPipelineOutputImage = ItkDisplacementFieldType::New(); // this is just an empty image for we have a SmartPointer we can pass around downstream. The actual image data will be grafted into this image.
m_ImportFilter = ImportFilterType::New();
}
template< int Dimensionality, class TPixel >
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::~DisplacementFieldNiftiToItkImageSinkComponent()
{
}
template< int Dimensionality, class TPixel >
int
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::Set( typename NiftiDisplacementFieldInterfaceType::Pointer other )
{
// Store pointer to the m_WarpedImageInterface for getting the result image after in has been generated (registration).
// TODO: sanity check that m_WarpedImageInterface was Null to detect if Set was called more than once erroneously.
this->m_DisplacementFieldInterface = other;
return 0;
}
template< int Dimensionality, class TPixel >
int
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::Set( typename ItkImageDomainInterfaceType::Pointer other )
{
// Store pointer to the m_ImageDomainInterface for getting the result image after in has been generated (registration).
// TODO: sanity check that m_ImageDomainInterface was Null to detect if Set was called more than once erroneously.
m_MiniPipelineOutputImage->SetRegions( other->GetItkImageDomainFixed()->GetLargestPossibleRegion() );
//this->m_ImageDomainInterface = other;
return 0;
}
template< int Dimensionality, class TPixel >
void
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::SetMiniPipelineOutput( itk::DataObject::Pointer NetworkBuilderOutput )
{
/** Tries to cast the NetworkBuilderOutput to an image (data object) and stores the result.
* The resulting output image will be grafted into when the sink component is connected to an other component.
* */
//
/*
this->m_NetworkBuilderOutputImage = dynamic_cast< ItkImageType * >(NetworkBuilderOutput.GetPointer());
if (this->m_NetworkBuilderOutputImage == nullptr)
{
throw std::runtime_error("DisplacementFieldNiftiToItkImageSinkComponent cannot cast the NetworkBuilder's Output to the required type");
}
this->m_MiniPipelineOutputImage = this->m_NetworkBuilderOutputImage;
*/
}
template< int Dimensionality, class TPixel >
typename itk::DataObject::Pointer
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::GetMiniPipelineOutput()
{
return this->m_MiniPipelineOutputImage.GetPointer();
}
template< int Dimensionality, class TPixel >
typename AnyFileWriter::Pointer
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::GetOutputFileWriter()
{
// Instanstiate an image file writer, decorated such that it can be implicitly cast to an AnyFileWriterType
return DecoratedWriterType::New().GetPointer();
}
template< int Dimensionality, class TPixel >
typename itk::DataObject::Pointer
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::GetInitializedOutput()
{
return ItkDisplacementFieldType::New().GetPointer();
}
template< int Dimensionality, class TPixel >
void
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::AfterRegistration()
{
auto displacementFieldNiftiImage = this->m_DisplacementFieldInterface->GetDisplacementFieldNiftiImage();
auto displacementFieldItkImage = NiftiToItkImage< ItkDisplacementFieldType, TPixel >::Convert( displacementFieldNiftiImage );
//auto displacementFieldItkImage = NiftiToItkImage< itk::Image< TPixel, Dimensionality >, TPixel >::Convert(displacementFieldNiftiImage);
this->m_MiniPipelineOutputImage->Graft( displacementFieldItkImage );
}
template< int Dimensionality, class TPixel >
bool
DisplacementFieldNiftiToItkImageSinkComponent< Dimensionality, TPixel >::MeetsCriterion( const ComponentBase::CriterionType & criterion )
{
bool hasUndefinedCriteria( false );
bool meetsCriteria( false );
auto status = CheckTemplateProperties( this->TemplateProperties(), criterion );
if( status == CriterionStatus::Satisfied )
{
return true;
}
else if( status == CriterionStatus::Failed )
{
return false;
} // else: CriterionStatus::Unknown
return meetsCriteria;
}
} //end namespace selx
......@@ -28,8 +28,37 @@
#include "itkImageIOBase.h"
#include "selxItkImageProperties.h"
// forward declaration of functions declared in ITK\Modules\IO\NIFTI\src\itkNiftiImageIO.cxx,
// since itk did not declare these in header files. Be aware that dynamic linking to itk might give issues.
namespace itk
{
int * UpperToLowerOrder(int dim);
int * LowerToUpperOrder(int dim);
int SymMatDim(int count);
}
namespace selx
{
// redirecting functions from itkNiftiImageIO. See forward declarations
inline int *
UpperToLowerOrder(int dim)
{
return ::itk::UpperToLowerOrder(dim);
}
inline int *
LowerToUpperOrder(int dim)
{
return ::itk::LowerToUpperOrder(dim);
}
inline int
SymMatDim(int count)
{
return ::itk::SymMatDim(count);
}
/** \class ItkToNiftiImage
* Convert an itk image to an nifti image object.
* Adapted from itkNiftiImageIO that is originally by Hans J. Johnson, The University of Iowa 2002
......
......@@ -23,104 +23,6 @@
namespace selx
{
// returns an ordering array for converting upper triangular symmetric matrix
// to lower triangular symmetric matrix
int *
UpperToLowerOrder( int dim )
{
int ** mat = new int *[ dim ];
for( int i = 0; i < dim; i++ )
{
mat[ i ] = new int[ dim ];
}
// fill in
int index( 0 );
for( int i = 0; i < dim; i++ )
{
for( int j = i; j < dim; j++ )
{
mat[ i ][ j ] = index;
mat[ j ][ i ] = index;
index++;
}
}
int * rval = new int[ index + 1 ];
int index2( 0 );
for( int i = 0; i < dim; i++ )
{
for( int j = 0; j <= i; j++, index2++ )
{
rval[ index2 ] = mat[ i ][ j ];
}
}
rval[ index2 ] = -1;
for( int i = 0; i < dim; i++ )
{
delete[] mat[ i ];
}
delete[] mat;
return rval;
}
// returns an ordering array for converting lower triangular symmetric matrix
// to upper triangular symmetric matrix
int *
LowerToUpperOrder( int dim )
{
int ** mat = new int *[ dim ];
for( int i = 0; i < dim; i++ )
{
mat[ i ] = new int[ dim ];
}
// fill in
int index( 0 );
for( int i = 0; i < dim; i++ )
{
for( int j = 0; j <= i; j++, index++ )
{
mat[ i ][ j ] = index;
mat[ j ][ i ] = index;
}
}
int * rval = new int[ index + 1 ];
int index2( 0 );
for( int i = 0; i < dim; i++ )
{
for( int j = i; j < dim; j++, index2++ )
{
rval[ index2 ] = mat[ i ][ j ];
}
}
rval[ index2 ] = -1;
for( int i = 0; i < dim; i++ )
{
delete[] mat[ i ];
}
delete[] mat;
return rval;
}
// compute the rank of the symmetric matrix from
// the count of the triangular matrix elements
int
SymMatDim( int count )
{
int dim = 0;
int row = 1;
while( count > 0 )
{
count -= row;
dim++;
row++;
}
return dim;
}
template< class ItkImageType, class NiftiPixelType >
std::shared_ptr< nifti_image >
......@@ -164,18 +66,15 @@ ItkToNiftiImage< ItkImageType, NiftiPixelType >::Convert( typename ItkImageType:
throw std::runtime_error( msg.str() );
}
typename ItkImageType::PixelContainerPointer pixelContainer = input->GetPixelContainer();
bool wasCopied = TransferImageData( pixelContainer->GetBufferPointer(), output );
if( !wasCopied )
{
// Take ownership of the pixelbuffer
pixelContainer->ContainerManageMemoryOff();
}
const bool wasCopied = TransferImageData(input->GetBufferPointer(), output );
std::shared_ptr< nifti_image > ptr( output, nifti_image_free );
return ptr;
// If the data was not copied to the nifti_image, it is shared between the ITK image and the nifti_image. Therefore, in that case,
// the ITK image is captured by the deleter of the shared_ptr, to ensure the lifetime of the ITK Image is extended to the end of
// the nifti_image lifetime. Note that in that case nifti_image_free should free all dynamically allocated memory of nifti_image
// except for its data, therefore its data pointer is in that case set to null.
return wasCopied ?
std::shared_ptr< nifti_image >(output, nifti_image_free):
std::shared_ptr< nifti_image >(output, [input](nifti_image* ptr) { ptr->data = nullptr; nifti_image_free(ptr); });
}
......
......@@ -39,6 +39,7 @@
#include "itkMetaDataObject.h"
#include "itkSpatialOrientationAdapter.h"
#include "itkImportImageFilter.h"
#include "itkShiftScaleImageFilter.h"
namespace selx
{
......@@ -86,7 +87,7 @@ NiftiToItkImage< ItkImageType, NiftiPixelType >
importImageFilter->SetSpacing( spacing );
importImageFilter->SetDirection( direction );
//importImageFilter->UpdateOutputInformation();
importImageFilter->SetImportPointer( dataFromNifti.buffer, dataFromNifti.numberOfElements, false );
importImageFilter->SetImportPointer( dataFromNifti.buffer, dataFromNifti.numberOfElements, true );
importImageFilter->UpdateOutputInformation();
auto resultImage = importImageFilter->GetOutput();
......@@ -122,6 +123,23 @@ RescaleFunction( TBuffer * buffer,
}
}
// Internal function to rescale pixel according to Rescale Slope/Intercept
template< typename TBuffer, unsigned int Dims >
void
RescaleFunction(itk::Vector<TBuffer, Dims> * buffer,
double slope,
double intercept,
size_t size)
{
for (size_t i = 0; i < size; i++)
{
itk::Vector<double, Dims> tmp;
tmp.CastFrom(buffer[i]);
tmp *= slope;
tmp += intercept;
buffer[i].CastFrom(tmp);
}
}
template< typename PixelType >
void
......@@ -168,7 +186,7 @@ NiftiToItkImage< ItkImageType, NiftiPixelType >
_size[i] = 1;
}
*/
unsigned int numComponents = ItkImageProperties< ItkImageType >::GetNumberOfComponents();
const unsigned int numComponents = ItkImageProperties< ItkImageType >::GetNumberOfComponents();
//
// special case for images of vector pixels
/*
......@@ -243,7 +261,6 @@ NiftiToItkImage< ItkImageType, NiftiPixelType >
data = _data;
}
//void* buffer = (void*) new char[NumBytes];
typename ItkImageType::PixelType * buffer = new typename ItkImageType::PixelType[ imageSizeInComponents ];
//
// if single or complex, nifti layout == itk layout
......@@ -310,6 +327,7 @@ NiftiToItkImage< ItkImageType, NiftiPixelType >
if( data != input_image->data )
{
free( data );
data = NULL;
}
//dumpdata(data);
//dumpdata(buffer);
......
......@@ -60,6 +60,26 @@ public:
virtual std::shared_ptr< nifti_image > GetWarpedNiftiImage() = 0;
};
template< class TPixel >
class NiftyregControlPointPositionImageInterface
{
public:
using Type = NiftyregControlPointPositionImageInterface< TPixel >;
using Pointer = std::shared_ptr< Type >;
virtual std::shared_ptr< nifti_image > GetControlPointPositionImage() = 0;
};
template< class TPixel >
class NiftyregDisplacementFieldImageInterface
{
public:
using Type = NiftyregDisplacementFieldImageInterface< TPixel >;
using Pointer = std::shared_ptr< Type >;
virtual std::shared_ptr< nifti_image > GetDisplacementFieldNiftiImage() = 0;
};
template< class TPixel >
struct Properties< NiftyregReferenceImageInterface< TPixel >>
{
......@@ -86,6 +106,24 @@ struct Properties< NiftyregWarpedImageInterface< TPixel >>
return { { keys::NameOfInterface, "NiftyregWarpedImageInterface" }, { keys::PixelType, PodString< TPixel >::Get() } };
}
};
template< class TPixel >
struct Properties< NiftyregControlPointPositionImageInterface< TPixel >>
{
static const std::map< std::string, std::string > Get()
{
return{ { keys::NameOfInterface, "NiftyregControlPointPositionImageInterface" }, { keys::PixelType, PodString< TPixel >::Get() } };
}
};
template< class TPixel >
struct Properties< NiftyregDisplacementFieldImageInterface< TPixel >>
{
static const std::map< std::string, std::string > Get()
{
return{ { keys::NameOfInterface, "NiftyregDisplacementFieldImageInterface" }, { keys::PixelType, PodString< TPixel >::Get() } };
}
};
}
#endif // #define selxNiftyregInterfaces_h
/*=========================================================================
*
* Copyright Leiden University Medical Center, Erasmus University Medical
* Center and contributors
*
* 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 selxNiftyregSplineToDisplacementFieldComponent_h
#define selxNiftyregSplineToDisplacementFieldComponent_h
#include "selxSuperElastixComponent.h"
#include "selxInterfaces.h"
#include "selxNiftyregInterfaces.h"
#include <string.h>
#include <array>
namespace selx
{
template< class TPixel >
class NiftyregSplineToDisplacementFieldComponent :
public SuperElastixComponent<
Accepting< NiftyregControlPointPositionImageInterface< TPixel >, NiftyregReferenceImageInterface< TPixel > >,
Providing< NiftyregDisplacementFieldImageInterface< TPixel >, ReconnectTransformInterface >
>
{