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
|
#include "thresholdconfig.h"
#include <algorithm>
#include <cmath>
#include <iostream>
#include "../structures/image2d.h"
#include "../util/rng.h"
#include "thresholdtools.h"
#include "sumthreshold.h"
#include "sumthresholdmissing.h"
namespace algorithms {
ThresholdConfig::ThresholdConfig() : _distribution(Gaussian) {}
void ThresholdConfig::InitializeLengthsDefault(unsigned count) {
if (count > 9 || count == 0) count = 9;
constexpr size_t lengths[9] = {1, 2, 4, 8, 16, 32, 64, 128, 256};
_horizontalOperations.clear();
_horizontalOperations.reserve(count);
_verticalOperations.clear();
_verticalOperations.reserve(count);
for (unsigned i = 0; i < count; ++i) {
_horizontalOperations.emplace_back(lengths[i]);
_verticalOperations.emplace_back(lengths[i]);
}
}
void ThresholdConfig::InitializeLengthsSingleSample() {
_horizontalOperations.emplace_back(1);
_verticalOperations.emplace_back(1);
}
void ThresholdConfig::InitializeThresholdsFromFirstThreshold(
num_t firstThreshold, Distribution noiseDistribution) {
constexpr num_t exp_base = 1.5;
for (ThresholdConfig::ThresholdOperation& op : _horizontalOperations) {
op.threshold =
firstThreshold * std::pow(exp_base, std::log2(op.length)) / op.length;
}
for (ThresholdConfig::ThresholdOperation& op : _verticalOperations) {
op.threshold =
firstThreshold * std::pow(exp_base, std::log2(op.length)) / op.length;
}
_distribution = noiseDistribution;
}
void ThresholdConfig::Execute(const Image2D* image, Mask2D* mask, bool additive,
num_t timeSensitivity,
num_t frequencySensitivity) const {
ExecuteWithMissing(image, mask, nullptr, additive, timeSensitivity,
frequencySensitivity);
}
void ThresholdConfig::ExecuteWithMissing(const Image2D* image, Mask2D* mask,
const Mask2D* missing, bool additive,
num_t timeSensitivity,
num_t frequencySensitivity) const {
num_t timeFactor = timeSensitivity;
num_t frequencyFactor = frequencySensitivity;
Image2D modified_image;
if (!image->AllFinite()) {
modified_image = image->MakeFiniteCopy();
image = &modified_image;
}
switch (_distribution) {
case Gaussian: {
num_t mean, stddev;
if (missing == nullptr)
ThresholdTools::WinsorizedMeanAndStdDev(image, mask, mean, stddev);
else
ThresholdTools::WinsorizedMeanAndStdDev(image, mask, missing, mean,
stddev);
if (stddev != 0.0L) {
timeFactor = stddev * timeSensitivity;
frequencyFactor = stddev * frequencySensitivity;
}
} break;
case Rayleigh: {
num_t mode;
if (missing == nullptr)
mode = ThresholdTools::WinsorizedMode(image, mask);
else
mode = ThresholdTools::WinsorizedMode(image, mask, missing);
if (mode != 0.0L) {
timeFactor = timeSensitivity * mode;
frequencyFactor = frequencySensitivity * mode;
}
} break;
}
if (!additive) mask->SetAll<false>();
Mask2D scratch(*mask);
const size_t operationCount =
std::max(_horizontalOperations.size(), _verticalOperations.size());
SumThreshold::VerticalScratch normalScratch(mask->Width(), mask->Height());
SumThresholdMissing::VerticalCache vMissingCache;
if (missing)
SumThresholdMissing::InitializeVertical(vMissingCache, *image, *missing);
for (size_t i = 0; i < operationCount; ++i) {
if (i < _horizontalOperations.size()) {
if (missing == nullptr) {
SumThreshold::HorizontalLarge(
image, mask, &scratch, _horizontalOperations[i].length,
_horizontalOperations[i].threshold * timeFactor);
} else {
SumThresholdMissing::Horizontal(
*image, *mask, *missing, scratch, _horizontalOperations[i].length,
_horizontalOperations[i].threshold * timeFactor);
}
}
if (i < _verticalOperations.size()) {
if (missing == nullptr)
SumThreshold::VerticalLarge(
image, mask, &scratch, &normalScratch,
_verticalOperations[i].length,
_verticalOperations[i].threshold * frequencyFactor);
else
SumThresholdMissing::VerticalStacked(
vMissingCache, *image, *mask, *missing, scratch,
_verticalOperations[i].length,
_verticalOperations[i].threshold * frequencyFactor);
}
}
}
} // namespace algorithms
|
