Luts.hpp 4.4 KB
Newer Older
David Frank's avatar
David Frank committed
1
2
3
#pragma once

#include "Blobs.h"
David Frank's avatar
David Frank committed
4
#include "BSplines.h"
David Frank's avatar
David Frank committed
5
#include "Logger.h"
David Frank's avatar
David Frank committed
6
#include "Timer.h"
David Frank's avatar
David Frank committed
7
8
9
10
11
12
13
14
15
16

#include <array>

namespace elsa
{
    namespace detail
    {
        template <typename data_t, index_t N>
        constexpr std::array<data_t, N> blob_lut(ProjectedBlob<data_t> blob)
        {
David Frank's avatar
David Frank committed
17
18
            Logger::get("blob_lut")->debug("Calculating lut");

David Frank's avatar
David Frank committed
19
20
21
22
23
24
25
26
27
28
29
30
31
            std::array<data_t, N> lut;

            auto t = static_cast<data_t>(0);
            const auto step = blob.radius() / N;

            for (std::size_t i = 0; i < N; ++i) {
                lut[i] = blob(t);
                t += step;
            }

            return lut;
        }

David Frank's avatar
David Frank committed
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
        template <typename data_t, index_t N>
        constexpr std::array<data_t, N> bspline_lut(ProjectedBSpline<data_t> bspline)
        {
            Logger::get("bspline_lut")->debug("Calculating lut");

            std::array<data_t, N> lut;

            auto t = static_cast<data_t>(0);
            const auto step = 2. / N;

            for (std::size_t i = 0; i < N; ++i) {
                lut[i] = bspline(t);
                t += step;
            }

            return lut;
        }

David Frank's avatar
David Frank committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
        template <typename data_t>
        data_t lerp(data_t a, SelfType_t<data_t> b, SelfType_t<data_t> t)
        {
            if ((a <= 0 && b >= 0) || (a >= 0 && b <= 0))
                return t * b + (1 - t) * a;

            if (t == 1)
                return b;

            const data_t x = a + t * (b - a);

            if ((t > 1) == (b > a))
                return b < x ? x : b;
            else
                return x < b ? x : b;
        }
    } // namespace detail

    template <typename data_t, std::size_t N>
    class Lut
    {
    public:
David Frank's avatar
David Frank committed
72
        constexpr Lut(std::array<data_t, N> data) : data_(std::move(data)) {}
David Frank's avatar
David Frank committed
73
74

        template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
David Frank's avatar
David Frank committed
75
        constexpr data_t operator()(T index) const
David Frank's avatar
David Frank committed
76
77
78
79
80
81
82
83
84
85
86
87
        {
            if (index < 0 || index > N) {
                return 0;
            }

            return data_[index];
        }

        /// TODO: Handle boundary conditions
        /// lerp(last, last+1, t), for some t > 0, yields f(last) / 2, as f(last + 1) = 0,
        /// this should be handled
        template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
David Frank's avatar
David Frank committed
88
        constexpr data_t operator()(T index) const
David Frank's avatar
David Frank committed
89
90
91
92
93
94
95
96
97
98
        {
            if (index < 0 || index > N) {
                return 0;
            }

            // Get the two closes indices
            const auto a = static_cast<std::size_t>(std::floor(index));
            const auto b = static_cast<std::size_t>(std::ceil(index));

            // Get the function values
David Frank's avatar
David Frank committed
99
100
101
102
            const auto fa = a == N ? 0 : data_[a];
            const auto fb = b == N ? 0 : data_[b];

            auto ret = detail::lerp(fa, fb, index - static_cast<data_t>(a));
David Frank's avatar
David Frank committed
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

            // Bilinear interpolation
            return detail::lerp(fa, fb, index - static_cast<data_t>(a));
        }

    private:
        std::array<data_t, N> data_;
    };

    // User defined deduction guide
    template <typename data_t, std::size_t N>
    Lut(std::array<data_t, N>) -> Lut<data_t, N>;

    template <typename data_t, index_t N>
    class ProjectedBlobLut
    {
    public:
        constexpr ProjectedBlobLut(data_t radius, SelfType_t<data_t> alpha,
                                   SelfType_t<data_t> order)
            : blob_(radius, alpha, order), lut_(detail::blob_lut<data_t, N>(blob_))
        {
        }

David Frank's avatar
David Frank committed
126
        constexpr data_t radius() const { return blob_.radius(); }
David Frank's avatar
David Frank committed
127

David Frank's avatar
David Frank committed
128
        constexpr data_t alpha() const { return blob_.alpha(); }
David Frank's avatar
David Frank committed
129

David Frank's avatar
David Frank committed
130
        constexpr data_t order() const { return blob_.order(); }
David Frank's avatar
David Frank committed
131

David Frank's avatar
David Frank committed
132
133
134
135
        constexpr data_t operator()(data_t distance) const
        {
            return lut_((distance / blob_.radius()) * N);
        }
David Frank's avatar
David Frank committed
136
137
138
139
140

    private:
        ProjectedBlob<data_t> blob_;
        Lut<data_t, N> lut_;
    };
David Frank's avatar
David Frank committed
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

    template <typename data_t, index_t N>
    class ProjectedBSplineLut
    {
    public:
        constexpr ProjectedBSplineLut(int dim, int degree)
            : bspline_(dim, degree), lut_(detail::bspline_lut<data_t, N>(bspline_))
        {
        }

        constexpr data_t order() const { return bspline_.order(); }

        constexpr data_t operator()(data_t distance) const
        {
            return lut_((std::abs(distance) / 2.) * N);
        }

    private:
        ProjectedBSpline<data_t> bspline_;
        Lut<data_t, N> lut_;
    };
David Frank's avatar
David Frank committed
162
} // namespace elsa