Skip to content

Commit

Permalink
ENH: TransformToStrainFilter Green-Lagrangian strain
Browse files Browse the repository at this point in the history
  • Loading branch information
@thewtex
thewtex committed Sep 30, 2016
1 parent 306e458 commit 6559ad7
Show file tree
Hide file tree
Showing 4 changed files with 81 additions and 55 deletions.
4 changes: 0 additions & 4 deletions include/itkStrainImageFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,18 +120,14 @@ StrainImageFilter< TInputImage, TOperatorValueType, TOutputValueType >
const GradientOutputPixelType gradientPixel = gradientIt.Get();
for( unsigned int j = 0; j < i; ++j )
{
std::cout << "gradient: " << gradientPixel[j] << std::endl;
outputPixel( i, j ) += gradientPixel[j] / static_cast< TOutputValueType >( 2 );
}
// j == i
std::cout << "gradient: " << gradientPixel[i] << std::endl;
outputPixel( i, i ) += gradientPixel[i];
for( unsigned int j = i + 1; j < ImageDimension; ++j )
{
std::cout << "gradient: " << gradientPixel[j] << std::endl;
outputPixel( i, j ) += gradientPixel[j] / static_cast< TOutputValueType >( 2 );
}
std::cout << "output: " << outputPixel << std::endl;
outputIt.Set( outputPixel );
}
}
Expand Down
68 changes: 31 additions & 37 deletions include/itkTransformToStrainFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,38 +77,32 @@ TransformToStrainFilter< TTransform, TOperatorValue, TOutputValue >
}
outputIt.Set( outputPixel );
}
//for( unsigned int i = 0; i < ImageDimension; ++i )
//{
//itk::ImageRegionConstIterator< GradientOutputImageType >
//gradientIt( reinterpret_cast< GradientOutputImageType* >(
//dynamic_cast< GradientOutputImageType* >(
//this->ProcessObject::GetOutput( i + 1 ) ) )
//, region );
//for( outputIt.GoToBegin(), gradientIt.GoToBegin();
//!gradientIt.IsAtEnd();
//++outputIt, ++gradientIt )
//{
//outputPixel = outputIt.Get();
//gradientPixel = gradientIt.Get();
//for( j = 0; j < i; ++j )
//{
//outputPixel( i, j ) += gradientPixel[j] / static_cast< TOutputValue >( 2 );
//}
//// j == i
//outputPixel( i, i ) += gradientPixel[i];
//for( j = i + 1; j < ImageDimension; ++j )
//{
//outputPixel( i, j ) += gradientPixel[j] / static_cast< TOutputValue >( 2 );
//}
//outputIt.Set( outputPixel );
//}
//}
//switch( m_StrainForm )
//{
//case INFINITESIMAL:
//break;
//// e_ij += 1/2 du_m/du_i du_m/du_j
//case GREENLAGRANGIAN:
switch( m_StrainForm )
{
case INFINITESIMAL:
break;
// e_ij += 1/2 du_m/du_i du_m/du_j
case GREENLAGRANGIAN:
for( outputIt.GoToBegin(); !outputIt.IsAtEnd(); ++outputIt )
{
const typename OutputImageType::IndexType index = outputIt.GetIndex();
typename OutputImageType::PointType point;
output->TransformIndexToPhysicalPoint( index, point );
typename TransformType::JacobianType jacobian;
input->ComputeJacobianWithRespectToPosition( point, jacobian );
typename OutputImageType::PixelType outputPixel = outputIt.Get();
for( unsigned int ii = 0; ii < ImageDimension; ++ii )
{
for( unsigned int jj = 0; jj < ImageDimension; ++jj )
{
for( unsigned int kk = 0; kk <= jj; ++kk )
{
//outputPixel( jj, kk ) += jacobian( ii, jj ) * jacobian( ii, kk ) / static_cast< TOutputValue >( 2 );
}
}
}
outputIt.Set( outputPixel );
}
//for( unsigned int i = 0; i < ImageDimension; ++i )
//{
//itk::ImageRegionConstIterator< GradientOutputImageType >
Expand All @@ -133,8 +127,8 @@ TransformToStrainFilter< TTransform, TOperatorValue, TOutputValue >
//outputIt.Set( outputPixel );
//}
//}
//break;
//// e_ij -= 1/2 du_m/du_i du_m/du_j
break;
// e_ij -= 1/2 du_m/du_i du_m/du_j
//case EULERIANALMANSI:
//for( unsigned int i = 0; i < ImageDimension; ++i )
//{
Expand All @@ -161,9 +155,9 @@ TransformToStrainFilter< TTransform, TOperatorValue, TOutputValue >
//}
//}
//break;
//default:
//itkExceptionMacro( << "Unknown strain form." );
//}
default:
itkExceptionMacro( << "Unknown strain form." );
}
}

} // end namespace itk
Expand Down
26 changes: 20 additions & 6 deletions test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,15 +60,29 @@ itk_add_test(NAME itkStrainImageFilterRecursiveGaussianTest
#${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTest.mha
#DATA{${ITK_DATA_ROOT}/Input/parametersBSpline.txt})

itk_add_test(NAME itkTransformToStrainFilterTest
itk_add_test(NAME itkTransformToStrainFilterInfinitesimalTest
COMMAND StrainTestDriver
--compare
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestAffine.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestAffineToDisplacementFieldStrain.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestInfinitesimalAffine.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestInfinitesimalAffineToDisplacementFieldStrain.mha
itkTransformToStrainFilterTest
"INFINITESIMAL"
"Affine"
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestAffine.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestAffineToDisplacementField.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestAffineToDisplacementFieldStrain.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestInfinitesimalAffine.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestInfinitesimalAffineToDisplacementField.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestInfinitesimalAffineToDisplacementFieldStrain.mha
)

itk_add_test(NAME itkTransformToStrainFilterGreenLagrangianTest
COMMAND StrainTestDriver
--compare
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestGreenLagrangianAffine.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestGreenLagrangianAffineToDisplacementFieldStrain.mha
itkTransformToStrainFilterTest
"GREENLAGRANGIAN"
"Affine"
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestGreenLagrangianAffine.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestGreenLagrangianAffineToDisplacementField.mha
${ITK_TEST_OUTPUT_DIR}/itkTransformToStrainFilterTestGreenLagrangianAffineToDisplacementFieldStrain.mha
)

38 changes: 30 additions & 8 deletions test/itkTransformToStrainFilterTest.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,21 +25,41 @@
int itkTransformToStrainFilterTest( int argc, char * argv [] )
{
/** Check command line arguments. */
if( argc < 5 )
if( argc < 6 )
{
std::cerr << "Usage: ";
std::cerr << argv[0] << "<transformName> <strainFieldFileName> <displacementField> <displacementFieldStrain> [bSplineParametersFile]" << std::endl;
std::cerr << argv[0] << "<strainForm> <transformName> <strainFieldFileName> <displacementField> <displacementFieldStrain> [bSplineParametersFile]" << std::endl;
return EXIT_FAILURE;
}

const std::string transformName = argv[1];
const std::string strainFieldFileName = argv[2];
const std::string displacementFieldFileName = argv[3];
const std::string displacementFieldStrainFileName = argv[4];

int strainForm = 0;
if( ! strcmp( argv[1], "INFINITESIMAL" ) )
{
strainForm = 0;
}
else if( ! strcmp( argv[1], "GREENLAGRANGIAN" ) )
{
strainForm = 1;
}
else if( ! strcmp( argv[1], "EULERIANALMANSI" ) )
{
strainForm = 2;
}
else
{
std::cerr << "Unknown strain form: " << argv[1] << std::endl;
return EXIT_FAILURE;
}

const std::string transformName = argv[2];
const std::string strainFieldFileName = argv[3];
const std::string displacementFieldFileName = argv[4];
const std::string displacementFieldStrainFileName = argv[5];
std::string bSplineParametersFile;
if( argc > 5 )
if( argc > 6 )
{
bSplineParametersFile = argv[5];
bSplineParametersFile = argv[6];
}

/** Typedefs. */
Expand All @@ -55,6 +75,7 @@ int itkTransformToStrainFilterTest( int argc, char * argv [] )

typedef itk::TransformToStrainFilter< TransformType, ScalarPixelType, ScalarPixelType > TransformToStrainFilterType;
TransformToStrainFilterType::Pointer transformToStrainFilter = TransformToStrainFilterType::New();
transformToStrainFilter->SetStrainForm( static_cast< TransformToStrainFilterType::StrainFormType >( strainForm ) );

// Create output information.
typedef TransformToStrainFilterType::SizeType SizeType;
Expand Down Expand Up @@ -196,6 +217,7 @@ int itkTransformToStrainFilterTest( int argc, char * argv [] )
typedef itk::StrainImageFilter< DisplacementFieldType, CoordRepresentationType > StrainImageFilterType;
StrainImageFilterType::Pointer strainImageFilter = StrainImageFilterType::New();
strainImageFilter->SetInput( transformToDisplacement->GetOutput() );
strainImageFilter->SetStrainForm( static_cast< StrainImageFilterType::StrainFormType >( strainForm ) );
writer->SetInput( strainImageFilter->GetOutput() );
writer->SetFileName( displacementFieldStrainFileName );

Expand Down

0 comments on commit 6559ad7

Please sign in to comment.