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
|
/*=========================================================================
*
* 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 itkBresenhamLine_hxx
#define itkBresenhamLine_hxx
#include "itkPoint.h"
#include "itkMath.h"
namespace itk
{
template <unsigned int VDimension>
auto
BresenhamLine<VDimension>::BuildLine(LType Direction, IdentifierType length) -> OffsetArray
{
// copied from the line iterator
/** Variables that drive the Bresenham-Algorithm */
// The dimension with the largest difference between start and end
unsigned int m_MainDirection;
// Accumulated error for the other dimensions
IndexType m_AccumulateError;
// Increment for the error for each step. Two times the difference between
// start and end
IndexType m_IncrementError;
// If enough is accumulated for a dimension, the index has to be
// incremented. Will be the number of pixels in the line
IndexType m_MaximalError;
// Direction of increment. -1 or 1
IndexType m_OverflowIncrement;
// After an overflow, the accumulated error is reduced again. Will be
// two times the number of pixels in the line
IndexType m_ReduceErrorAfterIncrement;
OffsetArray result(length);
IndexType m_CurrentImageIndex, LastIndex;
Direction.Normalize();
// we are going to start at 0
m_CurrentImageIndex.Fill(0);
constexpr IndexType StartIndex = { { 0 } };
for (unsigned int i = 0; i < VDimension; ++i)
{
LastIndex[i] = (IndexValueType)(length * Direction[i]);
}
// Find the dominant direction
IndexValueType maxDistance = 0;
unsigned int maxDistanceDimension = 0;
for (unsigned int i = 0; i < VDimension; ++i)
{
auto distance = static_cast<long>(itk::Math::abs(LastIndex[i]));
if (distance > maxDistance)
{
maxDistance = distance;
maxDistanceDimension = i;
}
m_IncrementError[i] = 2 * distance;
m_OverflowIncrement[i] = (LastIndex[i] < 0 ? -1 : 1);
}
m_MainDirection = maxDistanceDimension;
m_MaximalError.Fill(maxDistance);
m_ReduceErrorAfterIncrement.Fill(2 * maxDistance);
m_AccumulateError.Fill(0);
unsigned int steps = 1;
result[0] = m_CurrentImageIndex - StartIndex;
while (steps < length)
{
// This part is from ++ in LineConstIterator
// We need to modify m_AccumulateError, m_CurrentImageIndex, m_IsAtEnd
for (unsigned int i = 0; i < VDimension; ++i)
{
if (i == m_MainDirection)
{
m_CurrentImageIndex[i] += m_OverflowIncrement[i];
}
else
{
m_AccumulateError[i] += m_IncrementError[i];
if (m_AccumulateError[i] >= m_MaximalError[i])
{
m_CurrentImageIndex[i] += m_OverflowIncrement[i];
m_AccumulateError[i] -= m_ReduceErrorAfterIncrement[i];
}
}
}
result[steps] = m_CurrentImageIndex - StartIndex; // produce an offset
++steps;
}
return (result);
}
template <unsigned int VDimension>
auto
BresenhamLine<VDimension>::BuildLine(IndexType p0, IndexType p1) -> IndexArray
{
itk::Point<float, VDimension> point0;
itk::Point<float, VDimension> point1;
IdentifierType maxDistance = 0;
for (unsigned int i = 0; i < VDimension; ++i)
{
point0[i] = p0[i];
point1[i] = p1[i];
IdentifierType distance = itk::Math::abs(p0[i] - p1[i]) + 1;
if (distance > maxDistance)
{
maxDistance = distance;
}
}
OffsetArray offsets = this->BuildLine(point1 - point0, maxDistance + 1);
IndexArray indices;
indices.reserve(offsets.size()); // we might not have to use the last one
for (unsigned int i = 0; i < offsets.size(); ++i)
{
indices.push_back(p0 + offsets[i]);
if (indices.back() == p1)
{
break;
}
}
return indices;
}
} // namespace itk
#endif
|
