test_JosephsMethod.cpp 60.1 KB
Newer Older
Nikola Dinev's avatar
Nikola Dinev committed
1 2 3 4
#include <catch2/catch.hpp>

#include "JosephsMethod.h"
#include "Geometry.h"
5 6
#include "Logger.h"
#include "testHelpers.h"
Nikola Dinev's avatar
Nikola Dinev committed
7 8 9 10 11

using namespace elsa;

SCENARIO("Testing BinaryVoxelTraversal with only one ray")
{
12 13 14
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Nikola Dinev's avatar
Nikola Dinev committed
15 16 17 18 19 20 21 22 23
    IndexVector_t sizeDomain(2);
    sizeDomain << 5, 5;

    IndexVector_t sizeRange(2);
    sizeRange << 1, 1;

    auto domain = DataDescriptor(sizeDomain);
    auto range = DataDescriptor(sizeRange);

Jens Petit's avatar
Jens Petit committed
24 25
    Eigen::IOFormat CommaInitFmt(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ", ", "", "",
                                 " << ", ";");
Nikola Dinev's avatar
Nikola Dinev committed
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

    GIVEN("A JosephsMethod for 2D and a domain data with all ones")
    {
        std::vector<Geometry> geom;

        auto dataDomain = DataContainer(domain);
        dataDomain = 1;

        auto dataRange = DataContainer(range);
        dataRange = 0;

        WHEN("We have a single ray with 0 degrees")
        {
            geom.emplace_back(100, 5, 0, domain, range);
            auto op = JosephsMethod(domain, range, geom);

            THEN("A^t A x should be close to the original data")
            {
                auto AtAx = DataContainer(domain);

                op.apply(dataDomain, dataRange);

                REQUIRE(dataRange[0] == Approx(5));

                op.applyAdjoint(dataRange, AtAx);

                auto cmp = RealVector_t(sizeDomain.prod());
Jens Petit's avatar
Jens Petit committed
53 54 55
                cmp << 0, 0, 5, 0, 0, 0, 0, 5, 0, 0, 0, 0, 5, 0, 0, 0, 0, 5, 0, 0, 0, 0, 5, 0, 0;

                REQUIRE(DataContainer(domain, cmp) == AtAx);
Nikola Dinev's avatar
Nikola Dinev committed
56 57 58 59 60
            }
        }

        WHEN("We have a single ray with 180 degrees")
        {
61
            geom.emplace_back(100, 5, pi_t, domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
62 63 64 65 66 67 68 69 70 71 72 73 74
            auto op = JosephsMethod(domain, range, geom);

            THEN("A^t A x should be close to the original data")
            {
                auto AtAx = DataContainer(domain);

                op.apply(dataDomain, dataRange);

                REQUIRE(dataRange[0] == Approx(5));

                op.applyAdjoint(dataRange, AtAx);

                auto cmp = RealVector_t(sizeDomain.prod());
Jens Petit's avatar
Jens Petit committed
75
                cmp << 0, 0, 5, 0, 0, 0, 0, 5, 0, 0, 0, 0, 5, 0, 0, 0, 0, 5, 0, 0, 0, 0, 5, 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
76

Jens Petit's avatar
Jens Petit committed
77
                REQUIRE(isApprox(DataContainer(domain, cmp), AtAx));
Nikola Dinev's avatar
Nikola Dinev committed
78 79 80 81 82
            }
        }

        WHEN("We have a single ray with 90 degrees")
        {
83
            geom.emplace_back(100, 5, pi_t / 2, domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
84 85 86 87 88 89 90 91 92 93 94 95 96
            auto op = JosephsMethod(domain, range, geom);

            THEN("A^t A x should be close to the original data")
            {
                auto AtAx = DataContainer(domain);

                op.apply(dataDomain, dataRange);

                REQUIRE(dataRange[0] == Approx(5));

                op.applyAdjoint(dataRange, AtAx);

                auto cmp = RealVector_t(sizeDomain.prod());
Jens Petit's avatar
Jens Petit committed
97
                cmp << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
98

Jens Petit's avatar
Jens Petit committed
99
                REQUIRE(isApprox(DataContainer(domain, cmp), AtAx));
Nikola Dinev's avatar
Nikola Dinev committed
100 101 102 103 104
            }
        }

        WHEN("We have a single ray with 90 degrees")
        {
105
            geom.emplace_back(100, 5, 3 * pi_t / 2., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
106 107 108 109 110 111 112 113 114 115 116 117 118
            auto op = JosephsMethod(domain, range, geom);

            THEN("A^t A x should be close to the original data")
            {
                auto AtAx = DataContainer(domain);

                op.apply(dataDomain, dataRange);

                REQUIRE(dataRange[0] == Approx(5));

                op.applyAdjoint(dataRange, AtAx);

                auto cmp = RealVector_t(sizeDomain.prod());
Jens Petit's avatar
Jens Petit committed
119
                cmp << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
120

Jens Petit's avatar
Jens Petit committed
121
                REQUIRE(isApprox(DataContainer(domain, cmp), AtAx));
Nikola Dinev's avatar
Nikola Dinev committed
122 123 124 125 126
            }
        }

        WHEN("We have a single ray with 90 degrees")
        {
127
            geom.emplace_back(100, 5, 45 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
128 129 130 131 132 133 134 135 136 137
            auto op = JosephsMethod(domain, range, geom);

            THEN("A^t A x should be close to the original data")
            {
                auto AtAx = DataContainer(domain);

                op.apply(dataDomain, dataRange);
                op.applyAdjoint(dataRange, AtAx);

                auto cmp = RealVector_t(sizeDomain.prod());
Jens Petit's avatar
Jens Petit committed
138 139
                cmp << 10, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0,
                    10;
Nikola Dinev's avatar
Nikola Dinev committed
140

141
                REQUIRE(isApprox(DataContainer(domain, cmp), AtAx, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
142 143 144 145 146
            }
        }

        WHEN("We have a single ray with 90 degrees")
        {
147
            geom.emplace_back(100, 5, 225 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
148 149 150 151 152 153 154 155 156 157
            auto op = JosephsMethod(domain, range, geom);

            THEN("A^t A x should be close to the original data")
            {
                auto AtAx = DataContainer(domain);

                op.apply(dataDomain, dataRange);
                op.applyAdjoint(dataRange, AtAx);

                auto cmp = RealVector_t(sizeDomain.prod());
Jens Petit's avatar
Jens Petit committed
158 159
                cmp << 10, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0,
                    10;
Nikola Dinev's avatar
Nikola Dinev committed
160

161
                REQUIRE(isApprox(DataContainer(domain, cmp), AtAx, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
162 163 164 165 166 167 168
            }
        }
    }
}

SCENARIO("Testing JosephsMethod with only 4 ray")
{
169 170 171
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Nikola Dinev's avatar
Nikola Dinev committed
172 173 174 175 176 177 178 179 180
    IndexVector_t sizeDomain(2);
    sizeDomain << 5, 5;

    IndexVector_t sizeRange(2);
    sizeRange << 1, 4;

    auto domain = DataDescriptor(sizeDomain);
    auto range = DataDescriptor(sizeRange);

Jens Petit's avatar
Jens Petit committed
181 182
    Eigen::IOFormat CommaInitFmt(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ", ", "", "",
                                 " << ", ";");
Nikola Dinev's avatar
Nikola Dinev committed
183 184 185 186 187 188 189 190 191 192 193 194 195 196

    GIVEN("A JosephsMethod for 2D and a domain data with all ones")
    {
        std::vector<Geometry> geom;

        auto dataDomain = DataContainer(domain);
        dataDomain = 1;

        auto dataRange = DataContainer(range);
        dataRange = 0;

        WHEN("We have a single ray with 0, 90, 180, 270 degrees")
        {
            geom.emplace_back(100, 5, 0, domain, range);
197 198 199
            geom.emplace_back(100, 5, pi_t / 2, domain, range);
            geom.emplace_back(100, 5, pi_t, domain, range);
            geom.emplace_back(100, 5, 3 * pi_t / 2, domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
200 201 202 203 204 205 206 207 208 209
            auto op = JosephsMethod(domain, range, geom);

            THEN("A^t A x should be close to the original data")
            {
                auto AtAx = DataContainer(domain);

                op.apply(dataDomain, dataRange);
                op.applyAdjoint(dataRange, AtAx);

                auto cmp = RealVector_t(sizeDomain.prod());
Jens Petit's avatar
Jens Petit committed
210 211
                cmp << 0, 0, 10, 0, 0, 0, 0, 10, 0, 0, 10, 10, 20, 10, 10, 0, 0, 10, 0, 0, 0, 0, 10,
                    0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
212

Jens Petit's avatar
Jens Petit committed
213
                REQUIRE(isApprox(DataContainer(domain, cmp), AtAx));
Nikola Dinev's avatar
Nikola Dinev committed
214 215 216 217 218
            }
        }
    }
}

Jens Petit's avatar
Jens Petit committed
219 220
SCENARIO("Calls to functions of super class")
{
221 222 223
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
224 225
    GIVEN("A projector")
    {
Nikola Dinev's avatar
Nikola Dinev committed
226 227 228 229
        IndexVector_t volumeDims(2), sinoDims(2);
        const index_t volSize = 50;
        const index_t detectorSize = 50;
        const index_t numImgs = 50;
Jens Petit's avatar
Jens Petit committed
230 231
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
232 233 234
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
235
        volume = 0;
Nikola Dinev's avatar
Nikola Dinev committed
236
        DataContainer sino(sinoDescriptor);
237
        sino = 0;
Nikola Dinev's avatar
Nikola Dinev committed
238
        std::vector<Geometry> geom;
Jens Petit's avatar
Jens Petit committed
239
        for (std::size_t i = 0; i < numImgs; i++) {
240
            real_t angle = static_cast<real_t>(i * 2) * pi_t / 50;
Jens Petit's avatar
Jens Petit committed
241
            geom.emplace_back(20 * volSize, volSize, angle, volumeDescriptor, sinoDescriptor);
Nikola Dinev's avatar
Nikola Dinev committed
242 243
        }

Jens Petit's avatar
Jens Petit committed
244
        JosephsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
245

Jens Petit's avatar
Jens Petit committed
246 247
        WHEN("Projector is cloned")
        {
Nikola Dinev's avatar
Nikola Dinev committed
248 249 250 251
            auto opClone = op.clone();
            auto sinoClone = sino;
            auto volumeClone = volume;

Jens Petit's avatar
Jens Petit committed
252 253 254 255
            THEN("Results do not change (may still be slightly different due to summation being "
                 "performed in a different order)")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
256

Jens Petit's avatar
Jens Petit committed
257 258
                opClone->apply(volume, sinoClone);
                REQUIRE(isApprox(sino, sinoClone));
Nikola Dinev's avatar
Nikola Dinev committed
259

Jens Petit's avatar
Jens Petit committed
260 261
                op.applyAdjoint(sino, volume);
                opClone->applyAdjoint(sino, volumeClone);
Jens Petit's avatar
Jens Petit committed
262 263 264

                DataContainer resultsDifference = volume - volumeClone;
                REQUIRE(resultsDifference.squaredL2Norm() == Approx(0.0).margin(1e-5));
Nikola Dinev's avatar
Nikola Dinev committed
265 266 267 268 269
            }
        }
    }
}

Jens Petit's avatar
Jens Petit committed
270 271
SCENARIO("Output DataContainer is not zero initialized")
{
272 273 274
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
275 276
    GIVEN("A 2D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
277 278 279 280
        IndexVector_t volumeDims(2), sinoDims(2);
        const index_t volSize = 5;
        const index_t detectorSize = 1;
        const index_t numImgs = 1;
Jens Petit's avatar
Jens Petit committed
281 282
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
283 284 285 286 287
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;
Jens Petit's avatar
Jens Petit committed
288 289 290
        geom.emplace_back(20 * volSize, volSize, 0.0, volumeDescriptor, sinoDescriptor);
        JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                         JosephsMethod<>::Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
291

Jens Petit's avatar
Jens Petit committed
292 293
        WHEN("Sinogram conatainer is not zero initialized and we project through an empty volume")
        {
Nikola Dinev's avatar
Nikola Dinev committed
294 295 296
            volume = 0;
            sino = 1;

Jens Petit's avatar
Jens Petit committed
297 298 299
            THEN("Result is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
300 301

                DataContainer zero(sinoDescriptor);
302
                zero = 0;
303
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
304 305 306
            }
        }

Jens Petit's avatar
Jens Petit committed
307 308
        WHEN("Volume container is not zero initialized and we backproject from an empty sinogram")
        {
Nikola Dinev's avatar
Nikola Dinev committed
309 310 311
            sino = 0;
            volume = 1;

Jens Petit's avatar
Jens Petit committed
312 313 314
            THEN("Result is zero")
            {
                op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
315 316

                DataContainer zero(volumeDescriptor);
317
                zero = 0;
318
                REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
319 320 321 322
            }
        }
    }

Jens Petit's avatar
Jens Petit committed
323 324
    GIVEN("A 3D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
325 326 327 328
        IndexVector_t volumeDims(3), sinoDims(3);
        const index_t volSize = 3;
        const index_t detectorSize = 1;
        const index_t numImgs = 1;
Jens Petit's avatar
Jens Petit committed
329 330
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
331 332 333 334 335 336
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
337 338 339
        geom.emplace_back(volSize * 20, volSize, volumeDescriptor, sinoDescriptor, 0);
        JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                         JosephsMethod<>::Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
340

Jens Petit's avatar
Jens Petit committed
341 342
        WHEN("Sinogram conatainer is not zero initialized and we project through an empty volume")
        {
Nikola Dinev's avatar
Nikola Dinev committed
343 344 345
            volume = 0;
            sino = 1;

Jens Petit's avatar
Jens Petit committed
346 347 348
            THEN("Result is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
349
                DataContainer zero(sinoDescriptor);
350
                zero = 0;
351
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
352 353 354
            }
        }

Jens Petit's avatar
Jens Petit committed
355 356
        WHEN("Volume container is not zero initialized and we backproject from an empty sinogram")
        {
Nikola Dinev's avatar
Nikola Dinev committed
357 358 359
            sino = 0;
            volume = 1;

Jens Petit's avatar
Jens Petit committed
360 361 362
            THEN("Result is zero")
            {
                op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
363
                DataContainer zero(volumeDescriptor);
364
                zero = 0;
365
                REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
366 367 368 369 370
            }
        }
    }
}

Jens Petit's avatar
Jens Petit committed
371 372
SCENARIO("Rays not intersecting the bounding box are present")
{
373 374 375
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
376 377
    GIVEN("A 2D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
378 379 380 381
        IndexVector_t volumeDims(2), sinoDims(2);
        const index_t volSize = 5;
        const index_t detectorSize = 1;
        const index_t numImgs = 1;
Jens Petit's avatar
Jens Petit committed
382 383
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
384 385 386 387 388 389 390 391
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        volume = 1;
        sino = 1;
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
392 393 394 395
        WHEN("Tracing along a y-axis-aligned ray with a negative x-coordinate of origin")
        {
            geom.emplace_back(20 * volSize, volSize, 0.0, volumeDescriptor, sinoDescriptor, 0.0,
                              volSize);
Nikola Dinev's avatar
Nikola Dinev committed
396

Jens Petit's avatar
Jens Petit committed
397 398
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                             JosephsMethod<>::Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
399

Jens Petit's avatar
Jens Petit committed
400 401 402
            THEN("Result of forward projection is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
403
                DataContainer zero(sinoDescriptor);
404
                zero = 0;
405
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
406

Jens Petit's avatar
Jens Petit committed
407 408 409
                AND_THEN("Result of backprojection is zero")
                {
                    op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
410
                    DataContainer zero(volumeDescriptor);
411
                    zero = 0;
412
                    REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
413 414 415 416
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
417 418 419 420 421 422 423 424
        WHEN("Tracing along a y-axis-aligned ray with a x-coordinate of origin beyond the bounding "
             "box")
        {
            geom.emplace_back(20 * volSize, volSize, 0.0, volumeDescriptor, sinoDescriptor, 0.0,
                              -volSize);

            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                             JosephsMethod<>::Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
425

Jens Petit's avatar
Jens Petit committed
426 427 428
            THEN("Result of forward projection is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
429
                DataContainer zero(sinoDescriptor);
430
                zero = 0;
431
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
432

Jens Petit's avatar
Jens Petit committed
433 434 435
                AND_THEN("Result of backprojection is zero")
                {
                    op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
436
                    DataContainer zero(volumeDescriptor);
437
                    zero = 0;
438
                    REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
439 440 441 442
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
443 444
        WHEN("Tracing along a x-axis-aligned ray with a negative y-coordinate of origin")
        {
445 446
            geom.emplace_back(20 * volSize, volSize, pi_t / 2, volumeDescriptor, sinoDescriptor,
                              0.0, 0.0, volSize);
Nikola Dinev's avatar
Nikola Dinev committed
447

Jens Petit's avatar
Jens Petit committed
448 449
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                             JosephsMethod<>::Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
450

Jens Petit's avatar
Jens Petit committed
451 452 453
            THEN("Result of forward projection is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
454
                DataContainer zero(sinoDescriptor);
455
                zero = 0;
456
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
457

Jens Petit's avatar
Jens Petit committed
458 459 460
                AND_THEN("Result of backprojection is zero")
                {
                    op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
461
                    DataContainer zero(volumeDescriptor);
462
                    zero = 0;
463
                    REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
464 465 466 467
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
468 469 470
        WHEN("Tracing along a x-axis-aligned ray with a y-coordinate of origin beyond the bounding "
             "box")
        {
471 472
            geom.emplace_back(20 * volSize, volSize, pi_t / 2, volumeDescriptor, sinoDescriptor,
                              0.0, 0.0, -volSize);
Nikola Dinev's avatar
Nikola Dinev committed
473

Jens Petit's avatar
Jens Petit committed
474 475
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                             JosephsMethod<>::Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
476

Jens Petit's avatar
Jens Petit committed
477 478 479
            THEN("Result of forward projection is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
480
                DataContainer zero(sinoDescriptor);
481
                zero = 0;
482
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
483

Jens Petit's avatar
Jens Petit committed
484 485 486
                AND_THEN("Result of backprojection is zero")
                {
                    op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
487
                    DataContainer zero(volumeDescriptor);
488
                    zero = 0;
489
                    REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
490 491 492 493 494
                }
            }
        }
    }

Jens Petit's avatar
Jens Petit committed
495 496
    GIVEN("A 3D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
497 498 499 500
        IndexVector_t volumeDims(3), sinoDims(3);
        const index_t volSize = 5;
        const index_t detectorSize = 1;
        const index_t numImgs = 1;
Jens Petit's avatar
Jens Petit committed
501 502
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
503 504 505 506 507 508 509 510 511
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        volume = 1;
        sino = 1;
        std::vector<Geometry> geom;

        const index_t numCases = 9;
Jens Petit's avatar
Jens Petit committed
512
        real_t alpha[numCases] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
513 514 515
        real_t beta[numCases] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, pi_t / 2, pi_t / 2, pi_t / 2};
        real_t gamma[numCases] = {0.0,      0.0,      0.0,      pi_t / 2, pi_t / 2,
                                  pi_t / 2, pi_t / 2, pi_t / 2, pi_t / 2};
Jens Petit's avatar
Jens Petit committed
516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534
        real_t offsetx[numCases] = {volSize, 0.0, volSize, 0.0, 0.0, 0.0, volSize, 0.0, volSize};
        real_t offsety[numCases] = {0.0, volSize, volSize, volSize, 0.0, volSize, 0.0, 0.0, 0.0};
        real_t offsetz[numCases] = {0.0, 0.0, 0.0, 0.0, volSize, volSize, 0.0, volSize, volSize};
        std::string neg[numCases] = {"x", "y", "x and y", "y", "z", "y and z", "x", "z", "x and z"};
        std::string ali[numCases] = {"z", "z", "z", "x", "x", "x", "y", "y", "y"};

        for (int i = 0; i < numCases; i++) {
            WHEN("Tracing along a " + ali[i] + "-axis-aligned ray with negative " + neg[i]
                 + "-coodinate of origin")
            {
                geom.emplace_back(20 * volSize, volSize, volumeDescriptor, sinoDescriptor, gamma[i],
                                  beta[i], alpha[i], 0.0, 0.0, offsetx[i], offsety[i], offsetz[i]);

                JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                                 JosephsMethod<>::Interpolation::LINEAR);

                THEN("Result of forward projection is zero")
                {
                    op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
535
                    DataContainer zero(sinoDescriptor);
536
                    zero = 0;
537
                    REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
538

Jens Petit's avatar
Jens Petit committed
539 540 541
                    AND_THEN("Result of backprojection is zero")
                    {
                        op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
542
                        DataContainer zero(volumeDescriptor);
543
                        zero = 0;
544
                        REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
545 546 547 548 549 550 551 552 553
                    }
                }
            }
        }
    }
}

SCENARIO("Axis-aligned rays are present")
{
554 555 556
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
557 558
    GIVEN("A 2D setting with a single ray")
    {
Nikola Dinev's avatar
Nikola Dinev committed
559 560 561 562
        IndexVector_t volumeDims(2), sinoDims(2);
        const index_t volSize = 5;
        const index_t detectorSize = 1;
        const index_t numImgs = 1;
Jens Petit's avatar
Jens Petit committed
563 564
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
565 566 567 568 569 570 571
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

        const index_t numCases = 4;
572
        const real_t angles[numCases] = {0.0, pi_t / 2, pi_t, 3 * pi_t / 2};
Jens Petit's avatar
Jens Petit committed
573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591
        RealVector_t backProj[2];
        backProj[0].resize(volSize * volSize);
        backProj[1].resize(volSize * volSize);
        backProj[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;

        backProj[0] << 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0;

        for (index_t i = 0; i < numCases; i++) {
            WHEN("An axis-aligned ray with an angle of " + std::to_string(angles[i])
                 + " radians passes through the center of a pixel")
            {
                geom.emplace_back(volSize * 20, volSize, angles[i], volumeDescriptor,
                                  sinoDescriptor);

                JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                                 JosephsMethod<>::Interpolation::LINEAR);
                THEN("The result of projecting through a pixel is exactly the pixel value")
                {
                    for (index_t j = 0; j < volSize; j++) {
Nikola Dinev's avatar
Nikola Dinev committed
592
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
593 594
                        if (i % 2 == 0)
                            volume(volSize / 2, j) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
595
                        else
Jens Petit's avatar
Jens Petit committed
596 597 598
                            volume(j, volSize / 2) = 1;

                        op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
599 600 601
                        REQUIRE(sino[0] == Approx(1));
                    }

Jens Petit's avatar
Jens Petit committed
602 603 604 605 606
                    AND_THEN("The backprojection sets the values of all hit pixels to the detector "
                             "value")
                    {
                        op.applyAdjoint(sino, volume);
                        REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[i % 2])));
Nikola Dinev's avatar
Nikola Dinev committed
607 608 609 610 611
                    }
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632
        real_t offsetx[numCases] = {0.25, 0.0, 0.25, 0.0};
        real_t offsety[numCases] = {0.0, 0.25, 0.0, 0.25};

        backProj[0] << 0, 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0.25, 0.75,
            0, 0, 0, 0.25, 0.75, 0, 0;

        backProj[1] << 0, 0, 0, 0, 0, 0.25, 0.25, 0.25, 0.25, 0.25, 0.75, 0.75, 0.75, 0.75, 0.75, 0,
            0, 0, 0, 0, 0, 0, 0, 0, 0;

        for (index_t i = 0; i < numCases; i++) {
            WHEN("An axis-aligned ray with an angle of " + std::to_string(angles[i])
                 + " radians does not pass through the center of a pixel")
            {
                geom.emplace_back(volSize * 2000, volSize, angles[i], volumeDescriptor,
                                  sinoDescriptor, 0.0, -offsetx[i], -offsety[i]);
                JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                                 JosephsMethod<>::Interpolation::LINEAR);
                THEN("The result of projecting through a pixel is the interpolated value between "
                     "the two pixels closest to the ray")
                {
                    for (index_t j = 0; j < volSize; j++) {
Nikola Dinev's avatar
Nikola Dinev committed
633
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
634 635
                        if (i % 2 == 0)
                            volume(volSize / 2, j) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
636
                        else
Jens Petit's avatar
Jens Petit committed
637 638 639
                            volume(j, volSize / 2) = 1;

                        op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
640 641 642
                        REQUIRE(sino[0] == Approx(0.75));
                    }

Jens Petit's avatar
Jens Petit committed
643 644
                    AND_THEN("The backprojection yields the exact adjoint")
                    {
Nikola Dinev's avatar
Nikola Dinev committed
645
                        sino[0] = 1;
Jens Petit's avatar
Jens Petit committed
646 647
                        op.applyAdjoint(sino, volume);
                        REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[i % 2])));
Nikola Dinev's avatar
Nikola Dinev committed
648 649 650 651
                    }
                }
            }
        }
Jens Petit's avatar
Jens Petit committed
652 653 654 655 656 657 658 659 660 661 662 663 664

        backProj[0] << 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1;

        WHEN("A y-axis-aligned ray runs along the right volume boundary")
        {
            geom.emplace_back(volSize * 2000, volSize, 0.0, volumeDescriptor, sinoDescriptor, 0.0,
                              (volSize * 0.5));
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom);

            THEN("The result of projecting through a pixel is exactly the pixel's value (we mirror "
                 "values at the border for the purpose of interpolation)")
            {
                for (index_t j = 0; j < volSize; j++) {
665
                    volume = 0;
Jens Petit's avatar
Jens Petit committed
666 667 668
                    volume(volSize - 1, j) = 1;

                    op.apply(volume, sino);
669 670
                    REQUIRE(sino[0] == 1);
                }
Nikola Dinev's avatar
Nikola Dinev committed
671

Jens Petit's avatar
Jens Petit committed
672 673 674 675
                AND_THEN(
                    "The slow backprojection yields the exact adjoint, the fast backprojection "
                    "also yields the exact adjoint for a very distant x-ray source")
                {
Nikola Dinev's avatar
Nikola Dinev committed
676
                    sino[0] = 1;
Jens Petit's avatar
Jens Petit committed
677 678
                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[0])));
Nikola Dinev's avatar
Nikola Dinev committed
679 680 681
                }
            }
        }
Jens Petit's avatar
Jens Petit committed
682 683 684 685 686 687 688 689 690 691 692 693 694

        backProj[0] << 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0;

        WHEN("A y-axis-aligned ray runs along the left volume boundary")
        {
            geom.emplace_back(volSize * 2000, volSize, 0.0, volumeDescriptor, sinoDescriptor, 0.0,
                              -volSize / 2.0);
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom);

            THEN("The result of projecting through a pixel is exactly the pixel's value (we mirror "
                 "values at the border for the purpose of interpolation)")
            {
                for (index_t j = 0; j < volSize; j++) {
695
                    volume = 0;
Jens Petit's avatar
Jens Petit committed
696 697 698
                    volume(0, j) = 1;

                    op.apply(volume, sino);
699 700
                    REQUIRE(sino[0] == 1);
                }
Nikola Dinev's avatar
Nikola Dinev committed
701

Jens Petit's avatar
Jens Petit committed
702 703 704 705
                AND_THEN(
                    "The slow backprojection yields the exact adjoint, the fast backprojection "
                    "also yields the exact adjoint for a very distant x-ray source")
                {
Nikola Dinev's avatar
Nikola Dinev committed
706
                    sino[0] = 1;
Jens Petit's avatar
Jens Petit committed
707 708
                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[0])));
Nikola Dinev's avatar
Nikola Dinev committed
709 710 711 712 713
                }
            }
        }
    }

Jens Petit's avatar
Jens Petit committed
714 715
    GIVEN("A 3D setting with a single ray")
    {
Nikola Dinev's avatar
Nikola Dinev committed
716 717 718 719
        IndexVector_t volumeDims(3), sinoDims(3);
        const index_t volSize = 3;
        const index_t detectorSize = 1;
        const index_t numImgs = 1;
Jens Petit's avatar
Jens Petit committed
720 721
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
722 723 724 725 726 727 728
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

        const index_t numCases = 6;
729 730
        real_t beta[numCases] = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
        real_t gamma[numCases] = {0.0, pi_t, pi_t / 2, 3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
Jens Petit's avatar
Jens Petit committed
731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763
        std::string al[numCases] = {"z", "-z", "x", "-x", "y", "-y"};

        RealVector_t backProj[numCases];
        for (auto& backPr : backProj)
            backPr.resize(volSize * volSize * volSize);
        backProj[2] << 0, 0, 0, 0, 0, 0, 0, 0, 0,

            0, 1, 0, 0, 1, 0, 0, 1, 0,

            0, 0, 0, 0, 0, 0, 0, 0, 0;

        backProj[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0,

            0, 0, 0, 1, 1, 1, 0, 0, 0,

            0, 0, 0, 0, 0, 0, 0, 0, 0;

        backProj[0] << 0, 0, 0, 0, 1, 0, 0, 0, 0,

            0, 0, 0, 0, 1, 0, 0, 0, 0,

            0, 0, 0, 0, 1, 0, 0, 0, 0;

        for (index_t i = 0; i < numCases; i++) {
            WHEN("A " + al[i] + "-axis-aligned ray passes through the center of a pixel")
            {
                geom.emplace_back(volSize * 20, volSize, volumeDescriptor, sinoDescriptor, gamma[i],
                                  beta[i]);
                JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                                 JosephsMethod<>::Interpolation::LINEAR);
                THEN("The result of projecting through a voxel is exactly the voxel value")
                {
                    for (index_t j = 0; j < volSize; j++) {
Nikola Dinev's avatar
Nikola Dinev committed
764
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
765 766 767 768 769 770 771 772 773
                        if (i / 2 == 0)
                            volume(volSize / 2, volSize / 2, j) = 1;
                        else if (i / 2 == 1)
                            volume(j, volSize / 2, volSize / 2) = 1;
                        else if (i / 2 == 2)
                            volume(volSize / 2, j, volSize / 2) = 1;

                        op.apply(volume, sino);
                        REQUIRE(sino[0] == Approx(1.0));
Nikola Dinev's avatar
Nikola Dinev committed
774 775
                    }

Jens Petit's avatar
Jens Petit committed
776 777 778 779 780
                    AND_THEN("The backprojection sets the values of all hit voxels to the detector "
                             "value")
                    {
                        op.applyAdjoint(sino, volume);
                        REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[i / 2])));
Nikola Dinev's avatar
Nikola Dinev committed
781 782 783 784 785
                    }
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819
        real_t offsetx[numCases] = {0.25, 0.25, 0.0, 0.0, 0.0, 0.0};
        real_t offsety[numCases] = {0.0, 0.0, 0.25, 0.25, 0.0, 0.0};
        real_t offsetz[numCases] = {0.0, 0.0, 0.0, 0.0, 0.25, 0.25};

        backProj[2] << 0, 0.25, 0, 0, 0.25, 0, 0, 0.25, 0,

            0, 0.75, 0, 0, 0.75, 0, 0, 0.75, 0,

            0, 0, 0, 0, 0, 0, 0, 0, 0;

        backProj[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0,

            0.25, 0.25, 0.25, 0.75, 0.75, 0.75, 0, 0, 0,

            0, 0, 0, 0, 0, 0, 0, 0, 0;

        backProj[0] << 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0,

            0, 0, 0, 0.25, 0.75, 0, 0, 0, 0,

            0, 0, 0, 0.25, 0.75, 0, 0, 0, 0;

        for (index_t i = 0; i < numCases; i++) {
            WHEN("A " + al[i] + "-axis-aligned ray does not pass through the center of a voxel")
            {
                geom.emplace_back(volSize * 2000, volSize, volumeDescriptor, sinoDescriptor,
                                  gamma[i], beta[i], 0.0, 0.0, 0.0, -offsetx[i], -offsety[i],
                                  -offsetz[i]);
                JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                                 JosephsMethod<>::Interpolation::LINEAR);
                THEN("The result of projecting through a voxel is the interpolated value between "
                     "the four voxels nearest to the ray")
                {
                    for (index_t j = 0; j < volSize; j++) {
Nikola Dinev's avatar
Nikola Dinev committed
820
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
821 822 823 824 825 826 827 828 829
                        if (i / 2 == 0)
                            volume(volSize / 2, volSize / 2, j) = 1;
                        else if (i / 2 == 1)
                            volume(j, volSize / 2, volSize / 2) = 1;
                        else if (i / 2 == 2)
                            volume(volSize / 2, j, volSize / 2) = 1;

                        op.apply(volume, sino);
                        REQUIRE(sino[0] == Approx(0.75));
Nikola Dinev's avatar
Nikola Dinev committed
830 831
                    }

Jens Petit's avatar
Jens Petit committed
832 833 834
                    AND_THEN("The backprojection yields the exact adjoint")
                    {
                        sino[0] = 1;
Nikola Dinev's avatar
Nikola Dinev committed
835

Jens Petit's avatar
Jens Petit committed
836 837
                        op.applyAdjoint(sino, volume);
                        REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[i / 2])));
Nikola Dinev's avatar
Nikola Dinev committed
838 839 840 841 842
                    }
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
843 844
        offsetx[0] = volSize / 2.0;
        offsetx[1] = -(volSize / 2.0);
Nikola Dinev's avatar
Nikola Dinev committed
845 846
        offsetx[2] = 0.0;
        offsetx[3] = 0.0;
Jens Petit's avatar
Jens Petit committed
847 848
        offsetx[4] = -(volSize / 2.0);
        offsetx[5] = (volSize / 2.0);
Nikola Dinev's avatar
Nikola Dinev committed
849 850 851

        offsety[0] = 0.0;
        offsety[1] = 0.0;
Jens Petit's avatar
Jens Petit committed
852 853 854 855
        offsety[2] = volSize / 2.0;
        offsety[3] = -(volSize / 2.0);
        offsety[4] = -(volSize / 2.0);
        offsety[5] = (volSize / 2.0);
Nikola Dinev's avatar
Nikola Dinev committed
856 857 858

        al[0] = "left border";
        al[1] = "right border";
859 860
        al[2] = "bottom border";
        al[3] = "top border";
Nikola Dinev's avatar
Nikola Dinev committed
861 862
        al[4] = "top right edge";
        al[5] = "bottom left edge";
Jens Petit's avatar
Jens Petit committed
863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911

        backProj[0] << 0, 0, 0, 1, 0, 0, 0, 0, 0,

            0, 0, 0, 1, 0, 0, 0, 0, 0,

            0, 0, 0, 1, 0, 0, 0, 0, 0;

        backProj[1] << 0, 0, 0, 0, 0, 1, 0, 0, 0,

            0, 0, 0, 0, 0, 1, 0, 0, 0,

            0, 0, 0, 0, 0, 1, 0, 0, 0;

        backProj[2] << 0, 1, 0, 0, 0, 0, 0, 0, 0,

            0, 1, 0, 0, 0, 0, 0, 0, 0,

            0, 1, 0, 0, 0, 0, 0, 0, 0;

        backProj[3] << 0, 0, 0, 0, 0, 0, 0, 1, 0,

            0, 0, 0, 0, 0, 0, 0, 1, 0,

            0, 0, 0, 0, 0, 0, 0, 1, 0;

        backProj[4] << 0, 0, 0, 0, 0, 0, 0, 0, 1,

            0, 0, 0, 0, 0, 0, 0, 0, 1,

            0, 0, 0, 0, 0, 0, 0, 0, 1;

        backProj[5] << 1, 0, 0, 0, 0, 0, 0, 0, 0,

            1, 0, 0, 0, 0, 0, 0, 0, 0,

            1, 0, 0, 0, 0, 0, 0, 0, 0;

        for (index_t i = 0; i < numCases; i++) {
            WHEN("A z-axis-aligned ray runs along the " + al[i] + " of the volume")
            {
                // x-ray source must be very far from the volume center to make testing of the fast
                // backprojection simpler
                geom.emplace_back(volSize * 2000, volSize, volumeDescriptor, sinoDescriptor, 0.0,
                                  0.0, 0.0, 0.0, 0.0, -offsetx[i], -offsety[i]);
                JosephsMethod op(volumeDescriptor, sinoDescriptor, geom);
                THEN("The result of projecting through a voxel is exactly the voxel's value (we "
                     "mirror values at the border for the purpose of interpolation)")
                {
                    for (index_t j = 0; j < volSize; j++) {
912
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933
                        switch (i) {
                            case 0:
                                volume(0, volSize / 2, j) = 1;
                                break;
                            case 1:
                                volume(volSize - 1, volSize / 2, j) = 1;
                                break;
                            case 2:
                                volume(volSize / 2, 0, j) = 1;
                                break;
                            case 3:
                                volume(volSize / 2, volSize - 1, j) = 1;
                                break;
                            case 4:
                                volume(volSize - 1, volSize - 1, j) = 1;
                                break;
                            case 5:
                                volume(0, 0, j) = 1;
                                break;
                            default:
                                break;
934
                        }
Jens Petit's avatar
Jens Petit committed
935 936 937

                        op.apply(volume, sino);
                        REQUIRE(sino[0] == 1);
938
                    }
Nikola Dinev's avatar
Nikola Dinev committed
939

Jens Petit's avatar
Jens Petit committed
940 941 942
                    AND_THEN("The backprojection yields the exact adjoint")
                    {
                        sino[0] = 1;
943

Jens Petit's avatar
Jens Petit committed
944 945
                        op.applyAdjoint(sino, volume);
                        REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[i])));
Nikola Dinev's avatar
Nikola Dinev committed
946 947 948 949 950 951
                    }
                }
            }
        }
    }

Jens Petit's avatar
Jens Petit committed
952 953
    GIVEN("A 3D setting with multiple projection angles")
    {
Nikola Dinev's avatar
Nikola Dinev committed
954 955 956 957
        IndexVector_t volumeDims(3), sinoDims(3);
        const index_t volSize = 3;
        const index_t detectorSize = 1;
        const index_t numImgs = 6;
Jens Petit's avatar
Jens Petit committed
958 959
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
960 961 962 963 964 965
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
966 967
        WHEN("x-, y and z-axis-aligned rays are present")
        {
968 969
            real_t beta[numImgs] = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
            real_t gamma[numImgs] = {0.0, pi_t, pi_t / 2, 3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
Nikola Dinev's avatar
Nikola Dinev committed
970

Jens Petit's avatar
Jens Petit committed
971 972 973
            for (index_t i = 0; i < numImgs; i++)
                geom.emplace_back(volSize * 20, volSize, volumeDescriptor, sinoDescriptor, gamma[i],
                                  beta[i]);
Nikola Dinev's avatar
Nikola Dinev committed
974

Jens Petit's avatar
Jens Petit committed
975 976
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                             JosephsMethod<>::Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
977

Jens Petit's avatar
Jens Petit committed
978 979
            THEN("Values are accumulated correctly along each ray's path")
            {
Nikola Dinev's avatar
Nikola Dinev committed
980 981
                volume = 0;

Jens Petit's avatar
Jens Petit committed
982 983 984 985 986 987
                // set only values along the rays' path to one to make sure interpolation is done
                // correctly
                for (index_t i = 0; i < volSize; i++) {
                    volume(i, volSize / 2, volSize / 2) = 1;
                    volume(volSize / 2, i, volSize / 2) = 1;
                    volume(volSize / 2, volSize / 2, i) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
988 989
                }

Jens Petit's avatar
Jens Petit committed
990 991 992 993 994 995 996
                op.apply(volume, sino);
                for (index_t i = 0; i < numImgs; i++)
                    REQUIRE(sino[i] == Approx(3.0));

                AND_THEN("Backprojections yield the exact adjoint")
                {
                    RealVector_t cmp(volSize * volSize * volSize);
Nikola Dinev's avatar
Nikola Dinev committed
997

Jens Petit's avatar
Jens Petit committed
998
                    cmp << 0, 0, 0, 0, 6, 0, 0, 0, 0,
Nikola Dinev's avatar
Nikola Dinev committed
999

Jens Petit's avatar
Jens Petit committed
1000
                        0, 6, 0, 6, 18, 6, 0, 6, 0,
Nikola Dinev's avatar
Nikola Dinev committed
1001

Jens Petit's avatar
Jens Petit committed
1002
                        0, 0, 0, 0, 6, 0, 0, 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
1003

Jens Petit's avatar
Jens Petit committed
1004 1005
                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, cmp)));
Nikola Dinev's avatar
Nikola Dinev committed
1006 1007 1008 1009 1010 1011
                }
            }
        }
    }
}

