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