Add Rectangle and Circle phantom generator

elsa/generators/PhantomGenerator.cpp

@@ -2,12 +2,50 @@
 #include "EllipseGenerator.h"
 #include "Logger.h"
 #include "VolumeDescriptor.h"
+#include "CartesianIndices.h"
 
 #include <cmath>
 #include <stdexcept>
 
 namespace elsa
 {
+    template <typename data_t>
+    DataContainer<data_t> PhantomGenerator<data_t>::createCirclePhantom(IndexVector_t volumesize,
+                                                                        data_t radius)
+    {
+        VolumeDescriptor dd(volumesize);
+        DataContainer<data_t> dc(dd);
+        dc = 0;
+
+        const Vector_t<data_t> sizef = volumesize.template cast<data_t>();
+        const auto center = (sizef.array() / 2).matrix();
+
+        for (auto pos : CartesianIndices(volumesize)) {
+            const Vector_t<data_t> p = pos.template cast<data_t>();
+            if ((p - center).norm() <= radius) {
+                dc(pos) = 1;
+            }
+        }
+
+        return dc;
+    }
+
+    template <typename data_t>
+    DataContainer<data_t> PhantomGenerator<data_t>::createRectanglePhantom(IndexVector_t volumesize,
+                                                                           IndexVector_t lower,
+                                                                           IndexVector_t upper)
+    {
+        VolumeDescriptor dd(volumesize);
+        DataContainer<data_t> dc(dd);
+        dc = 0;
+
+        for (auto pos : CartesianIndices(lower, upper)) {
+            dc(pos) = 1;
+        }
+
+        return dc;
+    }
+
     template <typename data_t>
     DataContainer<data_t> PhantomGenerator<data_t>::createModifiedSheppLogan(IndexVector_t sizes)
     {

elsa/generators/PhantomGenerator.h

@@ -21,18 +21,39 @@ namespace elsa
         /**
          * @brief Create a modified Shepp-Logan phantom in 2d or 3d (with enhanced contrast).
          *
-         * @param[in] sizes a 2d/3d vector indicating the requested size (has to be square!)
-         *
-         * @returns DataContainer of specified size containing the phantom.
-         *
          * The phantom specifications are adapted from Matlab (which in turn references A.K. Jain,
          * "Fundamentals of Digital Image Processing", p. 439, and P.A. Toft, "The Radon Transform,
          * Theory and Implementation", p. 199).
          *
          * Warning: the 3D version is currently very inefficient to compute (cubic algorithm).
+         *
+         * @param[in] sizes a 2d/3d vector indicating the requested size (has to be square!)
+         *
+         * @returns DataContainer of specified size containing the phantom.
          */
         static DataContainer<data_t> createModifiedSheppLogan(IndexVector_t sizes);
 
+        /**
+         * @brief Create a phantom with a simple n-dimensional rectangle  going from lower to upper.
+         * It is assumed that lower < upper.
+         *
+         * @param[in] volumesize size of the volume
+         * @param[in] lower the lower corner of the rectangle
+         * @param[in] upper the upper corner of the rectangle
+         */
+        static DataContainer<data_t> createRectanglePhantom(IndexVector_t volumesize,
+                                                            IndexVector_t lower,
+                                                            IndexVector_t upper);
+
+        /**
+         * @brief Create a phantom with a simple n-dimensional sphere centered in the middle with
+         * given raidus
+         *
+         * @param[in] volumesize size of the volume
+         * @param[in] radius the radius of the circle
+         */
+        static DataContainer<data_t> createCirclePhantom(IndexVector_t volumesize, data_t radius);
+
     private:
         /// scale sizes from [0,1] to the (square) phantom size, producing indices (integers)
         static index_t scale(const DataDescriptor& dd, data_t value);

elsa/generators/tests/test_PhantomGenerator.cpp

@@ -10,6 +10,7 @@
 #include "PhantomGenerator.h"
 #include "VolumeDescriptor.h"
 #include "testHelpers.h"
+#include "CartesianIndices.h"
 #include <array>
 #include <iostream>
 
@@ -18,6 +19,129 @@ using namespace doctest;
 
 RealVector_t get2dModifiedSheppLogan45x45();
 
+TEST_CASE_TEMPLATE("PhantomGenerator: Drawing a simple 2d rectangle", data_t, float, double)
+{
+    const IndexVector_t size({{16, 16}});
+
+    WHEN("Drawing a rectangle going from the lower left corner to the upper right corner")
+    {
+        const IndexVector_t lower({{0, 0}});
+        const IndexVector_t upper({{16, 16}});
+        const auto dc = PhantomGenerator<data_t>::createRectanglePhantom(size, lower, upper);
+
+        THEN("Everything is set to 1")
+        {
+            for (int i = 0; i < dc.getSize(); ++i) {
+                CHECK_EQ(dc[i], 1);
+            }
+        }
+    }
+
+    WHEN("Drawing a rectangle")
+    {
+        const IndexVector_t lower({{4, 4}});
+        const IndexVector_t upper({{12, 12}});
+        const auto dc = PhantomGenerator<data_t>::createRectanglePhantom(size, lower, upper);
+
+        THEN("The pixels inside the rectangle are set to 1")
+        {
+            for (auto pos : CartesianIndices(lower, upper)) {
+
+                CHECK_EQ(dc(pos), 1);
+            }
+        }
+    }
+}
+
+TEST_CASE_TEMPLATE("PhantomGenerator: Drawing a simple 3d rectangle", data_t, float, double)
+{
+    const IndexVector_t size({{16, 16, 16}});
+
+    WHEN("Drawing a rectangle going from the lower left corner to the upper right corner")
+    {
+        const IndexVector_t lower({{0, 0, 0}});
+        const IndexVector_t upper({{16, 16, 16}});
+        const auto dc = PhantomGenerator<data_t>::createRectanglePhantom(size, lower, upper);
+
+        THEN("Everything is set to 1")
+        {
+            for (int i = 0; i < dc.getSize(); ++i) {
+                CHECK_EQ(dc[i], 1);
+            }
+        }
+    }
+
+    WHEN("Drawing a rectangle")
+    {
+        const IndexVector_t lower({{4, 4, 4}});
+        const IndexVector_t upper({{12, 12, 12}});
+        const auto dc = PhantomGenerator<data_t>::createRectanglePhantom(size, lower, upper);
+
+        THEN("The pixels inside the rectangle are set to 1")
+        {
+            for (auto pos : CartesianIndices(lower, upper)) {
+
+                CHECK_EQ(dc(pos), 1);
+            }
+        }
+    }
+}
+
+#undef SUBCASE
+#define SUBCASE(...) DOCTEST_SUBCASE(std::string(__VA_ARGS__).c_str())
+
+TEST_CASE_TEMPLATE("PhantomGenerator: Drawing a simple 2d circle", data_t, float, double)
+{
+    const IndexVector_t size({{16, 16}});
+
+    for (int i = 1; i < 9; ++i) {
+        SUBCASE("    When: Drawing a circle of radius " + std::to_string(i))
+        {
+            const data_t radius = i;
+            const auto dc = PhantomGenerator<data_t>::createCirclePhantom(size, radius);
+
+            THEN("Everything the correct pixels are set to 1")
+            {
+                const auto center = (size.template cast<data_t>().array() / 2).matrix();
+                for (auto pos : CartesianIndices(size)) {
+                    auto dist_to_center = (pos.template cast<data_t>() - center).norm();
+                    if (dist_to_center <= radius) {
+                        CHECK_EQ(dc(pos), 1);
+                    } else {
+                        CHECK_EQ(dc(pos), 0);
+                    }
+                }
+            }
+        }
+    }
+}
+
+TEST_CASE_TEMPLATE("PhantomGenerator: Drawing a simple 3d circle", data_t, float, double)
+{
+    const IndexVector_t size({{16, 16, 16}});
+
+    for (int i = 1; i < 9; ++i) {
+        SUBCASE("    When: Drawing a circle of radius " + std::to_string(i))
+        {
+            const data_t radius = i;
+            const auto dc = PhantomGenerator<data_t>::createCirclePhantom(size, radius);
+
+            THEN("Everything the correct pixels are set to 1")
+            {
+                const auto center = (size.template cast<data_t>().array() / 2).matrix();
+                for (auto pos : CartesianIndices(size)) {
+                    auto dist_to_center = (pos.template cast<data_t>() - center).norm();
+                    if (dist_to_center <= radius) {
+                        CHECK_EQ(dc(pos), 1);
+                    } else {
+                        CHECK_EQ(dc(pos), 0);
+                    }
+                }
+            }
+        }
+    }
+}
+
 TEST_CASE_TEMPLATE("PhantomGenerator: Drawing a 2d Shepp-Logan phantom", data_t, float, double)
 {
 
@@ -59,7 +183,7 @@ TEST_CASE_TEMPLATE("PhantomGenerator: Drawing a 2d Shepp-Logan phantom", data_t,
                 auto ref = DataContainer(VolumeDescriptor(size), expected);
 
                 INFO("Computed phantom: ", dc);
-                INFO("Refernce phantom: ", ref);
+                INFO("Reference phantom: ", ref);
 
                 for (int i = 0; i < ref.getSize(); ++i) {
                     INFO("Error at position: ", i);