Commit d51b8d7b authored by Cyril Poupon's avatar Cyril Poupon
Browse files

splitted the vector field into 2 direct and inverse vector field in VectorFieldTransform3d

parent b6902e49
......@@ -28,11 +28,11 @@ gkg::BundleMapOperatorTransform3d< L, Compare >::
// sanity checks for string parameters
if ( ( stringParameters.size() != 1U ) &&
( stringParameters.size() != 2U ) )
( stringParameters.size() != 3U ) )
{
throw std::runtime_error(
"string parameters must have 1 or 2 arguments" );
"string parameters must have 1 or 3 arguments" );
}
......@@ -44,7 +44,7 @@ gkg::BundleMapOperatorTransform3d< L, Compare >::
}
// reading rigid or affine3D transform
// reading rigid or affine 3D transform
std::string fileNameTransform3d = stringParameters[ 0 ];
try
{
......@@ -96,25 +96,38 @@ gkg::BundleMapOperatorTransform3d< L, Compare >::
}
// reading rigid or affine3D transform
// reading vector field 3D transform
if ( stringParameters.size() == 2U )
{
std::string fileNameVectorFieldTransform3d = stringParameters[ 1 ];
std::string fileNameDirectVectorFieldTransform3d = stringParameters[ 1 ];
gkg::Volume< float > directVectorField;
gkg::Reader::getInstance().read( fileNameDirectVectorFieldTransform3d,
directVectorField );
gkg::Volume< float > vectorField;
gkg::Reader::getInstance().read( fileNameVectorFieldTransform3d,
vectorField );
if ( directVectorField.getSizeT() != 3 )
{
throw std::runtime_error( "direct vector field T size must be 3" );
}
std::string fileNameInverseVectorFieldTransform3d = stringParameters[ 2 ];
gkg::Volume< float > inverseVectorField;
gkg::Reader::getInstance().read( fileNameInverseVectorFieldTransform3d,
inverseVectorField );
if ( vectorField.getSizeT() != 6 )
if ( inverseVectorField.getSizeT() != 3 )
{
throw std::runtime_error( "vector field T size must be 6" );
throw std::runtime_error( "inverse vector field T size must be 3" );
}
_vectorFieldTransform3d.reset(
new gkg::VectorFieldTransform3d< float >( vectorField ) );
new gkg::VectorFieldTransform3d< float >(
directVectorField,
inverseVectorField ) );
}
......@@ -349,7 +362,8 @@ gkg::BundleMapOperatorTransform3d< L, Compare >::getStringParameterHelp()
{
return std::string( ". <P1>: 3D affine or rigid transform filename\n"
". <P2>: 3D vector field transform filename "
". <P2>: 3D direct vector field transform filename "
". <P3>: 3D inverse vector field transform filename "
"(optional)\n" );
}
......
......@@ -327,7 +327,7 @@ EddyCurrentProcess::correct( gkg::Process& process,
volume->getSizeT() - 1 :
volume->getSizeT() );
std::vector< gkg::RCPointer< gkg::NonLinearTransform3d< float > > >
nonLinearTransform3ds( transform3dCount);
nonLinearTransform3ds( transform3dCount );
std::vector< gkg::RCPointer<
gkg::ParameterizedAffineTransform3d< float > > >
......@@ -802,8 +802,7 @@ EddyCurrentProcess::correct( gkg::Process& process,
linearName += gkg::StringConverter::toString( referenceRank );
linearName += ".trm";
nonLinearTransform3ds[ count ]->writeTrm( linearName,
"" );
nonLinearTransform3ds[ count ]->writeTrm( linearName, "", "" );
++ count;
}
......@@ -829,7 +828,7 @@ EddyCurrentProcess::correct( gkg::Process& process,
linearName += referenceName;
linearName += ".trm";
nonLinearTransform3ds[ t ]->writeTrm( linearName, "" );
nonLinearTransform3ds[ t ]->writeTrm( linearName, "", "" );
}
......
......@@ -383,10 +383,11 @@ void FiberLabellingProcess::label( gkg::Process& process,
}
nonLinearTransform3d->readTrm( fileNameTransform3dBundleMapToAtlas[ 0 ],
"",
"" );
}
else if ( fileNameTransform3dBundleMapToAtlas.size() == 2U )
else if ( fileNameTransform3dBundleMapToAtlas.size() == 3U )
{
if ( verbose )
......@@ -398,7 +399,8 @@ void FiberLabellingProcess::label( gkg::Process& process,
}
nonLinearTransform3d->readTrm(
fileNameTransform3dBundleMapToAtlas[ 0 ],
fileNameTransform3dBundleMapToAtlas[ 1 ] );
fileNameTransform3dBundleMapToAtlas[ 1 ],
fileNameTransform3dBundleMapToAtlas[ 2 ] );
}
......@@ -1407,7 +1409,7 @@ void gkg::DwiFiberLabellingCommand::parse()
"Bundle map to atlas 3D transformation "
"(default=id)",
fileNameTransform3dBundleMapToAtlas,
0, 2 );
0, 3 );
application.addSingleOption( "-d",
"Maximum distance to bundle map configuration "
"file",
......@@ -1538,6 +1540,24 @@ void gkg::DwiFiberLabellingCommand::execute(
try
{
//
// sanity check
//
if ( !fileNameTransform3dBundleMapToAtlas.empty() )
{
if ( ( fileNameTransform3dBundleMapToAtlas.size() != 1U ) &&
( fileNameTransform3dBundleMapToAtlas.size() != 3U ) )
{
throw std::runtime_error(
"number of 3D transformation must be 1 or 3" );
}
}
FiberLabellingProcess
fiberLabellingProcess( fileNameInputBundleMaps,
fileNameAtlasBundleMap,
......
......@@ -1118,16 +1118,26 @@ DwiMotionCorrectionProcess::correct( gkg::Process& process,
linearName += gkg::StringConverter::toString( referenceRank );
linearName += "_linear.trm";
std::string diffeomorphicName = transform3dName;
diffeomorphicName += "_";
diffeomorphicName += gkg::StringConverter::toString( t );
diffeomorphicName += "_to_";
diffeomorphicName += gkg::StringConverter::toString(
std::string diffeomorphicDirectName = transform3dName;
diffeomorphicDirectName += "_";
diffeomorphicDirectName += gkg::StringConverter::toString( t );
diffeomorphicDirectName += "_to_";
diffeomorphicDirectName += gkg::StringConverter::toString(
referenceRank );
diffeomorphicName += "_diffeomorphic.trm";
diffeomorphicDirectName += "_diffeomorphic.direct.trm";
nonLinearTransform3ds[ count ]->writeTrm( linearName,
diffeomorphicName );
std::string diffeomorphicInverseName = transform3dName;
diffeomorphicInverseName += "_";
diffeomorphicInverseName += gkg::StringConverter::toString( t );
diffeomorphicInverseName += "_to_";
diffeomorphicInverseName += gkg::StringConverter::toString(
referenceRank );
diffeomorphicInverseName += "_diffeomorphic.inverse.trm";
nonLinearTransform3ds[ count ]->writeTrm(
linearName,
diffeomorphicDirectName,
diffeomorphicInverseName );
++ count;
......@@ -1154,14 +1164,23 @@ DwiMotionCorrectionProcess::correct( gkg::Process& process,
linearName += referenceName;
linearName += "_linear.trm";
std::string diffeomorphicName = transform3dName;
diffeomorphicName += "_";
diffeomorphicName += gkg::StringConverter::toString( t );
diffeomorphicName += "_to_";
diffeomorphicName += referenceName;
diffeomorphicName += "_diffeomorphic.trm";
nonLinearTransform3ds[ t ]->writeTrm( linearName, diffeomorphicName );
std::string diffeomorphicDirectName = transform3dName;
diffeomorphicDirectName += "_";
diffeomorphicDirectName += gkg::StringConverter::toString( t );
diffeomorphicDirectName += "_to_";
diffeomorphicDirectName += referenceName;
diffeomorphicDirectName += "_diffeomorphic.direct.trm";
std::string diffeomorphicInverseName = transform3dName;
diffeomorphicInverseName += "_";
diffeomorphicInverseName += gkg::StringConverter::toString( t );
diffeomorphicInverseName += "_to_";
diffeomorphicInverseName += referenceName;
diffeomorphicInverseName += "_diffeomorphic.inverse.trm";
nonLinearTransform3ds[ t ]->writeTrm( linearName,
diffeomorphicDirectName,
diffeomorphicInverseName );
}
......
......@@ -470,7 +470,8 @@ void DwiResamplingProcess::resample( gkg::Process& process,
nonLinearTransform3d = new gkg::NonLinearTransform3d< double >;
nonLinearTransform3d->readTrm(
dwiResamplingProcess.getFileNameTransforms()[ 0 ],
dwiResamplingProcess.getFileNameTransforms()[ 1 ] );
dwiResamplingProcess.getFileNameTransforms()[ 1 ],
dwiResamplingProcess.getFileNameTransforms()[ 2 ] );
transform3d = nonLinearTransform3d;
foundTransformText = "non-linear transformation detected";
......@@ -1018,6 +1019,15 @@ void gkg::DwiResampling3dCommand::execute(
}
if ( !( ( fileNameTransforms.size() == 1U ) ||
( fileNameTransforms.size() == 3U ) ) )
{
throw std::runtime_error(
"number of 3D transformation filenames must be 1 or 3" );
}
//
// launching resampling process
//
......
......@@ -495,7 +495,7 @@ if ( ( options.transformationType == 'affine_and_diffeomorphic' ) or
# converting NIFTI direct diffeomorphism to GIS
fileNameDirectNiiDiffeomorphism = fileNameAntsTransform + \
'1Warp.nii.gz'
'1InverseWarp.nii.gz'
fileNameDirectImaDiffeomorphism = os.path.join( intermediateDirectory,
'diffeomorphism_direct.ima')
gkg.executeCommand( { 'algorithm' : 'Nifti2GisConverter',
......@@ -511,7 +511,7 @@ if ( ( options.transformationType == 'affine_and_diffeomorphic' ) or
# converting NIFTI iinverse diffeomorphism to GIS
fileNameInverseNiiDiffeomorphism = fileNameAntsTransform + \
'1InverseWarp.nii.gz'
'1Warp.nii.gz'
fileNameInverseImaDiffeomorphism = os.path.join( intermediateDirectory,
'diffeomorphism_inverse.ima' )
gkg.executeCommand( { 'algorithm' : 'Nifti2GisConverter',
......
......@@ -22,10 +22,10 @@ gkg::MotionCorrection< T >::MotionCorrection(
int32_t resamplingOrder,
const T& background,
bool verbose )
: _registrationAlgorithm( registrationAlgorithm ),
_resamplingOrder( resamplingOrder ),
_background( background ),
_verbose( verbose )
: _registrationAlgorithm( registrationAlgorithm ),
_resamplingOrder( resamplingOrder ),
_background( background ),
_verbose( verbose )
{
}
......
......@@ -563,10 +563,13 @@ void gkg::RegistrationAlgorithm< T1, T2 >::match(
std::max( diffeomorphismSize.y,
diffeomorphismSize.z ) );
gkg::Volume< float >
directAndInverseVectorFields( diffeomorphismSize, 6 );
directAndInverseVectorFields.setResolution( diffeomorphismResolution );
directAndInverseVectorFields.fill( 0.0f );
gkg::Volume< float > directVectorField( diffeomorphismSize, 3 );
directVectorField.setResolution( diffeomorphismResolution );
directVectorField.fill( 0.0f );
gkg::Volume< float > inverseVectorField( diffeomorphismSize, 3 );
inverseVectorField.setResolution( diffeomorphismResolution );
inverseVectorField.fill( 0.0f );
gkg::Volume< float > directFlowFieldAlongX( diffeomorphismSize );
directFlowFieldAlongX.setResolution( diffeomorphismResolution );
......@@ -583,7 +586,8 @@ void gkg::RegistrationAlgorithm< T1, T2 >::match(
inverseFlowFieldAlongZ.setResolution( diffeomorphismResolution );
vectorFieldTransform3d->setDirectAndInverseVectorFields(
directAndInverseVectorFields );
directVectorField,
inverseVectorField );
int32_t c = 0;
......
......@@ -3,6 +3,7 @@
#include <gkg-processing-io/TypedVolumeReaderProcess.h>
#include <gkg-core-io/Reader_i.h>
#include <gkg-core-io/Writer_i.h>
#include <gkg-core-exception/Exception.h>
#include <fstream>
......@@ -222,13 +223,24 @@ void gkg::NonLinearTransform3d< T >::getInverse(
template < class T >
void gkg::NonLinearTransform3d< T >::readTrm(
const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearTransform )
const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearDirectTransform,
const std::string& fileNameNonLinearInverseTransform )
{
try
{
// sanity check
if ( fileNameNonLinearDirectTransform.empty() !=
fileNameNonLinearInverseTransform.empty() )
{
throw std::runtime_error(
"direct or inverse non linear transform is missing" );
}
// saving the linear part of the non linear 3D transformation
if ( !fileNameLinearTransform.empty() )
{
......@@ -279,48 +291,72 @@ void gkg::NonLinearTransform3d< T >::readTrm(
}
if ( !fileNameNonLinearTransform.empty() )
if ( !fileNameNonLinearDirectTransform.empty() )
{
gkg::Volume< T > directAndInverseVectorFields;
gkg::Volume< T > directVectorField;
gkg::TypedVolumeReaderProcess< T >
typedVolumeReaderProcess( directAndInverseVectorFields );
typedVolumeReaderProcess.execute( fileNameNonLinearTransform );
typedDirectVolumeReaderProcess( directVectorField );
typedDirectVolumeReaderProcess.execute(
fileNameNonLinearDirectTransform );
gkg::Volume< T > inverseVectorField;
gkg::TypedVolumeReaderProcess< T >
typedInverseVolumeReaderProcess( inverseVectorField );
typedInverseVolumeReaderProcess.execute(
fileNameNonLinearInverseTransform );
if ( _vectorFieldTransform3d.isNull() )
{
_vectorFieldTransform3d.reset( new gkg::VectorFieldTransform3d< T >(
directAndInverseVectorFields ) );
directVectorField,
inverseVectorField ) );
}
else
{
_vectorFieldTransform3d->setDirectAndInverseVectorFields(
directAndInverseVectorFields );
directVectorField,
inverseVectorField );
}
}
}
GKG_CATCH( "template < class T > "
"void gkg::NonLinearTransform3d< T >::readTrm( "
"const std::string& fileNameLinearTransform, "
"const std::string& fileNameNonLinearTransform )" );
"const std::string& fileNameNonLinearDirectTransform, "
"const std::string& fileNameNonLinearInverseTransform )" );
}
template < class T >
void gkg::NonLinearTransform3d< T >::writeTrm(
const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearTransform ) const
const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearDirectTransform,
const std::string& fileNameNonLinearInverseTransform ) const
{
try
{
// sanity check
if ( fileNameNonLinearDirectTransform.empty() !=
fileNameNonLinearInverseTransform.empty() )
{
throw std::runtime_error(
"direct or inverse non linear transform is missing" );
}
// saving the linear part of the non linear 3D transformation
if ( !_linearTransform3d.isNull() && !fileNameLinearTransform.empty() )
{
......@@ -333,14 +369,19 @@ void gkg::NonLinearTransform3d< T >::writeTrm(
// saving the vector field
if ( !_vectorFieldTransform3d.isNull() &&
!fileNameNonLinearTransform.empty() )
!fileNameNonLinearDirectTransform.empty() &&
!fileNameNonLinearInverseTransform.empty() )
{
gkg::Volume< T > directAndInverseVectorFields;
gkg::Volume< T > directVectorField;
gkg::Volume< T > inverseVectorField;
_vectorFieldTransform3d->getDirectAndInverseVectorFields(
directAndInverseVectorFields );
gkg::Writer::getInstance().write( fileNameNonLinearTransform,
directAndInverseVectorFields );
directVectorField,
inverseVectorField );
gkg::Writer::getInstance().write( fileNameNonLinearDirectTransform,
directVectorField );
gkg::Writer::getInstance().write( fileNameNonLinearInverseTransform,
inverseVectorField );
}
......@@ -348,20 +389,32 @@ void gkg::NonLinearTransform3d< T >::writeTrm(
GKG_CATCH( "template < class T > "
"void gkg::NonLinearTransform3d< T >::writeTrm( "
"const std::string& fileNameLinearTransform, "
"const std::string& fileNameNonLinearTransform ) const" );
"const std::string& fileNameNonLinearDirectTransform, "
"const std::string& fileNameNonLinearInverseTransform ) const" );
}
template < class T >
void gkg::NonLinearTransform3d< T >::writeInverseTrm(
const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearTransform ) const
const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearDirectTransform,
const std::string& fileNameNonLinearInverseTransform ) const
{
try
{
// sanity check
if ( fileNameNonLinearDirectTransform.empty() !=
fileNameNonLinearInverseTransform.empty() )
{
throw std::runtime_error(
"direct or inverse non linear transform is missing" );
}
// saving the linear part of the non linear 3D transformation
if ( !_linearTransform3d.isNull() && !fileNameLinearTransform.empty() )
{
......@@ -424,48 +477,20 @@ void gkg::NonLinearTransform3d< T >::writeInverseTrm(
// saving the vector field
if ( !_vectorFieldTransform3d.isNull() &&
!fileNameNonLinearTransform.empty() )
!fileNameNonLinearDirectTransform.empty() &&
!fileNameNonLinearInverseTransform.empty() )
{
gkg::Volume< T > inverseAndDirectVectorFields;
gkg::Volume< T > directVectorField;
gkg::Volume< T > inverseVectorField;
_vectorFieldTransform3d->getDirectAndInverseVectorFields(
inverseAndDirectVectorFields );
int32_t sizeX = inverseAndDirectVectorFields.getSizeX();
int32_t sizeY = inverseAndDirectVectorFields.getSizeY();
int32_t sizeZ = inverseAndDirectVectorFields.getSizeZ();
int32_t x = 0;
int32_t y = 0;
int32_t z = 0;
for ( z = 0; z < sizeZ; z++ )
{
for ( y = 0; y < sizeY; y++ )
{
for ( x = 0; x < sizeX; x++ )
{
std::swap( inverseAndDirectVectorFields( x, y, z, 0 ),
inverseAndDirectVectorFields( x, y, z, 3 ) );
std::swap( inverseAndDirectVectorFields( x, y, z, 1 ),
inverseAndDirectVectorFields( x, y, z, 4 ) );
std::swap( inverseAndDirectVectorFields( x, y, z, 2 ),
inverseAndDirectVectorFields( x, y, z, 5 ) );
}
}
}
directVectorField,
inverseVectorField );
gkg::Writer::getInstance().write( fileNameNonLinearTransform,
inverseAndDirectVectorFields );
gkg::Writer::getInstance().write( fileNameNonLinearDirectTransform,
inverseVectorField );
gkg::Writer::getInstance().write( fileNameNonLinearInverseTransform,
directVectorField );
}
......@@ -473,7 +498,8 @@ void gkg::NonLinearTransform3d< T >::writeInverseTrm(
GKG_CATCH( "template < class T > "
"void gkg::NonLinearTransform3d< T >::writeTrm( "
"const std::string& fileNameLinearTransform, "
"const std::string& fileNameNonLinearTransform ) const" );
"const std::string& fileNameNonLinearDirectTransform, "
"const std::string& fileNameNonLinearInverseTransform ) const" );
}
......
......@@ -46,13 +46,16 @@ class NonLinearTransform3d : public Transform3d< T >
Vector3d< T >& from ) const;
void readTrm( const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearTransform );
const std::string& fileNameNonLinearDirectTransform,
const std::string& fileNameNonLinearInverseTransform );
void writeTrm( const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearTransform ) const;
void writeInverseTrm( const std::string& fileNameLinearTransform,
const std::string& fileNameNonLinearTransform ) const;
const std::string& fileNameNonLinearDirectTransform,
const std::string& fileNameNonLinearInverseTransform ) const;
void writeInverseTrm(
const std::string& fileNameLinearTransform,