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 #ifndef itkVoronoiSegmentationImageFilter_hxx
19 #define itkVoronoiSegmentationImageFilter_hxx
20 #include "itkVoronoiSegmentationImageFilter.h"
21
22 #include "itkImageRegionIteratorWithIndex.h"
23
24 namespace itk
25 {
26
27 template< typename TInputImage, typename TOutputImage, typename TBinaryPriorImage >
28 void
29 VoronoiSegmentationImageFilter< TInputImage, TOutputImage, TBinaryPriorImage >
SetMeanPercentError(double x)30 ::SetMeanPercentError(double x)
31 {
32 m_MeanPercentError = x;
33 m_MeanTolerance = x * m_Mean;
34 }
35
36 template< typename TInputImage, typename TOutputImage, typename TBinaryPriorImage >
37 void
38 VoronoiSegmentationImageFilter< TInputImage, TOutputImage, TBinaryPriorImage >
SetSTDPercentError(double x)39 ::SetSTDPercentError(double x)
40 {
41 m_STDPercentError = x;
42 m_STDTolerance = x * m_STD;
43 }
44
45 template< typename TInputImage, typename TOutputImage, typename TBinaryPriorImage >
46 bool
47 VoronoiSegmentationImageFilter< TInputImage, TOutputImage, TBinaryPriorImage >
TestHomogeneity(IndexList & Plist)48 ::TestHomogeneity(IndexList & Plist)
49 {
50 auto num = static_cast< int >( Plist.size() );
51 int i;
52 double getp;
53 double addp = 0;
54 double addpp = 0;
55
56 const InputImageType * inputImage = this->GetInput();
57
58 for ( i = 0; i < num; i++ )
59 {
60 getp = (double)( inputImage->GetPixel(Plist[i]) );
61 addp = addp + getp;
62 addpp = addpp + getp * getp;
63 }
64
65 double savemean, saveSTD;
66 if ( num > 1 )
67 {
68 savemean = addp / num;
69 saveSTD = std::sqrt( ( addpp - ( addp * addp ) / ( num ) ) / ( num - 1 ) );
70 }
71 else
72 {
73 savemean = 0;
74 saveSTD = -1;
75 }
76
77 // // jvm - Mahalanobis distance
78 // if (savevar > 0 && std::fabs(savemean - m_Mean) / m_Var < 2.5)
79 // return true;
80 // else
81 // return false;
82
83 savemean -= m_Mean;
84 saveSTD -= m_STD;
85 if ( ( savemean > -m_MeanTolerance ) && ( savemean < m_MeanTolerance )
86 && ( saveSTD < m_STDTolerance ) )
87 {
88 return true;
89 }
90 else
91 {
92 return false;
93 }
94 }
95
96 template< typename TInputImage, typename TOutputImage, typename TBinaryPriorImage >
97 void
98 VoronoiSegmentationImageFilter< TInputImage, TOutputImage, TBinaryPriorImage >
TakeAPrior(const BinaryObjectImage * aprior)99 ::TakeAPrior(const BinaryObjectImage *aprior)
100 {
101 RegionType region = this->GetInput()->GetRequestedRegion();
102
103 itk::ImageRegionConstIteratorWithIndex< BinaryObjectImage > ait(aprior, region);
104 itk::ImageRegionConstIteratorWithIndex< InputImageType > iit(this->GetInput(), region);
105
106 this->m_Size = this->GetInput()->GetRequestedRegion().GetSize();
107
108 int num = 0;
109 float addp = 0;
110 float addpp = 0;
111 float currp;
112
113 unsigned int i, j;
114 unsigned int minx = 0, miny = 0, maxx = 0, maxy = 0;
115 bool status = false;
116 for ( i = 0; i < this->m_Size[1]; i++ )
117 {
118 for ( j = 0; j < this->m_Size[0]; j++ )
119 {
120 if ( ( status == 0 ) && ( ait.Get() ) )
121 {
122 miny = i;
123 minx = j;
124 maxy = i;
125 maxx = j;
126 status = true;
127 }
128 else if ( ( status == 1 ) && ( ait.Get() ) )
129 {
130 maxy = i;
131 if ( minx > j ) { minx = j; }
132 if ( maxx < j ) { maxx = j; }
133 }
134 ++ait;
135 }
136 }
137
138 float addb = 0;
139 float addbb = 0;
140 int numb = 0;
141
142 ait.GoToBegin();
143 iit.GoToBegin();
144 for ( i = 0; i < miny; i++ )
145 {
146 for ( j = 0; j < this->m_Size[0]; j++ )
147 {
148 ++ait;
149 ++iit;
150 }
151 }
152 for ( i = miny; i <= maxy; i++ )
153 {
154 for ( j = 0; j < minx; j++ )
155 {
156 ++ait;
157 ++iit;
158 }
159 for ( j = minx; j <= maxx; j++ )
160 {
161 if ( ait.Get() )
162 {
163 num++;
164 currp = (float)( iit.Get() );
165 addp += currp;
166 addpp += currp * currp;
167 }
168 else
169 {
170 numb++;
171 currp = (float)( iit.Get() );
172 addb += currp;
173 addbb += currp * currp;
174 }
175 ++ait; ++iit;
176 }
177 for ( j = maxx + 1; j < this->m_Size[0]; j++ )
178 {
179 ++ait;
180 ++iit;
181 }
182 }
183
184 m_Mean = addp / num;
185 m_STD = std::sqrt( ( addpp - ( addp * addp ) / num ) / ( num - 1 ) );
186 float b_Mean = addb / numb;
187
188 if ( this->GetUseBackgroundInAPrior() )
189 {
190 m_MeanTolerance = std::fabs(m_Mean - b_Mean) * this->GetMeanDeviation();
191 }
192 else
193 {
194 m_MeanTolerance = m_Mean * m_MeanPercentError;
195 }
196 m_STDTolerance = m_STD * m_STDPercentError;
197 }
198
199 template< typename TInputImage, typename TOutputImage, typename TBinaryPriorImage >
200 void
201 VoronoiSegmentationImageFilter< TInputImage, TOutputImage, TBinaryPriorImage >
PrintSelf(std::ostream & os,Indent indent) const202 ::PrintSelf(std::ostream & os, Indent indent) const
203 {
204 Superclass::PrintSelf(os, indent);
205
206 os << indent << "Mean = " << m_Mean << std::endl;
207 os << indent << "MeanTolerance = " << m_MeanTolerance << std::endl;
208 os << indent << "MeanPercentError = " << m_MeanPercentError << std::endl;
209 os << indent << "STD = " << m_STD << std::endl;
210 os << indent << "STDTolerance = " << m_STDTolerance << std::endl;
211 os << indent << "STDPercentError = " << m_STDPercentError << std::endl;
212 }
213 } //end namespace
214
215 #endif
216