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 itkNeighborhoodAlgorithm_hxx
19 #define itkNeighborhoodAlgorithm_hxx
20 #include "itkNeighborhoodAlgorithm.h"
21 #include "itkImageRegionIterator.h"
22 #include "itkImageRegion.h"
23 #include "itkConstSliceIterator.h"
24
25 namespace itk
26 {
27 namespace NeighborhoodAlgorithm
28 {
29 template< typename TImage >
30 typename ImageBoundaryFacesCalculator< TImage >::FaceListType
31 ImageBoundaryFacesCalculator< TImage >
operator ()(const TImage * img,RegionType regionToProcess,RadiusType radius)32 ::operator()(const TImage *img, RegionType regionToProcess, RadiusType radius)
33 {
34 unsigned int j, i;
35 // Analyze the regionToProcess to determine if any of its faces are
36 // along a buffer boundary (we have no data in the buffer for pixels
37 // that are outside the boundary, but within the neighborhood radius and will
38 // have to treat them differently). We also determine the size of the non-
39 // boundary region that will be processed. For instance, given a 2D image
40 // and regionTOProcess (size = 5x5),
41
42 FaceListType faceList;
43 // The portion of the regionToProcess that is outside of the image bufferedRegion
44 // doesn't make sense. If the regionToProcess is completely outside of
45 // the image bufferedRegion, return a empty region list.
46 if( !regionToProcess.Crop( img->GetBufferedRegion() ) )
47 {
48 return faceList;
49 }
50
51 const IndexType bStart = img->GetBufferedRegion().GetIndex();
52 const SizeType bSize = img->GetBufferedRegion().GetSize();
53 const IndexType rStart = regionToProcess.GetIndex();
54 const SizeType rSize = regionToProcess.GetSize();
55
56 IndexValueType overlapLow, overlapHigh;
57 IndexType fStart; // Boundary, "face"
58 SizeType fSize; // region data.
59 RegionType fRegion;
60 SizeType nbSize = regionToProcess.GetSize(); // Non-boundary region
61 IndexType nbStart = regionToProcess.GetIndex(); // data.
62 RegionType nbRegion;
63
64 IndexType vrStart = rStart; //start index of variable processed region which has considered the boundary region in last direction
65 SizeType vrSize = rSize; //size of variable processed region which has considered the boundary region in last direction
66
67 for ( i = 0; i < ImageDimension; ++i )
68 {
69 overlapLow = static_cast< IndexValueType >( ( rStart[i] - radius[i] ) - bStart[i] );
70 // image buffered region should be more than twice of the radius,
71 // otherwise there would be overlap between two boundary regions.
72 // in the case, we reduce upper boundary size to remove overlap.
73 if ( bSize[i] > 2 * radius[i] )
74 {
75 overlapHigh = static_cast< IndexValueType >( ( bStart[i] + bSize[i] ) - ( rStart[i] + rSize[i] + radius[i] ) );
76 }
77 else
78 {
79 overlapHigh = static_cast< IndexValueType >( ( bStart[i] + radius[i] ) - ( rStart[i] + rSize[i] ) );
80 }
81
82 if ( overlapLow < 0 ) // out of bounds condition, define
83 // a region of
84 { // iteration along this face
85 for ( j = 0; j < ImageDimension; ++j ) // define the starting index
86 { // and size of the face region
87 fStart[j] = vrStart[j];
88 if ( j == i )
89 {
90 // Boundary region cannot be outside the region to process
91 if( -overlapLow > static_cast< IndexValueType >( rSize[i] ) )
92 {
93 overlapLow = - static_cast< IndexValueType >( rSize[i] );
94 }
95 fSize[j] = -overlapLow;
96 vrSize[j] += overlapLow; //change start and size in this direction
97 vrStart[j] -= overlapLow;//to ensure no duplicate pixels at corners
98 }
99 else
100 {
101 fSize[j] = vrSize[j];
102 }
103
104 // Boundary region cannot be outside the region to process
105 if ( fSize[j] > rSize[j] )
106 {
107 fSize[j] = rSize[j];
108 }
109 }
110 // avoid unsigned overflow if the non-boundary region is too small to
111 // process
112 if ( fSize[i] > nbSize[i] )
113 {
114 nbSize[i] = 0;
115 }
116 else
117 {
118 nbSize[i] -= fSize[i];
119 }
120 nbStart[i] += -overlapLow;
121 fRegion.SetIndex(fStart);
122 fRegion.SetSize(fSize);
123 faceList.push_back(fRegion);
124 }
125 if ( overlapHigh < 0 )
126 {
127 for ( j = 0; j < ImageDimension; ++j )
128 {
129 if ( j == i )
130 {
131 // Boundary region cannot be outside the region to process
132 if( -overlapHigh > static_cast< IndexValueType >( rSize[i] ) )
133 {
134 overlapHigh = - static_cast< IndexValueType >( rSize[i] );
135 }
136 fStart[j] = rStart[j] + static_cast< IndexValueType >( rSize[j] ) + overlapHigh;
137 fSize[j] = -overlapHigh;
138 vrSize[j] += overlapHigh; //change size in this direction
139 }
140 else
141 {
142 fStart[j] = vrStart[j];
143 fSize[j] = vrSize[j];
144 }
145 }
146 // avoid unsigned overflow if the non-boundary region is too small to
147 // process
148 if ( fSize[i] > nbSize[i] )
149 {
150 nbSize[i] = 0;
151 }
152 else
153 {
154 nbSize[i] -= fSize[i];
155 }
156
157 fRegion.SetIndex(fStart);
158 fRegion.SetSize(fSize);
159 faceList.push_back(fRegion);
160 }
161 }
162 nbRegion.SetSize(nbSize);
163 nbRegion.SetIndex(nbStart);
164
165 faceList.push_front(nbRegion);
166 return faceList;
167 }
168
169 template< typename TImage >
170 typename CalculateOutputWrapOffsetModifiers< TImage >::OffsetType
171 CalculateOutputWrapOffsetModifiers< TImage >
operator ()(TImage * input,TImage * output) const172 ::operator()(TImage *input, TImage *output) const
173 {
174 OffsetType ans;
175
176 const Size< TImage::ImageDimension > isz = input->GetBufferedRegion().GetSize();
177 const Size< TImage::ImageDimension > osz = output->GetBufferedRegion().GetSize();
178
179 for ( int i = 0; i < TImage::ImageDimension; ++i )
180 {
181 ans[i] = osz[i] - isz[i];
182 }
183 return ans;
184 }
185 } // end namespace NeighborhoodAlgorithm
186 } // end namespace itk
187
188 #endif
189