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 itkPeriodicBoundaryCondition_hxx
19 #define itkPeriodicBoundaryCondition_hxx
20
21 #include "itkConstNeighborhoodIterator.h"
22 #include "itkPeriodicBoundaryCondition.h"
23
24 namespace itk
25 {
26 template< typename TInputImage, typename TOutputImage >
27 typename PeriodicBoundaryCondition< TInputImage, TOutputImage >::OutputPixelType
28 PeriodicBoundaryCondition< TInputImage, TOutputImage >
operator ()(const OffsetType & point_index,const OffsetType & boundary_offset,const NeighborhoodType * data) const29 ::operator()(const OffsetType & point_index, const OffsetType & boundary_offset,
30 const NeighborhoodType *data) const
31 {
32 // This is guaranteed to be called with an object that is using
33 // PeriodicBoundaryCondition
34 const auto * iterator = reinterpret_cast< const ConstNeighborhoodIterator< TInputImage, Self > * >( data );
35 typename TInputImage::PixelType * ptr;
36 int linear_index = 0;
37 unsigned int i;
38
39 // Find the pointer of the closest boundary pixel
40
41 // Return the value of the pixel at the closest boundary point.
42 for ( i = 0; i < ImageDimension; ++i )
43 {
44 linear_index += ( point_index[i] + boundary_offset[i] ) * data->GetStride(i);
45 }
46 // (data->operator[](linear_index)) is guaranteed to be a pointer to
47 // TInputImage::PixelType except for VectorImage, in which case, it will be a
48 // pointer to TInputImage::InternalPixelType.
49 ptr = reinterpret_cast< PixelType * >( ( data->operator[](linear_index) ) );
50
51 // Wrap the pointer around the image in the necessary dimensions. If we have
52 // reached this point, we can assume that we are on the edge of the BUFFERED
53 // region of the image. Boundary conditions are only invoked if touching the
54 // actual memory boundary.
55
56 // These are the step sizes for increments in each dimension of the image.
57 const typename TInputImage::OffsetValueType * offset_table =
58 iterator->GetImagePointer()->GetOffsetTable();
59
60 for ( i = 0; i < ImageDimension; ++i )
61 {
62 if ( boundary_offset[i] != 0 )
63 { // If the neighborhood overlaps on the low edge, then wrap from the
64 // high edge of the image.
65 if ( point_index[i] < static_cast< OffsetValueType >( iterator->GetRadius(i) ) )
66 {
67 ptr += iterator->GetImagePointer()->GetBufferedRegion().GetSize()[i]
68 * offset_table[i] - boundary_offset[i] * offset_table[i];
69 }
70 else // wrap from the low side of the image
71 {
72 ptr -= iterator->GetImagePointer()->GetBufferedRegion().GetSize()[i]
73 * offset_table[i] + boundary_offset[i] * offset_table[i];
74 }
75 }
76 }
77
78 return static_cast< OutputPixelType >( *ptr );
79 }
80
81 template< typename TInputImage, typename TOutputImage >
82 typename PeriodicBoundaryCondition< TInputImage, TOutputImage >::OutputPixelType
83 PeriodicBoundaryCondition< TInputImage, TOutputImage >
operator ()(const OffsetType & point_index,const OffsetType & boundary_offset,const NeighborhoodType * data,const NeighborhoodAccessorFunctorType & neighborhoodAccessorFunctor) const84 ::operator()(const OffsetType & point_index, const OffsetType & boundary_offset,
85 const NeighborhoodType *data,
86 const NeighborhoodAccessorFunctorType & neighborhoodAccessorFunctor) const
87 {
88 // This is guaranteed to be called with an object that is using
89 // PeriodicBoundaryCondition
90 const auto * iterator = reinterpret_cast< const ConstNeighborhoodIterator< TInputImage, Self > * >( data );
91 typename TInputImage::InternalPixelType * ptr;
92 int linear_index = 0;
93 unsigned int i;
94
95 // Find the pointer of the closest boundary pixel
96 // std::cout << "Boundary offset = " << boundary_offset << std::endl;
97 // std::cout << "point index = " << point_index << std::endl;
98
99 // Return the value of the pixel at the closest boundary point.
100 for ( i = 0; i < ImageDimension; ++i )
101 {
102 linear_index += ( point_index[i] + boundary_offset[i] ) * data->GetStride(i);
103 }
104 ptr = data->operator[](linear_index);
105
106 // Wrap the pointer around the image in the necessary dimensions. If we have
107 // reached this point, we can assume that we are on the edge of the BUFFERED
108 // region of the image. Boundary conditions are only invoked if touching the
109 // actual memory boundary.
110
111 // These are the step sizes for increments in each dimension of the image.
112 const typename TInputImage::OffsetValueType * offset_table =
113 iterator->GetImagePointer()->GetOffsetTable();
114
115 for ( i = 0; i < ImageDimension; ++i )
116 {
117 if ( boundary_offset[i] != 0 )
118 { // If the neighborhood overlaps on the low edge, then wrap from the
119 // high edge of the image.
120 if ( point_index[i] < static_cast< OffsetValueType >( iterator->GetRadius(i) ) )
121 {
122 ptr += iterator->GetImagePointer()->GetBufferedRegion().GetSize()[i]
123 * offset_table[i] - boundary_offset[i] * offset_table[i];
124 }
125 else // wrap from the low side of the image
126 {
127 ptr -= iterator->GetImagePointer()->GetBufferedRegion().GetSize()[i]
128 * offset_table[i] + boundary_offset[i] * offset_table[i];
129 }
130 }
131 }
132
133 return static_cast< OutputPixelType >( neighborhoodAccessorFunctor.Get(ptr) );
134 }
135
136
137 template< typename TInputImage, typename TOutputImage >
138 typename PeriodicBoundaryCondition< TInputImage, TOutputImage >::RegionType
139 PeriodicBoundaryCondition< TInputImage, TOutputImage >
GetInputRequestedRegion(const RegionType & inputLargestPossibleRegion,const RegionType & outputRequestedRegion) const140 ::GetInputRequestedRegion( const RegionType & inputLargestPossibleRegion,
141 const RegionType & outputRequestedRegion ) const
142 {
143 IndexType imageIndex = inputLargestPossibleRegion.GetIndex();
144 SizeType imageSize = inputLargestPossibleRegion.GetSize();
145
146 IndexType outputIndex = outputRequestedRegion.GetIndex();
147 SizeType outputSize = outputRequestedRegion.GetSize();
148
149 IndexType inputRequestedIndex;
150 SizeType inputRequestedSize;
151
152 for ( unsigned int i = 0; i < ImageDimension; ++i )
153 {
154 // Check for image boundary overlap in the requested region
155 IndexValueType lowIndex =
156 ( ( outputIndex[i] - imageIndex[i] ) % static_cast< IndexValueType >( imageSize[i] ) );
157 if ( lowIndex < 0 )
158 {
159 lowIndex += static_cast< IndexValueType >( imageSize[i] );
160 }
161 IndexValueType highIndex = lowIndex + static_cast< IndexValueType >( outputSize[i] );
162
163 bool overlap = ( highIndex >= static_cast< IndexValueType >( imageSize[i] ) );
164
165 if ( overlap )
166 {
167 // Request the totality of the image in this dimension
168 inputRequestedIndex[i] = imageIndex[i];
169 inputRequestedSize[i] = imageSize[i];
170 }
171 else
172 {
173 // Remap the requested portion in this dimension into the image region.
174 inputRequestedIndex[i] = lowIndex;
175 inputRequestedSize[i] = outputSize[i];
176 }
177 }
178 RegionType inputRequestedRegion( inputRequestedIndex, inputRequestedSize );
179
180 return inputRequestedRegion;
181 }
182
183
184 template< typename TInputImage, typename TOutputImage >
185 typename PeriodicBoundaryCondition< TInputImage, TOutputImage >::OutputPixelType
186 PeriodicBoundaryCondition< TInputImage, TOutputImage >
GetPixel(const IndexType & index,const TInputImage * image) const187 ::GetPixel( const IndexType & index, const TInputImage * image ) const
188 {
189 RegionType imageRegion = image->GetLargestPossibleRegion();
190 IndexType imageIndex = imageRegion.GetIndex();
191 SizeType imageSize = imageRegion.GetSize();
192
193 IndexType lookupIndex;
194
195 for ( unsigned int i = 0; i < ImageDimension; ++i )
196 {
197 IndexValueType modIndex = ( ( index[i] - imageIndex[i] ) %
198 static_cast< IndexValueType >( imageSize[i] ) );
199
200 if ( modIndex < 0 )
201 {
202 modIndex += static_cast< IndexValueType >( imageSize[i] );
203 }
204
205 lookupIndex[i] = modIndex + imageIndex[i];
206 }
207
208 return static_cast< OutputPixelType >( image->GetPixel( lookupIndex ) );
209 }
210
211 } // end namespace itk
212
213 #endif
214