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 itkRenyiEntropyThresholdCalculator_hxx
20 #define itkRenyiEntropyThresholdCalculator_hxx
21
22 #include "itkRenyiEntropyThresholdCalculator.h"
23 #include "itkProgressReporter.h"
24 #include "itkMath.h"
25
26 namespace itk
27 {
28
29 template<typename THistogram, typename TOutput>
30 void
31 RenyiEntropyThresholdCalculator<THistogram, TOutput>
GenerateData()32 ::GenerateData()
33 {
34 const HistogramType * histogram = this->GetInput();
35
36 TotalAbsoluteFrequencyType total = histogram->GetTotalFrequency();
37 if( total == NumericTraits< TotalAbsoluteFrequencyType >::ZeroValue() )
38 {
39 itkExceptionMacro(<< "Histogram is empty");
40 }
41 m_Size = histogram->GetSize( 0 );
42 ProgressReporter progress( this, 0, m_Size );
43 if( m_Size == 1 )
44 {
45 this->GetOutput()->Set( static_cast<OutputType>(histogram->GetMeasurement( 0, 0 )) );
46 return;
47 }
48
49 const double tolerance = itk::Math::eps;
50
51 InstanceIdentifier ih;
52
53 std::vector<double> norm_histo(m_Size); /* normalized histogram */
54 std::vector<double> P1(m_Size); /* cumulative normalized histogram */
55 std::vector<double> P2(m_Size);
56
57 for( ih = 0; ih < m_Size; ih++ )
58 {
59 norm_histo[ih] = static_cast< double >( histogram->GetFrequency(ih, 0) ) / static_cast< double >( total );
60 }
61
62 P1[0] = norm_histo[0];
63 P2[0] = 1.0 - P1[0];
64 for( ih = 1; ih < m_Size; ih++ )
65 {
66 P1[ih] = P1[ih-1] + norm_histo[ih];
67 P2[ih] = 1.0 - P1[ih];
68 }
69
70 // Determine the first non-zero bin
71 m_FirstBin = 0;
72 for( ih = 0; ih < m_Size; ih++ )
73 {
74 if( !( std::abs( P1[ih] ) < tolerance ) )
75 {
76 m_FirstBin = ih;
77 break;
78 }
79 }
80
81 // Determine the last non-zero bin
82 m_LastBin = static_cast< InstanceIdentifier >( m_Size - 1 );
83 for( ih = m_Size - 1; ih >= m_FirstBin; ih-- )
84 {
85 if( !( std::abs( P2[ih] ) < tolerance ) )
86 {
87 m_LastBin = ih;
88 break;
89 }
90 }
91
92 InstanceIdentifier t_star2 = this->MaxEntropyThresholding( histogram, norm_histo, P1, P2 );
93 InstanceIdentifier t_star1 = this->MaxEntropyThresholding2( histogram, norm_histo, P1, P2 );
94 InstanceIdentifier t_star3 = this->MaxEntropyThresholding3( histogram, norm_histo, P1, P2 );
95
96 InstanceIdentifier tmp_var;
97
98 // Sort t_star values
99 if( t_star2 < t_star1 )
100 {
101 tmp_var = t_star1;
102 t_star1 = t_star2;
103 t_star2 = tmp_var;
104 }
105 if( t_star3 < t_star2 )
106 {
107 tmp_var = t_star2;
108 t_star2 = t_star3;
109 t_star3 = tmp_var;
110 }
111 if( t_star2 < t_star1 )
112 {
113 tmp_var = t_star1;
114 t_star1 = t_star2;
115 t_star2 = tmp_var;
116 }
117
118 double beta1 = 0.;
119 double beta2 = 0.;
120 double beta3 = 0.;
121
122 // Adjust beta values.
123 // Note that t_star1, t_star2, t_star3 are unsigned.
124 if( std::abs( static_cast< double >( t_star1 ) - static_cast< double >( t_star2 ) ) <= 5. )
125 {
126 if( std::abs( static_cast< double >( t_star2 ) - static_cast< double >( t_star3 ) ) <= 5. )
127 {
128 beta1 = 1.;
129 beta2 = 2.;
130 beta3 = 1.;
131 }
132 else
133 {
134 beta1 = 0.;
135 beta2 = 1.;
136 beta3 = 3.;
137 }
138 }
139 else
140 {
141 if( std::abs( static_cast< double >( t_star2 ) - static_cast< double >( t_star3 ) ) <= 5. )
142 {
143 beta1 = 3.;
144 beta2 = 1.;
145 beta3 = 0.;
146 }
147 else
148 {
149 beta1 = 1.;
150 beta2 = 2.;
151 beta3 = 1.;
152 }
153 }
154
155 itkAssertInDebugAndIgnoreInReleaseMacro( t_star1 < m_Size );
156 itkAssertInDebugAndIgnoreInReleaseMacro( t_star2 < m_Size );
157 itkAssertInDebugAndIgnoreInReleaseMacro( t_star3 < m_Size );
158
159 double omega = P1[t_star3] - P1[t_star1];
160
161 // Determine the optimal threshold value
162 double realOptThreshold = static_cast< double >( t_star1 ) * ( P1[t_star1]+ 0.25 * omega * beta1 ) +
163 static_cast< double >( t_star2 ) * 0.25 * omega * beta2 +
164 static_cast< double >( t_star3 ) * ( P2[t_star3] + 0.25 * omega * beta3 );
165
166 auto opt_threshold = static_cast< InstanceIdentifier >( realOptThreshold );
167
168 this->GetOutput()->Set( static_cast<OutputType>( histogram->GetMeasurement( opt_threshold, 0 ) ) );
169 }
170
171 template<typename THistogram, typename TOutput>
172 typename RenyiEntropyThresholdCalculator<THistogram, TOutput>::InstanceIdentifier
173 RenyiEntropyThresholdCalculator<THistogram, TOutput>
MaxEntropyThresholding(const HistogramType * histogram,const std::vector<double> & normHisto,const std::vector<double> & P1,const std::vector<double> & P2)174 ::MaxEntropyThresholding( const HistogramType* histogram,
175 const std::vector< double >& normHisto,
176 const std::vector< double >& P1,
177 const std::vector< double >& P2 )
178 {
179
180 // Calculate the total entropy each gray-level and fin the threshold that
181 // maximizes it
182
183 InstanceIdentifier threshold = 0; // was MIN_INT in original code, but if an empty image is processed it gives an error later on.
184 double max_ent = NumericTraits< double >::min(); // max entropy
185
186 for( InstanceIdentifier it = m_FirstBin; it <= m_LastBin; it++ )
187 {
188 // Entropy of the background pixels
189 double ent_back = 0.0;
190 for( InstanceIdentifier ih = 0; ih <= it; ih++ )
191 {
192 if( histogram->GetFrequency(ih, 0) != NumericTraits< AbsoluteFrequencyType >::ZeroValue() )
193 {
194 double x = ( normHisto[ih] / P1[it] );
195 ent_back -= x * std::log ( x );
196 }
197 }
198
199 // Entropy of the object pixels
200 double ent_obj = 0.0;
201 for( InstanceIdentifier ih = it + 1; ih < m_Size; ih++ )
202 {
203 if( histogram->GetFrequency(ih, 0) != NumericTraits< AbsoluteFrequencyType >::ZeroValue() )
204 {
205 double x = ( normHisto[ih] / P2[it] );
206 ent_obj -= x * std::log( x );
207 }
208 }
209
210 // Total entropy
211 double tot_ent = ent_back + ent_obj;
212
213 // IJ.log(""+max_ent+" "+tot_ent);
214
215 if( max_ent < tot_ent )
216 {
217 max_ent = tot_ent;
218 threshold = it;
219 }
220 }
221 return threshold;
222 }
223
224 template<typename THistogram, typename TOutput>
225 typename RenyiEntropyThresholdCalculator<THistogram, TOutput>::InstanceIdentifier
226 RenyiEntropyThresholdCalculator<THistogram, TOutput>
MaxEntropyThresholding2(const HistogramType * itkNotUsed (histogram),const std::vector<double> & normHisto,const std::vector<double> & P1,const std::vector<double> & P2)227 ::MaxEntropyThresholding2( const HistogramType* itkNotUsed( histogram ),
228 const std::vector< double >& normHisto,
229 const std::vector< double >& P1,
230 const std::vector< double >& P2 )
231 {
232
233 InstanceIdentifier threshold = 0; //was MIN_INT in original code, but if an empty image is processed it gives an error later on.
234 double max_ent = NumericTraits< double >::min();
235 double alpha = 0.5;
236 double term = 1.0 / ( 1.0 - alpha );
237
238 for( InstanceIdentifier it = m_FirstBin; it <= m_LastBin; it++ )
239 {
240 // Entropy of the background pixels
241 double ent_back = 0.0;
242 for( InstanceIdentifier ih = 0; ih <= it; ih++ )
243 {
244 ent_back += std::sqrt( normHisto[ih] / P1[it] );
245 }
246
247 // Entropy of the object pixel
248 double ent_obj = 0.0;
249 for( InstanceIdentifier ih = it + 1; ih < m_Size; ih++ )
250 {
251 ent_obj += std::sqrt( normHisto[ih] / P2[it] );
252 }
253
254 // Total entropy
255 double product = ent_back * ent_obj;
256 double tot_ent = 0.;
257
258 if( product > 0.0 )
259 {
260 tot_ent = term * std::log( ent_back * ent_obj );
261 }
262
263 if( tot_ent > max_ent )
264 {
265 max_ent = tot_ent;
266 threshold = it;
267 }
268 }
269
270 return threshold;
271 }
272
273 template<typename THistogram, typename TOutput>
274 typename RenyiEntropyThresholdCalculator<THistogram, TOutput>::InstanceIdentifier
275 RenyiEntropyThresholdCalculator<THistogram, TOutput>
MaxEntropyThresholding3(const HistogramType * itkNotUsed (histogram),const std::vector<double> & normHisto,const std::vector<double> & P1,const std::vector<double> & P2)276 ::MaxEntropyThresholding3( const HistogramType* itkNotUsed( histogram ),
277 const std::vector< double >& normHisto,
278 const std::vector< double >& P1,
279 const std::vector< double >& P2 )
280 {
281 InstanceIdentifier threshold = 0; //was MIN_INT in original code, but if an empty image is processed it gives an error later on.
282 double max_ent = 0.0;
283 double alpha = 2.0;
284 double term = 1.0 / ( 1.0 - alpha );
285 for( InstanceIdentifier it = m_FirstBin; it <= m_LastBin; it++ )
286 {
287 // Entropy of the background pixels
288 double ent_back = 0.0;
289 for( InstanceIdentifier ih = 0; ih <= it; ih++ )
290 {
291 double x = normHisto[ih] / P1[it];
292 ent_back += x * x;
293 }
294
295 // Entropy of the object pixels
296 double ent_obj = 0.0;
297 for( InstanceIdentifier ih = it + 1; ih < m_Size; ih++ )
298 {
299 double x = normHisto[ih] / P2[it];
300 ent_obj += x * x;
301 }
302
303 // Total entropy
304 double tot_ent = 0.0;
305 double product = ent_back * ent_obj;
306 if( product > 0.0 )
307 {
308 tot_ent = term * std::log( product );
309 }
310
311 if( tot_ent > max_ent )
312 {
313 max_ent = tot_ent;
314 threshold = it;
315 }
316 }
317
318 return threshold;
319 }
320
321 template<typename THistogram, typename TOutput>
322 void
323 RenyiEntropyThresholdCalculator<THistogram, TOutput>
PrintSelf(std::ostream & os,Indent indent) const324 ::PrintSelf( std::ostream& os, Indent indent ) const
325 {
326 Superclass::PrintSelf(os, indent);
327
328 os << indent << "FirstBin: " << static_cast< typename itk::NumericTraits<
329 InstanceIdentifier >::PrintType >( m_FirstBin ) << std::endl;
330 os << indent << "LastBin: " << static_cast< typename itk::NumericTraits<
331 InstanceIdentifier >::PrintType >( m_LastBin ) << std::endl;
332 os << indent << "Size: " << static_cast< typename itk::NumericTraits<
333 SizeValueType >::PrintType >( m_Size ) << std::endl;
334 }
335 } // end namespace itk
336
337 #endif
338