Search code examples
c++matrixeigen

Container for Eigen Matrices

I am looking for a way to store matrix-like structures of sensor measurements, that are again matrices or vectors themselves. To give a more concrete example: Lets take an image, where a single pixel's RGB value can be represented as a Eigen::Vector3f. As such an entire image would become something like

typedef Eigen::Array<Eigen::Vector3f, Eigen::Dynamic, Eigen::Dynamic> Image.

Adding one more level on top we could have

typedef Eigen::Array<Image, Eigen::Dynamic, 1> VideoRecording

or

typedef Eigen::Array<Image, Eigen::Dynamic, Eigen::Dynamic Image> CamaraArrayCapture.

Now, I don't need to do a lot of math on those outer storage containers, though some top-level multiplication and addition operator would be convenient for normalization purposes. However, what would be essential for me is block operations, e.g. for accessing an images Blue channel or retrieving all images in a VideoRecording from a certain time. Last but not least I'd like the container to be iterable.

While thinking this through I came up with the following ideas:

  1. Using plain a Eigen::Array or Eigen::Matrix as shown in the introduction and implement a few free functions returning something like of Eigen::Array<Eigen::Map<Eigen::Matrix<float,Eigen::Dynamic,1>>> for nested block operations. While this might be the fastest to implement I don't really like this, since nested Eigen Matrices are officially not really supported...

  2. Use the same as above, but instead of having free functions, wrap it in a custom class, so that the a library's user won't expect unsupported operation like VideoRecording x VideoRecording multiplication to work.

  3. Use an Eigen::Tensor from Eigen-unsupported for data storage. It does very powerful slicing operations, but I would loose all semantics about the included data. Also resizing a dataset would become quite inefficient size, since it would require to copy all data from all past images.

  4. There is an Eigen::Image class, which stores an RGB image as a 3 x (row*cols) matrix, where each column represents a pixel. Judging from the fact, that the class does not even have its own entry in Eigen's Wiki and the few Google hits I found I guess this class isn't widely known and/or used. Also it would still leave the problem on how to represent a structure of multiple images.

  5. When you tool is a hammer, everything looks like a nail - Forget about using an Eigen type for the outer structure and use something else, e.g. by implementing a class with an std::vector member. As this implies to re-build all block operation from scratch this feels a bit like re-inventing the wheel to me...

While I currently favor option two, I must admit that I am not really happy with either of these approaches... Since I am for sure not the first one asking how to represent an Image using Eigen types I am wondering if I am missing something and what would be the best practice to implement something this?

Solution

  • After thinking this through for some time and playing around with various approaches I ended up wrapping an Eigen::Tensor in a way that I hid one rank to the outside and returned and Eigen::Map to a particular "pixel". The basic idea looks as follows:

    template <Eigen::Index Rank, typename Scalar, Eigen::Index MeasurementRows>
class MeasurementTensor {
    public:
        typedef Eigen::Tensor<Scalar, Rank+1, Eigen::ColMajor> tensor_type;
        typedef Eigen::Map<Eigen::Matrix<Scalar, MeasurementRows, 1>> measurement_type;

        template<typename... IndexTypes>
        MeasurementTensor(Eigen::Index firstDimension, IndexTypes... otherDimensions)
        : _tensor(MeasurementRows, firstDimension, otherDimensions...) {
            _tensor.setZero();
        }

        template<typename... IndexTypes>
        inline measurement_type operator()(Eigen::Index firstDimension, Eigen::Index secondIndex, IndexTypes... otherDimensions) {
            return measurement_type(parent::tensor().data() + getScalarIndex(0, firstDimension, secondIndex, otherDimensions...));
        }

        inline measurement_type operator[](Eigen::Index linearIndex) const {
            return measurement_type(parent::tensor().data() + MeasurementRows * linearIndex);
        }

    private:
        template<typename... IndexTypes>
        Eigen::Index getScalarIndex(Eigen::Index firstIndex, Eigen::Index secondIndex, IndexTypes... otherIndices) const {
            return _tensor.dimensions().IndexOfColMajor(
                Eigen::array<Eigen::Index, Rank+1>{{firstIndex, secondIndex, otherIndices...}}
            );
        }

        tensor_type _tensor;
};

    /// Specialization for dynamic measurement sizes
template <Eigen::Index Rank, typename Scalar>
class MeasurementTensor<Rank, Scalar, Eigen::Dynamic> {
    ...
};

    On top I used some typedefs to keep my semantics, e.g.:

    typedef MeasurementTensor<2,float,3> RBGImage;
typedef MeasurementTensor<3,float,3> RBGImageSequence;

    Unfortunately Eigen's slicing operations on tensors don't seem to be as straightforward to use as its Matrices block object. So far I did not find a way to access a single value to an Eigen::TensorSlicingOp without assigning it to a second Tensor first. Thus slicing operations always require allocation of new memory and copying all values... However, in most cases it is possible to use a same class as above that wraps an Eigen::TensorMap to work an a subset of the Tensor.

    Since I don't do online processing I also decided to use an std::vector<RBGImage> for more efficient collection of raw data and only copy all data together in one continuous piece of memory once the recording has been completed.