uniform.hh 30 KB
Newer Older
Philip Trettner's avatar
Philip Trettner committed
1
2
#pragma once

Philip Trettner's avatar
Philip Trettner committed
3
#include <cstddef>
Philip Trettner's avatar
Philip Trettner committed
4
5
#include <initializer_list>

Philip Trettner's avatar
Philip Trettner committed
6
#include <typed-geometry/detail/scalar_traits.hh>
Philip Trettner's avatar
Philip Trettner committed
7

Philip Trettner's avatar
Philip Trettner committed
8
#include <typed-geometry/types/color.hh>
Philip Trettner's avatar
Philip Trettner committed
9
#include <typed-geometry/types/pos.hh>
10
#include <typed-geometry/types/range.hh>
Philip Trettner's avatar
Philip Trettner committed
11
#include <typed-geometry/types/scalars/scalars.hh>
12

Philip Trettner's avatar
Philip Trettner committed
13
#include <typed-geometry/types/objects/aabb.hh>
Philip Trettner's avatar
Philip Trettner committed
14
#include <typed-geometry/types/objects/box.hh>
15
#include <typed-geometry/types/objects/capsule.hh>
Philip Trettner's avatar
Philip Trettner committed
16
#include <typed-geometry/types/objects/cone.hh>
17
#include <typed-geometry/types/objects/cylinder.hh>
Aaron Grabowy's avatar
Aaron Grabowy committed
18
#include <typed-geometry/types/objects/ellipse.hh>
19
#include <typed-geometry/types/objects/hemisphere.hh>
Aaron Grabowy's avatar
Aaron Grabowy committed
20
#include <typed-geometry/types/objects/pyramid.hh>
Aaron Grabowy's avatar
Aaron Grabowy committed
21
#include <typed-geometry/types/objects/segment.hh>
Philip Trettner's avatar
Philip Trettner committed
22
23
#include <typed-geometry/types/objects/sphere.hh>
#include <typed-geometry/types/objects/triangle.hh>
24

Philip Trettner's avatar
Philip Trettner committed
25
26
27
#include <typed-geometry/functions/basic/minmax.hh>
#include <typed-geometry/functions/basic/mix.hh>
#include <typed-geometry/functions/vector/math.hh>
Philip Trettner's avatar
Philip Trettner committed
28

Philip Trettner's avatar
Philip Trettner committed
29
#include "random.hh"
30

Philip Trettner's avatar
Philip Trettner committed
31
/*
32
33
 * uniform<T>(rng)    - uniform sample for T's where no parameter is needed
 * uniform(rng)       - returns an object that implicitly converts to T (and calls uniform<T>(rng))
34
 * uniform(rng, a, b) - uniform between a..b (inclusive!)
Philip Trettner's avatar
Philip Trettner committed
35
 * uniform(rng, obj)  - uniform sample from obj
36
 *
37
38
39
40
 * Note:
 *  to "uniformly sample" vec3(0)..vec3(1) independent per component, use uniform(rng, tg::aabb3(0, 1))
 *  to "uniformly sample" vec3(0)..vec3(1) as a segment, use uniform(rng, tg::segment3({0,0,0}, {1,1,1}))
 *
41
 * uniform_vec(rng, ...) same as uniform but returns uniform(rng, ...) as a vector
Philip Trettner's avatar
Philip Trettner committed
42
43
 */

44
45
// TODO: uniform int/uint distribution might need some improvement but is currently faster than the stdlib versions

Philip Trettner's avatar
Philip Trettner committed
46
47
namespace tg
{
Philip Trettner's avatar
Philip Trettner committed
48
49
50
51
52
53
54
namespace detail
{
template <class T>
struct sampler
{
    static_assert(sizeof(T) >= 0, "sampling of this type not supported without parameters");
};
55
56
57
58
59
60
61
62
63
64
65
66

// NOTE: this class should not be saved
template <class Rng>
struct uniform_sampler
{
    explicit uniform_sampler(Rng& rng) : rng(rng) {}

