24.09., 9:00 - 11:00: Due to updates GitLab will be unavailable for some minutes between 09:00 and 11:00.

test_SiddonsMethod.cpp 55.5 KB
Newer Older
Nikola Dinev's avatar
Nikola Dinev committed
1 2 3 4
#include "catch2/catch.hpp"

#include "SiddonsMethod.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 SiddonsMethod projector 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

    GIVEN("A SiddonsMethod 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 = SiddonsMethod(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
50
                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
51

Jens Petit's avatar
Jens Petit committed
52
                REQUIRE(DataContainer(domain, cmp) == AtAx);
Nikola Dinev's avatar
Nikola Dinev committed
53 54 55 56 57 58
                REQUIRE(dataRange[0] == Approx(5));
            }
        }

        WHEN("We have a single ray with 180 degrees")
        {
59
            geom.emplace_back(100, 5, 180 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
60 61 62 63 64 65 66 67 68 69
            auto op = SiddonsMethod(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
70
                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
71

Jens Petit's avatar
Jens Petit committed
72
                REQUIRE(DataContainer(domain, cmp) == AtAx);
Nikola Dinev's avatar
Nikola Dinev committed
73 74 75 76 77 78
                REQUIRE(dataRange[0] == Approx(5));
            }
        }

        WHEN("We have a single ray with 90 degrees")
        {
79
            geom.emplace_back(100, 5, 90 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
80 81 82 83 84 85 86 87 88 89
            auto op = SiddonsMethod(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
90
                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
91

Jens Petit's avatar
Jens Petit committed
92
                REQUIRE(DataContainer(domain, cmp) == AtAx);
Nikola Dinev's avatar
Nikola Dinev committed
93 94 95 96 97 98
                REQUIRE(dataRange[0] == Approx(5));
            }
        }

        WHEN("We have a single ray with 270 degrees")
        {
99
            geom.emplace_back(100, 5, 270 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
100 101 102 103 104 105 106 107 108 109
            auto op = SiddonsMethod(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
110
                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
111

Jens Petit's avatar
Jens Petit committed
112
                REQUIRE(DataContainer(domain, cmp) == AtAx);
Nikola Dinev's avatar
Nikola Dinev committed
113 114 115 116
                REQUIRE(dataRange[0] == Approx(5));
            }
        }

Jens Petit's avatar
Jens Petit committed
117 118
        // FIXME This does not yield the desired result/if fixed the overall results in a
        // reconstruction is bad
Nikola Dinev's avatar
Nikola Dinev committed
119 120
        /*WHEN("We have a single ray with 45 degrees")
        {
121
            geom.emplace_back(100, 5, 45 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
            auto op = SiddonsMethod(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);

                std::cout << dataRange.getData().format(CommaInitFmt) << "\n";
                std::cout << AtAx.getData().format(CommaInitFmt) << "\n";

                auto cmp = RealVector_t(sizeDomain.prod());
                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;

                REQUIRE(cmp.isApprox(AtAx.getData(), 1e-2));
                REQUIRE(dataRange[0] == Approx(7.071));
            }
        }

        WHEN("We have a single ray with 225 degrees")
        {
148
            geom.emplace_back(100, 5, 225 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
            auto op = SiddonsMethod(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);

                std::cout << dataRange.getData().format(CommaInitFmt) << "\n";
                std::cout << AtAx.getData().format(CommaInitFmt) << "\n";

                auto cmp = RealVector_t(sizeDomain.prod());
                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;

                REQUIRE(cmp.isApprox(AtAx.getData()));
                REQUIRE(dataRange[0] == Approx(7.071));
            }
        }*/

Jens Petit's avatar
Jens Petit committed
173 174
        // TODO fix this direction, currently it does not work correctly. Consider changing
        // geometry, to mirror stuff
Nikola Dinev's avatar
Nikola Dinev committed
175 176
        /*WHEN("We have a single ray with 135 degrees")
        {
177
            geom.emplace_back(100, 5, 135 * pi_t / 180., domain, range);
Nikola Dinev's avatar
Nikola Dinev committed
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
            auto op = SiddonsMethod(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);

                std::cout << dataRange.getData().format(CommaInitFmt) << "\n";
                std::cout << AtAx.getData().format(CommaInitFmt) << "\n";

                auto cmp = RealVector_t(sizeDomain.prod());
                cmp << 0, 0, 0, 0, 10,
                        0, 0, 0, 10, 0,
                        0, 0, 10, 0, 0,
                        0, 10, 0, 0, 0,
                        10, 0, 0, 0, 0;

                REQUIRE(cmp.isApprox(AtAx.getData()));
                REQUIRE(dataRange[0] == Approx(7.071));
            }
        }*/
    }
}

Jens Petit's avatar
Jens Petit committed
204 205
SCENARIO("Calls to functions of super class")
{
206 207 208
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
209 210
    GIVEN("A projector")
    {
Nikola Dinev's avatar
Nikola Dinev committed
211 212 213 214
        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
215 216
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
217 218 219 220 221 222 223 224
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);

        RealVector_t randomStuff(volumeDescriptor.getNumberOfCoefficients());
        randomStuff.setRandom();
        DataContainer volume(volumeDescriptor, randomStuff);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;
Jens Petit's avatar
Jens Petit committed
225
        for (std::size_t i = 0; i < numImgs; i++) {
226
            real_t angle = static_cast<real_t>(i * 2) * pi_t / 50;
Jens Petit's avatar
Jens Petit committed
227
            geom.emplace_back(20 * volSize, volSize, angle, volumeDescriptor, sinoDescriptor);
Nikola Dinev's avatar
Nikola Dinev committed
228
        }
Jens Petit's avatar
Jens Petit committed
229
        SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
230

Jens Petit's avatar
Jens Petit committed
231 232
        WHEN("Projector is cloned")
        {
Nikola Dinev's avatar
Nikola Dinev committed
233 234 235 236
            auto opClone = op.clone();
            auto sinoClone = DataContainer(sino.getDataDescriptor());
            auto volumeClone = DataContainer(volume.getDataDescriptor());

Jens Petit's avatar
Jens Petit committed
237 238 239 240 241 242
            THEN("Results do not change (may still be slightly different due to summation being "
                 "performed in a different order)")
            {
                op.apply(volume, sino);
                opClone->apply(volume, sinoClone);
                REQUIRE(isApprox(sino, sinoClone));
Nikola Dinev's avatar
Nikola Dinev committed
243

Jens Petit's avatar
Jens Petit committed
244 245
                op.applyAdjoint(sino, volume);
                opClone->applyAdjoint(sino, volumeClone);
Jens Petit's avatar
Jens Petit committed
246 247 248

                DataContainer resultsDifference = volume - volumeClone;
                REQUIRE(resultsDifference.squaredL2Norm() == Approx(0.0).margin(1e-5));
Nikola Dinev's avatar
Nikola Dinev committed
249 250 251 252 253
            }
        }
    }
}

