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 itkDifferenceOfGaussiansGradientImageFilter_hxx
19 #define itkDifferenceOfGaussiansGradientImageFilter_hxx
20
21 #include <cmath>
22 #include "itkProgressReporter.h"
23 #include "itkDifferenceOfGaussiansGradientImageFilter.h"
24 #include "itkImageRegionIterator.h"
25
26 namespace itk
27 {
28 template< typename TInputImage, typename TDataType >
29 DifferenceOfGaussiansGradientImageFilter< TInputImage, TDataType >
DifferenceOfGaussiansGradientImageFilter()30 ::DifferenceOfGaussiansGradientImageFilter()
31 {
32 itkDebugMacro(<< "DifferenceOfGaussiansGradientImageFilter::DifferenceOfGaussiansGradientImageFilter() called");
33
34 m_Width = 2;
35 }
36
37 template< typename TInputImage, typename TDataType >
38 void
39 DifferenceOfGaussiansGradientImageFilter< TInputImage, TDataType >
GenerateData()40 ::GenerateData()
41 {
42 itkDebugMacro(<< "DifferenceOfGaussiansGradientImageFilter::GenerateData() called");
43
44 // Get the input and output pointers
45 typename Superclass::InputImagePointer inputPtr =
46 const_cast< TInputImage * >( this->GetInput(0) );
47 typename Superclass::OutputImagePointer outputPtr = this->GetOutput(0);
48
49 // Make sure we're getting everything
50 inputPtr->SetRequestedRegionToLargestPossibleRegion();
51
52 // How big is the input image?
53 typename TInputImage::SizeType size = inputPtr->GetLargestPossibleRegion().GetSize();
54
55 // Create a region object native to the output image type
56 OutputImageRegionType outputRegion;
57
58 // Resize the output region
59 outputRegion.SetSize(size);
60
61 // Set the largest legal region size (i.e. the size of the whole image)
62 // to what we just defined
63 outputPtr->SetBufferedRegion(outputRegion);
64 outputPtr->Allocate();
65
66 // Create a progress reporter
67 ProgressReporter progress( this, 0, outputPtr->GetRequestedRegion().GetNumberOfPixels() );
68
69 // Create an iterator that will walk the output region
70 using OutputIterator = ImageRegionIterator< TOutputImage >;
71
72 OutputIterator outIt = OutputIterator( outputPtr,
73 outputPtr->GetRequestedRegion() );
74
75 // Define a few indices that will be used to translate from an input pixel
76 // to an output pixel
77 typename TOutputImage::IndexType outputIndex;
78 typename TOutputImage::IndexType upperIndex;
79 typename TOutputImage::IndexType lowerIndex;
80
81 // walk the output image, and sample the input image
82 for (; !outIt.IsAtEnd(); ++outIt )
83 {
84 // determine the index of the output pixel
85 outputIndex = outIt.GetIndex();
86
87 // is the current index an acceptable distance from the edges
88 // of the image?
89 bool isValidGrad = true;
90
91 for ( unsigned int i = 0; i < NDimensions; ++i )
92 {
93 const auto width = static_cast< int >( m_Width );
94 const int sizeDifference =
95 static_cast< int >( size.m_InternalArray[i] ) - width;
96
97 if ( !( ( outputIndex[i] < sizeDifference )
98 && ( outputIndex[i] >= width ) ) )
99 {
100 isValidGrad = false;
101 }
102 }
103
104 if ( isValidGrad )
105 {
106 // We're in a safe position, so calculate the gradient for
107 // each dimension
108 for ( unsigned int i = 0; i < NDimensions; i++ )
109 {
110 // Build the indices for each pixel
111 for ( unsigned int j = 0; j < NDimensions; j++ )
112 {
113 if ( j == i )
114 {
115 upperIndex[j] = outputIndex[j] + static_cast< IndexValueType >( m_Width );
116 lowerIndex[j] = outputIndex[j] - static_cast< IndexValueType >( m_Width );
117 }
118 else
119 {
120 upperIndex[j] = outputIndex[j];
121 lowerIndex[j] = outputIndex[j];
122 }
123 }
124 // Remember, output type is a covariant vector of TDataType
125 outputPtr->GetPixel(outputIndex)[i] =
126 inputPtr->GetPixel(upperIndex) - inputPtr->GetPixel(lowerIndex);
127 }
128 }
129 else // We're not in a safe position, gradient is zero
130 {
131 for ( unsigned int i = 0; i < NDimensions; ++i )
132 {
133 outputPtr->GetPixel(outputIndex)[i] = 0.0;
134 }
135 }
136 progress.CompletedPixel();
137 }
138
139 itkDebugMacro(<< "DifferenceOfGaussiansGradientImageFilter::GenerateData() finished");
140 }
141
142 template< typename TInputImage, typename TDataType >
143 void
144 DifferenceOfGaussiansGradientImageFilter< TInputImage, TDataType >
PrintSelf(std::ostream & os,Indent indent) const145 ::PrintSelf(std::ostream & os, Indent indent) const
146 {
147 Superclass::PrintSelf(os, indent);
148
149 os << indent << "Width is " << m_Width << std::endl;
150 }
151 } // end namespace
152
153 #endif
154