    template <class T>
    operator T() &&
    {
        return sampler<T>::uniform(rng);
    }
67
68
69
70
71
    template <class T>
    operator T() &
    {
        static_assert(always_false<T>, "the result of uniform(rng) must not be stored!");
    }
72
73
74
75

private:
    Rng& rng;
};
Philip Trettner's avatar
Philip Trettner committed
76
77
78
}

template <class T, class Rng>
79
[[nodiscard]] constexpr T uniform(Rng& rng)
Philip Trettner's avatar
Philip Trettner committed
80
{
Philip Trettner's avatar
Philip Trettner committed
81
    return detail::sampler<T>::uniform(rng);
Philip Trettner's avatar
Philip Trettner committed
82
83
}

84
85
86
87
88
89
template <class Rng>
[[nodiscard]] constexpr detail::uniform_sampler<Rng> uniform(Rng& rng)
{
    return detail::uniform_sampler<Rng>{rng};
}

90
91
// ======== Uniform between two numbers ========

Philip Trettner's avatar
Philip Trettner committed
92
template <class Rng>
93
[[nodiscard]] constexpr f32 uniform(Rng& rng, f32 a, f32 b)
Philip Trettner's avatar
Philip Trettner committed
94
95
96
97
{
    return mix(a, b, detail::uniform01<f32>(rng));
}
template <class Rng>
98
[[nodiscard]] constexpr f64 uniform(Rng& rng, f64 a, f64 b)
Philip Trettner's avatar
Philip Trettner committed
99
100
101
102
{
    return mix(a, b, detail::uniform01<f64>(rng));
}
template <class Rng>
103
[[nodiscard]] constexpr i32 uniform(Rng& rng, i32 a, i32 b_inc)
Philip Trettner's avatar
Philip Trettner committed
104
{
Philip Trettner's avatar
Philip Trettner committed
105
    TG_CONTRACT(a <= b_inc);
Philip Trettner's avatar
Philip Trettner committed
106
    i32 r = 0;
107
    auto fa = f32(a);
108
    auto fb = f32(b_inc) + 1;
109
110
111
    do
    {
        r = tg::ifloor(uniform(rng, fa, fb));
112
    } while (r > b_inc);
113
    return r;
Philip Trettner's avatar
Philip Trettner committed
114
115
}
template <class Rng>
116
[[nodiscard]] constexpr long uniform(Rng& rng, long a, long b_inc)
Philip Trettner's avatar
Philip Trettner committed
117
{
Philip Trettner's avatar
Philip Trettner committed
118
    TG_CONTRACT(a <= b_inc);
119
    long r = 0;
120
    auto fa = f64(a);
121
    auto fb = f64(b_inc) + 1;
122
123
    do
    {
124
125
126
127
128
129
130
        r = (long)tg::ifloor(uniform(rng, fa, fb));
    } while (r > b_inc);
    return r;
}
template <class Rng>
[[nodiscard]] constexpr i64 uniform(Rng& rng, long long a, long long b_inc)
{
Philip Trettner's avatar
Philip Trettner committed
131
    TG_CONTRACT(a <= b_inc);
132
133
134
135
136
137
    long long r = 0;
    auto fa = f64(a);
    auto fb = f64(b_inc) + 1;
    do
    {
        r = (long long)tg::ifloor(uniform(rng, fa, fb));
138
    } while (r > b_inc);
139
    return r;
Philip Trettner's avatar
Philip Trettner committed
140
141
}
template <class Rng>
142
[[nodiscard]] constexpr u32 uniform(Rng& rng, u32 a, u32 b_inc)
Philip Trettner's avatar
Philip Trettner committed
143
{
Philip Trettner's avatar
Philip Trettner committed
144
    u32 r = 0;
145
    auto fa = f32(a);
146
    auto fb = f32(b_inc) + 1;
147
148
149
    do
    {
        r = u32(tg::ifloor(uniform(rng, fa, fb)));
150
    } while (r > b_inc);
151
    return r;
Philip Trettner's avatar
Philip Trettner committed
152
153
}
template <class Rng>
154
[[nodiscard]] constexpr unsigned long uniform(Rng& rng, unsigned long a, unsigned long b_inc)
Philip Trettner's avatar
Philip Trettner committed
155
{
Philip Trettner's avatar
Philip Trettner committed
156
    TG_CONTRACT(a <= b_inc);
157
    unsigned long r = 0;
158
    auto fa = f64(a);
159
    auto fb = f64(b_inc) + 1;
160
161
    do
    {
162
        r = (unsigned long)(tg::ifloor(uniform(rng, fa, fb)));
163
    } while (r > b_inc);
164
    return r;
Philip Trettner's avatar
Philip Trettner committed
165
}
166
167
template <class Rng>
[[nodiscard]] constexpr unsigned long long uniform(Rng& rng, unsigned long long a, unsigned long long b_inc)
Philip Trettner's avatar
Philip Trettner committed
168
{
Philip Trettner's avatar
Philip Trettner committed
169
    TG_CONTRACT(a <= b_inc);
170
    unsigned long long r = 0;
171
    auto fa = f64(a);
172
    auto fb = f64(b_inc) + 1;
173
174
    do
    {
175
        r = (unsigned long long)(tg::ifloor(uniform(rng, fa, fb)));
176
    } while (r > b_inc);
177
    return r;
Philip Trettner's avatar
Philip Trettner committed
178
}
Philip Trettner's avatar
Philip Trettner committed
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
template <class Rng>
[[nodiscard]] constexpr char uniform(Rng& rng, char a, char b_inc)
{
    return char(uniform(rng, i32(a), i32(b_inc)));
}
template <class Rng>
[[nodiscard]] constexpr i8 uniform(Rng& rng, i8 a, i8 b_inc)
{
    return i8(uniform(rng, i32(a), i32(b_inc)));
}
template <class Rng>
[[nodiscard]] constexpr i16 uniform(Rng& rng, i16 a, i16 b_inc)
{
    return i16(uniform(rng, i32(a), i32(b_inc)));
}
template <class Rng>
[[nodiscard]] constexpr u8 uniform(Rng& rng, u8 a, u8 b_inc)
{
    return u8(uniform(rng, i32(a), i32(b_inc)));
}
template <class Rng>
[[nodiscard]] constexpr u16 uniform(Rng& rng, u16 a, u16 b_inc)
{
    return u16(uniform(rng, i32(a), i32(b_inc)));
}
Philip Trettner's avatar
Philip Trettner committed
204

205
template <class T, class Rng>
206
[[nodiscard]] constexpr angle_t<T> uniform(Rng& rng, angle_t<T> a, angle_t<T> b)
207
208
209
210
{
    return mix(a, b, detail::uniform01<T>(rng));
}

211
212
// ======== Uniform selection from a list ========

213
template <class Rng, class T>
214
[[nodiscard]] constexpr T uniform(Rng& rng, std::initializer_list<T> list)
215
216
217
218
219
{
    TG_CONTRACT(list.size() > 0 && "cannot pick from an empty list");
    return list.begin()[uniform(rng, u64(0), u64(list.size() - 1))];
}

220
template <class Rng, class Container>
221
[[nodiscard]] constexpr auto uniform(Rng& rng, Container& c) -> decltype(c[c.size()])
222
223
{
    TG_CONTRACT(c.size() > 0 && "cannot pick from an empty container");
224
    return c[uniform(rng, u64(0), u64(c.size() - 1))];
225
226
}

227
228
// ======== Uniform point on multiple objects ========

229
230
template <class Rng, class ObjectA, class... ObjectRest>
constexpr auto uniform_by_length(Rng& rng, ObjectA const& obj, ObjectRest const&... rest) -> typename object_traits<ObjectA>::pos_t
231
232
233
234
235
{
    using ScalarT = typename object_traits<ObjectA>::scalar_t;

    static_assert(object_traits<ObjectA>::object_dimension == 1, "objects must be 1D (but might be embedded in another dimension)");
    static_assert((... && (object_traits<ObjectRest>::object_dimension == 1)), "objects must be 1D (but might be embedded in another dimension)");
Philip Trettner's avatar
Philip Trettner committed
236
237
    static_assert((... && (object_traits<ObjectA>::domain_dimension == object_traits<ObjectRest>::domain_dimension)), "objects must be in the same "
                                                                                                                      "domain");
238
239
240
241
    static_assert((... && std::is_same_v<typename object_traits<ObjectRest>::scalar_t, ScalarT>), "objects must have same scalar type");
    static_assert(!std::is_integral_v<ScalarT>, "sampling from integer objects not supported (yet)");

    ScalarT lengthSums[1 + sizeof...(rest)];
Aaron Grabowy's avatar
Aaron Grabowy committed
242
    auto lengthSum = ScalarT(0);
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
    { // compute cumulative areas
        auto i = 0;
        lengthSums[i++] = lengthSum += length(obj);
        ((lengthSums[i++] = lengthSum += length(rest)), ...);
    }

    auto part = uniform(rng, ScalarT(0), lengthSum);

    // special case for the first object
    if (part <= lengthSums[0])
        return uniform(rng, obj);

    typename object_traits<ObjectA>::pos_t result;
    { // otherwise, sample from rest
        auto i = 1;
        // this code finds the interval where part is within lengthSums, samples from the associated object (stored in result), then stops by short circuiting
        auto ok = ((part <= lengthSums[i++] ? (result = uniform(rng, rest), true) : false) || ...);
        TG_ASSERT(ok);
    }
    return result;
}

265
266
template <class Rng, class ObjectA, class... ObjectRest>
constexpr auto uniform_by_area(Rng& rng, ObjectA const& obj, ObjectRest const&... rest) -> typename object_traits<ObjectA>::pos_t
267
268
269
270
271
{
    using ScalarT = typename object_traits<ObjectA>::scalar_t;

    static_assert(object_traits<ObjectA>::object_dimension == 2, "objects must be 2D (but might be embedded in another dimension)");
    static_assert((... && (object_traits<ObjectRest>::object_dimension == 2)), "objects must be 2D (but might be embedded in another dimension)");
Philip Trettner's avatar
Philip Trettner committed
272
273
    static_assert((... && (object_traits<ObjectA>::domain_dimension == object_traits<ObjectRest>::domain_dimension)), "objects must be in the same "
                                                                                                                      "domain");
274
275
276
277
    static_assert((... && std::is_same_v<typename object_traits<ObjectRest>::scalar_t, ScalarT>), "objects must have same scalar type");
    static_assert(!std::is_integral_v<ScalarT>, "sampling from integer objects not supported (yet)");

    ScalarT areaSums[1 + sizeof...(rest)];
Aaron Grabowy's avatar
Aaron Grabowy committed
278
    auto areaSum = ScalarT(0);
279
280
    { // compute cumulative areas
        auto i = 0;
Aaron Grabowy's avatar
Aaron Grabowy committed
281
282
        areaSums[i++] = areaSum += area_of(obj);
        ((areaSums[i++] = areaSum += area_of(rest)), ...);
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
    }

    auto part = uniform(rng, ScalarT(0), areaSum);

    // special case for the first object
    if (part <= areaSums[0])
        return uniform(rng, obj);

    typename object_traits<ObjectA>::pos_t result;
    { // otherwise, sample from rest
        auto i = 1;
        // this code finds the interval where part is within areaSums, samples from the associated object (stored in result), then stops by short circuiting
        auto ok = ((part <= areaSums[i++] ? (result = uniform(rng, rest), true) : false) || ...);
        TG_ASSERT(ok);
    }
    return result;
}

301
302
template <class Rng, class ObjectA, class... ObjectRest>
constexpr auto uniform_by_volume(Rng& rng, ObjectA const& obj, ObjectRest const&... rest) -> typename object_traits<ObjectA>::pos_t
303
304
305
306
307
{
    using ScalarT = typename object_traits<ObjectA>::scalar_t;

    static_assert(object_traits<ObjectA>::object_dimension == 3, "objects must be 3D (but might be embedded in another dimension)");
    static_assert((... && (object_traits<ObjectRest>::object_dimension == 3)), "objects must be 3D (but might be embedded in another dimension)");
Philip Trettner's avatar
Philip Trettner committed
308
309
    static_assert((... && (object_traits<ObjectA>::domain_dimension == object_traits<ObjectRest>::domain_dimension)), "objects must be in the same "
                                                                                                                      "domain");
310
311
312
313
    static_assert((... && std::is_same_v<typename object_traits<ObjectRest>::scalar_t, ScalarT>), "objects must have same scalar type");
    static_assert(!std::is_integral_v<ScalarT>, "sampling from integer objects not supported (yet)");

    ScalarT volumeSums[1 + sizeof...(rest)];
Aaron Grabowy's avatar
Aaron Grabowy committed
314
    auto volumeSum = ScalarT(0);
315
316
    { // compute cumulative volumes
        auto i = 0;
Aaron Grabowy's avatar
Aaron Grabowy committed
317
318
        volumeSums[i++] = volumeSum += volume_of(obj);
        ((volumeSums[i++] = volumeSum += volume_of(rest)), ...);
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
    }

    auto part = uniform(rng, ScalarT(0), volumeSum);

    // special case for the first object
    if (part <= volumeSums[0])
        return uniform(rng, obj);

    typename object_traits<ObjectA>::pos_t result;
    { // otherwise, sample from rest
        auto i = 1;
        // this code finds the interval where part is within volumeSums, samples from the associated object (stored in result), then stops by short circuiting
        auto ok = ((part <= volumeSums[i++] ? (result = uniform(rng, rest), true) : false) || ...);
        TG_ASSERT(ok);
    }
    return result;
}

// ======== Uniform point on an object ========

339
template <class Rng>
340
[[nodiscard]] constexpr int uniform(Rng& rng, range1 const& b)
341
342
343
344
345
{
    return uniform(rng, b.min, b.max - 1);
}

template <class Rng>
346
[[nodiscard]] constexpr comp<2, int> uniform(Rng& rng, range2 const& b)
347
348
349
350
351
352
{
    return {uniform(rng, b.min.comp0, b.max.comp0 - 1), //
            uniform(rng, b.min.comp1, b.max.comp1 - 1)};
}

template <class Rng>
353
[[nodiscard]] constexpr comp<3, int> uniform(Rng& rng, range3 const& b)
354
355
356
357
358
359
{
    return {uniform(rng, b.min.comp0, b.max.comp0 - 1), //
            uniform(rng, b.min.comp1, b.max.comp1 - 1), //
            uniform(rng, b.min.comp2, b.max.comp2 - 1)};
}

Philip Trettner's avatar
Philip Trettner committed
360
template <class ScalarT, class Rng>
361
[[nodiscard]] constexpr pos<1, ScalarT> uniform(Rng& rng, aabb<1, ScalarT> const& b)
Philip Trettner's avatar
Philip Trettner committed
362
{
363
    return {uniform(rng, b.min.x, b.max.x)};
Philip Trettner's avatar
Philip Trettner committed
364
365
}
template <class ScalarT, class Rng>
366
[[nodiscard]] constexpr pos<2, ScalarT> uniform(Rng& rng, aabb<2, ScalarT> const& b)
Philip Trettner's avatar
Philip Trettner committed
367
368
369
370
371
{
    return {uniform(rng, b.min.x, b.max.x), //
            uniform(rng, b.min.y, b.max.y)};
}
template <class ScalarT, class Rng>
372
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, aabb<3, ScalarT> const& b)
Philip Trettner's avatar
Philip Trettner committed
373
374
375
376
377
378
{
    return {uniform(rng, b.min.x, b.max.x), //
            uniform(rng, b.min.y, b.max.y), //
            uniform(rng, b.min.z, b.max.z)};
}
template <class ScalarT, class Rng>
379
[[nodiscard]] constexpr pos<4, ScalarT> uniform(Rng& rng, aabb<4, ScalarT> const& b)
Philip Trettner's avatar
Philip Trettner committed
380
381
382
383
384
385
{
    return {uniform(rng, b.min.x, b.max.x), //
            uniform(rng, b.min.y, b.max.y), //
            uniform(rng, b.min.z, b.max.z), //
            uniform(rng, b.min.w, b.max.w)};
}
386
template <class ScalarT, class Rng>
387
[[nodiscard]] constexpr pos<1, ScalarT> uniform(Rng& rng, aabb_boundary<1, ScalarT> const& b)
388
{
389
    return uniform<bool>(rng) ? b.min.x : b.max.x;
390
391
}
template <class ScalarT, class Rng>
392
[[nodiscard]] constexpr pos<2, ScalarT> uniform(Rng& rng, aabb_boundary<2, ScalarT> const& b)
393
394
395
396
{
    auto extends = b.max - b.min;
    if (uniform(rng, ScalarT(0), extends.x + extends.y) < extends.x)
        return {uniform(rng, b.min.x, b.max.x), //
397
                uniform<bool>(rng) ? b.min.y : b.max.y};
398

399
    return {uniform<bool>(rng) ? b.min.x : b.max.x, //
400
401
402
            uniform(rng, b.min.y, b.max.y)};
}
template <class ScalarT, class Rng>
403
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, aabb_boundary<3, ScalarT> const& b)
404
405
406
407
408
409
{
    auto extends = b.max - b.min;
    auto areaX = extends.y * extends.z;
    auto areaY = extends.x * extends.z;
    auto areaZ = extends.x * extends.y;

Aaron Grabowy's avatar
Aaron Grabowy committed
410
    auto res = uniform(rng, solid_of(b)); // Sample a random point inside the aabb
411
412
413
414
415
416
417
    // Project to one of the sides, proportional to their area
    auto part = uniform(rng, ScalarT(0), areaX + areaY + areaZ);
    int i = part < areaX ? 0 : part < areaX + areaY ? 1 : 2;
    res[i] = uniform<bool>(rng) ? b.min[i] : b.max[i];
    return res;
}
template <class ScalarT, class Rng>
418
[[nodiscard]] constexpr pos<4, ScalarT> uniform(Rng& rng, aabb_boundary<4, ScalarT> const& b)
419
420
421
422
423
424
425
{
    auto extends = b.max - b.min;
    auto volX = extends.y * extends.z * extends.w;
    auto volY = extends.x * extends.z * extends.w;
    auto volZ = extends.x * extends.y * extends.w;
    auto volW = extends.x * extends.y * extends.z;

Aaron Grabowy's avatar
Aaron Grabowy committed
426
    auto res = uniform(rng, solid_of(b)); // Sample a random point inside the aabb
427
428
429
430
431
432
    // Project to one of the borders, proportional to their volume
    auto part = uniform(rng, ScalarT(0), volX + volY + volZ + volW);
    int i = part < volX + volY ? (part < volX ? 0 : 1) : (part < volX + volY + volZ ? 2 : 3);
    res[i] = uniform<bool>(rng) ? b.min[i] : b.max[i];
    return res;
}
Philip Trettner's avatar
Philip Trettner committed
433

