ITK  6.0.0
Insight Toolkit
itkBoxUtilities.h
Go to the documentation of this file.
1/*=========================================================================
2 *
3 * Copyright NumFOCUS
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 * https://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 itkBoxUtilities_h
19#define itkBoxUtilities_h
20
21#include "itkProgressReporter.h"
27#include "itkOffset.h"
30#include <algorithm> // For min.
31
32/*
33 *
34 * This code was contributed in the Insight Journal paper:
35 * "Efficient implementation of kernel filtering"
36 * by Beare R., Lehmann G
37 * https://doi.org/10.54294/igq8fn
38 *
39 */
40
42{
43
44template <typename TIterator>
45inline TIterator *
46setConnectivityEarlyBox(TIterator * it, bool fullyConnected = false)
47{
48 // activate the "previous" neighbours
49 typename TIterator::OffsetType offset;
50 it->ClearActiveList();
51 if (!fullyConnected)
52 {
53 // only activate the neighbors that are face connected
54 // to the current pixel. do not include the center pixel
55 offset.Fill(0);
56 for (unsigned int d = 0; d < TIterator::Dimension; ++d)
57 {
58 offset[d] = -1;
59 it->ActivateOffset(offset);
60 offset[d] = 0;
61 }
62 }
63 else
64 {
65 // activate all neighbors that are face+edge+vertex
66 // connected to the current pixel. do not include the center pixel
67 unsigned int centerIndex = it->GetCenterNeighborhoodIndex();
68 for (unsigned int d = 0; d < centerIndex; ++d)
69 {
70 offset = it->GetOffset(d);
71 // check for positives in any dimension
72 bool keep = true;
73 for (unsigned int i = 0; i < TIterator::Dimension; ++i)
74 {
75 if (offset[i] > 0)
76 {
77 keep = false;
78 break;
79 }
80 }
81 if (keep)
82 {
83 it->ActivateOffset(offset);
84 }
85 }
86 offset.Fill(0);
87 it->DeactivateOffset(offset);
88 }
89 return it;
90}
91
92} // namespace itk_impl_details
93
94namespace itk
95{
96
97template <typename TInputImage, typename TOutputImage>
98void
99BoxAccumulateFunction(const TInputImage * inputImage,
100 const TOutputImage * outputImage,
101 typename TInputImage::RegionType inputRegion,
102 typename TOutputImage::RegionType outputRegion)
103{
104 // type alias
105 using InputImageType = TInputImage;
106 using OffsetType = typename TInputImage::OffsetType;
107 using OutputImageType = TOutputImage;
108 using OutputPixelType = typename TOutputImage::PixelType;
109
110 using InputIterator = ImageRegionConstIterator<TInputImage>;
111
112 using NOutputIterator = ShapedNeighborhoodIterator<TOutputImage>;
113 InputIterator inIt(inputImage, inputRegion);
114 auto kernelRadius = TInputImage::SizeType::Filled(1);
115
116 NOutputIterator noutIt(kernelRadius, outputImage, outputRegion);
117 // this iterator is fully connected
119
121 oBC.SetConstant(OutputPixelType{});
122 noutIt.OverrideBoundaryCondition(&oBC);
123 // This uses several iterators. An alternative and probably better
124 // approach would be to copy the input to the output and convolve
125 // with the following weights (in 2D)
126 // -(dim - 1) 1
127 // 1 1
128 // The result of each convolution needs to get written back to the
129 // image being convolved so that the accumulation propagates
130 // This should be implementable with neighborhood operators.
131
132 std::vector<int> weights;
133 typename NOutputIterator::ConstIterator sIt;
134 for (auto idxIt = noutIt.GetActiveIndexList().begin(); idxIt != noutIt.GetActiveIndexList().end(); ++idxIt)
135 {
136 OffsetType offset = noutIt.GetOffset(*idxIt);
137 int w = -1;
138 for (unsigned int k = 0; k < InputImageType::ImageDimension; ++k)
139 {
140 if (offset[k] != 0)
141 {
142 w *= offset[k];
143 }
144 }
145 // std::cout << offset << " " << w << std::endl;
146 weights.push_back(w);
147 }
148
149 for (inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt)
150 {
151 OutputPixelType sum = 0;
152 int k;
153 for (k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd(); ++sIt, ++k)
154 {
155 sum += sIt.Get() * weights[k];
156 }
157 noutIt.SetCenterPixel(sum + inIt.Get());
158 }
159}
160
161// a function to generate corners of arbitrary dimension box
162template <typename TImage>
163std::vector<typename TImage::OffsetType>
164CornerOffsets(const TImage * im)
165{
166 using NIterator = ShapedNeighborhoodIterator<TImage>;
167 auto unitradius = TImage::SizeType::Filled(1);
168 NIterator n1(unitradius, im, im->GetRequestedRegion());
169 unsigned int centerIndex = n1.GetCenterNeighborhoodIndex();
170 typename NIterator::OffsetType offset;
171 std::vector<typename TImage::OffsetType> result;
172 for (unsigned int d = 0; d < centerIndex * 2 + 1; ++d)
173 {
174 offset = n1.GetOffset(d);
175 // check whether this is a corner - corners have no zeros
176 bool corner = true;
177 for (unsigned int k = 0; k < TImage::ImageDimension; ++k)
178 {
179 if (offset[k] == 0)
180 {
181 corner = false;
182 break;
183 }
184 }
185 if (corner)
186 {
187 result.push_back(offset);
188 }
189 }
190 return (result);
191}
192
193template <typename TInputImage, typename TOutputImage>
194void
195BoxMeanCalculatorFunction(const TInputImage * accImage,
196 TOutputImage * outputImage,
197 typename TInputImage::RegionType inputRegion,
198 typename TOutputImage::RegionType outputRegion,
199 typename TInputImage::SizeType radius)
200{
201 // type alias
202 using InputImageType = TInputImage;
203 using RegionType = typename TInputImage::RegionType;
204 using SizeType = typename TInputImage::SizeType;
205 using IndexType = typename TInputImage::IndexType;
206 using OffsetType = typename TInputImage::OffsetType;
207 using OutputImageType = TOutputImage;
208 using OutputPixelType = typename TOutputImage::PixelType;
209 // use the face generator for speed
211 using FaceListType = typename FaceCalculatorType::FaceListType;
212 FaceCalculatorType faceCalculator;
213
215
216 // this process is actually slightly asymmetric because we need to
217 // subtract rectangles that are next to our kernel, not overlapping it
218 SizeType kernelSize;
219 SizeType internalRadius;
220 SizeType regionLimit;
221
222 IndexType regionStart = inputRegion.GetIndex();
223 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
224 {
225 kernelSize[i] = radius[i] * 2 + 1;
226 internalRadius[i] = radius[i] + 1;
227 regionLimit[i] = inputRegion.GetSize()[i] + regionStart[i] - 1;
228 }
229
230 using AccPixType = typename NumericTraits<OutputPixelType>::RealType;
231 // get a set of offsets to corners for a unit hypercube in this image
232 std::vector<OffsetType> unitCorners = CornerOffsets<TInputImage>(accImage);
233 std::vector<OffsetType> realCorners;
234 std::vector<AccPixType> weights;
235 // now compute the weights
236 for (unsigned int k = 0; k < unitCorners.size(); ++k)
237 {
238 int prod = 1;
239 OffsetType thisCorner;
240 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
241 {
242 prod *= unitCorners[k][i];
243 if (unitCorners[k][i] > 0)
244 {
245 thisCorner[i] = radius[i];
246 }
247 else
248 {
249 thisCorner[i] = -(static_cast<OffsetValueType>(radius[i]) + 1);
250 }
251 }
252 weights.push_back((AccPixType)prod);
253 realCorners.push_back(thisCorner);
254 }
255
256 FaceListType faceList = faceCalculator(accImage, outputRegion, internalRadius);
257 // start with the body region
258 for (const auto & face : faceList)
259 {
260 if (&face == &faceList.front())
261 {
262 // this is the body region. This is meant to be an optimized
263 // version that doesn't use neighborhood regions
264 // compute the various offsets
265 AccPixType pixelscount = 1;
266 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
267 {
268 pixelscount *= (AccPixType)(2 * radius[i] + 1);
269 }
270
271 using OutputIteratorType = ImageRegionIterator<OutputImageType>;
272 using InputIteratorType = ImageRegionConstIterator<InputImageType>;
273
274 using CornerItVecType = std::vector<InputIteratorType>;
275 CornerItVecType cornerItVec;
276 // set up the iterators for each corner
277 for (unsigned int k = 0; k < realCorners.size(); ++k)
278 {
279 typename InputImageType::RegionType tReg = face;
280 tReg.SetIndex(tReg.GetIndex() + realCorners[k]);
281 InputIteratorType tempIt(accImage, tReg);
282 tempIt.GoToBegin();
283 cornerItVec.push_back(tempIt);
284 }
285 // set up the output iterator
286 OutputIteratorType oIt(outputImage, face);
287 // now do the work
288 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
289 {
290 AccPixType sum = 0;
291 // check each corner
292 for (unsigned int k = 0; k < cornerItVec.size(); ++k)
293 {
294 sum += weights[k] * cornerItVec[k].Get();
295 // increment each corner iterator
296 ++(cornerItVec[k]);
297 }
298 oIt.Set(static_cast<OutputPixelType>(sum / pixelscount));
299 }
300 }
301 else
302 {
303 // now we need to deal with the border regions
304 using OutputIteratorType = ImageRegionIteratorWithIndex<OutputImageType>;
305 OutputIteratorType oIt(outputImage, face);
306 // now do the work
307 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
308 {
309 // figure out the number of pixels in the box by creating an
310 // equivalent region and cropping - this could probably be
311 // included in the loop below.
312 RegionType currentKernelRegion;
313 currentKernelRegion.SetSize(kernelSize);
314 // compute the region's index
315 IndexType kernelRegionIdx = oIt.GetIndex();
316 IndexType centIndex = kernelRegionIdx;
317 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
318 {
319 kernelRegionIdx[i] -= radius[i];
320 }
321 currentKernelRegion.SetIndex(kernelRegionIdx);
322 currentKernelRegion.Crop(inputRegion);
323 OffsetValueType edgepixelscount = currentKernelRegion.GetNumberOfPixels();
324 AccPixType sum = 0;
325 // rules are : for each corner,
326 // for each dimension
327 // if dimension offset is positive -> this is
328 // a leading edge. Crop if outside the input
329 // region
330 // if dimension offset is negative -> this is
331 // a trailing edge. Ignore if it is outside
332 // image region
333 for (unsigned int k = 0; k < realCorners.size(); ++k)
334 {
335 IndexType thisCorner = centIndex + realCorners[k];
336 bool includeCorner = true;
337 for (unsigned int j = 0; j < TInputImage::ImageDimension; ++j)
338 {
339 if (unitCorners[k][j] > 0)
340 {
341 // leading edge - crop it
342 thisCorner[j] = std::min(thisCorner[j], static_cast<OffsetValueType>(regionLimit[j]));
343 }
344 else
345 {
346 // trailing edge - check bounds
347 if (thisCorner[j] < regionStart[j])
348 {
349 includeCorner = false;
350 break;
351 }
352 }
353 }
354 if (includeCorner)
355 {
356 sum += accImage->GetPixel(thisCorner) * weights[k];
357 }
358 }
359
360 oIt.Set(static_cast<OutputPixelType>(sum / (AccPixType)edgepixelscount));
361 }
362 }
363 }
364}
365
366template <typename TInputImage, typename TOutputImage>
367void
368BoxSigmaCalculatorFunction(const TInputImage * accImage,
369 TOutputImage * outputImage,
370 typename TInputImage::RegionType inputRegion,
371 typename TOutputImage::RegionType outputRegion,
372 typename TInputImage::SizeType radius)
373{
374 // type alias
375 using InputImageType = TInputImage;
376 using RegionType = typename TInputImage::RegionType;
377 using SizeType = typename TInputImage::SizeType;
378 using IndexType = typename TInputImage::IndexType;
379 using OffsetType = typename TInputImage::OffsetType;
380 using OutputImageType = TOutputImage;
381 using OutputPixelType = typename TOutputImage::PixelType;
382 using InputPixelType = typename TInputImage::PixelType;
383 // use the face generator for speed
385 using FaceListType = typename FaceCalculatorType::FaceListType;
386 FaceCalculatorType faceCalculator;
387
389
390 // this process is actually slightly asymmetric because we need to
391 // subtract rectangles that are next to our kernel, not overlapping it
392 SizeType kernelSize;
393 SizeType internalRadius;
394 SizeType regionLimit;
395 IndexType regionStart = inputRegion.GetIndex();
396 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
397 {
398 kernelSize[i] = radius[i] * 2 + 1;
399 internalRadius[i] = radius[i] + 1;
400 regionLimit[i] = inputRegion.GetSize()[i] + regionStart[i] - 1;
401 }
402
403 using AccPixType = typename NumericTraits<OutputPixelType>::RealType;
404 // get a set of offsets to corners for a unit hypercube in this image
405 std::vector<OffsetType> unitCorners = CornerOffsets<TInputImage>(accImage);
406 std::vector<OffsetType> realCorners;
407 std::vector<AccPixType> weights;
408 // now compute the weights
409 for (unsigned int k = 0; k < unitCorners.size(); ++k)
410 {
411 int prod = 1;
412 OffsetType thisCorner;
413 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
414 {
415 prod *= unitCorners[k][i];
416 if (unitCorners[k][i] > 0)
417 {
418 thisCorner[i] = radius[i];
419 }
420 else
421 {
422 thisCorner[i] = -(static_cast<OffsetValueType>(radius[i]) + 1);
423 }
424 }
425 weights.push_back((AccPixType)prod);
426 realCorners.push_back(thisCorner);
427 }
428
429 FaceListType faceList = faceCalculator(accImage, outputRegion, internalRadius);
430 // start with the body region
431 for (const auto & face : faceList)
432 {
433 if (&face == &faceList.front())
434 {
435 // this is the body region. This is meant to be an optimized
436 // version that doesn't use neighborhood regions
437 // compute the various offsets
438 AccPixType pixelscount = 1;
439 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
440 {
441 pixelscount *= (AccPixType)(2 * radius[i] + 1);
442 }
443
444 using OutputIteratorType = ImageRegionIterator<OutputImageType>;
445 using InputIteratorType = ImageRegionConstIterator<InputImageType>;
446
447 using CornerItVecType = std::vector<InputIteratorType>;
448 CornerItVecType cornerItVec;
449 // set up the iterators for each corner
450 for (unsigned int k = 0; k < realCorners.size(); ++k)
451 {
452 typename InputImageType::RegionType tReg = face;
453 tReg.SetIndex(tReg.GetIndex() + realCorners[k]);
454 InputIteratorType tempIt(accImage, tReg);
455 tempIt.GoToBegin();
456 cornerItVec.push_back(tempIt);
457 }
458 // set up the output iterator
459 OutputIteratorType oIt(outputImage, face);
460 // now do the work
461 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
462 {
463 AccPixType sum = 0;
464 AccPixType squareSum = 0;
465 // check each corner
466 for (unsigned int k = 0; k < cornerItVec.size(); ++k)
467 {
468 const InputPixelType & i = cornerItVec[k].Get();
469 sum += weights[k] * i[0];
470 squareSum += weights[k] * i[1];
471 // increment each corner iterator
472 ++(cornerItVec[k]);
473 }
474
475 oIt.Set(static_cast<OutputPixelType>(std::sqrt((squareSum - sum * sum / pixelscount) / (pixelscount - 1))));
476 }
477 }
478 else
479 {
480 // now we need to deal with the border regions
481 using OutputIteratorType = ImageRegionIteratorWithIndex<OutputImageType>;
482 OutputIteratorType oIt(outputImage, face);
483 // now do the work
484 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
485 {
486 // figure out the number of pixels in the box by creating an
487 // equivalent region and cropping - this could probably be
488 // included in the loop below.
489 RegionType currentKernelRegion;
490 currentKernelRegion.SetSize(kernelSize);
491 // compute the region's index
492 IndexType kernelRegionIdx = oIt.GetIndex();
493 IndexType centIndex = kernelRegionIdx;
494 for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
495 {
496 kernelRegionIdx[i] -= radius[i];
497 }
498 currentKernelRegion.SetIndex(kernelRegionIdx);
499 currentKernelRegion.Crop(inputRegion);
500 SizeValueType edgepixelscount = currentKernelRegion.GetNumberOfPixels();
501 AccPixType sum = 0;
502 AccPixType squareSum = 0;
503 // rules are : for each corner,
504 // for each dimension
505 // if dimension offset is positive -> this is
506 // a leading edge. Crop if outside the input
507 // region
508 // if dimension offset is negative -> this is
509 // a trailing edge. Ignore if it is outside
510 // image region
511 for (unsigned int k = 0; k < realCorners.size(); ++k)
512 {
513 IndexType thisCorner = centIndex + realCorners[k];
514 bool includeCorner = true;
515 for (unsigned int j = 0; j < TInputImage::ImageDimension; ++j)
516 {
517 if (unitCorners[k][j] > 0)
518 {
519 // leading edge - crop it
520 thisCorner[j] = std::min(thisCorner[j], static_cast<OffsetValueType>(regionLimit[j]));
521 }
522 else
523 {
524 // trailing edge - check bounds
525 if (thisCorner[j] < regionStart[j])
526 {
527 includeCorner = false;
528 break;
529 }
530 }
531 }
532 if (includeCorner)
533 {
534 const InputPixelType & i = accImage->GetPixel(thisCorner);
535 sum += weights[k] * i[0];
536 squareSum += weights[k] * i[1];
537 }
538 }
539
540 oIt.Set(
541 static_cast<OutputPixelType>(std::sqrt((squareSum - sum * sum / edgepixelscount) / (edgepixelscount - 1))));
542 }
543 }
544 }
545}
546
547template <typename TInputImage, typename TOutputImage>
548void
549BoxSquareAccumulateFunction(const TInputImage * inputImage,
550 TOutputImage * outputImage,
551 typename TInputImage::RegionType inputRegion,
552 typename TOutputImage::RegionType outputRegion)
553{
554 // type alias
555 using InputImageType = TInputImage;
556 using OffsetType = typename TInputImage::OffsetType;
557 using OutputImageType = TOutputImage;
558 using OutputPixelType = typename TOutputImage::PixelType;
559 using ValueType = typename OutputPixelType::ValueType;
560 using InputPixelType = typename TInputImage::PixelType;
561
562 using InputIterator = ImageRegionConstIterator<TInputImage>;
563
564 using NOutputIterator = ShapedNeighborhoodIterator<TOutputImage>;
565 InputIterator inIt(inputImage, inputRegion);
566 auto kernelRadius = TInputImage::SizeType::Filled(1);
567
568 NOutputIterator noutIt(kernelRadius, outputImage, outputRegion);
569 // this iterator is fully connected
571
573 oBC.SetConstant(OutputPixelType{});
574 noutIt.OverrideBoundaryCondition(&oBC);
575 // This uses several iterators. An alternative and probably better
576 // approach would be to copy the input to the output and convolve
577 // with the following weights (in 2D)
578 // -(dim - 1) 1
579 // 1 1
580 // The result of each convolution needs to get written back to the
581 // image being convolved so that the accumulation propagates
582 // This should be implementable with neighborhood operators.
583
584 std::vector<int> weights;
585 typename NOutputIterator::ConstIterator sIt;
586 for (auto idxIt = noutIt.GetActiveIndexList().begin(); idxIt != noutIt.GetActiveIndexList().end(); ++idxIt)
587 {
588 OffsetType offset = noutIt.GetOffset(*idxIt);
589 int w = -1;
590 for (unsigned int k = 0; k < InputImageType::ImageDimension; ++k)
591 {
592 if (offset[k] != 0)
593 {
594 w *= offset[k];
595 }
596 }
597 weights.push_back(w);
598 }
599
600 for (inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt)
601 {
602 ValueType sum = 0;
603 ValueType squareSum = 0;
604 int k;
605 for (k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd(); ++sIt, ++k)
606 {
607 const OutputPixelType & v = sIt.Get();
608 sum += v[0] * weights[k];
609 squareSum += v[1] * weights[k];
610 }
611 OutputPixelType o;
612 const InputPixelType & i = inIt.Get();
613 o[0] = sum + i;
614 o[1] = squareSum + i * i;
615 noutIt.SetCenterPixel(o);
616 }
617}
618} // namespace itk
619
620#endif
This boundary condition returns a constant value for out-of-bounds image pixels.
void SetConstant(const OutputPixelType &c)
A multi-dimensional iterator templated over image type that walks a region of pixels.
A multi-dimensional iterator templated over image type that walks pixels within a region and is speci...
A multi-dimensional iterator templated over image type that walks a region of pixels.
void SetSize(const SizeType &size)
void SetIndex(const IndexType &index)
Define additional traits for native types such as int or float.
A neighborhood iterator which can take on an arbitrary shape.
A function object that determines a neighborhood of values at an image boundary according to a Neuman...
TIterator * setConnectivityEarlyBox(TIterator *it, bool fullyConnected=false)
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
std::vector< typename TImage::OffsetType > CornerOffsets(const TImage *im)
void BoxSigmaCalculatorFunction(const TInputImage *accImage, TOutputImage *outputImage, typename TInputImage::RegionType inputRegion, typename TOutputImage::RegionType outputRegion, typename TInputImage::SizeType radius)
void BoxMeanCalculatorFunction(const TInputImage *accImage, TOutputImage *outputImage, typename TInputImage::RegionType inputRegion, typename TOutputImage::RegionType outputRegion, typename TInputImage::SizeType radius)
unsigned long SizeValueType
Definition: itkIntTypes.h:86
void BoxAccumulateFunction(const TInputImage *inputImage, const TOutputImage *outputImage, typename TInputImage::RegionType inputRegion, typename TOutputImage::RegionType outputRegion)
long OffsetValueType
Definition: itkIntTypes.h:97
void BoxSquareAccumulateFunction(const TInputImage *inputImage, TOutputImage *outputImage, typename TInputImage::RegionType inputRegion, typename TOutputImage::RegionType outputRegion)
const IndexValueType * GetIndex() const
Definition: itkIndex.h:230
void SetIndex(const IndexValueType val[VDimension])
Definition: itkIndex.h:240
Splits an image into a main region and several "face" regions which are used to handle computations o...