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
|
#include "siroperator.h"
namespace algorithms {
template <typename MaskLikeA, typename MaskLikeB>
void SIROperator::operateHorizontallyMissing(MaskLikeA& mask,
const MaskLikeB& missing,
num_t eta) {
const unsigned width = mask.Width(), maxWSize = width + 1;
std::unique_ptr<num_t[]> values(new num_t[width]), w(new num_t[maxWSize]);
std::unique_ptr<unsigned[]> minIndices(new unsigned[maxWSize]),
maxIndices(new unsigned[maxWSize]);
for (unsigned row = 0; row < mask.Height(); ++row) {
unsigned nAvailable = 0;
for (unsigned i = 0; i < width; ++i) {
if (!missing.Value(i, row)) {
if (mask.Value(i, row))
values[nAvailable] = eta;
else
values[nAvailable] = eta - 1.0;
++nAvailable;
}
}
if (nAvailable != 0) {
const unsigned wSize = nAvailable + 1;
w[0] = 0.0;
unsigned currentMinIndex = 0;
minIndices[0] = 0;
for (unsigned i = 1; i != wSize; ++i) {
w[i] = w[i - 1] + values[i - 1];
if (w[i] < w[currentMinIndex]) {
currentMinIndex = i;
}
minIndices[i] = currentMinIndex;
}
// Calculate the maximum suffixes
unsigned currentMaxIndex = wSize - 1;
for (unsigned i = nAvailable - 1; i != 0; --i) {
maxIndices[i] = currentMaxIndex;
if (w[i] > w[currentMaxIndex]) {
currentMaxIndex = i;
}
}
maxIndices[0] = currentMaxIndex;
// See if max sequence exceeds limit.
nAvailable = 0;
for (unsigned i = 0; i != width; ++i) {
if (!missing.Value(i, row)) {
const num_t maxW =
w[maxIndices[nAvailable]] - w[minIndices[nAvailable]];
mask.SetValue(i, row, (maxW >= 0.0));
++nAvailable;
}
}
}
}
}
template void SIROperator::operateHorizontallyMissing(Mask2D& mask,
const Mask2D& missing,
num_t eta);
template void SIROperator::operateHorizontallyMissing(
XYSwappedMask2D<Mask2D>& mask, const XYSwappedMask2D<const Mask2D>& missing,
num_t eta);
template <typename MaskLikeA, typename MaskLikeB>
void SIROperator::operateHorizontallyMissing(MaskLikeA& mask,
const MaskLikeB& missing,
num_t eta, num_t penalty) {
const size_t width = mask.Width();
const size_t maxWSize = width + 1;
const std::unique_ptr<num_t[]> values(new num_t[width]);
const std::unique_ptr<num_t[]> w(new num_t[maxWSize]);
const std::unique_ptr<size_t[]> minIndices(new size_t[maxWSize]);
const std::unique_ptr<size_t[]> maxIndices(new size_t[maxWSize]);
const num_t penaltyValue = (eta - 1.0) * penalty;
for (size_t row = 0; row < mask.Height(); ++row) {
for (size_t i = 0; i != width; ++i) {
if (missing.Value(i, row))
values[i] = penaltyValue;
else if (mask.Value(i, row))
values[i] = eta;
else
values[i] = eta - 1.0;
}
w[0] = 0.0;
size_t currentMinIndex = 0;
minIndices[0] = 0;
for (size_t i = 1; i != width + 1; ++i) {
w[i] = w[i - 1] + values[i - 1];
if (w[i] < w[currentMinIndex]) {
currentMinIndex = i;
}
minIndices[i] = currentMinIndex;
}
// Calculate the maximum suffixes
size_t currentMaxIndex = width;
for (size_t i = width - 1; i != 0; --i) {
maxIndices[i] = currentMaxIndex;
if (w[i] > w[currentMaxIndex]) {
currentMaxIndex = i;
}
}
maxIndices[0] = currentMaxIndex;
// See if max sequence exceeds limit.
for (size_t i = 0; i != width; ++i) {
if (!missing.Value(i, row)) {
const num_t maxW = w[maxIndices[i]] - w[minIndices[i]];
mask.SetValue(i, row, maxW >= 0.0);
}
}
}
}
template void SIROperator::operateHorizontallyMissing(Mask2D& mask,
const Mask2D& missing,
num_t eta, num_t penalty);
template void SIROperator::operateHorizontallyMissing(
XYSwappedMask2D<Mask2D>& mask, const XYSwappedMask2D<const Mask2D>& missing,
num_t eta, num_t penalty);
} // namespace algorithms
|
