1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkNeighborhoodAlgorithm_hxx
#define itkNeighborhoodAlgorithm_hxx
#include "itkImageRegionIterator.h"
#include "itkImageRegion.h"
#include "itkConstSliceIterator.h"
#include <algorithm> // For min.
namespace itk
{
namespace NeighborhoodAlgorithm
{
template <typename TImage>
auto
ImageBoundaryFacesCalculator<TImage>::Compute(const TImage & img, RegionType regionToProcess, RadiusType radius)
-> Result
{
// Analyze the regionToProcess to determine if any of its faces are
// along a buffer boundary (we have no data in the buffer for pixels
// that are outside the boundary, but within the neighborhood radius and will
// have to treat them differently). We also determine the size of the non-
// boundary region that will be processed. For instance, given a 2D image
// and regionTOProcess (size = 5x5),
Result result;
const RegionType & bufferedRegion = img.GetBufferedRegion();
// The portion of the regionToProcess that is outside of the image bufferedRegion
// doesn't make sense. If the regionToProcess is completely outside of
// the image bufferedRegion, return an empty (default-constructed) result.
if (!regionToProcess.Crop(bufferedRegion))
{
return result;
}
const IndexType bStart = bufferedRegion.GetIndex();
const SizeType bSize = bufferedRegion.GetSize();
const IndexType rStart = regionToProcess.GetIndex();
const SizeType rSize = regionToProcess.GetSize();
SizeType nbSize = rSize; // Non-boundary region
IndexType nbStart = rStart; // data.
IndexType vrStart =
rStart; // start index of variable processed region which has considered the boundary region in last direction
SizeType vrSize =
rSize; // size of variable processed region which has considered the boundary region in last direction
for (unsigned int i = 0; i < ImageDimension; ++i)
{
IndexType fStart; // Boundary, "face"
SizeType fSize; // region data.
auto overlapLow = static_cast<IndexValueType>((rStart[i] - radius[i]) - bStart[i]);
// image buffered region should be more than twice of the radius,
// otherwise there would be overlap between two boundary regions.
// in the case, we reduce upper boundary size to remove overlap.
IndexValueType overlapHigh;
if (bSize[i] > 2 * radius[i])
{
overlapHigh = static_cast<IndexValueType>((bStart[i] + bSize[i]) - (rStart[i] + rSize[i] + radius[i]));
}
else
{
overlapHigh = static_cast<IndexValueType>((bStart[i] + radius[i]) - (rStart[i] + rSize[i]));
}
if (overlapLow < 0) // out of bounds condition, define
// a region of
{ // iteration along this face
for (unsigned int j = 0; j < ImageDimension; ++j) // define the starting index
{ // and size of the face region
fStart[j] = vrStart[j];
if (j == i)
{
// Boundary region cannot be outside the region to process
if (-overlapLow > static_cast<IndexValueType>(rSize[i]))
{
overlapLow = -static_cast<IndexValueType>(rSize[i]);
}
fSize[j] = -overlapLow;
vrSize[j] += overlapLow; // change start and size in this direction
vrStart[j] -= overlapLow; // to ensure no duplicate pixels at corners
}
else
{
fSize[j] = vrSize[j];
}
// Boundary region cannot be outside the region to process
fSize[j] = std::min(fSize[j], rSize[j]);
}
// avoid unsigned overflow if the non-boundary region is too small to
// process
if (fSize[i] > nbSize[i])
{
nbSize[i] = 0;
}
else
{
nbSize[i] -= fSize[i];
}
nbStart[i] += -overlapLow;
result.m_BoundaryFaces.emplace_back(fStart, fSize);
}
if (overlapHigh < 0)
{
for (unsigned int j = 0; j < ImageDimension; ++j)
{
if (j == i)
{
// Boundary region cannot be outside the region to process
if (-overlapHigh > static_cast<IndexValueType>(rSize[i]))
{
overlapHigh = -static_cast<IndexValueType>(rSize[i]);
}
fStart[j] = rStart[j] + static_cast<IndexValueType>(rSize[j]) + overlapHigh;
fSize[j] = -overlapHigh;
vrSize[j] += overlapHigh; // change size in this direction
}
else
{
fStart[j] = vrStart[j];
fSize[j] = vrSize[j];
}
}
// avoid unsigned overflow if the non-boundary region is too small to
// process
if (fSize[i] > nbSize[i])
{
nbSize[i] = 0;
}
else
{
nbSize[i] -= fSize[i];
}
result.m_BoundaryFaces.emplace_back(fStart, fSize);
}
}
result.m_NonBoundaryRegion.SetSize(nbSize);
result.m_NonBoundaryRegion.SetIndex(nbStart);
return result;
}
template <typename TImage>
auto
ImageBoundaryFacesCalculator<TImage>::operator()(const TImage * img, RegionType regionToProcess, RadiusType radius)
-> FaceListType
{
const auto result = Compute(*img, regionToProcess, radius);
if (result == Result{})
{
return FaceListType{};
}
else
{
FaceListType faceList = std::move(result.m_BoundaryFaces);
faceList.push_front(result.m_NonBoundaryRegion);
return faceList;
}
}
template <typename TImage>
auto
CalculateOutputWrapOffsetModifiers<TImage>::operator()(TImage * input, TImage * output) const -> OffsetType
{
OffsetType ans;
const Size<TImage::ImageDimension> isz = input->GetBufferedRegion().GetSize();
const Size<TImage::ImageDimension> osz = output->GetBufferedRegion().GetSize();
for (int i = 0; i < TImage::ImageDimension; ++i)
{
ans[i] = osz[i] - isz[i];
}
return ans;
}
} // end namespace NeighborhoodAlgorithm
} // end namespace itk
#endif
|