Jens Petit's avatar
Jens Petit committed
254 255
SCENARIO("Output DataContainer is not zero initialized")
{
256 257 258
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
259 260
    GIVEN("A 2D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
261 262 263 264
        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
265 266
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
267 268 269 270 271
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;
Jens Petit's avatar
Jens Petit committed
272 273
        geom.emplace_back(20 * volSize, volSize, 0.0, volumeDescriptor, sinoDescriptor);
        SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
274

Jens Petit's avatar
Jens Petit committed
275 276
        WHEN("Sinogram conatainer is not zero initialized and we project through an empty volume")
        {
Nikola Dinev's avatar
Nikola Dinev committed
277 278 279
            volume = 0;
            sino = 1;

Jens Petit's avatar
Jens Petit committed
280 281 282
            THEN("Result is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
283
                DataContainer zero(sinoDescriptor);
284
                zero = 0;
285
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
286 287 288
            }
        }

Jens Petit's avatar
Jens Petit committed
289 290
        WHEN("Volume container is not zero initialized and we backproject from an empty sinogram")
        {
Nikola Dinev's avatar
Nikola Dinev committed
291 292 293
            sino = 0;
            volume = 1;

Jens Petit's avatar
Jens Petit committed
294 295 296
            THEN("Result is zero")
            {
                op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
297
                DataContainer zero(volumeDescriptor);
298
                zero = 0;
299
                REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
300 301 302 303
            }
        }
    }

Jens Petit's avatar
Jens Petit committed
304 305
    GIVEN("A 3D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
306 307 308 309
        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
310 311
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
312 313 314 315 316 317
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
318 319
        geom.emplace_back(volSize * 20, volSize, volumeDescriptor, sinoDescriptor, 0);
        SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
320

Jens Petit's avatar
Jens Petit committed
321 322
        WHEN("Sinogram conatainer is not zero initialized and we project through an empty volume")
        {
Nikola Dinev's avatar
Nikola Dinev committed
323 324 325
            volume = 0;
            sino = 1;

Jens Petit's avatar
Jens Petit committed
326 327 328
            THEN("Result is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
329
                DataContainer zero(sinoDescriptor);
330
                zero = 0;
331
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
332 333 334
            }
        }

Jens Petit's avatar
Jens Petit committed
335 336
        WHEN("Volume container is not zero initialized and we backproject from an empty sinogram")
        {
Nikola Dinev's avatar
Nikola Dinev committed
337 338 339
            sino = 0;
            volume = 1;

Jens Petit's avatar
Jens Petit committed
340 341 342
            THEN("Result is zero")
            {
                op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
343
                DataContainer zero(volumeDescriptor);
344
                zero = 0;
345
                REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
346 347 348 349 350
            }
        }
    }
}

