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 "itkMath.h"
20 #include "itkNumericTraits.h"
21 #include "itkLevelSetContainer.h"
22 #include "itkLevelSetEquationChanAndVeseExternalTerm.h"
23 #include "itkSinRegularizedHeavisideStepFunction.h"
24 #include "itkBinaryImageToLevelSetImageAdaptor.h"
25 #include "itkTestingMacros.h"
26
itkLevelSetEquationChanAndVeseExternalTermTest(int argc,char * argv[])27 int itkLevelSetEquationChanAndVeseExternalTermTest( int argc, char* argv[] )
28 {
29
30 if( argc < 2 )
31 {
32 std::cerr << "Missing Arguments" << std::endl;
33 std::cerr << "Program " << itkNameOfTestExecutableMacro(argv) << std::endl;
34 return EXIT_FAILURE;
35 }
36
37 constexpr unsigned int Dimension = 2;
38
39 using InputPixelType = unsigned short;
40 using InputImageType = itk::Image< InputPixelType, Dimension >;
41 using IdentifierType = itk::IdentifierType;
42
43 using InputPixelType = unsigned short;
44 using InputImageType = itk::Image< InputPixelType, Dimension >;
45
46 using PixelType = float;
47 using SparseLevelSetType =
48 itk::WhitakerSparseLevelSetImage< PixelType, Dimension >;
49 using BinaryToSparseAdaptorType =
50 itk::BinaryImageToLevelSetImageAdaptor< InputImageType, SparseLevelSetType >;
51
52 using LevelSetContainerType = itk::LevelSetContainer< IdentifierType, SparseLevelSetType >;
53 using ChanAndVeseExternalTermType =
54 itk::LevelSetEquationChanAndVeseExternalTerm< InputImageType, LevelSetContainerType >;
55
56 using IdListType = std::list< IdentifierType >;
57 using IdListImageType = itk::Image< IdListType, Dimension >;
58 using CacheImageType = itk::Image< short, Dimension >;
59 using DomainMapImageFilterType =
60 itk::LevelSetDomainMapImageFilter< IdListImageType, CacheImageType >;
61
62 using ChanAndVeseExternalTermType =
63 itk::LevelSetEquationChanAndVeseExternalTerm< InputImageType, LevelSetContainerType >;
64
65 using LevelSetOutputRealType = SparseLevelSetType::OutputRealType;
66 using HeavisideFunctionBaseType =
67 itk::SinRegularizedHeavisideStepFunction< LevelSetOutputRealType, LevelSetOutputRealType >;
68 using InputImageIteratorType = itk::ImageRegionIteratorWithIndex< InputImageType >;
69
70 // load binary mask
71 InputImageType::SizeType size;
72 size.Fill( 50 );
73
74 InputImageType::PointType origin;
75 origin[0] = 0.0;
76 origin[1] = 0.0;
77
78 InputImageType::SpacingType spacing;
79 spacing[0] = 1.0;
80 spacing[1] = 1.0;
81
82 InputImageType::IndexType index;
83 index.Fill( 0 );
84
85 InputImageType::RegionType region;
86 region.SetIndex( index );
87 region.SetSize( size );
88
89 // Binary initialization
90 InputImageType::Pointer binary = InputImageType::New();
91 binary->SetRegions( region );
92 binary->SetSpacing( spacing );
93 binary->SetOrigin( origin );
94 binary->Allocate();
95 binary->FillBuffer( itk::NumericTraits<InputPixelType>::ZeroValue() );
96
97 index.Fill( 10 );
98 size.Fill( 30 );
99
100 region.SetIndex( index );
101 region.SetSize( size );
102
103 InputImageIteratorType iIt( binary, region );
104 iIt.GoToBegin();
105 while( !iIt.IsAtEnd() )
106 {
107 iIt.Set( itk::NumericTraits<InputPixelType>::OneValue() );
108 ++iIt;
109 }
110
111 // Convert binary mask to sparse level set
112 BinaryToSparseAdaptorType::Pointer adaptor = BinaryToSparseAdaptorType::New();
113 adaptor->SetInputImage( binary );
114 adaptor->Initialize();
115 std::cout << "Finished converting to sparse format" << std::endl;
116
117 SparseLevelSetType::Pointer level_set = adaptor->GetModifiableLevelSet();
118
119 IdListType list_ids;
120 list_ids.push_back( 1 );
121
122 IdListImageType::Pointer id_image = IdListImageType::New();
123 id_image->SetRegions( binary->GetLargestPossibleRegion() );
124 id_image->Allocate();
125 id_image->FillBuffer( list_ids );
126
127 DomainMapImageFilterType::Pointer domainMapFilter = DomainMapImageFilterType::New();
128 domainMapFilter->SetInput( id_image );
129 domainMapFilter->Update();
130 std::cout << "Domain map computed" << std::endl;
131
132 // Define the Heaviside function
133 HeavisideFunctionBaseType::Pointer heaviside = HeavisideFunctionBaseType::New();
134 heaviside->SetEpsilon( 1.0 );
135
136 // Insert the levelsets in a levelset container
137 LevelSetContainerType::Pointer lscontainer = LevelSetContainerType::New();
138 lscontainer->SetHeaviside( heaviside );
139 lscontainer->SetDomainMapFilter( domainMapFilter );
140
141 const bool LevelSetNotYetAdded = lscontainer->AddLevelSet( 0, level_set, false );
142 if ( !LevelSetNotYetAdded )
143 {
144 return EXIT_FAILURE;
145 }
146
147 // Create ChanAndVese External term for phi_{1}
148 ChanAndVeseExternalTermType::Pointer cvExternalTerm0 = ChanAndVeseExternalTermType::New();
149 cvExternalTerm0->SetInput( binary );
150 cvExternalTerm0->SetCoefficient( 1.0 );
151 cvExternalTerm0->SetCurrentLevelSetId( 0 );
152 cvExternalTerm0->SetLevelSetContainer( lscontainer );
153 std::cout << "CV External term created" << std::endl;
154
155 // Initialize the ChanAndVese term here
156 cvExternalTerm0->InitializeParameters();
157 InputImageIteratorType it( binary, binary->GetLargestPossibleRegion() );
158 it.GoToBegin();
159
160 while( !it.IsAtEnd() )
161 {
162 cvExternalTerm0->Initialize( it.GetIndex() );
163 ++it;
164 }
165
166 cvExternalTerm0->Update();
167
168 if( itk::Math::NotAlmostEquals( cvExternalTerm0->GetMean(),
169 itk::NumericTraits< ChanAndVeseExternalTermType::InputPixelRealType >::OneValue() ) )
170 {
171 return EXIT_FAILURE;
172 }
173
174 return EXIT_SUCCESS;
175 }
176