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
itkSingleLevelSetWhitakerImage2DTest(int argc,char * argv[])30 int itkSingleLevelSetWhitakerImage2DTest( 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 PixelType = float;
46
47 using SparseLevelSetType =
48 itk::WhitakerSparseLevelSetImage< PixelType, Dimension >;
49 using BinaryToSparseAdaptorType =
50 itk::BinaryImageToLevelSetImageAdaptor< InputImageType, SparseLevelSetType >;
51
52 using IdentifierType = itk::IdentifierType;
53
54 using LevelSetContainerType =
55 itk::LevelSetContainer< IdentifierType, SparseLevelSetType >;
56
57 using ChanAndVeseInternalTermType =
58 itk::LevelSetEquationChanAndVeseInternalTerm< InputImageType, LevelSetContainerType >;
59 using ChanAndVeseExternalTermType =
60 itk::LevelSetEquationChanAndVeseExternalTerm< InputImageType, LevelSetContainerType >;
61 using TermContainerType =
62 itk::LevelSetEquationTermContainer< InputImageType, LevelSetContainerType >;
63
64 using EquationContainerType = itk::LevelSetEquationContainer<TermContainerType>;
65
66 using LevelSetEvolutionType =
67 itk::LevelSetEvolution< EquationContainerType, SparseLevelSetType >;
68
69 using LevelSetOutputRealType = SparseLevelSetType::OutputRealType;
70 using HeavisideFunctionBaseType =
71 itk::SinRegularizedHeavisideStepFunction< LevelSetOutputRealType, LevelSetOutputRealType >;
72 using InputIteratorType = itk::ImageRegionIteratorWithIndex< InputImageType >;
73
74 // load binary mask
75 ReaderType::Pointer reader = ReaderType::New();
76 reader->SetFileName( argv[1] );
77 reader->Update();
78 InputImageType::Pointer input = reader->GetOutput();
79
80 // Binary initialization
81 InputImageType::Pointer binary = InputImageType::New();
82 binary->SetRegions( input->GetLargestPossibleRegion() );
83 binary->CopyInformation( input );
84 binary->Allocate();
85 binary->FillBuffer( itk::NumericTraits<InputPixelType>::ZeroValue() );
86
87 InputImageType::RegionType region;
88 InputImageType::IndexType index;
89 InputImageType::SizeType size;
90
91 index.Fill( 10 );
92 size.Fill( 30 );
93
94 region.SetIndex( index );
95 region.SetSize( size );
96
97 InputIteratorType iIt( binary, region );
98 iIt.GoToBegin();
99 while( !iIt.IsAtEnd() )
100 {
101 iIt.Set( itk::NumericTraits<InputPixelType>::OneValue() );
102 ++iIt;
103 }
104
105 // Convert binary mask to sparse level set
106 BinaryToSparseAdaptorType::Pointer adaptor = BinaryToSparseAdaptorType::New();
107 adaptor->SetInputImage( binary );
108 adaptor->Initialize();
109 std::cout << "Finished converting to sparse format" << std::endl;
110
111 SparseLevelSetType::Pointer level_set = adaptor->GetModifiableLevelSet();
112
113 // Define the Heaviside function
114 HeavisideFunctionBaseType::Pointer heaviside = HeavisideFunctionBaseType::New();
115 heaviside->SetEpsilon( 1.0 );
116
117 // Insert the levelsets in a levelset container
118 LevelSetContainerType::Pointer lscontainer = LevelSetContainerType::New();
119 lscontainer->SetHeaviside( heaviside );
120
121 bool LevelSetNotYetAdded = lscontainer->AddLevelSet( 0, level_set, false );
122 if ( !LevelSetNotYetAdded )
123 {
124 return EXIT_FAILURE;
125 }
126 std::cout << "Level set container created" << std::endl;
127
128 // **************** CREATE ALL TERMS ****************
129
130 // -----------------------------
131 // *** 1st Level Set phi ***
132
133 // Create ChanAndVese internal term for phi_{1}
134 ChanAndVeseInternalTermType::Pointer cvInternalTerm0 = ChanAndVeseInternalTermType::New();
135 cvInternalTerm0->SetInput( input );
136 cvInternalTerm0->SetCoefficient( 1.0 );
137 std::cout << "LevelSet 1: CV internal term created" << std::endl;
138
139 // Create ChanAndVese external term for phi_{1}
140 ChanAndVeseExternalTermType::Pointer cvExternalTerm0 = ChanAndVeseExternalTermType::New();
141 cvExternalTerm0->SetInput( input );
142 cvExternalTerm0->SetCoefficient( 1.0 );
143 std::cout << "LevelSet 1: CV external term created" << std::endl;
144
145 // **************** CREATE ALL EQUATIONS ****************
146
147 // Create Term Container
148 TermContainerType::Pointer termContainer0 = TermContainerType::New();
149 termContainer0->SetInput( input );
150 termContainer0->SetCurrentLevelSetId( 0 );
151 termContainer0->SetLevelSetContainer( lscontainer );
152
153 termContainer0->AddTerm( 0, cvInternalTerm0 );
154 termContainer0->AddTerm( 1, cvExternalTerm0 );
155 std::cout << "Term container 0 created" << std::endl;
156
157 EquationContainerType::Pointer equationContainer = EquationContainerType::New();
158 equationContainer->AddEquation( 0, termContainer0 );
159 equationContainer->SetLevelSetContainer( lscontainer );
160
161 using StoppingCriterionType = itk::LevelSetEvolutionNumberOfIterationsStoppingCriterion<LevelSetContainerType>;
162 StoppingCriterionType::Pointer criterion = StoppingCriterionType::New();
163 criterion->SetNumberOfIterations( std::stoi( argv[2]) );
164
165 LevelSetEvolutionType::Pointer evolution = LevelSetEvolutionType::New();
166 evolution->SetEquationContainer( equationContainer );
167 evolution->SetStoppingCriterion( criterion );
168 evolution->SetLevelSetContainer( lscontainer );
169 evolution->SetNumberOfWorkUnits( 1 );
170
171 try
172 {
173 evolution->Update();
174 }
175 catch ( itk::ExceptionObject& err )
176 {
177 std::cerr << err << std::endl;
178 return EXIT_FAILURE;
179 }
180
181 using OutputImageType = itk::Image< signed char, Dimension >;
182 OutputImageType::Pointer outputImage = OutputImageType::New();
183 outputImage->SetRegions( input->GetLargestPossibleRegion() );
184 outputImage->CopyInformation( input );
185 outputImage->Allocate();
186 outputImage->FillBuffer( 0 );
187
188 using OutputIteratorType = itk::ImageRegionIteratorWithIndex< OutputImageType >;
189 OutputIteratorType oIt( outputImage, outputImage->GetLargestPossibleRegion() );
190 oIt.GoToBegin();
191
192 OutputImageType::IndexType idx;
193
194 while( !oIt.IsAtEnd() )
195 {
196 idx = oIt.GetIndex();
197 oIt.Set( level_set->GetLabelMap()->GetPixel(idx) );
198 ++oIt;
199 }
200
201 using OutputWriterType = itk::ImageFileWriter< OutputImageType >;
202 OutputWriterType::Pointer writer = OutputWriterType::New();
203 writer->SetFileName( argv[3] );
204 writer->SetInput( outputImage );
205
206 try
207 {
208 writer->Update();
209 }
210 catch ( itk::ExceptionObject& err )
211 {
212 std::cout << err << std::endl;
213 }
214
215 if ( evolution->GetNumberOfWorkUnits() != 1 )
216 {
217 return EXIT_FAILURE;
218 }
219
220 return EXIT_SUCCESS;
221 }
222