Aaron Grabowy's avatar
Aaron Grabowy committed
434
template <int D, class ScalarT, class Rng>
435
[[nodiscard]] constexpr pos<D, ScalarT> uniform(Rng& rng, segment<D, ScalarT> const& s)
Aaron Grabowy's avatar
Aaron Grabowy committed
436
437
438
439
{
    return mix(s.pos0, s.pos1, detail::uniform01<ScalarT>(rng));
}

Aaron Grabowy's avatar
Aaron Grabowy committed
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
template <int D, class ScalarT, class Rng, class = enable_if<is_floating_point<ScalarT>>>
[[nodiscard]] constexpr pos<D, ScalarT> uniform(Rng& rng, triangle<D, ScalarT> const& t)
{
    auto e0 = t.pos1 - t.pos0;
    auto e1 = t.pos2 - t.pos0;
    auto u0 = uniform(rng, ScalarT(0), ScalarT(1));
    auto u1 = uniform(rng, ScalarT(0), ScalarT(1));
    if (u0 + u1 > ScalarT(1))
    {
        u0 = ScalarT(1) - u0;
        u1 = ScalarT(1) - u1;
    }
    return t.pos0 + u0 * e0 + u1 * e1;
}

455
456
template <int ObjectD, class ScalarT, int DomainD, class TraitsT, class Rng>
[[nodiscard]] constexpr pos<DomainD, ScalarT> uniform(Rng& rng, box<ObjectD, ScalarT, DomainD, TraitsT> const& b)
457
{
458
    return b.center + b.half_extents * uniform_vec(rng, aabb<ObjectD, ScalarT, TraitsT>::minus_one_to_one);
459
}
Philip Trettner's avatar
Philip Trettner committed
460