Jens Petit's avatar
Jens Petit committed
351 352
SCENARIO("Rays not intersecting the bounding box are present")
{
353 354 355
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
356 357
    GIVEN("A 2D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
358 359 360 361
        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
362 363
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
364 365 366 367 368 369 370 371
        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
372 373 374 375
        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
376

Jens Petit's avatar
Jens Petit committed
377
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
378

Jens Petit's avatar
Jens Petit committed
379 380 381
            THEN("Result of forward projection is zero")
            {
                op.apply(volume, sino);
Nikola Dinev's avatar
Nikola Dinev committed
382
                DataContainer zero(sinoDescriptor);
383
                zero = 0;
384
                REQUIRE(isApprox(sino, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
385

Jens Petit's avatar
Jens Petit committed
386 387 388
                AND_THEN("Result of backprojection is zero")
                {
                    op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
389
                    DataContainer zero(volumeDescriptor);
390
                    zero = 0;
391
                    REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
392 393 394 395
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
396 397 398 399 400
        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);
Nikola Dinev's avatar
Nikola Dinev committed
401

Jens Petit's avatar
Jens Petit committed
402
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
403

Jens Petit's avatar
Jens Petit committed
404 405 406 407
            THEN("Result of forward projection is zero")
            {
                op.apply(volume, sino);
                DataContainer zero(sinoDescriptor);
408
                zero = 0;
409
                REQUIRE(isApprox(sino, zero, epsilon));
Jens Petit's avatar
Jens Petit committed
410 411 412 413 414

                AND_THEN("Result of backprojection is zero")
                {
                    op.applyAdjoint(sino, volume);
                    DataContainer zero(volumeDescriptor);
415
                    zero = 0;
416
                    REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
417 418 419 420
                }
            }
        }

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

Jens Petit's avatar
Jens Petit committed
426
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
427

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

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

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

Jens Petit's avatar
Jens Petit committed
451
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
452

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

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

Jens Petit's avatar
Jens Petit committed
471 472
    GIVEN("A 3D setting")
    {
Nikola Dinev's avatar
Nikola Dinev committed
473 474 475 476
        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
477 478
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
479 480 481 482 483 484 485 486 487
        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
488
        real_t alpha[numCases] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
489 490 491
        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
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511
        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]);

                SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);

                THEN("Result of forward projection is zero")
                {
                    op.apply(volume, sino);
                    DataContainer zero(sinoDescriptor);
512
                    zero = 0;
513
                    REQUIRE(isApprox(sino, zero, epsilon));
Jens Petit's avatar
Jens Petit committed
514 515 516 517 518

                    AND_THEN("Result of backprojection is zero")
                    {
                        op.applyAdjoint(sino, volume);
                        DataContainer zero(volumeDescriptor);
519
                        zero = 0;
520
                        REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
521 522 523 524 525 526 527 528 529
                    }
                }
            }
        }
    }
}

SCENARIO("Axis-aligned rays are present")
{
530 531 532
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
533 534
    GIVEN("A 2D setting with a single ray")
    {
Nikola Dinev's avatar
Nikola Dinev committed
535 536 537 538
        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
539 540
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
541 542 543 544 545 546 547
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

        const index_t numCases = 4;
548
        const real_t angles[numCases] = {0.0, pi_t / 2, pi_t, 3 * pi_t / 2};
Jens Petit's avatar
Jens Petit committed
549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565
        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);
                SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
                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
566
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
567 568
                        if (i % 2 == 0)
                            volume(volSize / 2, j) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
569
                        else
Jens Petit's avatar
Jens Petit committed
570 571 572 573
                            volume(j, volSize / 2) = 1;

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

Jens Petit's avatar
Jens Petit committed
576 577 578 579
                    AND_THEN("The backprojection sets the values of all hit pixels to the detector "
                             "value")
                    {
                        op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
580

Jens Petit's avatar
Jens Petit committed
581
                        REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[i % 2])));
