Merge pull request #89 from ami-iit/add_slices_eigen

Added possibility to select slice when converting MultiDimensionalArray to an EigenMatrix

.github/workflows/ci.yml

@@ -6,7 +6,7 @@ on:
   workflow_dispatch:
 
 env:
-  Catch2_TAG: v3.0.1
+  Catch2_TAG: v3.7.1
 
 jobs:
   build:

CMakeLists.txt

@@ -3,7 +3,7 @@
 # BSD-2-Clause license (https://opensource.org/licenses/BSD-2-Clause).
 
 cmake_minimum_required(VERSION 3.10)
-project(matioCpp VERSION 0.2.6 LANGUAGES CXX)
+project(matioCpp VERSION 0.3.0 LANGUAGES CXX)
 
 # Defines the CMAKE_INSTALL_LIBDIR, CMAKE_INSTALL_BINDIR and many other useful macros.
 # See https://cmake.org/cmake/help/latest/module/GNUInstallDirs.html

README.md

@@ -17,7 +17,7 @@ It can be used for reading and writing binary MATLAB `.mat` files from C++, with
 
 The depencies are [``CMake``](https://cmake.org/) (minimum version 3.10) and [``matio``](https://github.com/tbeu/matio). While we suggest to follow the build instructions provided in the [``matio`` home page](https://github.com/tbeu/matio), it can also installed from common package managers:
 - Linux: ``sudo apt install libmatio-dev``
-- Linux, macOS, Windows, via [``conda-forge``](https://conda-forge.org/): ``mamba install -c conda-forge libmatio``
+- Linux, macOS, Windows, via [``conda-forge``](https://conda-forge.org/): ``conda install -c conda-forge libmatio``
 
 [`Eigen`](https://eigen.tuxfamily.org/index.php) is an optional dependency. If available, some conversions are defined.
 
@@ -153,6 +153,44 @@ Eigen::Vector3i eigenVec;
 eigenVec << 2, 4, 6;                                                            
 auto toMatioEigenVec = matioCpp::make_variable("testEigen", eigenVec);          
 ```
+It is also possible to slice a ``MultiDimensionalArray`` into an Eigen matrix:
+```c++
+std::vector<float> tensor(12);
+for (size_t i = 0; i < 12; ++i)
+{
+    tensor[i] = i + 1.0;
+}
+
+matioCpp::MultiDimensionalArray<float> matioCppMatrix2("matrix", { 2, 2, 3 }, tensor.data());
+
+/*
+So we have a tensor of the type
+| 1  3 | | 5 7 | |  9 11 |
+| 2  4 | | 6 8 | | 10 12 |
+*/
+
+Eigen::MatrixXf slice1 = matioCpp::to_eigen(matioCppMatrix2, { -1, -1, 0 }; //Equivalent to the Matlab operation matioCppMatrix2(:,:,1)
+/* 
+Obtain
+| 1  3 |
+| 2  4 |
+*/
+
+Eigen::MatrixXf slice2 = matioCpp::to_eigen(matioCppMatrix2, { 1, -1, -1 }; //Equivalent to the Matlab operation matioCppMatrix2(2,:,:)
+/*
+Obtain
+| 2 6 10|
+| 4 8 12|
+*/
+
+Eigen::MatrixXf slice3 = matioCpp::to_eigen(matioCppMatrix2, { -1, 0, 0 }; //Equivalent to the Matlab operation matioCppMatrix2(:,1,1)
+/*
+Obtain
+| 1 |
+| 2 |
+*/
+```
+In the slice, the value `-1` means that the entire dimension is taken.
 
 ``matioCpp`` also exploits [``visit_struct``](https://github.com/garbageslam/visit_struct) to parse C++ structs into ``matioCpp`` structs. Example:
 ```c++

include/matioCpp/EigenConversions.h

@@ -20,6 +20,12 @@
 namespace matioCpp
 {
 
+template <typename type>
+using EigenMapWithStride = Eigen::Map<Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>;
+
+template <typename type>
+using ConstEigenMapWithStride = Eigen::Map<const Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>;
+
 /**
  * @brief Conversion from a MultiDimensionalArray to an Eigen matrix
  * @param input The MultiDimensionalArray
@@ -34,7 +40,31 @@ inline Eigen::Map<Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>> to_eigen(
  * @return A const map from the internal data of the MultiDimensionalArray
  */
 template <typename type>
-inline const Eigen::Map<Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>> to_eigen(const MultiDimensionalArray<type>& input);
+inline Eigen::Map<const Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>> to_eigen(const MultiDimensionalArray<type>& input);
+
+/**
+ * @brief Conversion from a MultiDimensionalArray to an Eigen matrix
+ * @param input The MultiDimensionalArray
+ * @param slice The slice to extract from the MultiDimensionalArray. Use a negative value to extract the whole dimension.
+ *              If not provided, it is assumed that the input is a matrix.
+ *              At most 2 slices are allowed, one for the rows and one for the columns of the output matrix.
+ *              If only one dimension is selected, the output is a column vector.
+ * @return A map from the internal data of the MultiDimensionalArray
+ */
+template <typename type>
+inline EigenMapWithStride<type> to_eigen(MultiDimensionalArray<type>& input, const std::vector<int>& slice);
+
+/**
+ * @brief Conversion from a const MultiDimensionalArray to an Eigen matrix
+ * @param input The MultiDimensionalArray
+ * @param slice The slice to extract from the MultiDimensionalArray. Use a negative value to extract the whole dimension.
+ *              If not provided, it is assumed that the input is a matrix.
+ *              At most 2 slices are allowed, one for the rows and one for the columns of the output matrix.
+ *              If only one dimension is selected, the output is a column vector.
+ * @return A const map from the internal data of the MultiDimensionalArray
+ */
+template <typename type>
+inline ConstEigenMapWithStride<type> to_eigen(const MultiDimensionalArray<type>& input, const std::vector<int>& slice);
 
 /**
  * @brief Conversion from a Vector to an Eigen vector

include/matioCpp/impl/EigenConversions.tpp

@@ -10,6 +10,91 @@
 
 #include <cassert>
 
+namespace matioCpp
+{
+    struct SlicingInfo
+    {
+        std::pair<size_t, size_t> dimensions{0,0};
+        size_t offset{ 0 };
+        size_t innerStride{ 0 };
+        size_t outerStride{ 0 };
+        bool valid{ false };
+    };
+
+    template <typename type>
+    SlicingInfo computeSlicingInfo(const matioCpp::MultiDimensionalArray<type>& input, const std::vector<int>& slice)
+    {
+        std::string errorPrefix = "[ERROR][matioCpp::to_eigen] ";
+
+        using index_type = typename matioCpp::MultiDimensionalArray<type>::index_type;
+        SlicingInfo info;
+
+        const auto& dimensions = input.dimensions();
+
+        if (slice.size() != dimensions.size())
+        {
+            std::cerr << errorPrefix << "The number of slices must be equal to the number of dimensions of the input MultiDimensionalArray" << std::endl;
+            assert(false);
+            return info;
+        }
+
+        std::vector<index_type> startingIndex(slice.size(), 0);
+        index_type previousDimensionsFactorial = 1;
+
+        for (size_t i = 0; i < slice.size(); ++i)
+        {
+            if (slice[i] >= 0)
+            {
+                if (slice[i] >= dimensions(i))
+                {
+                    std::cerr << errorPrefix << "The slice is larger than the dimension of the input MultiDimensionalArray" << std::endl;
+                    assert(false);
+                    return SlicingInfo();
+                }
+                startingIndex[i] = static_cast<size_t>(slice[i]);
+            }
+            else
+            {
+                if (info.dimensions.first == 0)
+                {
+                    info.dimensions.first = dimensions(i);
+                    info.innerStride = previousDimensionsFactorial;
+                }
+                else if (info.dimensions.second == 0)
+                {
+                    info.dimensions.second = dimensions(i);
+                    info.outerStride = previousDimensionsFactorial;
+                }
+                else
+                {
+                    std::cerr << errorPrefix << "Only at most two free dimensions are allowed" << std::endl;
+                    assert(false);
+                    return SlicingInfo();
+                }
+            }
+            previousDimensionsFactorial *= dimensions(i);
+        }
+
+        if (info.dimensions.first == 0) //Single element
+        {
+            info.dimensions.first = 1;
+            info.innerStride = 1;
+        }
+
+        if (info.dimensions.second == 0) //Vector case
+        {
+            info.dimensions.second = 1;
+            info.outerStride = info.dimensions.first * info.innerStride;
+        }
+
+        info.offset = input.rawIndexFromIndices(startingIndex);
+
+        info.valid = true;
+        return info;
+
+    }
+}
+
 template <typename type>
 inline Eigen::Map<Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>> matioCpp::to_eigen(matioCpp::MultiDimensionalArray<type>& input)
 {
@@ -19,13 +104,44 @@ inline Eigen::Map<Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>> matioCpp:
 }
 
 template <typename type>
-inline const Eigen::Map<Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>> matioCpp::to_eigen(const matioCpp::MultiDimensionalArray<type>& input)
+inline Eigen::Map<const Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>> matioCpp::to_eigen(const matioCpp::MultiDimensionalArray<type>& input)
 {
     assert(input.isValid());
     assert(input.dimensions().size() == 2);
     return Eigen::Map<const Eigen::Matrix<type, Eigen::Dynamic, Eigen::Dynamic>>(input.data(), input.dimensions()(0), input.dimensions()(1));
 }
 
+template <typename type>
+inline matioCpp::EigenMapWithStride<type> matioCpp::to_eigen(matioCpp::MultiDimensionalArray<type>& input, const std::vector<int>& slice)
+{
+    assert(input.isValid());
+
+    auto slicingInfo = computeSlicingInfo(input, slice);
+    Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> stride(slicingInfo.outerStride, slicingInfo.innerStride);
+    if (!slicingInfo.valid)
+    {
+        return matioCpp::EigenMapWithStride<type>(nullptr, 0, 0, stride);
+    }
+    return matioCpp::EigenMapWithStride<type>(input.data() + slicingInfo.offset, slicingInfo.dimensions.first, slicingInfo.dimensions.second, stride);
+}
+
+template <typename type>
+inline matioCpp::ConstEigenMapWithStride<type> matioCpp::to_eigen(const matioCpp::MultiDimensionalArray<type>& input, const std::vector<int>& slice)
+{
+    assert(input.isValid());
+
+    auto slicingInfo = computeSlicingInfo(input, slice);
+    Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> stride(slicingInfo.outerStride, slicingInfo.innerStride);
+
+
+    if (!slicingInfo.valid)
+    {
+        return ConstEigenMapWithStride<type>(nullptr, 0, 0, stride);
+    }
+
+    return ConstEigenMapWithStride<type>(input.data() + slicingInfo.offset, slicingInfo.dimensions.first, slicingInfo.dimensions.second, stride);
+}
+
 template <typename type>
 inline Eigen::Map<Eigen::Matrix<type, Eigen::Dynamic, 1>> matioCpp::to_eigen(matioCpp::Vector<type>& input)
 {

test/ExogenousConversionsUnitTest.cpp

@@ -85,6 +85,56 @@ TEST_CASE("Eigen Conversions")
 
         Eigen::MatrixXf toEigenMatrix = matioCpp::to_eigen(matioCppMatrix);
         checkSameMatrix(toEigenMatrix, matioCppMatrix);
+
+        const auto& constMatioCppMatrix = matioCppMatrix;
+        Eigen::MatrixXf toEigenMatrixConst = matioCpp::to_eigen(constMatioCppMatrix);
+        checkSameMatrix(toEigenMatrixConst, constMatioCppMatrix);
+
+        std::vector<float> tensor(12);
+        for (size_t i = 0; i < 12; ++i)
+        {
+            tensor[i] = i + 1.0;
+        }
+
+        matioCpp::MultiDimensionalArray<float> matioCppMatrix2("matrix", { 2, 2, 3 }, tensor.data());
+
+        /*
+        So we have a tensor of the type
+        | 1  3 | | 5 7 | |  9 11 |
+        | 2  4 | | 6 8 | | 10 12 |
+        */
+
+        Eigen::MatrixXf expectedSlice(2, 2);
+        expectedSlice << 1.0, 3.0,
+                         2.0, 4.0;
+        REQUIRE(expectedSlice.isApprox(matioCpp::to_eigen(matioCppMatrix2, { -1, -1, 0 }), 1e-5));
+
+
+        expectedSlice.resize(2,3);
+        expectedSlice << 2.0, 6.0, 10.0,
+                         4.0, 8.0, 12.0;
+        REQUIRE(expectedSlice.isApprox(matioCpp::to_eigen(matioCppMatrix2, {1, -1, -1 }), 1e-5));
+
+        expectedSlice.resize(2, 3);
+        expectedSlice << 3.0, 7.0, 11.0,
+                         4.0, 8.0, 12.0;
+        REQUIRE(expectedSlice.isApprox(matioCpp::to_eigen(matioCppMatrix2, { -1, 1, -1 }), 1e-5));
+
+        expectedSlice.resize(2, 1);
+        expectedSlice << 11.0,
+                         12.0;
+        REQUIRE(expectedSlice.isApprox(matioCpp::to_eigen(matioCppMatrix2, { -1, 1, 2 }), 1e-5));
+
+        expectedSlice.resize(1, 1);
+        expectedSlice << 8.0;
+        REQUIRE(expectedSlice.isApprox(matioCpp::to_eigen(matioCppMatrix2, { 1, 1, 1 }), 1e-5));
+
+        const auto& constMatioCppMatrix2 = matioCppMatrix2;
+        expectedSlice.resize(2, 3);
+        expectedSlice << 1.0, 5.0, 9.0,
+                         3.0, 7.0, 11.0;
+        REQUIRE(expectedSlice.isApprox(matioCpp::to_eigen(constMatioCppMatrix2, { 0, -1, -1 }), 1e-5));
+
     }
 
     SECTION("From Eigen")