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