Nikola Dinev's avatar
Nikola Dinev committed
582 583 584 585 586
                    }
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
587 588 589 590 591 592 593 594 595
        WHEN("A y-axis-aligned ray runs along a voxel boundary")
        {
            geom.emplace_back(volSize * 2000, volSize, 0.0, volumeDescriptor, sinoDescriptor, 0.0,
                              -0.5);
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
            THEN("The result of projecting through a pixel is the value of the pixel with the "
                 "higher index")
            {
                for (index_t j = 0; j < volSize; j++) {
Nikola Dinev's avatar
Nikola Dinev committed
596
                    volume = 0;
Jens Petit's avatar
Jens Petit committed
597 598 599
                    volume(volSize / 2, j) = 1;

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

Jens Petit's avatar
Jens Petit committed
603 604
                AND_THEN("The backprojection yields the exact adjoint")
                {
Nikola Dinev's avatar
Nikola Dinev committed
605
                    sino[0] = 1;
Jens Petit's avatar
Jens Petit committed
606 607
                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[0])));
Nikola Dinev's avatar
Nikola Dinev committed
608 609 610 611
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
612 613 614 615 616 617 618 619 620 621
        WHEN("A y-axis-aligned ray runs along the right volume boundary")
        {
            // For Siddon's values in the range [0,boxMax) are considered, a ray running along
            // boxMax should be ignored
            geom.emplace_back(volSize * 2000, volSize, 0.0, volumeDescriptor, sinoDescriptor, 0.0,
                              (volSize * 0.5));
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);

            THEN("The result of projecting is zero")
            {
Nikola Dinev's avatar
Nikola Dinev committed
622
                volume = 0;
Jens Petit's avatar
Jens Petit committed
623 624
                op.apply(volume, sino);
                REQUIRE(sino[0] == 0.0);
Nikola Dinev's avatar
Nikola Dinev committed
625

Jens Petit's avatar
Jens Petit committed
626 627
                AND_THEN("The result of backprojection is also zero")
                {
Nikola Dinev's avatar
Nikola Dinev committed
628 629
                    sino[0] = 1;

Jens Petit's avatar
Jens Petit committed
630 631
                    op.applyAdjoint(sino, volume);
                    DataContainer zero(volumeDescriptor);
632
                    zero = 0;
633
                    REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
634 635 636
                }
            }
        }
Jens Petit's avatar
Jens Petit committed
637

Nikola Dinev's avatar
Nikola Dinev committed
638
        /**
Jens Petit's avatar
Jens Petit committed
639 640
         * CPU version exhibits slightly different behaviour, ignoring all ray running along a
         * bounding box border
Nikola Dinev's avatar
Nikola Dinev committed
641 642 643 644 645 646
         */
        // 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;
Jens Petit's avatar
Jens Petit committed
647

Nikola Dinev's avatar
Nikola Dinev committed
648 649 650
        // 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);
        //     SiddonsMethod op(volumeDescriptor,sinoDescriptor,geom);
Jens Petit's avatar
Jens Petit committed
651
        //     THEN("The result of projecting through a pixel is exactly the pixel's value") {
Nikola Dinev's avatar
Nikola Dinev committed
652 653 654
        //         for (index_t j=0; j<volSize;j++) {
        //             volume = 0;
        //             volume(0,j) = 1;
Jens Petit's avatar
Jens Petit committed
655

Nikola Dinev's avatar
Nikola Dinev committed
656 657 658 659 660 661 662 663 664 665 666 667 668
        //             op.apply(volume,sino);
        //             REQUIRE(sino[0] == 1);
        //         }

        //         AND_THEN("The backprojection yields the exact adjoint") {
        //             sino[0] = 1;
        //             op.applyAdjoint(sino,volume);
        //             REQUIRE(volume.getData().isApprox(backProj[0]));
        //         }
        //     }
        // }
    }

Jens Petit's avatar
Jens Petit committed
669 670
    GIVEN("A 3D setting with a single ray")
    {
Nikola Dinev's avatar
Nikola Dinev committed
671 672 673 674
        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
675 676
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
677 678 679 680 681 682 683
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

        const index_t numCases = 6;
684 685
        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
686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718
        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]);
                SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
                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
719
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
720 721 722 723 724 725 726 727 728
                        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] == 1);
Nikola Dinev's avatar
Nikola Dinev committed
729 730
                    }

Jens Petit's avatar
Jens Petit committed
731 732 733 734 735
                    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
736 737 738 739 740 741 742 743
                    }
                }
            }
        }

        real_t offsetx[numCases];
        real_t offsety[numCases];

Jens Petit's avatar
Jens Petit committed
744 745
        offsetx[0] = volSize / 2.0;
        offsetx[3] = -(volSize / 2.0);
Nikola Dinev's avatar
Nikola Dinev committed
746 747
        offsetx[1] = 0.0;
        offsetx[4] = 0.0;
Jens Petit's avatar
Jens Petit committed
748 749
        offsetx[5] = -(volSize / 2.0);
        offsetx[2] = (volSize / 2.0);
Nikola Dinev's avatar
Nikola Dinev committed
750 751 752

        offsety[0] = 0.0;
        offsety[3] = 0.0;
Jens Petit's avatar
Jens Petit committed
753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774
        offsety[1] = volSize / 2.0;
        offsety[4] = -(volSize / 2.0);
        offsety[5] = -(volSize / 2.0);
        offsety[2] = (volSize / 2.0);

        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, 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[2] << 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;
Nikola Dinev's avatar
Nikola Dinev committed
775 776

        al[0] = "left border";
777 778 779 780 781
        al[1] = "bottom border";
        al[2] = "bottom left edge";
        al[3] = "right border";
        al[4] = "top border";
        al[5] = "top right edge";
