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 "itkImageFileReader.h"
20 #include "itkShiSparseLevelSetImage.h"
21 #include "itkLevelSetContainer.h"
22 #include "itkLevelSetEquationChanAndVeseInternalTerm.h"
23 #include "itkLevelSetEquationChanAndVeseExternalTerm.h"
24 #include "itkLevelSetEquationTermContainer.h"
25 #include "itkLevelSetEquationContainer.h"
26 #include "itkLevelSetEvolution.h"
27 #include "itkLevelSetEvolutionNumberOfIterationsStoppingCriterion.h"
28 #include "itkBinaryImageToLevelSetImageAdaptor.h"
29 #include "itkAtanRegularizedHeavisideStepFunction.h"
30 #include "itkLevelSetDomainMapImageFilter.h"
31 #include "itkTestingMacros.h"
32
itkMultiLevelSetShiImageSubset2DTest(int argc,char * argv[])33 int itkMultiLevelSetShiImageSubset2DTest( int argc, char* argv[] )
34 {
35 if( argc < 1 )
36 {
37 std::cerr << "Missing Arguments: " << itkNameOfTestExecutableMacro(argv) << std::endl;
38 return EXIT_FAILURE;
39 }
40
41 constexpr unsigned int Dimension = 2;
42
43 using InputPixelType = unsigned short;
44 using InputImageType = itk::Image< InputPixelType, Dimension >;
45 using InputIteratorType = itk::ImageRegionIteratorWithIndex< InputImageType >;
46
47 using PixelType = float;
48 using LevelSetType = itk::ShiSparseLevelSetImage< Dimension >;
49
50 using LevelSetOutputRealType = LevelSetType::OutputRealType;
51
52 using IdentifierType = itk::IdentifierType;
53 using LevelSetContainerType = itk::LevelSetContainer< IdentifierType, LevelSetType >;
54
55 using ChanAndVeseInternalTermType =
56 itk::LevelSetEquationChanAndVeseInternalTerm< InputImageType, LevelSetContainerType >;
57 using ChanAndVeseExternalTermType =
58 itk::LevelSetEquationChanAndVeseExternalTerm< InputImageType, LevelSetContainerType >;
59 using TermContainerType =
60 itk::LevelSetEquationTermContainer< InputImageType, LevelSetContainerType >;
61
62 using EquationContainerType = itk::LevelSetEquationContainer< TermContainerType >;
63 using LevelSetEvolutionType = itk::LevelSetEvolution< EquationContainerType, LevelSetType >;
64
65 using HeavisideFunctionBaseType = itk::AtanRegularizedHeavisideStepFunction<
66 LevelSetOutputRealType, LevelSetOutputRealType >;
67
68 using BinaryImageToLevelSetType = itk::BinaryImageToLevelSetImageAdaptor< InputImageType, LevelSetType >;
69
70 using IdListType = std::list< IdentifierType >;
71 using IdListImageType = itk::Image< IdListType, Dimension >;
72 using CacheImageType = itk::Image< short, Dimension >;
73 using IdIteratorType = itk::ImageRegionIteratorWithIndex< IdListImageType >;
74 using DomainMapImageFilterType =
75 itk::LevelSetDomainMapImageFilter< IdListImageType, CacheImageType >;
76
77 // load binary input
78 InputImageType::SizeType size;
79 size.Fill( 1000 );
80
81 InputImageType::PointType origin;
82 origin[0] = 0.0;
83 origin[1] = 0.0;
84
85 InputImageType::SpacingType spacing;
86 spacing[0] = 1.0;
87 spacing[1] = 1.0;
88
89 InputImageType::IndexType index;
90 index.Fill( 0 );
91
92 InputImageType::RegionType region;
93 region.SetIndex( index );
94 region.SetSize( size );
95
96 // Input initialization
97 InputImageType::Pointer input = InputImageType::New();
98 input->SetRegions( region );
99 input->SetSpacing( spacing );
100 input->SetOrigin( origin );
101 input->Allocate();
102 input->FillBuffer( itk::NumericTraits<InputPixelType>::ZeroValue() );
103
104 index.Fill( 910 );
105 size.Fill( 80 );
106
107 region.SetIndex( index );
108 region.SetSize( size );
109
110 InputIteratorType iIt( input, region );
111 iIt.GoToBegin();
112 while( !iIt.IsAtEnd() )
113 {
114 iIt.Set( 100*itk::NumericTraits<InputPixelType>::OneValue() );
115 ++iIt;
116 }
117
118 index.Fill( 0 );
119 size.Fill( 100 );
120 origin.Fill( 900.0 );
121 region.SetIndex( index );
122 region.SetSize( size );
123
124 // Binary initialization
125 InputImageType::Pointer binary = InputImageType::New();
126 binary->SetRegions( region );
127 binary->SetSpacing( spacing );
128 binary->SetOrigin( origin );
129 binary->Allocate();
130 binary->FillBuffer( itk::NumericTraits<InputPixelType>::ZeroValue() );
131
132 index.Fill( 30 );
133 size.Fill( 40 );
134 region.SetIndex( index );
135 region.SetSize( size );
136
137 InputIteratorType bIt( binary, region );
138 bIt.GoToBegin();
139 while( !bIt.IsAtEnd() )
140 {
141 bIt.Set( itk::NumericTraits<InputPixelType>::OneValue() );
142 ++bIt;
143 }
144
145 // Convert binary mask to dense level set
146 BinaryImageToLevelSetType::Pointer adaptor1 = BinaryImageToLevelSetType::New();
147 adaptor1->SetInputImage( binary );
148 adaptor1->Initialize();
149 LevelSetType::Pointer levelSet1 = adaptor1->GetModifiableLevelSet();
150
151 input->TransformPhysicalPointToIndex( binary->GetOrigin(), index );
152 InputImageType::OffsetType offset;
153 for( unsigned int i = 0; i < Dimension; i++ )
154 {
155 offset[i] = index[i];
156 }
157 levelSet1->SetDomainOffset( offset );
158
159 IdListType listIds;
160 listIds.clear();
161
162 index.Fill( 900 );
163 size.Fill( 100 );
164 region.SetIndex( index );
165 region.SetSize( size );
166
167 IdListImageType::Pointer idImage = IdListImageType::New();
168 idImage->SetRegions( input->GetLargestPossibleRegion() );
169 idImage->Allocate();
170 idImage->FillBuffer( listIds );
171
172 listIds.push_back( 1 );
173 IdIteratorType it( idImage, region );
174 it.GoToBegin();
175 while( !it.IsAtEnd() )
176 {
177 it.Set( listIds );
178 ++it;
179 }
180
181 DomainMapImageFilterType::Pointer domainMapFilter = DomainMapImageFilterType::New();
182 domainMapFilter->SetInput( idImage );
183 domainMapFilter->Update();
184 std::cout << "Domain map computed" << std::endl;
185
186 // Define the Heaviside function
187 HeavisideFunctionBaseType::Pointer heaviside = HeavisideFunctionBaseType::New();
188 heaviside->SetEpsilon( 1.0 );
189
190 // Insert the levelsets in a levelset container
191 LevelSetContainerType::Pointer lscontainer = LevelSetContainerType::New();
192 lscontainer->SetHeaviside( heaviside );
193 lscontainer->SetDomainMapFilter( domainMapFilter );
194
195 bool levelSetNotYetAdded = lscontainer->AddLevelSet( 0, levelSet1, false );
196 if ( !levelSetNotYetAdded )
197 {
198 return EXIT_FAILURE;
199 }
200 std::cout << "Level set container created" << std::endl;
201
202 // Create ChanAndVese internal term for phi_{1}
203 ChanAndVeseInternalTermType::Pointer cvInternalTerm0 = ChanAndVeseInternalTermType::New();
204 cvInternalTerm0->SetInput( input );
205 cvInternalTerm0->SetCoefficient( 1.0 );
206 std::cout << "LevelSet 0: CV internal term created" << std::endl;
207
208 // Create ChanAndVese external term for phi_{1}
209 ChanAndVeseExternalTermType::Pointer cvExternalTerm0 = ChanAndVeseExternalTermType::New();
210 cvExternalTerm0->SetInput( input );
211 cvExternalTerm0->SetCoefficient( 1.0 );
212 std::cout << "LevelSet 0: CV external term created" << std::endl;
213
214
215 // Create Term Container
216 TermContainerType::Pointer termContainer0 = TermContainerType::New();
217 termContainer0->SetInput( input );
218 termContainer0->SetCurrentLevelSetId( 0 );
219 termContainer0->SetLevelSetContainer( lscontainer );
220
221 termContainer0->AddTerm( 0, cvInternalTerm0 );
222 termContainer0->AddTerm( 1, cvExternalTerm0 );
223 std::cout << "Term container 0 created" << std::endl;
224
225
226 EquationContainerType::Pointer equationContainer = EquationContainerType::New();
227 equationContainer->SetLevelSetContainer( lscontainer );
228 equationContainer->AddEquation( 0, termContainer0 );
229 std::cout << "Equation container created" << std::endl;
230
231 using StoppingCriterionType = itk::LevelSetEvolutionNumberOfIterationsStoppingCriterion<LevelSetContainerType>;
232 StoppingCriterionType::Pointer criterion = StoppingCriterionType::New();
233 criterion->SetNumberOfIterations( 500 );
234 std::cout << "Stopping criterion created" << std::endl;
235
236 LevelSetEvolutionType::Pointer evolution = LevelSetEvolutionType::New();
237 evolution->SetEquationContainer( equationContainer );
238 evolution->SetStoppingCriterion( criterion );
239 evolution->SetLevelSetContainer( lscontainer );
240
241 try
242 {
243 evolution->Update();
244 }
245 catch ( itk::ExceptionObject& err )
246 {
247 std::cerr << err << std::endl;
248 return EXIT_FAILURE;
249 }
250
251 PixelType mean = cvInternalTerm0->GetMean();
252 if ( ( mean < 90 ) || ( mean > 105 ) )
253 {
254 std::cerr << "( ( mean < 95 ) || ( mean > 105 ) )" <<std::endl;
255 std::cerr << "mean = " <<mean <<std::endl;
256 return EXIT_FAILURE;
257 }
258
259 mean = cvExternalTerm0->GetMean();
260 if ( ( mean > 20.0 ) )
261 {
262 std::cerr << "( ( mean < 0 ) || ( mean > 5 ) )" <<std::endl;
263 std::cerr << "External mean = " <<mean <<std::endl;
264 return EXIT_FAILURE;
265 }
266
267 return EXIT_SUCCESS;
268 }
269