Merge branch 'master' into shearlet_transform

# Conflicts:
#	elsa/core/CMakeLists.txt

elsa/core/CMakeLists.txt

@@ -5,6 +5,7 @@ set(MODULE_HEADERS
     Cloneable.h
     Utilities/Badge.hpp
     Utilities/DataContainerFormatter.hpp
+    Utilities/FormatConfig.h
     Utilities/Math.hpp
     Utilities/Statistics.hpp
     Utilities/TypeCasts.hpp

elsa/core/DataContainer.cpp

 #include "DataContainer.h"
 #include "DataContainerFormatter.hpp"
+#include "FormatConfig.h"
 #include "DataHandlerCPU.h"
 #include "DataHandlerMapCPU.h"
 #include "BlockDescriptor.h"
@@ -509,9 +510,9 @@ namespace elsa
     }
 
     template <typename data_t>
-    void DataContainer<data_t>::format(std::ostream& os) const
+    void DataContainer<data_t>::format(std::ostream& os, format_config cfg) const
     {
-        DataContainerFormatter<data_t> fmt;
+        DataContainerFormatter<data_t> fmt{cfg};
         fmt.format(os, *this);
     }
 
@@ -604,6 +605,58 @@ namespace elsa
         return concatenated;
     }
 
+    template <typename data_t>
+    DataContainer<data_t> fftShift2D(DataContainer<data_t> dc)
+    {
+        assert(dc.getDataDescriptor().getNumberOfDimensions() == 2
+               && "DataContainer::fftShift2D: currently only supporting 2D signals");
+
+        const DataDescriptor& dataDescriptor = dc.getDataDescriptor();
+        IndexVector_t numOfCoeffsPerDim = dataDescriptor.getNumberOfCoefficientsPerDimension();
+        index_t m = numOfCoeffsPerDim[0];
+        index_t n = numOfCoeffsPerDim[1];
+
+        index_t firstShift = m / 2;
+        index_t secondShift = n / 2;
+
+        DataContainer<data_t> copyDC(dataDescriptor);
+
+        for (index_t i = 0; i < m; ++i) {
+            for (index_t j = 0; j < n; ++j) {
+                copyDC((i + firstShift) % m, (j + secondShift) % n) = dc(i, j);
+            }
+        }
+
+        return copyDC;
+    }
+
+    template <typename data_t>
+    DataContainer<data_t> ifftShift2D(DataContainer<data_t> dc)
+    {
+        assert(dc.getDataDescriptor().getNumberOfDimensions() == 2
+               && "DataContainer::ifftShift2D: currently only supporting 2D signals");
+
+        const DataDescriptor& dataDescriptor = dc.getDataDescriptor();
+        IndexVector_t numOfCoeffsPerDim = dataDescriptor.getNumberOfCoefficientsPerDimension();
+        index_t m = numOfCoeffsPerDim[0];
+        index_t n = numOfCoeffsPerDim[1];
+
+        index_t firstShift = -m / 2;
+        index_t secondShift = -n / 2;
+
+        DataContainer<data_t> copyDC(dataDescriptor);
+
+        for (index_t i = 0; i < m; ++i) {
+            for (index_t j = 0; j < n; ++j) {
+                index_t leftIndex = (((i + firstShift) % m) + m) % m;
+                index_t rightIndex = (((j + secondShift) % n) + n) % n;
+                copyDC(leftIndex, rightIndex) = dc(i, j);
+            }
+        }
+
+        return copyDC;
+    }
+
     // ------------------------------------------
     // explicit template instantiation
     template class DataContainer<float>;
@@ -623,4 +676,18 @@ namespace elsa
         concatenate<std::complex<double>>(const DataContainer<std::complex<double>>&,
                                           const DataContainer<std::complex<double>>&);
 