Jens Petit's avatar
Jens Petit committed
782 783 784 785 786 787 788 789 790 791 792 793

        for (index_t i = 0; i < numCases / 2; 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 op
                // backprojection simpler
                geom.emplace_back(volSize * 2000, volSize, volumeDescriptor, sinoDescriptor, 0.0,
                                  0.0, 0.0, 0.0, 0.0, -offsetx[i], -offsety[i]);
                SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
                THEN("The result of projecting through a voxel is exactly the voxel's value")
                {
                    for (index_t j = 0; j < volSize; j++) {
794
                        volume = 0;
Jens Petit's avatar
Jens Petit committed
795 796 797 798 799 800 801 802 803 804 805 806 807
                        switch (i) {
                            case 0:
                                volume(0, volSize / 2, j) = 1;
                                break;
                            case 1:
                                volume(volSize / 2, 0, j) = 1;
                                break;
                                break;
                            case 2:
                                volume(0, 0, j) = 1;
                                break;
                            default:
                                break;
808
                        }
Jens Petit's avatar
Jens Petit committed
809 810 811

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

Jens Petit's avatar
Jens Petit committed
814 815 816 817
                    AND_THEN("The backprojection yields the exact adjoints")
                    {
                        sino[0] = 1;
                        op.applyAdjoint(sino, volume);
Nikola Dinev's avatar
Nikola Dinev committed
818

Jens Petit's avatar
Jens Petit committed
819
                        REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, backProj[i])));
820 821 822 823
                    }
                }
            }
        }
Nikola Dinev's avatar
Nikola Dinev committed
824

Jens Petit's avatar
Jens Petit committed
825 826 827 828 829 830 831 832 833 834
        for (index_t i = 3; 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 op
                // backprojection simpler
                geom.emplace_back(volSize * 2000, volSize, volumeDescriptor, sinoDescriptor, 0.0,
                                  0.0, 0.0, 0.0, 0.0, -offsetx[i], -offsety[i]);
                SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
                THEN("The result of projecting is zero")
                {
Nikola Dinev's avatar
Nikola Dinev committed
835 836
                    volume = 1;

Jens Petit's avatar
Jens Petit committed
837 838
                    op.apply(volume, sino);
                    REQUIRE(sino[0] == 0);
Nikola Dinev's avatar
Nikola Dinev committed
839

Jens Petit's avatar
Jens Petit committed
840 841 842 843 844 845
                    AND_THEN("The result of backprojection is also zero")
                    {
                        sino[0] = 1;
                        op.applyAdjoint(sino, volume);

                        DataContainer zero(volumeDescriptor);
846
                        zero = 0;
847
                        REQUIRE(isApprox(volume, zero, epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
848 849 850 851 852 853
                    }
                }
            }
        }
    }

Jens Petit's avatar
Jens Petit committed
854 855
    GIVEN("A 2D setting with multiple projection angles")
    {
Nikola Dinev's avatar
Nikola Dinev committed
856 857 858 859
        IndexVector_t volumeDims(2), sinoDims(2);
        const index_t volSize = 5;
        const index_t detectorSize = 1;
        const index_t numImgs = 4;
Jens Petit's avatar
Jens Petit committed
860 861
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
862 863 864 865 866 867
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
868 869 870
        WHEN("Both x- and y-axis-aligned rays are present")
        {
            geom.emplace_back(20 * volSize, volSize, 0, volumeDescriptor, sinoDescriptor);
871
            geom.emplace_back(20 * volSize, volSize, 90 * pi_t / 180., volumeDescriptor,
Jens Petit's avatar
Jens Petit committed
872
                              sinoDescriptor);
873
            geom.emplace_back(20 * volSize, volSize, 180 * pi_t / 180., volumeDescriptor,
Jens Petit's avatar
Jens Petit committed
874
                              sinoDescriptor);
875
            geom.emplace_back(20 * volSize, volSize, 270 * pi_t / 180., volumeDescriptor,
Jens Petit's avatar
Jens Petit committed
876
                              sinoDescriptor);
Nikola Dinev's avatar
Nikola Dinev committed
877

Jens Petit's avatar
Jens Petit committed
878
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
879

Jens Petit's avatar
Jens Petit committed
880 881
            THEN("Values are accumulated correctly along each ray's path")
            {
Nikola Dinev's avatar
Nikola Dinev committed
882 883
                volume = 0;

Jens Petit's avatar
Jens Petit committed
884 885 886 887 888
                // set only values along the rays' path to one to make sure interpolation is dones
                // correctly
                for (index_t i = 0; i < volSize; i++) {
                    volume(i, volSize / 2) = 1;
                    volume(volSize / 2, i) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
889
                }
Jens Petit's avatar
Jens Petit committed
890 891 892 893 894 895 896 897 898 899 900 901 902 903

                op.apply(volume, sino);
                for (index_t i = 0; i < numImgs; i++)
                    REQUIRE(sino[i] == Approx(5.0));

                AND_THEN("Backprojection yields the exact adjoint")
                {
                    RealVector_t cmp(volSize * volSize);

                    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;

                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, cmp)));
Nikola Dinev's avatar
Nikola Dinev committed
904 905 906 907 908
                }
            }
        }
    }

