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 itkBilateralImageFilter_hxx
19 #define itkBilateralImageFilter_hxx
20
21 #include "itkBilateralImageFilter.h"
22 #include "itkImageRegionIterator.h"
23 #include "itkGaussianImageSource.h"
24 #include "itkNeighborhoodAlgorithm.h"
25 #include "itkZeroFluxNeumannBoundaryCondition.h"
26 #include "itkProgressReporter.h"
27 #include "itkStatisticsImageFilter.h"
28
29 namespace itk
30 {
31 template< typename TInputImage, typename TOutputImage >
32 BilateralImageFilter< TInputImage, TOutputImage >
BilateralImageFilter()33 ::BilateralImageFilter()
34 {
35 this->m_Radius.Fill(1);
36 this->m_AutomaticKernelSize = true;
37 this->m_DomainSigma.Fill(4.0);
38 this->m_RangeSigma = 50.0;
39 this->m_FilterDimensionality = ImageDimension;
40 this->m_NumberOfRangeGaussianSamples = 100;
41 this->m_DynamicRange = 0.0;
42 this->m_DynamicRangeUsed = 0.0;
43 this->m_DomainMu = 2.5; // keep small to keep kernels small
44 this->m_RangeMu = 4.0; // can be bigger then DomainMu since we only
45 // index into a single table
46 this->DynamicMultiThreadingOn();
47 }
48
49 template< typename TInputImage, typename TOutputImage >
50 void
51 BilateralImageFilter< TInputImage, TOutputImage >
SetRadius(const SizeValueType i)52 ::SetRadius(const SizeValueType i)
53 {
54 m_Radius.Fill(i);
55 }
56
57 template< typename TInputImage, typename TOutputImage >
58 void
59 BilateralImageFilter< TInputImage, TOutputImage >
GenerateInputRequestedRegion()60 ::GenerateInputRequestedRegion()
61 {
62 // call the superclass' implementation of this method. this should
63 // copy the output requested region to the input requested region
64 Superclass::GenerateInputRequestedRegion();
65
66 // get pointers to the input and output
67 typename Superclass::InputImagePointer inputPtr =
68 const_cast< TInputImage * >( this->GetInput() );
69
70 if ( !inputPtr )
71 {
72 return;
73 }
74
75 // Pad the image by 2.5*sigma in all directions
76 typename TInputImage::SizeType radius;
77 unsigned int i;
78
79 if ( m_AutomaticKernelSize )
80 {
81 for ( i = 0; i < ImageDimension; i++ )
82 {
83 radius[i] =
84 ( typename TInputImage::SizeType::SizeValueType )
85 std::ceil(m_DomainMu * m_DomainSigma[i] / this->GetInput()->GetSpacing()[i]);
86 }
87 }
88 else
89 {
90 for ( i = 0; i < ImageDimension; i++ )
91 {
92 radius[i] = m_Radius[i];
93 }
94 }
95
96 // get a copy of the input requested region (should equal the output
97 // requested region)
98 typename TInputImage::RegionType inputRequestedRegion;
99 inputRequestedRegion = inputPtr->GetRequestedRegion();
100
101 // pad the input requested region by the operator radius
102 inputRequestedRegion.PadByRadius(radius);
103
104 // crop the input requested region at the input's largest possible region
105 if ( inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() ) )
106 {
107 inputPtr->SetRequestedRegion(inputRequestedRegion);
108 return;
109 }
110 else
111 {
112 // Couldn't crop the region (requested region is outside the largest
113 // possible region). Throw an exception.
114
115 // store what we tried to request (prior to trying to crop)
116 inputPtr->SetRequestedRegion(inputRequestedRegion);
117
118 // build an exception
119 InvalidRequestedRegionError e(__FILE__, __LINE__);
120 e.SetLocation(ITK_LOCATION);
121 e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
122 e.SetDataObject(inputPtr);
123 throw e;
124 }
125 }
126
127 template< typename TInputImage, typename TOutputImage >
128 void
129 BilateralImageFilter< TInputImage, TOutputImage >
BeforeThreadedGenerateData()130 ::BeforeThreadedGenerateData()
131 {
132 // Build a small image of the N-dimensional Gaussian used for domain filter
133 //
134 // Gaussian image size will be (2*std::ceil(2.5*sigma)+1) x
135 // (2*std::ceil(2.5*sigma)+1)
136 unsigned int i;
137
138 typename InputImageType::SizeType radius;
139 typename InputImageType::SizeType domainKernelSize;
140
141 const InputImageType *inputImage = this->GetInput();
142
143 const typename InputImageType::SpacingType inputSpacing = inputImage->GetSpacing();
144 const typename InputImageType::PointType inputOrigin = inputImage->GetOrigin();
145
146 if ( m_AutomaticKernelSize )
147 {
148 for ( i = 0; i < ImageDimension; i++ )
149 {
150 radius[i] =
151 ( typename TInputImage::SizeType::SizeValueType )
152 std::ceil(m_DomainMu * m_DomainSigma[i] / inputSpacing[i]);
153 domainKernelSize[i] = 2 * radius[i] + 1;
154 }
155 }
156 else
157 {
158 for ( i = 0; i < ImageDimension; i++ )
159 {
160 radius[i] = m_Radius[i];
161 domainKernelSize[i] = 2 * radius[i] + 1;
162 }
163 }
164
165 typename GaussianImageSource< GaussianImageType >::Pointer gaussianImage;
166 typename GaussianImageSource< GaussianImageType >::ArrayType mean;
167 typename GaussianImageSource< GaussianImageType >::ArrayType sigma;
168
169 gaussianImage = GaussianImageSource< GaussianImageType >::New();
170 gaussianImage->SetSize( domainKernelSize );
171 gaussianImage->SetSpacing(inputSpacing);
172 gaussianImage->SetOrigin(inputOrigin);
173 gaussianImage->SetScale(1.0);
174 gaussianImage->SetNormalized(true);
175
176 for ( i = 0; i < ImageDimension; i++ )
177 {
178 mean[i] = inputSpacing[i] * radius[i] + inputOrigin[i]; // center pixel pos
179 sigma[i] = m_DomainSigma[i];
180 }
181 gaussianImage->SetSigma(sigma);
182 gaussianImage->SetMean(mean);
183
184 gaussianImage->Update();
185
186 // copy this small Gaussian image into a neighborhood
187 m_GaussianKernel.SetRadius(radius);
188
189 KernelIteratorType kernel_it;
190 ImageRegionIterator< GaussianImageType > git =
191 ImageRegionIterator< GaussianImageType >( gaussianImage->GetOutput(),
192 gaussianImage->GetOutput()
193 ->GetBufferedRegion() );
194 double norm = 0.0;
195 for ( git.GoToBegin(); !git.IsAtEnd(); ++git )
196 {
197 norm += git.Get();
198 }
199 for ( git.GoToBegin(), kernel_it = m_GaussianKernel.Begin(); !git.IsAtEnd();
200 ++git, ++kernel_it )
201 {
202 *kernel_it = git.Get() / norm;
203 }
204
205 // Build a lookup table for the range gaussian
206
207 auto localInput = TInputImage::New();
208 localInput->Graft( this->GetInput() );
209
210 // First, determine the min and max intensity range
211 typename StatisticsImageFilter< TInputImage >::Pointer statistics =
212 StatisticsImageFilter< TInputImage >::New();
213
214 statistics->SetInput(localInput);
215 statistics->Update();
216
217 // Now create the lookup table whose domain runs from 0.0 to
218 // (max-min) and range is gaussian evaluated at
219 // that point
220 double rangeVariance = m_RangeSigma * m_RangeSigma;
221
222 // denominator (normalization factor) for Gaussian used for range
223 double rangeGaussianDenom;
224 rangeGaussianDenom = m_RangeSigma * std::sqrt(2.0 * itk::Math::pi);
225
226 // Maximum delta for the dynamic range
227 double tableDelta;
228 double v;
229
230 m_DynamicRange = ( static_cast< double >( statistics->GetMaximum() )
231 - static_cast< double >( statistics->GetMinimum() ) );
232
233 m_DynamicRangeUsed = m_RangeMu * m_RangeSigma;
234
235 tableDelta = m_DynamicRangeUsed
236 / static_cast< double >( m_NumberOfRangeGaussianSamples );
237
238 // Finally, build the table
239 m_RangeGaussianTable.resize(m_NumberOfRangeGaussianSamples);
240 for ( i = 0, v = 0.0; i < m_NumberOfRangeGaussianSamples;
241 ++i, v += tableDelta )
242 {
243 m_RangeGaussianTable[i] = std::exp(-0.5 * v * v / rangeVariance) / rangeGaussianDenom;
244 }
245 }
246
247 template< typename TInputImage, typename TOutputImage >
248 void
249 BilateralImageFilter< TInputImage, TOutputImage >
DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread)250 ::DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread)
251 {
252 typename TInputImage::ConstPointer input = this->GetInput();
253 typename TOutputImage::Pointer output = this->GetOutput();
254 typename TInputImage::IndexValueType i;
255 const double rangeDistanceThreshold = m_DynamicRangeUsed;
256
257 ZeroFluxNeumannBoundaryCondition< TInputImage > BC;
258
259 // Find the boundary "faces"
260 typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< InputImageType >::FaceListType faceList;
261 NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< InputImageType > fC;
262 faceList = fC( this->GetInput(), outputRegionForThread,
263 m_GaussianKernel.GetRadius() );
264
265 typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< InputImageType >::FaceListType::iterator fit;
266
267 OutputPixelRealType centerPixel;
268 OutputPixelRealType val, tableArg, normFactor, rangeGaussian,
269 rangeDistance, pixel, gaussianProduct;
270
271 const double distanceToTableIndex =
272 static_cast< double >( m_NumberOfRangeGaussianSamples ) / m_DynamicRangeUsed;
273
274 // Process all the faces, the NeighborhoodIterator will deteremine
275 // whether a specified region needs to use the boundary conditions or
276 // not.
277 NeighborhoodIteratorType b_iter;
278 ImageRegionIterator< OutputImageType > o_iter;
279 KernelConstIteratorType k_it;
280 KernelConstIteratorType kernelEnd = m_GaussianKernel.End();
281
282 for ( fit = faceList.begin(); fit != faceList.end(); ++fit )
283 {
284 // walk the boundary face and the corresponding section of the output
285 b_iter = NeighborhoodIteratorType(m_GaussianKernel.GetRadius(),
286 this->GetInput(), *fit);
287 b_iter.OverrideBoundaryCondition(&BC);
288 o_iter = ImageRegionIterator< OutputImageType >(this->GetOutput(), *fit);
289
290 while ( !b_iter.IsAtEnd() )
291 {
292 // Setup
293 centerPixel = static_cast< OutputPixelRealType >( b_iter.GetCenterPixel() );
294 val = 0.0;
295 normFactor = 0.0;
296
297 // Walk the neighborhood of the input and the kernel
298 for ( i = 0, k_it = m_GaussianKernel.Begin(); k_it < kernelEnd;
299 ++k_it, ++i )
300 {
301 // range distance between neighborhood pixel and neighborhood center
302 pixel = static_cast< OutputPixelRealType >( b_iter.GetPixel(i) );
303 rangeDistance = pixel - centerPixel;
304 // flip sign if needed
305 if ( rangeDistance < 0.0 )
306 {
307 rangeDistance *= -1.0;
308 }
309
310 // if the range distance is close enough, then use the pixel
311 if ( rangeDistance < rangeDistanceThreshold )
312 {
313 // look up the range gaussian in a table
314 tableArg = rangeDistance * distanceToTableIndex;
315 rangeGaussian = m_RangeGaussianTable[Math::Floor < SizeValueType > ( tableArg )];
316
317 // normalization factor so filter integrates to one
318 // (product of the domain and the range gaussian)
319 gaussianProduct = ( *k_it ) * rangeGaussian;
320 normFactor += gaussianProduct;
321
322 // Input Image * Domain Gaussian * Range Gaussian
323 val += pixel * gaussianProduct;
324 }
325 }
326 // normalize the value
327 val /= normFactor;
328
329 // store the filtered value
330 o_iter.Set( static_cast< OutputPixelType >( val ) );
331
332 ++b_iter;
333 ++o_iter;
334 }
335 }
336 }
337
338 template< typename TInputImage, typename TOutputImage >
339 void
340 BilateralImageFilter< TInputImage, TOutputImage >
PrintSelf(std::ostream & os,Indent indent) const341 ::PrintSelf(std::ostream & os, Indent indent) const
342 {
343 Superclass::PrintSelf(os, indent);
344
345 os << indent << "DomainSigma: " << m_DomainSigma << std::endl;
346 os << indent << "RangeSigma: " << m_RangeSigma << std::endl;
347 os << indent << "FilterDimensionality: " << m_FilterDimensionality << std::endl;
348 os << indent << "NumberOfRangeGaussianSamples: " << m_NumberOfRangeGaussianSamples << std::endl;
349 os << indent << "Input dynamic range: " << m_DynamicRange << std::endl;
350 os << indent << "Amount of dynamic range used: " << m_DynamicRangeUsed << std::endl;
351 os << indent << "AutomaticKernelSize: " << m_AutomaticKernelSize << std::endl;
352 os << indent << "Radius: " << m_Radius << std::endl;
353 }
354 } // end namespace itk
355
356 #endif
357