Commit 60bd6976 authored by francois_budin's avatar francois_budin
Browse files

BUG:measurement frame with HField transform was not taken into account;...

BUG:measurement frame with HField transform was not taken into account; problem with transforms when multithreading;improvement of interpolators implementation

git-svn-id: http://svn.slicer.org/Slicer4/trunk@13073 3bd1e089-480b-0410-8dfb-8563597acbee
parent 79c105c6
......@@ -144,15 +144,13 @@ This work is part of the National Alliance for Medical Image Computing (NAMIC),
<element>a</element>
<default>a</default>
</string-enumeration>
<string-enumeration>
<boolean>
<name>space</name>
<longflag>--space</longflag>
<description>Space Orientation: RAS/LPS (warning: if the transform is a Transform Node in Slicer3, use LPS)</description>
<label>Transform Space Orientation</label>
<element>RAS</element>
<element>LPS</element>
<default>LPS</default>
</string-enumeration>
<longflag>--spaceChange</longflag>
<description>Space Orientation between transform and image is different (RAS/LPS) (warning: if the transform is a Transform Node in Slicer3, do not select)</description>
<label>Space Orientation inconsistency (between transform and image)</label>
<default>false</default>
</boolean>
</parameters>
<parameters advanced="true">
<label>Rigid/Affine Parameters</label>
......@@ -172,6 +170,15 @@ This work is part of the National Alliance for Medical Image Computing (NAMIC),
<label>Centered Transform</label>
<default>false</default>
</boolean>
<string-enumeration>
<name>imageCenter</name>
<longflag>--image_center</longflag>
<description>Image to use to center the transform (used only if "Centered Transform" is selected)</description>
<label>Image Center</label>
<element>input</element>
<element>output</element>
<default>input</default>
</string-enumeration>
<boolean>
<name>inverseITKTransformation</name>
<flag>-b</flag>
......@@ -183,14 +190,17 @@ This work is part of the National Alliance for Medical Image Computing (NAMIC),
</parameters>
<parameters advanced="true">
<label>Affine Transform Parameters</label>
<boolean>
<label>Tensor Transform Type</label>
<string-enumeration>
<name>ppd</name>
<longflag>--ppd</longflag>
<description>Uses the Preservation of the Principal Direction to compute the affine transform</description>
<label>Preservation of Principal Direction</label>
<default>false</default>
</boolean>
<longflag>--transform_tensor_method</longflag>
<flag>-T</flag>
<description>Chooses between 2 methods to transform the tensors: Finite Strain (FS), faster but less accurate, or Preservation of the Principal Direction (PPD)</description>
<label>Finite Strain (FS) or \nPreservation of the Principal Direction (PPD)</label>
<element>PPD</element>
<element>FS</element>
<default>PPD</default>
</string-enumeration>
</parameters>
......
......@@ -35,7 +35,7 @@ DiffusionTensor3DAffineTransform< TData >
}*/
this->m_TransformMatrix = transform->GetMatrix() ;
this->m_Translation = transform->GetTranslation() ;
this->Modified() ;
}
template< class TData >
......@@ -56,13 +56,14 @@ DiffusionTensor3DAffineTransform< TData >
::SetMatrix4x4( MatrixTransform4x4Type matrix )
{
for( int i = 0 ; i < 3 ; i++ )
{
{
for( int j = 0 ; j < 3 ; j++ )
{
{
this->m_TransformMatrix[ i ][ j ] = matrix[ i ][ j ] ;
}
this->m_Translation[ i ] = matrix[ i ][ 3 ] ;
}
this->m_Translation[ i ] = matrix[ i ][ 3 ] ;
}
this->Modified() ;
}
......
......@@ -25,18 +25,18 @@ namespace itk
*
* Implementation of blockwise spline interpolation for diffusion tensor images
*/
template< class TData >
template< class TData , class TCoordRep = double >
class DiffusionTensor3DBSplineInterpolateImageFunction :
public DiffusionTensor3DInterpolateImageFunctionReimplementation< TData >
public DiffusionTensor3DInterpolateImageFunctionReimplementation< TData , TCoordRep >
{
public:
typedef TData DataType ;
typedef DiffusionTensor3DBSplineInterpolateImageFunction Self ;
typedef DiffusionTensor3DInterpolateImageFunctionReimplementation< DataType > Superclass ;
typedef DiffusionTensor3DInterpolateImageFunctionReimplementation< DataType , TCoordRep > Superclass ;
typedef typename Superclass::ImageType ImageType ;
typedef SmartPointer< Self > Pointer ;
typedef SmartPointer< const Self > ConstPointer ;
typedef BSplineInterpolateImageFunction< ImageType , double , double >
typedef BSplineInterpolateImageFunction< ImageType , TCoordRep , double >
BSplineInterpolateFunction ;
itkNewMacro( Self ) ;
......
......@@ -19,16 +19,16 @@
namespace itk
{
template< class TData >
DiffusionTensor3DBSplineInterpolateImageFunction< TData >
template< class TData , class TCoordRep >
DiffusionTensor3DBSplineInterpolateImageFunction< TData , TCoordRep >
::DiffusionTensor3DBSplineInterpolateImageFunction()
{
m_SplineOrder = 1 ;
}
template< class TData >
template< class TData , class TCoordRep >
void
DiffusionTensor3DBSplineInterpolateImageFunction< TData >
DiffusionTensor3DBSplineInterpolateImageFunction< TData , TCoordRep >
::AllocateInterpolator()
{
for( int i = 0 ; i < 6 ; i++ )
......
......@@ -40,7 +40,7 @@ public :
///Get a Symmetric matrix representing the tensor
MatrixType GetTensor2Matrix() ;
///Set the Tensor from a symmetric matrix
void SetTensorFromMatrix( MatrixType matrix ) ;
template< class C > void SetTensorFromMatrix( Matrix< C , 3 , 3 > matrix ) ;
///Cast the component values of the tensor
template< class C > operator DiffusionTensor3DExtended< C > const() ;
......
......@@ -47,15 +47,16 @@ DiffusionTensor3DExtended< T >
}
template< class T >
template< class C >
void
DiffusionTensor3DExtended< T >
::SetTensorFromMatrix( MatrixType matrix )
::SetTensorFromMatrix( Matrix< C , 3 , 3 > matrix )
{
for( int i = 0 ; i < 3 ; i++ )
{
for( int j = i ; j < 3 ; j++ )
{
( *this )( i , j ) = matrix[ i ][ j ] ;
( *this )( i , j ) = static_cast< T >( matrix[ i ][ j ] ) ;
}
}
}
......
......@@ -56,15 +56,18 @@ DiffusionTensor3DFSAffineTransform< TData >
::PreCompute()
{
InternalMatrixTransformType m_RotationMatrix ;
this->latestTime = Object::GetMTime() ;
m_RotationMatrix=ComputeRotationMatrixFromTransformationMatrix() ;
m_RotationMatrix = ComputeRotationMatrixFromTransformationMatrix() ;
InternalMatrixTransformType MeasurementFrameTranspose = this->m_MeasurementFrame.GetTranspose() ;
InternalMatrixTransformType RotationMatrixTranspose = m_RotationMatrix.GetTranspose() ;
//Instead of computing the inverse transform (which go from the input image to the output image)
//and compute the rotation matrix from the inverse, we compute the rotation matrix from the original
//transform and takes its inverse (=transpose since it is a rotation). Both rotation matrices are equals,
//therefore we avoid computing the inverse of our original transform since it is not necessary.
this->m_TransformT = MeasurementFrameTranspose * m_RotationMatrix ;
this->m_Transform = RotationMatrixTranspose * this->m_MeasurementFrame ;
this->ComputeOffset() ;
this->latestTime = Object::GetMTime() ;
}
......
......@@ -19,9 +19,9 @@
#include <itkOrientedImage.h>
#include <itkPoint.h>
//#include <itkSemaphore.h>
#include <itkNumericTraits.h>
#include "define.h"
//#include <itkNumericTraits.h>
//#include "define.h"
#include <itkImageFunction.h>
namespace itk
{
......@@ -31,8 +31,12 @@ namespace itk
* Virtual class to implement diffusion tensor interpolation classes
*
*/
template< class TData >
class DiffusionTensor3DInterpolateImageFunction : public Object
template< class TData , class TCoordRep = double >
class DiffusionTensor3DInterpolateImageFunction :
public ImageFunction< OrientedImage< DiffusionTensor3D < TData > , 3 > ,
DiffusionTensor3D < TData > ,
TCoordRep
>
{
public :
typedef TData TensorType ;
......@@ -44,19 +48,60 @@ public :
typedef SmartPointer< Self > Pointer ;
typedef SmartPointer< const Self > ConstPointer ;
typedef typename TensorDataType::RealValueType TensorRealType ;
///Set the input image
itkSetObjectMacro( InputImage , DiffusionImageType ) ;
///Evaluate the tensor value at a given position
virtual TensorDataType Evaluate( const PointType &point ) = 0 ;
void SetDefaultPixelValue( TensorRealType defaultPixelValue ) ;
itkGetMacro( DefaultPixelValue , TensorRealType ) ;
typedef ImageFunction< OrientedImage< DiffusionTensor3D < TData > , 3 > ,
DiffusionTensor3D < TData > ,
TCoordRep
> Superclass ;
typedef typename Superclass::ContinuousIndexType ContinuousIndexType ;
typedef typename Superclass::IndexType IndexType ;
/////Copied from itkInterpolateImageFunction.h
/** Interpolate the image at a point position
*
* Returns the interpolated image intensity at a
* specified point position. No bounds checking is done.
* The point is assume to lie within the image buffer.
*
* ImageFunction::IsInsideBuffer() can be used to check bounds before
* calling the method. */
virtual TensorDataType Evaluate( const PointType& point ) const
{
ContinuousIndexType index;
this->GetInputImage()->TransformPhysicalPointToContinuousIndex( point, index ) ;
return ( this->EvaluateAtContinuousIndex( index ) ) ;
}
/** Interpolate the image at a continuous index position
*
* Returns the interpolated image intensity at a
* specified index position. No bounds checking is done.
* The point is assume to lie within the image buffer.
*
* Subclasses must override this method.
*
* ImageFunction::IsInsideBuffer() can be used to check bounds before
* calling the method. */
virtual TensorDataType EvaluateAtContinuousIndex( const ContinuousIndexType & index ) const = 0 ;
/** Interpolate the image at an index position.
*
* Simply returns the image value at the
* specified index position. No bounds checking is done.
* The point is assume to lie within the image buffer.
*
* ImageFunction::IsInsideBuffer() can be used to check bounds before
* calling the method. */
virtual TensorDataType EvaluateAtIndex( const IndexType & index ) const
{
return this->GetInputImage()->GetPixel( index ) ;
}
// void SetDefaultPixelValue( TensorRealType defaultPixelValue ) ;
// itkGetMacro( DefaultPixelValue , TensorRealType ) ;
protected:
DiffusionTensor3DInterpolateImageFunction() ;
DiffusionImageTypePointer m_InputImage ;
unsigned long latestTime ;
TensorRealType m_DefaultPixelValue ;
TensorDataType m_DefaultPixel ;
// TensorRealType m_DefaultPixelValue ;
// TensorDataType m_DefaultPixel ;
};
}//end namespace itk
......
......@@ -19,15 +19,20 @@
namespace itk
{
template< class TData >
DiffusionTensor3DInterpolateImageFunction< TData >
template< class TData , class TCoordRep >
DiffusionTensor3DInterpolateImageFunction< TData , TCoordRep >
::DiffusionTensor3DInterpolateImageFunction()
{
m_InputImage = 0 ;
// m_InputImage = 0 ;
latestTime = 0 ;
SetDefaultPixelValue( ZERO ) ;
// SetDefaultPixelValue( ZERO ) ;
}
/*
template< class TData >
void
DiffusionTensor3DInterpolateImageFunction< TData >
......@@ -40,7 +45,7 @@ DiffusionTensor3DInterpolateImageFunction< TData >
m_DefaultPixel( i , i ) *= static_cast< TData >( this->m_DefaultPixelValue ) ;
}
}
*/
}//end namespace itk
......
......@@ -16,21 +16,22 @@
#include "itkDiffusionTensor3DInterpolateImageFunction.h"
#include <itkOrientedImage.h>
#include <itkImageRegionIteratorWithIndex.h>
//#include <itkImageRegionIteratorWithIndex.h>
#include <itkInterpolateImageFunction.h>
#include <itkMutexLock.h>
#include <itkSemaphore.h>
#include "itkSeparateComponentsOfADiffusionTensorImage.h"
//#include <itkMutexLock.h>
//#include <itkSemaphore.h>
namespace itk
{
struct RegionType
/*struct RegionType
{
ImageRegion< 3 > itkRegion ;
bool Done ;
Index< 3 > PositionInImage ;
bool Stop ;
};
};*/
/**
* \class DiffusionTensor3DInterpolateImageFunctionReimplementation
......@@ -38,14 +39,14 @@ struct RegionType
* Abstract class allowing to implement blockwise interpolation for diffusion tensor images
*/
template< class TData >
template< class TData , class TCoordRep = double >
class DiffusionTensor3DInterpolateImageFunctionReimplementation :
public DiffusionTensor3DInterpolateImageFunction< TData >
public DiffusionTensor3DInterpolateImageFunction< TData , TCoordRep >
{
public :
typedef TData DataType ;
typedef DiffusionTensor3DInterpolateImageFunctionReimplementation Self ;
typedef DiffusionTensor3DInterpolateImageFunction< DataType > Superclass ;
typedef DiffusionTensor3DInterpolateImageFunction< DataType , TCoordRep > Superclass ;
typedef typename Superclass::TensorDataType TensorDataType ;
typedef typename Superclass::DiffusionImageType DiffusionImageType ;
typedef typename Superclass::DiffusionImageTypePointer DiffusionImageTypePointer ;
......@@ -59,30 +60,35 @@ public :
typedef InterpolateImageFunction< ImageType , double > InterpolateImageFunctionType ;
typedef typename DiffusionImageType::RegionType itkRegionType ;
typedef typename DiffusionImageType::SizeType SizeType ;
typedef typename Superclass::ContinuousIndexType ContinuousIndexType ;
/** Evaluate the interpolated tensor at a position
*/
TensorDataType Evaluate( const PointType &point ) ;
virtual void SetInputImage( DiffusionImageTypePointer inputImage ) ;
//TensorDataType Evaluate( const PointType &point ) ;
TensorDataType EvaluateAtContinuousIndex( const ContinuousIndexType & index ) const ;
virtual void SetInputImage( const DiffusionImageType *inputImage ) ;
itkSetMacro( NumberOfThreads , int ) ;
protected:
DiffusionTensor3DInterpolateImageFunctionReimplementation() ;
virtual void AllocateInterpolator() = 0 ;
void SeparateImages() ;
void AllocateImages() ;
bool DivideRegion( int currentThread ) ;
int RegionToDivide() ;
// void SeparateImages() ;
// void AllocateImages() ;
// bool DivideRegion( int currentThread ) ;
// int RegionToDivide() ;
typename InterpolateImageFunctionType::Pointer m_Interpol[ 6 ] ;
ImagePointer m_Image[ 6 ] ;
Semaphore::Pointer m_Threads ;
int m_SplitAxis ;
bool m_SeparationDone ;
bool m_CannotSplit ;
MutexLock::Pointer m_Lock ;
MutexLock::Pointer m_LockNewThreadDetected ;
std::vector< RegionType > m_ListRegions ;
int m_NbThread ;
MutexLock::Pointer m_CheckRegionsDone ;
bool m_ExceptionThrown ;
bool m_AllocateInterpolatorsDone ;
ImagePointer m_ImageVec[ 6 ] ;
int m_NumberOfThreads ;
// Semaphore::Pointer m_Threads ;
// int m_SplitAxis ;
// bool m_SeparationDone ;
// bool m_CannotSplit ;
// MutexLock::Pointer m_Lock ;
// MutexLock::Pointer m_LockNewThreadDetected ;
// std::vector< RegionType > m_ListRegions ;
// int m_NbThread ;
// MutexLock::Pointer m_CheckRegionsDone ;
// bool m_ExceptionThrown ;
// SizeType m_Size ;
// bool m_AllocateInterpolatorsDone ;
};
}//end namespace itk
......
......@@ -20,305 +20,57 @@
namespace itk
{
template< class TData >
DiffusionTensor3DInterpolateImageFunctionReimplementation< TData >
template< class TData , class TCoordRep >
DiffusionTensor3DInterpolateImageFunctionReimplementation< TData , TCoordRep >
::DiffusionTensor3DInterpolateImageFunctionReimplementation()
{
m_Threads = Semaphore::New() ;
m_SeparationDone = false ;
m_CannotSplit = false ;
m_Lock = MutexLock::New() ;
m_LockNewThreadDetected = MutexLock::New() ;
m_NbThread = 0 ;
m_CheckRegionsDone = MutexLock::New() ;
m_ExceptionThrown = false ;
m_AllocateInterpolatorsDone = false ;
m_NumberOfThreads = 0 ;
}
//allocate the memory for the 6 images containing the 6 separated components of the DTI
//This function is run only by one single thread that blocks all the other threads if any
template< class TData >
template< class TData , class TCoordRep >
void
DiffusionTensor3DInterpolateImageFunctionReimplementation< TData >
::AllocateImages()
DiffusionTensor3DInterpolateImageFunctionReimplementation< TData , TCoordRep >
::SetInputImage( const DiffusionImageType *inputImage )
{
itkRegionType itkRegion ;
itkRegion = this->m_InputImage->GetLargestPossibleRegion() ;
SizeType size = itkRegion.GetSize() ;
typename DiffusionImageType::PointType origin = this->m_InputImage->GetOrigin() ;
typename DiffusionImageType::DirectionType direction = this->m_InputImage->GetDirection() ;
typename DiffusionImageType::SpacingType spacing = this->m_InputImage->GetSpacing() ;
//Read the DTI input image size, spacing, origin and direction and create 6 scalar images with the same properties
for( int i = 0 ; i < 6 ; i++ )
{
m_Image[ i ] = ImageType::New() ;
m_Image[ i ]->SetRegions( size ) ;
m_Image[ i ]->SetOrigin( origin ) ;
m_Image[ i ]->SetDirection( direction ) ;
m_Image[ i ]->SetSpacing( spacing ) ;
m_Image[ i ]->Allocate() ;
}
/* //allocate the interpolator in the derivated class. This is only a virtual class that doesn't know the interpolator type
DiffusionTensor3DInterpolateImageFunction< DataType >::SetInputImage( inputImage ) ;//separateFilter->GetOutput( 0 ) ) ;
if( !inputImage )
{
return ;
}
typedef SeparateComponentsOfADiffusionTensorImage< TData , TData > SeparateType ;
typename SeparateType::Pointer separateFilter = SeparateType::New() ;
separateFilter->SetInput( inputImage ) ;
separateFilter->SetNumberOfThreads( this->m_NumberOfThreads ) ;
separateFilter->Update() ;
AllocateInterpolator() ;
for( int i = 0 ; i < 6 ; i++ )
{
m_Interpol[ i ]->SetInputImage( m_Image[ i ] ) ;
}*/
// Copied from itkImageSource.txx and adapted
// split on the outermost dimension available
m_SplitAxis = this->m_InputImage->GetImageDimension() - 1 ;
while( size[ m_SplitAxis ] == 1)
{
--m_SplitAxis;
if( m_SplitAxis < 0 )
{ // cannot split
itkDebugMacro(" Cannot Split" ) ;
m_CannotSplit = true ;
}
}
m_NbThread = 0 ;
//Initialize the list of region with a first region that contains the whole image
m_ListRegions.clear() ;
RegionType region ;
region.itkRegion = itkRegion ;
region.Done = false ;
Index< 3 > position ;
position.Fill( 0 ) ;
region.PositionInImage = position ;
region.Stop = false ;
m_ListRegions.push_back( region ) ;
this->latestTime = Object::GetMTime() ;
m_ExceptionThrown = false ;
}
//Checks the region that has the most voxel not copied into the separated images and returns its ID
template< class TData >
int
DiffusionTensor3DInterpolateImageFunctionReimplementation< TData >
::RegionToDivide()
{
int largestRegion = -1 ;//returns -1 if all the regions have been processed completely
if( m_ListRegions.size() != 0 ) //at least one region exist
{
//Find largest region
long largestNumber = 0 ;
for( int i = 0 ; i < (signed)m_ListRegions.size() ; i++ )
{
if( !m_ListRegions[ i ].Done )
{
SizeType size = m_ListRegions[ i ].itkRegion.GetSize() ;
Index< 3 > start = m_ListRegions[ i ].itkRegion.GetIndex() ;//Where the region starts in the image
Index< 3 > position = m_ListRegions[ i ].PositionInImage ;//Where we are in the image
long currentPos = position[ 0 ] - start[ 0 ] + size[ 0 ] * ( position[ 1 ] - start[ 1 ] + size[ 1 ] * ( position[ 2 ] - start[ 2 ] ) ) ;//The number of voxels copied in this region
long currentSize = size[ 0 ] * size[ 1 ] * size[ 2 ] ;//The size of the current region
long currentNotDone = currentSize - currentPos ;//The number of voxels that are not copied yet
if( largestNumber < currentNotDone )
{
largestRegion = i ;
largestNumber = currentNotDone ;
}
}
}
m_Interpol[ i ]->SetInputImage( separateFilter->GetOutput( i ) ) ;
}
return largestRegion ;
}
//Divide a region in 2 regions
template< class TData >
bool
DiffusionTensor3DInterpolateImageFunctionReimplementation< TData >
::DivideRegion( int currentThread )
{
//Compute how much has been completed in the current region
itkRegionType itkRegion = m_ListRegions[ currentThread ].itkRegion ;
SizeType oldSize = itkRegion.GetSize() ;
typename ImageType::IndexType start = itkRegion.GetIndex() ;//position in the image of the origin of the current region
int portionDone ;
Index< 3 > position ;
position = m_ListRegions[ currentThread ].PositionInImage ;//what voxel has to be copied next (its position in the image)
portionDone = position[ m_SplitAxis ] - start[ m_SplitAxis ] ;//number of line/plane already copied along the split axis
//Split region
start[ m_SplitAxis ] += portionDone ;
oldSize[ m_SplitAxis ] -= portionDone ;//size of the region that still needs to be copied
itkRegion.SetIndex( start ) ;
SizeType size ;
size = oldSize ;
size[ m_SplitAxis ] = oldSize[ m_SplitAxis ] / 2 + oldSize[ m_SplitAxis ]%2 ;//size of the first new region
itkRegion.SetSize( size ) ;
m_ListRegions[ currentThread ].itkRegion = itkRegion ;
m_ListRegions[ currentThread ].Done = false ;