Jens Petit's avatar
Jens Petit committed
1012 1013
SCENARIO("Projection under an angle")
{
1014 1015 1016
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
1017 1018
    GIVEN("A 2D setting with a single ray")
    {
Nikola Dinev's avatar
Nikola Dinev committed
1019 1020 1021 1022
        IndexVector_t volumeDims(2), sinoDims(2);
        const index_t volSize = 4;
        const index_t detectorSize = 1;
        const index_t numImgs = 1;
Jens Petit's avatar
Jens Petit committed
1023 1024
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
1025 1026 1027 1028 1029 1030
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
1031 1032 1033 1034
        WHEN("Projecting under an angle of 30 degrees and ray goes through center of volume")
        {
            // In this case the ray enters and exits the volume through the borders along the main
            // direction Weighting for all interpolated values should be the same
1035
            geom.emplace_back(volSize * 20, volSize, -pi_t / 6, volumeDescriptor, sinoDescriptor);
Jens Petit's avatar
Jens Petit committed
1036 1037 1038
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                             JosephsMethod<>::Interpolation::LINEAR);

1039
            real_t weight = static_cast<real_t>(2 / std::sqrt(3.f));
Jens Petit's avatar
Jens Petit committed
1040 1041
            THEN("Ray intersects the correct pixels")
            {
Nikola Dinev's avatar
Nikola Dinev committed
1042
                volume = 1;
Jens Petit's avatar
Jens Petit committed
1043 1044 1045 1046
                volume(3, 0) = 0;
                volume(1, 1) = 0;
                volume(1, 2) = 0;
                volume(1, 3) = 0;
Nikola Dinev's avatar
Nikola Dinev committed
1047

Jens Petit's avatar
Jens Petit committed
1048 1049 1050 1051
                volume(2, 0) = 0;
                volume(2, 1) = 0;
                volume(2, 2) = 0;
                volume(0, 3) = 0;
Nikola Dinev's avatar
Nikola Dinev committed
1052 1053 1054

                op.apply(volume, sino);
                DataContainer zero(sinoDescriptor);
1055
                zero = 0;
1056
                CHECK(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
1057

Jens Petit's avatar
Jens Petit committed
1058 1059 1060 1061 1062 1063
                AND_THEN("The correct weighting is applied")
                {
                    volume(3, 0) = 1;
                    volume(1, 1) = 1;
                    volume(1, 2) = 1;
                    volume(1, 3) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
1064

Jens Petit's avatar
Jens Petit committed
1065 1066
                    op.apply(volume, sino);
                    REQUIRE(sino[0] == Approx(2 * weight));
Nikola Dinev's avatar
Nikola Dinev committed
1067 1068 1069

                    sino[0] = 1;

Jens Petit's avatar
Jens Petit committed
1070
                    RealVector_t opExpected(volSize * volSize);
1071 1072 1073 1074 1075 1076
                    opExpected << 0, 0, (3 - std::sqrt(3.f)) / 2, (std::sqrt(3.f) - 1) / 2, 0,
                        (std::sqrt(3.f) - 1) / (2 * std::sqrt(3.f)),
                        (std::sqrt(3.f) + 1) / (2 * std::sqrt(3.f)), 0, 0,
                        (std::sqrt(3.f) + 1) / (2 * std::sqrt(3.f)),
                        (std::sqrt(3.f) - 1) / (2 * std::sqrt(3.f)), 0, (std::sqrt(3.f) - 1) / 2,
                        (3 - std::sqrt(3.f)) / 2, 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
1077

Jens Petit's avatar
Jens Petit committed
1078 1079 1080
                    opExpected *= weight;
                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, opExpected)));
Nikola Dinev's avatar
Nikola Dinev committed
1081 1082 1083 1084
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
1085 1086 1087 1088 1089
        WHEN("Projecting under an angle of 30 degrees and ray enters through the right border")
        {
            // In this case the ray exits through a border along the main ray direction, but enters
            // through a border not along the main direction First pixel should be weighted
            // differently
1090 1091
            geom.emplace_back(volSize * 20, volSize, -pi_t / 6, volumeDescriptor, sinoDescriptor,
                              0.0, std::sqrt(3.f));
Jens Petit's avatar
Jens Petit committed
1092 1093 1094 1095 1096
            JosephsMethod op(volumeDescriptor, sinoDescriptor, geom,
                             JosephsMethod<>::Interpolation::LINEAR);

            THEN("Ray intersects the correct pixels")
            {
Nikola Dinev's avatar
Nikola Dinev committed
1097
                volume = 1;
Jens Petit's avatar
Jens Petit committed
1098 1099 1100 1101 1102 1103
                volume(3, 1) = 0;
                volume(3, 2) = 0;
                volume(3, 3) = 0;
                volume(2, 3) = 0;
                volume(2, 2) = 0;

Nikola Dinev's avatar
Nikola Dinev committed
1104 1105
                op.apply(volume, sino);
                DataContainer zero(sinoDescriptor);
1106
                zero = 0;
1107
                CHECK(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
1108

Jens Petit's avatar
Jens Petit committed
1109 1110 1111 1112 1113
                AND_THEN("The correct weighting is applied")
                {
                    volume(3, 1) = 1;
                    volume(2, 2)