Aaron Grabowy's avatar
Aaron Grabowy committed
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
template <int ObjectD, class ScalarT, int DomainD, class Rng>
[[nodiscard]] constexpr pos<DomainD, ScalarT> uniform(Rng& rng, ellipse<ObjectD, ScalarT, DomainD> const& e)
{
    // Note: The transformation preserves the uniform distribution in the interior. This is not the case for the boundary_tag variant.
    return e.center + e.semi_axes * uniform_vec(rng, sphere<ObjectD, ScalarT>::unit);
}

template <class ScalarT, class Rng>
[[nodiscard]] constexpr pos<1, ScalarT> uniform(Rng& rng, ellipse_boundary<1, ScalarT> const& e)
{
    return e.center + (uniform<bool>(rng) ? e.semi_axes[0][0] : -e.semi_axes[0][0]);
}
template <int ObjectD, class ScalarT, int DomainD, class Rng>
[[nodiscard]] constexpr pos<DomainD, ScalarT> uniform(Rng& rng, ellipse_boundary<ObjectD, ScalarT, DomainD> const& e)
{
    // based on https://math.stackexchange.com/a/982833
    auto ls = vec<ObjectD, ScalarT>();
    auto coeffs = comp<ObjectD, ScalarT>();
    for (auto i = 0; i < ObjectD; ++i)
        ls[i] = length(e.semi_axes[i]);

    if constexpr (ObjectD == 2)
        coeffs = {ls.x, ls.y};
    else if constexpr (ObjectD == 3)
        coeffs = {ls.y * ls.z, ls.x * ls.z, ls.x * ls.y};
    else if constexpr (ObjectD == 4)
        coeffs = {ls.y * ls.z * ls.w, ls.x * ls.z * ls.w, ls.x * ls.y * ls.w, ls.x * ls.y * ls.z};
    else
        static_assert(always_false<ObjectD>, "dimension not supported");

    auto probMax = max_element(coeffs);

    // sample transformed sphere_boundary with rejection of some samples for uniform distribution
    while (true)
    {
        auto v = uniform_vec(rng, sphere_boundary<ObjectD, ScalarT>::unit);
        auto prob = length(coeffs * v); // == sqrt(pow2(coeffs.x * v.x) + pow2(coeffs.y * v.y) + ...)
        if (detail::uniform01<ScalarT>(rng) <= prob / probMax)
            return e.center + e.semi_axes * v;
    }
}

