-
Notifications
You must be signed in to change notification settings - Fork 3
/
hist-2d-dist.h
285 lines (245 loc) · 8.69 KB
/
hist-2d-dist.h
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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
// hist-2d-dist.h -- Sampling distribution based on a 2d histogram
//
// Copyright (C) 2010-2012 Miles Bader <miles@gnu.org>
//
// This source code is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 3, or (at
// your option) any later version. See the file COPYING for more details.
//
// Written by Miles Bader <miles@gnu.org>
//
#ifndef SNOGRAY_HIST_2D_DIST_H
#define SNOGRAY_HIST_2D_DIST_H
#include <vector>
#include <algorithm>
#include "compiler.h"
#include "hist-2d.h"
namespace snogray {
// A sampling distribution based on a 2d histogram. This is useful for
// doing cheap re-sampling based on an arbitrary 2d input set.
//
class Hist2dDist
{
public:
// This constructor allocates the necessary memory, but won't be
// usable until a histogram has been specified using Hist2dDist::calc.
//
Hist2dDist (unsigned w, unsigned h)
: width (w), height (h), size (w*h),
column_width (1.f / width), row_height (1.f / height),
whole_row_cumulative_sums (height),
individual_row_cumulative_sums (size)
{ }
// This constructor automatically copies the size from HIST, and
// calculates the PDF. No references to HIST is kept.
//
Hist2dDist (const Hist2d &hist)
: width (hist.width), height (hist.height), size (height * width),
column_width (1.f / width), row_height (1.f / height),
whole_row_cumulative_sums (height),
individual_row_cumulative_sums (size)
{ calc (hist); }
// Calculate the PDF based from the histogram HIST. HIST's size must
// be the same as the size this object was created with. No reference
// to HIST is kept.
//
void calc (const Hist2d &hist)
{
// Note, the use of double-precision floats here is intentional --
// HDR images can cause precision problems if single-precision
// floats are used.
unsigned row_offs;
// Find sum of entire input array.
//
double bin_sum = 0;
row_offs = 0;
for (unsigned row = 0; row < height; row++)
{
for (unsigned col = 0; col < width; col++)
bin_sum += double (hist.bins[row_offs + col]);
row_offs += width;
}
// Find cumulative sums of entire rows, normalized to the range 0-1
// (so the last row will always have a value of 1, except in the
// degenerate case where all bins are zero).
//
double inv_bin_sum = (bin_sum == 0) ? 0 : 1 / bin_sum;
double normalized_sum = 0;
row_offs = 0;
for (unsigned row = 0; row < height; row++)
{
for (unsigned col = 0; col < width; col++)
normalized_sum += double (hist.bins[row_offs + col]) * inv_bin_sum;
whole_row_cumulative_sums[row] = normalized_sum;
row_offs += width;
}
// Find cumulative sums within each row, normalized to the range 0-1
// (so for each row, the last column within the row will always have
// value 1, except in the degenerate case where all bins are zero).
//
row_offs = 0;
for (unsigned row = 0; row < height; row++)
{
double row_sum = 0;
for (unsigned col = 0; col < width; col++)
row_sum += double (hist.bins[row_offs + col]);
double inv_row_sum = (row_sum == 0) ? 0 : 1 / row_sum;
double normalized_row_sum = 0;
for (unsigned col = 0; col < width; col++)
{
normalized_row_sum
+= double (hist.bins[row_offs + col]) * inv_row_sum;
individual_row_cumulative_sums[row_offs + col] = normalized_row_sum;
}
row_offs += width;
}
}
// Return a sample of this distribution based on the random
// variables in PARAM. The PDF at the sample location is returned
// in _PDF.
//
// The returned UV coordinates should have roughly the same
// distribution as the input data (limited by the granularity of
// the histogram).
//
UV sample (const UV ¶m, float &_pdf) const
{
unsigned col, row, row_offs;
if (sample (param, col, row, row_offs))
{
_pdf = pdf (col, row, row_offs);
return UV (col * column_width + fmod (param.u, column_width),
row * row_height + fmod (param.v, row_height));
}
else
{
_pdf = 0;
return UV (0, 0);
}
}
// Return a sample of this distribution based on the random
// variables in PARAM.
//
// The returned UV coordinates should have roughly the same
// distribution as the input data (limited by the granularity of the
// histogram).
//
UV sample (const UV ¶m) const
{
unsigned col, row, row_offs;
if (sample (param, col, row, row_offs))
return UV (col * column_width + fmod (param.u, column_width),
row * row_height + fmod (param.v, row_height));
else
return UV (0,0);
}
// Return the PDF of this distribution at location POS.
//
float pdf (const UV &pos) const
{
unsigned col = clamp (int (pos.u * width), 0, int (width) - 1);
unsigned row = clamp (int (pos.v * height), 0, int (height) - 1);
return pdf (col, row, row * width);
}
const unsigned width, height, size;
const float column_width, row_height;
private:
template<typename Ptr>
unsigned find_pos_in_sorted_vec (float val, const Ptr &beg, unsigned len)
const
{
#if 0
unsigned pos = 0;
while (val > beg[pos])
pos++;
return pos;
#else
// XXX std::lower_bound is much faster than a linear search for
// large histograms, but slower for small ones; it'd be nice to
// parameterize this somehow...
//
return std::lower_bound (beg, beg + len, val) - beg;
#endif
}
// Sample the histogram and return the coordinates of the chosen
// bin in COL and ROW. The offset of the beginning of the row in
// INDIVIDUAL_ROW_CUMULATIVE_SUMS is also returned in ROW_OFFSET.
//
// Normally the function return value is true, but in the rare
// case where sampling is impossible because _all_ the data was
// zero, false is returned instead (and all other return values
// are undefined).
//
bool sample (const UV ¶m, unsigned &col, unsigned &row,
unsigned &row_offs)
const
{
float u = min (param.u, 1.f), v = min (param.v, 1.f);
// look in y dir.
//
row = find_pos_in_sorted_vec (v, whole_row_cumulative_sums.begin(), height);
// If sampling totally failed, return false (this should only happen
// if all the data in the source histogram was zero).
//
if (unlikely (row == height))
return false;
// XXX this multiply actually uses a lot of time; it's nicer to
// accumulate the row-offset while finding the right how, or maybe
// keep a vector of row offsets?
//
row_offs = row * width;
// look in x dir
//
col = find_pos_in_sorted_vec (
u, individual_row_cumulative_sums.begin() + row_offs, width);
return true;
}
// Return the PDF of this distribution for locations in the bin
// located at (COL, ROW), where ROW_OFFS is the offset in
// INDIVIDUAL_ROW_CUMULATIVE_SUMS of the beginning of the row.
//
float pdf (unsigned col, unsigned row, unsigned row_offs) const
{
// Probability of choosing this row.
//
// As WHOLE_ROW_CUMULATIVE_SUMS contains cumulative sums of
// whole-row probabilities, the probability of this row is the
// difference of this row's cumulative-sum entry minus the
// previous row's entry.
//
float row_prob = whole_row_cumulative_sums[row];
if (row != 0)
row_prob -= whole_row_cumulative_sums[row - 1];
// Probability of choosing this column in the row. Similarly to
// ROW_PROB, this as the difference of the entries for the current
// and previous columns in INDIVIDUAL_ROW_CUMULATIVE_SUMS.
//
float col_prob = individual_row_cumulative_sums[row_offs + col];
if (col != 0)
col_prob -= individual_row_cumulative_sums[row_offs + col - 1];
// Probability of choosing this bin, which is just the probability
// of choosing this row (ROW_PROB) multiplied by the probability
// of choosing this columin within the row (COL_PROB).
//
float bin_prob = row_prob * col_prob;
// PDF = probability of choosing a bin / bin area. Since we
// consider the "total area" to be 1, then the bin area is just
// 1 / the number of bins (which is SIZE).
//
return bin_prob * size;
}
// Cumulative sum of whole-row probabilities. Each entry is the
// probability of choosing that row or any row before it (so the
// last entry is always 1).
//
std::vector<float> whole_row_cumulative_sums;
// For each row, the cumulative sum of column probabilities for
// that row. Each entry is the probability of choosing that colum
// in the row (assuming the row is chosen) or any column before it
// (so the last entry for each row is always 1).
//
std::vector<float> individual_row_cumulative_sums;
};
}
#endif // SNOGRAY_HIST_2D_DIST_H