Jens Petit's avatar
Jens Petit committed
909 910
    GIVEN("A 3D setting with multiple projection angles")
    {
Nikola Dinev's avatar
Nikola Dinev committed
911 912 913 914
        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
915 916
        volumeDims << volSize, volSize, volSize;
        sinoDims << detectorSize, detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
917 918 919 920 921 922
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
923 924
        WHEN("x-, y and z-axis-aligned rays are present")
        {
925 926
            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
927

Jens Petit's avatar
Jens Petit committed
928 929 930
            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
931

Jens Petit's avatar
Jens Petit committed
932
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
933

Jens Petit's avatar
Jens Petit committed
934 935
            THEN("Values are accumulated correctly along each ray's path")
            {
Nikola Dinev's avatar
Nikola Dinev committed
936 937
                volume = 0;

Jens Petit's avatar
Jens Petit committed
938 939 940 941 942 943
                // set only values along the rays' path to one to make sure interpolation is dones
                // 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
944
                }
Jens Petit's avatar
Jens Petit committed
945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961

                op.apply(volume, sino);
                for (index_t i = 0; i < numImgs; i++)
                    REQUIRE(sino[i] == Approx(3.0));

                AND_THEN("Backprojection yields the exact adjoint")
                {
                    RealVector_t cmp(volSize * volSize * volSize);

                    cmp << 0, 0, 0, 0, 6, 0, 0, 0, 0,

                        0, 6, 0, 6, 18, 6, 0, 6, 0,

                        0, 0, 0, 0, 6, 0, 0, 0, 0;

                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, cmp)));
Nikola Dinev's avatar
Nikola Dinev committed
962 963 964 965 966 967
                }
            }
        }
    }
}