+    template DataContainer<float> fftShift2D<float>(DataContainer<float>);
+    template DataContainer<std::complex<float>>
+        fftShift2D<std::complex<float>>(DataContainer<std::complex<float>>);
+    template DataContainer<double> fftShift2D<double>(DataContainer<double>);
+    template DataContainer<std::complex<double>>
+        fftShift2D<std::complex<double>>(DataContainer<std::complex<double>>);
+
+    template DataContainer<float> ifftShift2D<float>(DataContainer<float>);
+    template DataContainer<std::complex<float>>
+        ifftShift2D<std::complex<float>>(DataContainer<std::complex<float>>);
+    template DataContainer<double> ifftShift2D<double>(DataContainer<double>);
+    template DataContainer<std::complex<double>>
+        ifftShift2D<std::complex<double>>(DataContainer<std::complex<double>>);
+
 } // namespace elsa

elsa/core/DataContainer.h

@@ -8,6 +8,7 @@
 #include "DataContainerIterator.h"
 #include "Error.h"
 #include "Expression.h"
+#include "FormatConfig.h"
 #include "TypeCasts.hpp"
 
 #ifdef ELSA_CUDA_VECTOR
@@ -20,7 +21,6 @@
 
 namespace elsa
 {
-
     /**
      * @brief class representing and storing a linearized n-dimensional signal
      *
@@ -444,15 +444,15 @@ namespace elsa
         friend constexpr auto evaluateOrReturn(Operand const& operand);
 
         /// write a pretty-formatted string representation to stream
-        void format(std::ostream& os) const;
+        void format(std::ostream& os, format_config cfg = {}) const;
 
         /**
          * @brief Factory function which returns GPU based DataContainers
          *
          * @return the GPU based DataContainer
          *
-         * Note that if this function is called on a container which is already GPU based, it will
-         * throw an exception.
+         * Note that if this function is called on a container which is already GPU based, it
+         * will throw an exception.
          */
         DataContainer loadToGPU();
 
@@ -516,6 +516,18 @@ namespace elsa
     DataContainer<data_t> concatenate(const DataContainer<data_t>& dc1,
                                       const DataContainer<data_t>& dc2);
 
+    /// Perform the FFT shift operation to the provided signal. Refer to
+    /// https://numpy.org/doc/stable/reference/generated/numpy.fft.fftshift.html for further
+    /// details.
+    template <typename data_t>
+    DataContainer<data_t> fftShift2D(DataContainer<data_t> dc);
+
+    /// Perform the IFFT shift operation to the provided signal. Refer to
+    /// https://numpy.org/doc/stable/reference/generated/numpy.fft.ifftshift.html for further
+    /// details.
+    template <typename data_t>
+    DataContainer<data_t> ifftShift2D(DataContainer<data_t> dc);
+
     /// User-defined template argument deduction guide for the expression based constructor
     template <typename Source>
     DataContainer(Source const& source) -> DataContainer<typename Source::data_t>;

elsa/core/Utilities/DataContainerFormatter.hpp

 #pragma once
 
-#include <elsaDefines.h>
+#include "elsaDefines.h"
+#include "DataContainer.h"
+#include "FormatConfig.h"
+
 #include <iomanip>
 #include <iostream>
 #include <iterator>
@@ -25,29 +28,7 @@ namespace elsa
     template <typename data_t>
     class DataContainerFormatter
     {
-
     public:
-        /**
-         * Formatting output configuration.
-         */
-        struct format_config {
-            /// for too many elements, abbreviate the output
-            bool summary_enabled = true;
-
-            /// number of summary items to display - this also triggers the summary enabling
-            /// if theres more than 2 * summary_items elements.
-            index_t summary_items = 6;
-
-            /// what is inserted between the summary items
-            std::string summary_elem = "...";
-
-            // what is inserted vertically between summary items
-            std::string summary_elem_vertical = "\u22EE";
-
-            // what's inserted between elements and newlines
-            std::string separator = ", ";
-        };
-
         /**
          * Create a formatter with default config.
          */
@@ -73,7 +54,7 @@ namespace elsa
         }
 
         /// stream a pretty format to is using the current configuration.
-        void format(std::ostream& os, const elsa::DataContainer<data_t>& dc)
+        void format(std::ostream& os, const DataContainer<data_t>& dc)
         {
             using namespace std::string_literals;
 
@@ -224,7 +205,7 @@ namespace elsa
          */
         template <typename T>
         std::function<std::ostream&(std::ostream& os, const T& elem)>
-            get_element_formatter(const elsa::DataContainer<T>& dc)
+            get_element_formatter(const DataContainer<T>& dc)
         {
             if constexpr (elsa::isComplex<T>) {
                 // format both components independently
@@ -271,9 +252,14 @@ namespace elsa
                     teststream.str("");
                     teststream.clear();
 
-                    auto&& elem = dc[idx];
+                    auto&& elem = config.suppress_close_to_zero
+                                          && std::abs(dc[idx]) < config.suppression_epsilon
+                                      ? static_cast<data_t>(0)
+                                      : dc[idx];
+
                     teststream << elem;
                     auto len = static_cast<int>(teststream.tellp());
+
                     if (len > maxlen) {
                         maxlen = len;
                     }
@@ -281,11 +267,14 @@ namespace elsa
 
                 auto streamflags = teststream.flags();
 
-                return [maxlen, streamflags](std::ostream & os, const T& elem) -> auto&
+                return [
+                    maxlen, streamflags, do_suppress = config.suppress_close_to_zero,
+                    eps = config.suppression_epsilon
+                ](std::ostream & os, const T& elem) -> auto&
                 {
                     os.flags(streamflags);
                     os << std::setw(maxlen);
-                    os << elem;
+                    os << (do_suppress && std::abs(elem) < eps ? static_cast<data_t>(0) : elem);
                     return os;
                 };
             } else {

elsa/core/Utilities/FormatConfig.h

+#pragma once
+
+#include "elsaDefines.h"
+
+#include <string>
+
+namespace elsa
+{
+    /**
+     * Formatting output configuration.
+     */
+    struct format_config {
+        /// for too many elements, abbreviate the output
+        bool summary_enabled = true;
+
+        /// number of summary items to display - this also triggers the summary enabling
+        /// if theres more than 2 * summary_items elements.
+        index_t summary_items = 6;
+
+        /// what is inserted between the summary items
+        std::string summary_elem = "...";
+
+        /// what is inserted vertically between summary items
+        std::string summary_elem_vertical = "\u22EE";
+
+        /// what's inserted between elements and newlines
+        std::string separator = ", ";
+
+        /// if a value is smaller than some epsilon (`suppression_epsilon`), just print 0
+        /// instead
+        bool suppress_close_to_zero = false;
+
+        /// epsilon value for suppressing small numbers
+        real_t suppression_epsilon = static_cast<real_t>(0.0000001);
+    };
+} // namespace elsa

elsa/core/tests/CMakeLists.txt

@@ -26,3 +26,4 @@ ELSA_DOCTEST(ExpressionTemplates)
 ELSA_DOCTEST(DataHandlers)
 ELSA_DOCTEST(DataHandlerMap)
 ELSA_DOCTEST(DataContainer)
+ELSA_DOCTEST(DataContainerFormatter)

elsa/core/tests/test_DataContainer.cpp

@@ -1303,6 +1303,205 @@ TEST_CASE_TEMPLATE("DataContainer: Slice a DataContainer", data_t, float, double
     }
 }
 
+TEST_CASE_TEMPLATE("DataContainer: FFT shift and IFFT shift a DataContainer", data_t, float, double,
+                   std::complex<float>, std::complex<double>)
+{
+    GIVEN("a one-element 2D data container")
+    {
+        DataContainer<data_t> dc(VolumeDescriptor{{1, 1}});
+        dc[0] = 8;
+        WHEN("running the FFT shift operation to the container")
+        {
+            DataContainer<data_t> fftShiftedDC = fftShift2D(dc);
+            THEN("the data descriptors match")
+            {
+                REQUIRE_EQ(dc.getDataDescriptor(), fftShiftedDC.getDataDescriptor());
+            }
+            THEN("the data containers match") { REQUIRE_UNARY(fftShiftedDC == dc); }
+        }
+
+        WHEN("running the IFFT shift operation to the container")
+        {
+            DataContainer<data_t> ifftShiftedDC = ifftShift2D(dc);
+            THEN("the data descriptors match")
+            {
+                REQUIRE_EQ(dc.getDataDescriptor(), ifftShiftedDC.getDataDescriptor());
+            }
+            THEN("the data containers match") { REQUIRE_UNARY(ifftShiftedDC == dc); }
+        }
+    }
+
+    GIVEN("a 3x3 2D data container")
+    {
+        DataContainer<data_t> dc(VolumeDescriptor{{3, 3}});
+        dc(0, 0) = 0;
+        dc(0, 1) = 1;
+        dc(0, 2) = 2;
+        dc(1, 0) = 3;
+        dc(1, 1) = 4;
+        dc(1, 2) = -4;
+        dc(2, 0) = -3;
+        dc(2, 1) = -2;
+        dc(2, 2) = -1;
+
+        DataContainer<data_t> expectedFFTShiftDC(VolumeDescriptor{{3, 3}});
+        expectedFFTShiftDC(0, 0) = -1;
+        expectedFFTShiftDC(0, 1) = -3;
+        expectedFFTShiftDC(0, 2) = -2;
+        expectedFFTShiftDC(1, 0) = 2;
+        expectedFFTShiftDC(1, 1) = 0;
+        expectedFFTShiftDC(1, 2) = 1;
+        expectedFFTShiftDC(2, 0) = -4;
+        expectedFFTShiftDC(2, 1) = 3;
+        expectedFFTShiftDC(2, 2) = 4;
+
+        WHEN("running the FFT shift operation to the container")
+        {
+            DataContainer<data_t> fftShiftedDC = fftShift2D(dc);
+            THEN("the data descriptors match")
+            {
+                REQUIRE_EQ(fftShiftedDC.getDataDescriptor(),
+                           expectedFFTShiftDC.getDataDescriptor());
+            }
+            THEN("the data containers match") { REQUIRE_UNARY(fftShiftedDC == expectedFFTShiftDC); }
+        }
+
+        DataContainer<data_t> expectedIFFTShiftDC(VolumeDescriptor{{3, 3}});
+        expectedIFFTShiftDC(0, 0) = 4;
+        expectedIFFTShiftDC(0, 1) = -4;
+        expectedIFFTShiftDC(0, 2) = 3;
+        expectedIFFTShiftDC(1, 0) = -2;
+        expectedIFFTShiftDC(1, 1) = -1;
+        expectedIFFTShiftDC(1, 2) = -3;
+        expectedIFFTShiftDC(2, 0) = 1;
+        expectedIFFTShiftDC(2, 1) = 2;
+        expectedIFFTShiftDC(2, 2) = 0;
+
+        WHEN("running the IFFT shift operation to the container")
+        {
+            DataContainer<data_t> ifftShiftedDC = ifftShift2D(dc);
+            THEN("the data descriptors match")
+            {
+                REQUIRE_EQ(ifftShiftedDC.getDataDescriptor(),
+                           expectedIFFTShiftDC.getDataDescriptor());
+            }
+            THEN("the data containers match")
+            {
+                REQUIRE_UNARY(ifftShiftedDC == expectedIFFTShiftDC);
+            }
+        }
+    }
+
+    GIVEN("a 5x5 2D data container")
+    {
+        DataContainer<data_t> dc(VolumeDescriptor{{5, 5}});
+        dc(0, 0) = 28;
+        dc(0, 1) = 1;
+        dc(0, 2) = 5;
+        dc(0, 3) = -18;
+        dc(0, 4) = 8;
+        dc(1, 0) = 5;
+        dc(1, 1) = 6;
+        dc(1, 2) = 50;
+        dc(1, 3) = -8;
+        dc(1, 4) = 9;
+        dc(2, 0) = 8;
+        dc(2, 1) = 9;
+        dc(2, 2) = 10;
+        dc(2, 3) = 11;
+        dc(2, 4) = 12;
+        dc(3, 0) = -12;
+        dc(3, 1) = -41;
+        dc(3, 2) = -10;
+        dc(3, 3) = -9;
+        dc(3, 4) = -8;
+        dc(4, 0) = -70;
+        dc(4, 1) = -6;
+        dc(4, 2) = 22;
+        dc(4, 3) = -10;
+        dc(4, 4) = -3;
+
+        DataContainer<data_t> expectedFFTShiftDC(VolumeDescriptor{{5, 5}});
+        expectedFFTShiftDC(0, 0) = -9;
+        expectedFFTShiftDC(0, 1) = -8;
+        expectedFFTShiftDC(0, 2) = -12;
+        expectedFFTShiftDC(0, 3) = -41;
+        expectedFFTShiftDC(0, 4) = -10;
+        expectedFFTShiftDC(1, 0) = -10;
+        expectedFFTShiftDC(1, 1) = -3;
+        expectedFFTShiftDC(1, 2) = -70;
+        expectedFFTShiftDC(1, 3) = -6;
+        expectedFFTShiftDC(1, 4) = 22;
+        expectedFFTShiftDC(2, 0) = -18;
+        expectedFFTShiftDC(2, 1) = 8;
+        expectedFFTShiftDC(2, 2) = 28;
+        expectedFFTShiftDC(2, 3) = 1;
+        expectedFFTShiftDC(2, 4) = 5;
+        expectedFFTShiftDC(3, 0) = -8;
+        expectedFFTShiftDC(3, 1) = 9;
+        expectedFFTShiftDC(3, 2) = 5;
+        expectedFFTShiftDC(3, 3) = 6;
+        expectedFFTShiftDC(3, 4) = 50;
+        expectedFFTShiftDC(4, 0) = 11;
+        expectedFFTShiftDC(4, 1) = 12;
+        expectedFFTShiftDC(4, 2) = 8;
+        expectedFFTShiftDC(4, 3) = 9;
+        expectedFFTShiftDC(4, 4) = 10;
+
+        WHEN("running the FFT shift operation to the container")
+        {
+            DataContainer<data_t> fftShiftedDC = fftShift2D(dc);
+            THEN("the data descriptors match")
+            {
+                REQUIRE_EQ(fftShiftedDC.getDataDescriptor(),
+                           expectedFFTShiftDC.getDataDescriptor());
+            }
+            THEN("the data containers match") { REQUIRE_UNARY(fftShiftedDC == expectedFFTShiftDC); }
+        }
+
+        DataContainer<data_t> expectedIFFTShiftDC(VolumeDescriptor{{5, 5}});
+        expectedIFFTShiftDC(0, 0) = 10;
+        expectedIFFTShiftDC(0, 1) = 11;
+        expectedIFFTShiftDC(0, 2) = 12;
+        expectedIFFTShiftDC(0, 3) = 8;
+        expectedIFFTShiftDC(0, 4) = 9;
+        expectedIFFTShiftDC(1, 0) = -10;
+        expectedIFFTShiftDC(1, 1) = -9;
+        expectedIFFTShiftDC(1, 2) = -8;
+        expectedIFFTShiftDC(1, 3) = -12;
+        expectedIFFTShiftDC(1, 4) = -41;
+        expectedIFFTShiftDC(2, 0) = 22;
+        expectedIFFTShiftDC(2, 1) = -10;
+        expectedIFFTShiftDC(2, 2) = -3;
+        expectedIFFTShiftDC(2, 3) = -70;
+        expectedIFFTShiftDC(2, 4) = -6;
+        expectedIFFTShiftDC(3, 0) = 5;
+        expectedIFFTShiftDC(3, 1) = -18;
+        expectedIFFTShiftDC(3, 2) = 8;
+        expectedIFFTShiftDC(3, 3) = 28;
+        expectedIFFTShiftDC(3, 4) = 1;
+        expectedIFFTShiftDC(4, 0) = 50;
+        expectedIFFTShiftDC(4, 1) = -8;
+        expectedIFFTShiftDC(4, 2) = 9;
+        expectedIFFTShiftDC(4, 3) = 5;
+        expectedIFFTShiftDC(4, 4) = 6;
+
+        WHEN("running the IFFT shift operation to the container")
+        {
+            DataContainer<data_t> ifftShiftedDC = ifftShift2D(dc);
+            THEN("the data descriptors match")
+            {
+                REQUIRE_EQ(ifftShiftedDC.getDataDescriptor(),
+                           expectedIFFTShiftDC.getDataDescriptor());
+            }
+            THEN("the data containers match")
+            {
+                REQUIRE_UNARY(ifftShiftedDC == expectedIFFTShiftDC);
+            }
+        }
+    }
+}
+
 // "instantiate" the test templates for CPU types
 TEST_CASE_TEMPLATE_APPLY(datacontainer_construction, CPUTypeTuple);
 TEST_CASE_TEMPLATE_APPLY(datacontainer_reduction, CPUTypeTuple);

elsa/core/tests/test_DataContainerFormatter.cpp

+#include "doctest/doctest.h"
+#include "DataContainerFormatter.hpp"
+#include "VolumeDescriptor.h"
+
+using namespace elsa;
+using namespace doctest;
+
+TEST_SUITE_BEGIN("core");
+
+TEST_CASE("DataContainerFormatter: default config")
+{
+    using data_t = float;
+
+    std::stringstream buffer;
+    DataContainerFormatter<data_t> formatter;
+
+    GIVEN("A 1D DataContainer with arbitrary values")
+    {
+        IndexVector_t numCoeff(1);
+        numCoeff << 10;
+
+        DataContainer<data_t> dc(VolumeDescriptor{numCoeff});
+        dc = 1;
+
+        THEN("Formatting writes correctly to the stream")
+        {
+            formatter.format(buffer, dc);
+
+            auto resultString = buffer.str();
+            REQUIRE_EQ(resultString,
+                       "DataContainer<dims=1, shape=(10)>\n[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]");
+        }
+    }
+
+    GIVEN("A larger 1D DataContainer with arbitrary values")
+    {
+        IndexVector_t numCoeff(1);
+        numCoeff << 15;
+
+        DataContainer<data_t> dc(VolumeDescriptor{numCoeff});
+        dc = 1;
+
+        THEN("Formatting writes abbreviated to the stream")
+        {
+            formatter.format(buffer, dc);
+
+            auto resultString = buffer.str();
+            REQUIRE_EQ(
+                resultString,
+                "DataContainer<dims=1, shape=(15)>\n[1, 1, 1, 1, 1, 1, ..., 1, 1, 1, 1, 1, 1]");
+        }
+    }
+
+    GIVEN("A 2D DataContainer with arbitrary values")
+    {
+        IndexVector_t numCoeff(2);
+        numCoeff << 5, 3;
+
+        DataContainer<data_t> dc(VolumeDescriptor{numCoeff});
+        dc = 1;
+
+        THEN("Formatting writes correctly to the stream")
+        {
+            formatter.format(buffer, dc);
+
+            auto resultString = buffer.str();
+            REQUIRE_EQ(resultString,
+                       R"(DataContainer<dims=2, shape=(5 3)>
+[[1, 1, 1],
+ [1, 1, 1],
+ [1, 1, 1],
+ [1, 1, 1],
+ [1, 1, 1]])");
+        }
+    }
+}
+
+TEST_SUITE_END();

elsa/generators/tests/CMakeLists.txt

@@ -19,3 +19,4 @@ ELSA_DOCTEST(PhantomGenerator)
 ELSA_DOCTEST(CircleTrajectoryGenerator)
 ELSA_DOCTEST(LimitedAngleTrajectoryGenerator)
 ELSA_DOCTEST(SphereTrajectoryGenerator)
+ELSA_DOCTEST(NoiseGenerators)

elsa/generators/tests/test_NoiseGenerators.cpp

+#include "doctest/doctest.h"
+#include "NoiseGenerators.h"
+#include "VolumeDescriptor.h"
+#include "testHelpers.h"
+
+using namespace elsa;
+using namespace doctest;
+
+TEST_SUITE_BEGIN("generators");
+
+TEST_CASE_TEMPLATE("Noise generators:", data_t, float, double)
+{