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
|
/*=========================================================================
*
* 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 itkApproximateSignedDistanceMapImageFilter_hxx
#define itkApproximateSignedDistanceMapImageFilter_hxx
#include "itkNumericTraits.h"
#include "itkImageScanlineIterator.h"
#include "itkProgressAccumulator.h"
namespace itk
{
template <typename TInputImage, typename TOutputImage>
ApproximateSignedDistanceMapImageFilter<TInputImage, TOutputImage>::ApproximateSignedDistanceMapImageFilter()
: m_IsoContourFilter(IsoContourType::New())
, m_ChamferFilter(ChamferType::New())
, m_InsideValue(NumericTraits<InputPixelType>::min())
, m_OutsideValue(NumericTraits<InputPixelType>::max())
{}
template <typename TInputImage, typename TOutputImage>
void
ApproximateSignedDistanceMapImageFilter<TInputImage, TOutputImage>::GenerateData()
{
ThreadIdType numberOfWorkUnits = this->GetNumberOfWorkUnits();
OutputImagePointer output = this->GetOutput();
using OutputRegionType = typename OutputImageType::RegionType;
OutputRegionType oRegion = output->GetRequestedRegion();
// Calculate the largest possible distance in the output image.
// this maximum is the distance from one corner of the image to the other.
OutputSizeType outputSize = oRegion.GetSize();
OutputSizeValueType maximumDistance = 0;
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
maximumDistance += outputSize[i] * outputSize[i];
}
// Cast to double and back because there's no sqrt defined for unsigned long,
// which is the general SizeValueType.
maximumDistance = static_cast<OutputSizeValueType>(std::sqrt(static_cast<double>(maximumDistance)));
// Allocate the output
this->AllocateOutputs();
// Create a process accumulator for tracking the progress of this minipipeline
auto progress = ProgressAccumulator::New();
progress->SetMiniPipelineFilter(this);
progress->RegisterInternalFilter(m_IsoContourFilter, 0.5f);
progress->RegisterInternalFilter(m_ChamferFilter, 0.5f);
// Set up the isocontour filter
m_IsoContourFilter->SetInput(this->GetInput());
m_IsoContourFilter->SetFarValue(maximumDistance + 1);
m_IsoContourFilter->SetNumberOfWorkUnits(numberOfWorkUnits);
typename IsoContourType::PixelRealType levelSetValue =
(static_cast<typename IsoContourType::PixelRealType>(m_InsideValue) +
static_cast<typename IsoContourType::PixelRealType>(m_OutsideValue)) /
2.0;
m_IsoContourFilter->SetLevelSetValue(levelSetValue);
// Set up the chamfer filter
m_ChamferFilter->SetInput(m_IsoContourFilter->GetOutput());
m_ChamferFilter->SetMaximumDistance(maximumDistance);
m_ChamferFilter->SetNumberOfWorkUnits(numberOfWorkUnits);
// Graft our output to the chamfer filter to force the proper regions
// to be generated
m_ChamferFilter->GraftOutput(output);
// Create the distance map
m_ChamferFilter->Update();
// Graft the output of the chamfer filter back onto this filter's
// output. this is needed to get the appropriate regions passed
// back.
this->GraftOutput(m_ChamferFilter->GetOutput());
// Recall that we set the isocontour value to halfway between the inside and
// outside value. The above filters assume that regions "inside" objects are
// those with values *less* than the isocontour. This assumption is violated
// if the "inside" intensity value is greater than the "outside" value.
// (E.g. in the case that we're computing the distance from a mask where the
// background is zero and the objects are colored 255.)
// In this case, the distance will be calculated negative, so we need to
// flip the sign of the output image.
if (m_InsideValue > m_OutsideValue)
{
for (ImageScanlineIterator ot(output, oRegion); !ot.IsAtEnd(); ot.NextLine())
{
while (!ot.IsAtEndOfLine())
{
ot.Set(ot.Get() * -1);
++ot;
}
}
}
}
template <typename TInputImage, typename TOutputImage>
void
ApproximateSignedDistanceMapImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "Inside intensity value: " << m_InsideValue << std::endl;
os << indent << "Outside intensity value: " << m_OutsideValue << std::endl;
os << indent << "IsoContourDistanceImageFilter (used internally): " << m_IsoContourFilter << std::endl;
os << indent << "FastChamferDistanceImageFilter (used internally): " << m_ChamferFilter << std::endl;
}
} // end of namespace itk
#endif
|