Jens Petit's avatar
Jens Petit committed
968 969
SCENARIO("Projection under an angle")
{
970 971 972
    // Turn logger of
    Logger::setLevel(Logger::LogLevel::OFF);

Jens Petit's avatar
Jens Petit committed
973 974
    GIVEN("A 2D setting with a single ray")
    {
Nikola Dinev's avatar
Nikola Dinev committed
975 976 977 978
        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
979 980
        volumeDims << volSize, volSize;
        sinoDims << detectorSize, numImgs;
Nikola Dinev's avatar
Nikola Dinev committed
981 982 983 984 985 986
        DataDescriptor volumeDescriptor(volumeDims);
        DataDescriptor sinoDescriptor(sinoDims);
        DataContainer volume(volumeDescriptor);
        DataContainer sino(sinoDescriptor);
        std::vector<Geometry> geom;

Jens Petit's avatar
Jens Petit committed
987 988 989 990
        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
991
            geom.emplace_back(volSize * 20, volSize, -pi_t / 6, volumeDescriptor, sinoDescriptor);
Jens Petit's avatar
Jens Petit committed
992
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
993

Jens Petit's avatar
Jens Petit committed
994 995
            THEN("Ray intersects the correct pixels")
            {
Nikola Dinev's avatar
Nikola Dinev committed
996
                volume = 1;
Jens Petit's avatar
Jens Petit committed
997 998 999
                volume(3, 0) = 0;
                volume(2, 0) = 0;
                volume(2, 1) = 0;
Nikola Dinev's avatar
Nikola Dinev committed
1000

Jens Petit's avatar
Jens Petit committed
1001 1002 1003
                volume(1, 2) = 0;
                volume(1, 3) = 0;
                volume(0, 3) = 0;
Nikola Dinev's avatar
Nikola Dinev committed
1004
                // volume(2,2 also hit due to numerical errors)
Jens Petit's avatar
Jens Petit committed
1005
                volume(2, 2) = 0;
Nikola Dinev's avatar
Nikola Dinev committed
1006 1007 1008

                op.apply(volume, sino);
                DataContainer zero(sinoDescriptor);
1009
                zero = 0;
1010
                CHECK(std::abs(sino[0]) <= Approx(0.0001f).epsilon(epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
1011

Jens Petit's avatar
Jens Petit committed
1012 1013 1014 1015 1016
                AND_THEN("The correct weighting is applied")
                {
                    volume(3, 0) = 1;
                    volume(2, 0) = 2;
                    volume(2, 1) = 3;
Nikola Dinev's avatar
Nikola Dinev committed
1017

Jens Petit's avatar
Jens Petit committed
1018
                    op.apply(volume, sino);
1019
                    REQUIRE(sino[0] == Approx(2 * std::sqrt(3.f) + 2));
Nikola Dinev's avatar
Nikola Dinev committed
1020

Jens Petit's avatar
Jens Petit committed
1021
                    // on the other side of the center
Nikola Dinev's avatar
Nikola Dinev committed
1022
                    volume = 0;
Jens Petit's avatar
Jens Petit committed
1023 1024 1025 1026 1027
                    volume(1, 2) = 3;
                    volume(1, 3) = 2;
                    volume(0, 3) = 1;

                    op.apply(volume, sino);
1028
                    REQUIRE(sino[0] == Approx(2 * std::sqrt(3.f) + 2));
Nikola Dinev's avatar
Nikola Dinev committed
1029 1030 1031

                    sino[0] = 1;

Jens Petit's avatar
Jens Petit committed
1032
                    RealVector_t expected(volSize * volSize);
1033 1034 1035
                    expected << 0, 0, 2 - 2 / std::sqrt(3.f), 4 / std::sqrt(3.f) - 2, 0, 0,
                        2 / std::sqrt(3.f), 0, 0, 2 / std::sqrt(3.f), 0, 0, 4 / std::sqrt(3.f) - 2,
                        2 - 2 / std::sqrt(3.f), 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
1036

Jens Petit's avatar
Jens Petit committed
1037 1038
                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, expected)));
Nikola Dinev's avatar
Nikola Dinev committed
1039 1040 1041 1042
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
1043 1044 1045 1046
        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
1047 1048
            geom.emplace_back(volSize * 20, volSize, -pi_t / 6, volumeDescriptor, sinoDescriptor,
                              0.0, std::sqrt(3.f));
Jens Petit's avatar
Jens Petit committed
1049
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
1050

Jens Petit's avatar
Jens Petit committed
1051 1052
            THEN("Ray intersects the correct pixels")
            {
Nikola Dinev's avatar
Nikola Dinev committed
1053
                volume = 1;
Jens Petit's avatar
Jens Petit committed
1054 1055 1056 1057 1058
                volume(3, 1) = 0;
                volume(3, 2) = 0;
                volume(3, 3) = 0;
                volume(2, 3) = 0;

Nikola Dinev's avatar
Nikola Dinev committed
1059
                op.apply(volume, sino);
Jens Petit's avatar
Jens Petit committed
1060
                DataContainer zero(sinoDescriptor);
1061
                zero = 0;
1062
                REQUIRE(isApprox(sino, zero, epsilon));
Jens Petit's avatar
Jens Petit committed
1063 1064 1065 1066 1067 1068 1069

                AND_THEN("The correct weighting is applied")
                {
                    volume(3, 1) = 4;
                    volume(3, 2) = 3;
                    volume(3, 3) = 2;
                    volume(2, 3) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
1070

Jens Petit's avatar
Jens Petit committed
1071
                    op.apply(volume, sino);
1072
                    REQUIRE(sino[0] == Approx(14 - 4 * std::sqrt(3.f)));
Nikola Dinev's avatar
Nikola Dinev committed
1073 1074 1075

                    sino[0] = 1;

Jens Petit's avatar
Jens Petit committed
1076
                    RealVector_t expected(volSize * volSize);
1077 1078
                    expected << 0, 0, 0, 0, 0, 0, 0, 4 - 2 * std::sqrt(3.f), 0, 0, 0,
                        2 / std::sqrt(3.f), 0, 0, 2 - 2 / std::sqrt(3.f), 4 / std::sqrt(3.f) - 2;
Nikola Dinev's avatar
Nikola Dinev committed
1079

Jens Petit's avatar
Jens Petit committed
1080
                    op.applyAdjoint(sino, volume);
1081
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, expected), epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
1082 1083 1084 1085
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
1086 1087 1088 1089
        WHEN("Projecting under an angle of 30 degrees and ray exits through the left border")
        {
            // In this case the ray enters through a border along the main ray direction, but exits
            // through a border not along the main direction
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
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
1093

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

Nikola Dinev's avatar
Nikola Dinev committed
1102
                op.apply(volume, sino);
Jens Petit's avatar
Jens Petit committed
1103
                DataContainer zero(sinoDescriptor);
1104
                zero = 0;
1105
                REQUIRE(isApprox(sino, zero, epsilon));
Jens Petit's avatar
Jens Petit committed
1106 1107 1108 1109 1110 1111 1112

                AND_THEN("The correct weighting is applied")
                {
                    volume(1, 0) = 1;
                    volume(0, 0) = 2;
                    volume(0, 1) = 3;
                    volume(0, 2) = 4;
Nikola Dinev's avatar
Nikola Dinev committed
1113

Jens Petit's avatar
Jens Petit committed
1114
                    op.apply(volume, sino);
1115
                    REQUIRE(sino[0] == Approx(14 - 4 * std::sqrt(3.f)));
Nikola Dinev's avatar
Nikola Dinev committed
1116 1117 1118

                    sino[0] = 1;

Jens Petit's avatar
Jens Petit committed
1119
                    RealVector_t expected(volSize * volSize);
1120 1121
                    expected << 4 / std::sqrt(3.f) - 2, 2 - 2 / std::sqrt(3.f), 0, 0,
                        2 / std::sqrt(3.f), 0, 0, 0, 4 - 2 * std::sqrt(3.f), 0, 0, 0, 0, 0, 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
1122

Jens Petit's avatar
Jens Petit committed
1123 1124
                    op.applyAdjoint(sino, volume);
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, expected)));
Nikola Dinev's avatar
Nikola Dinev committed
1125 1126 1127 1128
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
1129 1130
        WHEN("Projecting under an angle of 30 degrees and ray only intersects a single pixel")
        {
1131 1132
            geom.emplace_back(volSize * 20, volSize, -pi_t / 6, volumeDescriptor, sinoDescriptor,
                              0.0, -2 - std::sqrt(3.f) / 2);
Jens Petit's avatar
Jens Petit committed
1133
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
1134

Jens Petit's avatar
Jens Petit committed
1135 1136
            THEN("Ray intersects the correct pixels")
            {
Nikola Dinev's avatar
Nikola Dinev committed
1137
                volume = 1;
Jens Petit's avatar
Jens Petit committed
1138 1139
                volume(0, 0) = 0;

Nikola Dinev's avatar
Nikola Dinev committed
1140
                op.apply(volume, sino);
Jens Petit's avatar
Jens Petit committed
1141
                DataContainer zero(sinoDescriptor);
1142
                zero = 0;
1143
                REQUIRE(isApprox(sino, zero, epsilon));
Jens Petit's avatar
Jens Petit committed
1144 1145 1146 1147

                AND_THEN("The correct weighting is applied")
                {
                    volume(0, 0) = 1;
Nikola Dinev's avatar
Nikola Dinev committed
1148

Jens Petit's avatar
Jens Petit committed
1149
                    op.apply(volume, sino);
1150
                    REQUIRE(sino[0] == Approx(1 / std::sqrt(3.f)));
Nikola Dinev's avatar
Nikola Dinev committed
1151 1152 1153

                    sino[0] = 1;

Jens Petit's avatar
Jens Petit committed
1154
                    RealVector_t expected(volSize * volSize);
1155
                    expected << 1 / std::sqrt(3.f), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
Nikola Dinev's avatar
Nikola Dinev committed
1156

Jens Petit's avatar
Jens Petit committed
1157
                    op.applyAdjoint(sino, volume);
1158
                    REQUIRE(isApprox(volume, DataContainer(volumeDescriptor, expected), epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
1159 1160 1161 1162
                }
            }
        }

Jens Petit's avatar
Jens Petit committed
1163 1164 1165 1166
        WHEN("Projecting under an angle of 120 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
1167 1168
            geom.emplace_back(volSize * 20, volSize, -2 * pi_t / 3, volumeDescriptor,
                              sinoDescriptor);
Jens Petit's avatar
Jens Petit committed
1169
            SiddonsMethod op(volumeDescriptor, sinoDescriptor, geom);
Nikola Dinev's avatar
Nikola Dinev committed
1170

Jens Petit's avatar
Jens Petit committed
1171 1172
            THEN("Ray intersects the correct pixels")
            {
Nikola Dinev's avatar
Nikola Dinev committed
1173
                volume = 1;
Jens Petit's avatar
Jens Petit committed
1174 1175 1176 1177 1178 1179
                volume(0, 0) = 0;
                volume(0, 1) = 0;
                volume(1, 1) = 0;
                volume(2, 2) = 0;
                volume(3, 2) = 0;
                volume(3, 3) = 0;
Nikola Dinev's avatar
Nikola Dinev committed
1180
                // volume(1,2) hit due to numerical error
Jens Petit's avatar
Jens Petit committed
1181
                volume(1, 2) = 0;
Nikola Dinev's avatar
Nikola Dinev committed
1182 1183

                op.apply(volume, sino);
Jens Petit's avatar
Jens Petit committed
1184
                DataContainer zero(sinoDescriptor);
1185
                zero = 0;
1186
                CHECK(std::abs(sino[0]) <= Approx(0.0001f).epsilon(epsilon));
Nikola Dinev's avatar
Nikola Dinev committed
1187

Jens Petit's avatar
Jens Petit committed
1188 1189 1190 1191 1192
                AND_THEN("The correct weighting is applied")
                {
                    volume(0, 0) = 1;
                    volume(0, 1) = 2;
                    volume(1, 1) = 3;
Nikola Dinev's avatar
Nikola Dinev committed
1193

Jens Petit's avatar
Jens Petit committed
1194
                    op.apply(volume, sino);
1195
                    REQUIRE(sino[0] == Approx(2 * std::sqrt(3.f) + 2));
Nikola Dinev's avatar
Nikola Dinev committed
1196

Jens Petit's avatar
Jens Petit committed
1197
                    // on the other side of the center