1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 #include "itkTimeVaryingVelocityFieldIntegrationImageFilter.h"
20
itkTimeVaryingVelocityFieldIntegrationImageFilterTest(int,char * [])21 int itkTimeVaryingVelocityFieldIntegrationImageFilterTest( int, char* [] )
22 {
23 using VectorType = itk::Vector<double, 3>;
24 using DisplacementFieldType = itk::Image<VectorType, 3>;
25 using TimeVaryingVelocityFieldType = itk::Image<VectorType, 4>;
26
27 TimeVaryingVelocityFieldType::PointType origin;
28 origin.Fill( 0.0 );
29
30 TimeVaryingVelocityFieldType::SpacingType spacing;
31 spacing.Fill( 2.0 );
32
33 TimeVaryingVelocityFieldType::SizeType size;
34 size.Fill( 25 );
35
36 VectorType displacement1;
37 displacement1.Fill( 0.1 );
38
39 TimeVaryingVelocityFieldType::Pointer timeVaryingVelocityField =
40 TimeVaryingVelocityFieldType::New();
41
42 timeVaryingVelocityField->SetOrigin( origin );
43 timeVaryingVelocityField->SetSpacing( spacing );
44 timeVaryingVelocityField->SetRegions( size );
45 timeVaryingVelocityField->Allocate();
46 timeVaryingVelocityField->FillBuffer( displacement1 );
47
48
49 using IntegratorType = itk::TimeVaryingVelocityFieldIntegrationImageFilter
50 <TimeVaryingVelocityFieldType, DisplacementFieldType>;
51 IntegratorType::Pointer integrator = IntegratorType::New();
52 integrator->SetInput( timeVaryingVelocityField );
53 integrator->SetLowerTimeBound( 0.3 );
54 integrator->SetUpperTimeBound( 0.75 );
55 integrator->SetNumberOfIntegrationSteps( 10 );
56 integrator->Update();
57
58 integrator->Print( std::cout, 3 );
59
60 DisplacementFieldType::IndexType index;
61 index.Fill( 0 );
62 VectorType displacement;
63
64 IntegratorType::Pointer inverseIntegrator = IntegratorType::New();
65 inverseIntegrator->SetInput( timeVaryingVelocityField );
66 inverseIntegrator->SetLowerTimeBound( 1.0 );
67 inverseIntegrator->SetUpperTimeBound( 0.0 );
68 inverseIntegrator->SetNumberOfIntegrationSteps( 10 );
69 inverseIntegrator->Update();
70
71 // This integration should result in a constant image of value
72 // -( 0.1 * 1.0 - ( 0.1 * 0.0 ) ) = -0.1 with ~epsilon deviation
73 // due to numerical computations
74 const DisplacementFieldType * inverseField = inverseIntegrator->GetOutput();
75 displacement = inverseField->GetPixel( index );
76
77 std::cout << "Estimated inverse displacement vector: " << displacement << std::endl;
78 if( itk::Math::abs( displacement[0] + 0.101852 ) > 0.01 )
79 {
80 std::cerr << "Failed to produce the correct inverse integration." << std::endl;
81 return EXIT_FAILURE;
82 }
83
84 inverseIntegrator->Print( std::cout, 3 );
85
86 std::cout << inverseIntegrator->GetNameOfClass() << std::endl;
87
88 // This integration should result in a constant image of value
89 // 0.75 * 0.1 - 0.3 * 0.1 = 0.045 with ~epsilon deviation
90 // due to numerical computations
91 const DisplacementFieldType * displacementField = integrator->GetOutput();
92
93 displacement = displacementField->GetPixel( index );
94
95 std::cout << "Estimated forward displacement vector: " << displacement << std::endl;
96 if( itk::Math::abs( displacement[0] - 0.045 ) > 0.0001 )
97 {
98 std::cerr << "Failed to produce the correct forward integration." << std::endl;
99 return EXIT_FAILURE;
100 }
101
102 return EXIT_SUCCESS;
103 }
104