503
504
505
template <int D, class ScalarT, class Rng>
[[nodiscard]] constexpr pos<D, ScalarT> uniform(Rng& rng, sphere_boundary<D, ScalarT> const& s)
{
Aaron Grabowy's avatar
Aaron Grabowy committed
506
    auto ub = aabb<D, ScalarT>::minus_one_to_one;
507
508
509
510
511
512
513
514
515
516
517
518
    while (true)
    {
        auto p = uniform_vec(rng, ub);
        auto l = length_sqr(p);
        if (l > ScalarT(0) && l <= ScalarT(1))
            return s.center + p * (s.radius / sqrt(l));
    }
}

template <int D, class ScalarT, class Rng>
[[nodiscard]] constexpr pos<D, ScalarT> uniform(Rng& rng, sphere<D, ScalarT> const& b)
{
Aaron Grabowy's avatar
Aaron Grabowy committed
519
    auto ub = aabb<D, ScalarT>::minus_one_to_one;
520
521
522
523
524
525
526
527
528
    while (true)
    {
        auto p = uniform_vec(rng, ub);
        auto l = length_sqr(p);
        if (l <= ScalarT(1))
            return b.center + p * b.radius;
    }
}

529
template <class ScalarT, class Rng>
530
[[nodiscard]] constexpr pos<2, ScalarT> uniform(Rng& rng, sphere_boundary<2, ScalarT> const& c)
531
{
532
    return c.center + c.radius * uniform<dir<2, ScalarT>>(rng);
533
534
}
template <class ScalarT, class Rng>
535
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, sphere_boundary<2, ScalarT, 3> const& c)
536
{
537
    auto direction = uniform<dir<2, ScalarT>>(rng);
538
539
    auto x = any_normal(c.normal);
    auto y = cross(c.normal, x);
540
    return c.center + c.radius * (direction.x * x + direction.y * y);
541
542
543
}

