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 #ifndef itkFastMarchingBase_hxx
20 #define itkFastMarchingBase_hxx
21
22 #include "itkFastMarchingBase.h"
23
24 #include "itkProgressReporter.h"
25 #include "itkMath.h"
26 #include "itkMath.h"
27
28 namespace itk
29 {
30 // -----------------------------------------------------------------------------
31 template< typename TInput, typename TOutput >
32 FastMarchingBase< TInput, TOutput >::
FastMarchingBase()33 FastMarchingBase()
34 {
35 this->ProcessObject::SetNumberOfRequiredInputs(0);
36
37 m_TrialPoints = nullptr;
38 m_AlivePoints = nullptr;
39 m_ProcessedPoints = nullptr;
40 m_ForbiddenPoints = nullptr;
41
42 //m_Heap = PriorityQueueType::New();
43 m_SpeedConstant = 1.;
44 m_InverseSpeed = -1.;
45 m_NormalizationFactor = 1.;
46 m_TargetReachedValue = NumericTraits< OutputPixelType >::ZeroValue();
47 m_TopologyCheck = Nothing;
48 m_LargeValue = NumericTraits< OutputPixelType >::max();
49 m_TopologyValue = m_LargeValue;
50 m_CollectPoints = false;
51 }
52 // -----------------------------------------------------------------------------
53
54 // -----------------------------------------------------------------------------
55 template< typename TInput, typename TOutput >
56 void
57 FastMarchingBase< TInput, TOutput >::
PrintSelf(std::ostream & os,Indent indent) const58 PrintSelf( std::ostream & os, Indent indent ) const
59 {
60 Superclass::PrintSelf( os, indent );
61 os << indent << "Speed constant: " << m_SpeedConstant << std::endl;
62 os << indent << "Topology check: " << m_TopologyCheck << std::endl;
63 os << indent << "Normalization Factor: " << m_NormalizationFactor << std::endl;
64 }
65
66 // -----------------------------------------------------------------------------
67
68 // -----------------------------------------------------------------------------
69 template< typename TInput, typename TOutput >
70 void
71 FastMarchingBase< TInput, TOutput >::
Initialize(OutputDomainType * oDomain)72 Initialize( OutputDomainType* oDomain )
73 {
74 if( m_TrialPoints.IsNull() )
75 {
76 itkExceptionMacro( <<"No Trial Nodes" );
77 }
78 if( m_StoppingCriterion.IsNull() )
79 {
80 itkExceptionMacro( <<"No Stopping Criterion Set" );
81 }
82 if( m_NormalizationFactor < itk::Math::eps )
83 {
84 itkExceptionMacro( <<"Normalization Factor is null or negative" );
85 }
86 if( m_SpeedConstant < itk::Math::eps )
87 {
88 itkExceptionMacro( <<"SpeedConstant is null or negative" );
89 }
90 if( m_CollectPoints )
91 {
92 if( m_ProcessedPoints.IsNull() )
93 {
94 m_ProcessedPoints = NodePairContainerType::New();
95 }
96 }
97
98 // make sure the heap is empty
99 while ( !m_Heap.empty() )
100 {
101 m_Heap.pop();
102 }
103 /*
104 while ( !m_Heap->Empty() )
105 {
106 m_Heap->Pop();
107 }
108 */
109
110 this->InitializeOutput( oDomain );
111
112 // By setting the output domain to the stopping criterion, we enable funky
113 // criterion based on informations extracted from it
114 m_StoppingCriterion->SetDomain( oDomain );
115 }
116 // -----------------------------------------------------------------------------
117
118 // -----------------------------------------------------------------------------
119 template< typename TInput, typename TOutput >
120 void
121 FastMarchingBase< TInput, TOutput >::
GenerateData()122 GenerateData()
123 {
124 OutputDomainType* output = this->GetOutput();
125
126 Initialize( output );
127
128 OutputPixelType current_value = 0.;
129
130 ProgressReporter progress( this, 0, this->GetTotalNumberOfNodes() );
131
132 m_StoppingCriterion->Reinitialize();
133
134 try
135 {
136 //while( !m_Heap->Empty() )
137 while( !m_Heap.empty() )
138 {
139 //PriorityQueueElementType element = m_Heap->Peek();
140 //m_Heap->Pop();
141 //
142 //NodeType current_node = element.m_Element;
143 //OutputPixelType current_value = element.m_Priority;
144
145
146 NodePairType current_node_pair = m_Heap.top();
147 m_Heap.pop();
148
149 NodeType current_node = current_node_pair.GetNode();
150 current_value = this->GetOutputValue( output, current_node );
151
152 if( Math::ExactlyEquals(current_value, current_node_pair.GetValue()) )
153 {
154 // is this node already alive ?
155 if( this->GetLabelValueForGivenNode( current_node ) != Traits::Alive )
156 {
157 m_StoppingCriterion->SetCurrentNodePair( current_node_pair );
158
159 if( m_StoppingCriterion->IsSatisfied() )
160 {
161 break;
162 }
163
164 if( this->CheckTopology( output, current_node ) )
165 {
166 if ( m_CollectPoints )
167 {
168 m_ProcessedPoints->push_back( current_node_pair );
169 }
170
171 // set this node as alive
172 this->SetLabelValueForGivenNode( current_node, Traits::Alive );
173
174 // update its neighbors
175 this->UpdateNeighbors( output, current_node );
176 }
177 }
178 progress.CompletedPixel();
179 }
180 }
181 }
182 catch ( ProcessAborted & )
183 {
184 // User aborted filter execution Here we catch an exception thrown by the
185 // progress reporter and rethrow it with the correct line number and file
186 // name. We also invoke AbortEvent in case some observer was interested on
187 // it.
188 //
189 // RELEASE MEMORY!!!
190 while( !m_Heap.empty() )
191 {
192 m_Heap.pop();
193 }
194 /*while( !m_Heap->Empty() )
195 {
196 m_Heap->Pop();
197 }*/
198
199 throw ProcessAborted(__FILE__, __LINE__);
200 }
201
202 m_TargetReachedValue = current_value;
203
204 // let's release some useless memory...
205 while( !m_Heap.empty() )
206 {
207 m_Heap.pop();
208 }
209 /*while( !m_Heap->Empty() )
210 {
211 m_Heap->Pop();
212 }*/
213 }
214 // -----------------------------------------------------------------------------
215
216 } // end of namespace itk
217
218 #endif
219