template <class ScalarT, class Rng>
544
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, sphere<2, ScalarT, 3> const& d)
545
{
546
    auto direction = uniform(rng, sphere<2, ScalarT>(pos<2, ScalarT>::zero, d.radius));
547
548
    auto x = any_normal(d.normal);
    auto y = cross(d.normal, x);
549
    return d.center + direction.x * x + direction.y * y;
550
551
}

552
553
554
555
556
557
558
559
560
561
562
563
564
template <class ScalarT, class Rng>
[[nodiscard]] constexpr pos<2, ScalarT> uniform(Rng& rng, sphere<1, ScalarT, 2> const& s)
{
    auto v = perpendicular(s.normal) * s.radius;
    return mix(s.center - v, s.center + v, detail::uniform01<ScalarT>(rng));
}
template <class ScalarT, class Rng>
[[nodiscard]] constexpr pos<2, ScalarT> uniform(Rng& rng, sphere_boundary<1, ScalarT, 2> const& s)
{
    auto v = perpendicular(s.normal) * s.radius;
    return uniform(rng) ? s.center + v : s.center - v;
}

Aaron Grabowy's avatar
Aaron Grabowy committed
565
template <class ScalarT, class Rng>
566
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, cylinder<3, ScalarT> const& t)
567
{
568
    auto d = sphere<2, ScalarT, 3>(pos<3, ScalarT>::zero, t.radius, normalize(t.axis.pos1 - t.axis.pos0));
569
570
571
    return uniform(rng, t.axis) + vec<3, ScalarT>(uniform(rng, d));
}
template <class ScalarT, class Rng>
572
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, cylinder_boundary_no_caps<3, ScalarT> const& t)
Aaron Grabowy's avatar
Aaron Grabowy committed
573
{
574
    auto c = sphere_boundary<2, ScalarT, 3>(pos<3, ScalarT>::zero, t.radius, normalize(t.axis.pos1 - t.axis.pos0));
Aaron Grabowy's avatar
Aaron Grabowy committed
575
    return uniform(rng, t.axis) + vec<3, ScalarT>(uniform(rng, c));
Aaron Grabowy's avatar
Aaron Grabowy committed
576
577
}

Aaron Grabowy's avatar
Aaron Grabowy committed
578
template <class ScalarT, class Rng>
579
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, cylinder_boundary<3, ScalarT> const& c)
Aaron Grabowy's avatar
Aaron Grabowy committed
580
581
582
{
    auto x = c.axis.pos1 - c.axis.pos0;
    auto h = length(x);
Aaron Grabowy's avatar
Aaron Grabowy committed
583
    auto sideArea = ScalarT(2) * c.radius * h; // * Pi, but that does not matter here
584
    auto capArea = c.radius * c.radius;        // * Pi
Aaron Grabowy's avatar
Aaron Grabowy committed
585
    auto totalArea = ScalarT(2) * capArea + sideArea;
Aaron Grabowy's avatar
Aaron Grabowy committed
586
587
    auto part = detail::uniform01<ScalarT>(rng) * totalArea;
    if (part < sideArea) // Uniform sampling on cylinder side
588
        return uniform(rng, boundary_no_caps_of(c));
Aaron Grabowy's avatar
Aaron Grabowy committed
589
590

    // Otherwise sampling on one of the caps
591
    auto capDisk = sphere<2, ScalarT, 3>(part < sideArea + capArea ? c.axis.pos0 : c.axis.pos1, c.radius, normalize(x));
Aaron Grabowy's avatar
Aaron Grabowy committed
592
593
    return uniform(rng, capDisk);
}
Aaron Grabowy's avatar
Aaron Grabowy committed
594

Aaron Grabowy's avatar
Aaron Grabowy committed
595
template <class ScalarT, class Rng>
596
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, capsule<3, ScalarT> const& c)
597
598
599
{
    auto x = c.axis.pos1 - c.axis.pos0;
    auto h = length(x);
600
    auto tubeVolume = c.radius * c.radius * h;                                 // * Pi, but that does not matter here
601
602
603
604
    auto capVolume = ScalarT(2) / ScalarT(3) * c.radius * c.radius * c.radius; // * Pi
    auto totalVolume = ScalarT(2) * capVolume + tubeVolume;
    auto part = detail::uniform01<ScalarT>(rng) * totalVolume;
    if (part < tubeVolume) // Uniform sampling in capsule tube part
605
        return uniform(rng, cylinder<3, ScalarT>(c.axis, c.radius));
606
607
608
609
610
611
612
613

    // Otherwise sampling in one of the caps
    auto capHemi = hemisphere<3, ScalarT>();
    capHemi.radius = c.radius;
    capHemi.center = part < tubeVolume + capVolume ? c.axis.pos0 : c.axis.pos1;
    capHemi.normal = part < tubeVolume + capVolume ? -normalize(x) : normalize(x);
    return uniform(rng, capHemi);
}
614

615
template <class ScalarT, class Rng>
616
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, capsule_boundary<3, ScalarT> const& c)
Aaron Grabowy's avatar
Aaron Grabowy committed
617
618
619
{
    auto x = c.axis.pos1 - c.axis.pos0;
    auto h = length(x);
620
    auto sideArea = ScalarT(2) * c.radius * h; // * Pi, but that does not matter here
621
    auto capArea = c.radius * c.radius;        // * Pi
Aaron Grabowy's avatar
Aaron Grabowy committed
622
    auto totalArea = ScalarT(2) * capArea + sideArea;
Aaron Grabowy's avatar
Aaron Grabowy committed
623
624
    auto part = detail::uniform01<ScalarT>(rng) * totalArea;
    if (part < sideArea) // Uniform sampling on capsule side
625
        return uniform(rng, cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
Aaron Grabowy's avatar
Aaron Grabowy committed
626
627

    // Otherwise sampling on one of the caps
628
    auto capHemi = hemisphere_boundary_no_caps<3, ScalarT>();
Aaron Grabowy's avatar
Aaron Grabowy committed
629
630
631
    capHemi.radius = c.radius;
    capHemi.center = part < sideArea + capArea ? c.axis.pos0 : c.axis.pos1;
    capHemi.normal = part < sideArea + capArea ? -normalize(x) : normalize(x);
632
    return uniform(rng, capHemi);
Aaron Grabowy's avatar
Aaron Grabowy committed
633
634
}

635
template <int D, class ScalarT, class Rng>
636
[[nodiscard]] constexpr pos<D, ScalarT> uniform(Rng& rng, hemisphere<D, ScalarT> const& h)
637
{
638
    auto p = uniform(rng, sphere<D, ScalarT>(h.center, h.radius));
639
640
641
642
643
    auto v = p - h.center;
    if (dot(v, h.normal) >= ScalarT(0))
        return p;
    else
        return h.center - v;
644
645
}
template <int D, class ScalarT, class Rng>
646
[[nodiscard]] constexpr pos<D, ScalarT> uniform(Rng& rng, hemisphere_boundary_no_caps<D, ScalarT> const& h)
647
{
648
    auto p = uniform(rng, sphere_boundary<D, ScalarT>(h.center, h.radius));
Aaron Grabowy's avatar
Aaron Grabowy committed
649
    auto v = p - h.center;
Aaron Grabowy's avatar
Aaron Grabowy committed
650
    if (dot(v, h.normal) >= ScalarT(0))
Aaron Grabowy's avatar
Aaron Grabowy committed
651
652
653
654
        return p;
    else
        return h.center - v;
}
655
656
template <int D, class ScalarT, class Rng>
[[nodiscard]] constexpr pos<D, ScalarT> uniform(Rng& rng, hemisphere_boundary<D, ScalarT> const& h)
657
{
658
659
    ScalarT ratio;
    if constexpr (D == 2)
660
        ratio = ScalarT(2) / (ScalarT(2) + pi_scalar<ScalarT>); // round length = pi * r, flat length = 2 * r => ratio pi : 2
661
    if constexpr (D == 3)
662
        ratio = ScalarT(1) / ScalarT(3); // round area = 2 * pi * r^2, flat area = pi * r^2 => ratio 2 : 1
663

664
    if (detail::uniform01<ScalarT>(rng) >= ratio)
665
666
        return uniform(rng, boundary_no_caps_of(h));

667
    return uniform(rng, caps_of(h));
668
}
Aaron Grabowy's avatar
Aaron Grabowy committed
669

670
template <class ScalarT, class Rng>
671
672
673
674
675
676
677
678
679
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, cone_boundary<3, ScalarT> const& c)
{
    // base area = pi * r^^2, lateral area (mantle) = pi * r * sqrt(r^2 + h^2) => ratio r : sqrt(r^2 + h^2)
    auto r = c.base.radius;
    ScalarT ratio = r / (r + sqrt(pow2(r) + pow2(c.height)));
    if (detail::uniform01<ScalarT>(rng) >= ratio)
        return uniform(rng, boundary_no_caps_of(c));

    return uniform(rng, caps_of(c));
680
681
}

682
template <class ScalarT, class Rng>
683
[[nodiscard]] constexpr pos<3, ScalarT> uniform(Rng& rng, cone_boundary_no_caps<3, ScalarT> const& c)
684
{
Aaron Grabowy's avatar
Aaron Grabowy committed
685
    auto ub = aabb<2, ScalarT>::minus_one_to_one;
686
687
688
689
690
691
    while (true)
    {
        auto p = uniform_vec(rng, ub);
        auto l = length_sqr(p);
        if (l <= ScalarT(1))
        {
692
            p *= c.base.radius;
693
694
            auto x = any_normal(c.base.normal);
            auto y = cross(c.base.normal, x);
Aaron Grabowy's avatar
Aaron Grabowy committed
695
            return c.base.center + p.x * x + p.y * y + (ScalarT(1) - sqrt(l)) * c.base.normal * c.height;
696
697
698
699
        }
    }
}

700
701
702
// All solid pyramid variants
template <class BaseT, class Rng>
[[nodiscard]] constexpr pos<3, typename BaseT::scalar_t> uniform(Rng& rng, pyramid<BaseT> const& p)
703
{
704
    using ScalarT = typename BaseT::scalar_t;
Aaron Grabowy's avatar
Aaron Grabowy committed
705
    const auto n = normal_of(p.base);
706
    const auto h = ScalarT(1) - cbrt(detail::uniform01<ScalarT>(rng));
707
    const auto pBase = uniform(rng, p.base);
708
    return mix(pBase, centroid_of(p.base), h) + h * p.height * n;
709
710
}

711
// All boundary and boundary_no_caps pyramid variants except cones
712
713
template <class BaseT, class Rng>
[[nodiscard]] constexpr pos<3, typename BaseT::scalar_t> uniform(Rng& rng, pyramid_boundary<BaseT> const& py)
714
{
715
716
717
718
719
720
721
722
723
724
725
726
727
728
    const auto triangles = faces_of(boundary_no_caps_of(py));
    if constexpr (triangles.size() == 3)
        return uniform_by_area(rng, py.base, triangles[0], triangles[1], triangles[2]);
    else // verts.size() == 4
        return uniform_by_area(rng, py.base, triangles[0], triangles[1], triangles[2], triangles[3]);
}
template <class BaseT, class Rng>
[[nodiscard]] constexpr pos<3, typename BaseT::scalar_t> uniform(Rng& rng, pyramid_boundary_no_caps<BaseT> const& py)
{
    const auto triangles = faces_of(py);
    if constexpr (triangles.size() == 3)
        return uniform_by_area(rng, triangles[0], triangles[1], triangles[2]);
    else // verts.size() == 4
        return uniform_by_area(rng, triangles[0], triangles[1], triangles[2], triangles[3]);
729
730
}

731
template <int D, class ScalarT, class Rng, class = enable_if<is_floating_point<ScalarT>>>
732
[[deprecated("potentially misleading operation. use uniform_vec(rng, tg::aabb3(..)) or uniform_vec(rng, tg::segment3(..)) depending on your intended "
733
             "semantics")]] [[nodiscard]] constexpr vec<D, ScalarT>
734
uniform(Rng& rng, vec<D, ScalarT> const& a, vec<D, ScalarT> const& b)
Philip Trettner's avatar
Philip Trettner committed
735
736
737
{
    return mix(a, b, detail::uniform01<ScalarT>(rng));
}
738
template <int D, class ScalarT, class Rng, class = enable_if<is_floating_point<ScalarT>>>
739
[[deprecated("potentially misleading operation. use uniform(rng, tg::aabb3(..)) or uniform(rng, tg::segment3(..)) depending on your intended "
740
             "semantics")]] [[nodiscard]] constexpr pos<D, ScalarT>
741
uniform(Rng& rng, pos<D, ScalarT> const& a, pos<D, ScalarT> const& b)
Philip Trettner's avatar
Philip Trettner committed
742
743
744
{
    return mix(a, b, detail::uniform01<ScalarT>(rng));
}
Kersten Schuster's avatar
Kersten Schuster committed
745

746
template <class Rng, class... Args>
747
[[nodiscard]] constexpr auto uniform_vec(Rng& rng, Args const&... args) -> decltype(uniform(rng, args...) - decltype(uniform(rng, args...))::zero)
748
749
750
751
{
    return uniform(rng, args...) - decltype(uniform(rng, args...))::zero;
}

752
753
754
755
756
757
758
759
760
761
762
763
764
765
/// returns uniformly sampled barycentric coordinates
/// i.e. uniform(rng, triangle) has the same distribution as
///      triangle[uniform_barycoords(rng)]
template <class ScalarT = float, class Rng>
[[nodiscard]] constexpr comp<3, ScalarT> uniform_barycoords(Rng& rng)
{
    auto const a = detail::uniform01<ScalarT>(rng);
    auto const b = detail::uniform01<ScalarT>(rng);
    if (a + b <= ScalarT(1))
        return {a, b, 1 - a - b};
    else
        return {1 - a, 1 - b, a + b - 1};
}

Philip Trettner's avatar
Philip Trettner committed
766
767
768
769
770
771
772
773
774
775
776
777

namespace detail
{
template <>
struct sampler<bool>
{
    template <class Rng>
    constexpr static bool uniform(Rng& rng)
    {
        return rng() & 1;
    }
};
Philip Trettner's avatar
Philip Trettner committed
778
779
780
781
782
783
784
785
786
template <>
struct sampler<std::byte>
{
    template <class Rng>
    constexpr static std::byte uniform(Rng& rng)
    {
        return std::byte(rng() & 0xFF);
    }
};
787
template <class T>
788
struct sampler<angle_t<T>>
789
790
{
    template <class Rng>
791
    constexpr static angle_t<T> uniform(Rng& rng)
792
793
794
795
    {
        return tg::uniform(rng, tg::radians(T(0)), 2 * tg::pi<T>);
    }
};
Philip Trettner's avatar
Philip Trettner committed
796
797
798
799
800
801
template <class T>
struct sampler<color<3, T>>
{
    template <class Rng>
    constexpr static color<3, T> uniform(Rng& rng)
    {
Aaron Grabowy's avatar
Aaron Grabowy committed
802
        return color<3, T>(tg::uniform(rng, T(0), T(1)), tg::uniform(rng, T(0), T(1)), tg::uniform(rng, T(0), T(1)));
Philip Trettner's avatar
Philip Trettner committed
803
804
805
806
807
808
809
810
    }
};
template <class T>
struct sampler<color<4, T>>
{
    template <class Rng>
    constexpr static color<4, T> uniform(Rng& rng)
    {
Aaron Grabowy's avatar
Aaron Grabowy committed
811
        return color<4, T>(tg::uniform(rng, T(0), T(1)), tg::uniform(rng, T(0), T(1)), tg::uniform(rng, T(0), T(1)), tg::uniform(rng, T(0), T(1)));
Philip Trettner's avatar
Philip Trettner committed
812
813
    }
};
Philip Trettner's avatar
Philip Trettner committed
814
815
816
817
818
819
template <int D, class ScalarT>
struct sampler<dir<D, ScalarT>>
{
    template <class Rng>
    constexpr static dir<D, ScalarT> uniform(Rng& rng)
    {
Aaron Grabowy's avatar
Aaron Grabowy committed
820
        auto ub = aabb<D, ScalarT>::minus_one_to_one;
Philip Trettner's avatar
Philip Trettner committed
821
822
823
        while (true)
        {
            auto p = uniform_vec(rng, ub);
824
            auto l = length_sqr(p);
Philip Trettner's avatar
Philip Trettner committed
825
            if (l > ScalarT(0) && l <= ScalarT(1))
826
                return tg::dir<D, ScalarT>(p / sqrt(l));
Philip Trettner's avatar
Philip Trettner committed
827
828
829
830
831
        }
    }
};
}

Philip Trettner's avatar
Philip Trettner committed
832
} // namespace tg