Transforming from spherical coordinates to Cartesian coordinates using Eigen

Question

I need to transform the coordinates from spherical to Cartesian space using the Eigen C++ Library. The following code serves the purpose:

const int size = 1000;
    Eigen::Array<std::pair<float, float>, Eigen::Dynamic, 1> direction(size);
    for(int i=0; i<direction.size();i++)
    {
            direction(i).first = (i+10)%360; // some value for this example (denoting the azimuth angle)
            direction(i).second = (i+20)%360; // some value for this example (denoting the elevation angle)

     }

    Eigen::MatrixX<T1> transformedMatrix(3, direction.size());
    for(int i=0; i<transformedMatrix.cols(); i++)
    {
        const T1 azimuthAngle = direction(i).first*M_PI/180;    //converting to radians
        const T1 elevationAngle = direction(i).second*M_PI/180; //converting to radians

        transformedMatrix(0,i) = std::cos(azimuthAngle)*std::cos(elevationAngle);
        transformedMatrix(1,i) = std::sin(azimuthAngle)*std::cos(elevationAngle);
        transformedMatrix(2,i) = std::sin(elevationAngle);
    }

I would like to know a better implementation is possible to improve the speed. I know that Eigen has supporting functions for Geometrical transformations. But I am yet to see a clear example to implement the same. Is it also possible to vectorize the code to improve the performance?

Note: This is a duplicate posting. I think that the question will be more relevant on this site.

Answer 1

There's not really much to review here. While I have some experience with Eigen, I have no idea what SSPL is. I'm going to assume SSPL:MatrixX is basically Eigen::Matrix3Xf.

const int size = 1000;

This should probably use constexpr rather than const.

for(int i=0; i<direction.size();i++)

It's been a while since I've used Eigen, but I believe you have a bug in this for prologue. If I recall, the return type of size() for Array types is not int. I think it is actually std::ptrdiff_t (or possibly std::size_t), but it's user-customizable. If it is std::ptrdiff_t (for example), and std::ptrdiff_t is larger than int (as it is on some platforms, I think including 64-bit Windows), then you will get UB if the values get cut off.

The way to fix this is to use decltype:

for (decltype(direction.size()) i = 0; i < direction.size(); ++i)

Now the first loop is just generating test data, so let's skip down to the next loop.

const T1 azimuthAngle = direction(i).first*M_PI/180;    //converting to radians

What's infinitely better than:

auto y = /* expression with x */; // convert x to foo

is:

auto y = convert_to_foo(x);

In other words, since you're converting to radians, you should have:

constexpr auto to_radians(float v) noexcept
{
    return (v * pi<float>) / 180.0f;
}

Now, in addition, M_PI is not actually portable. If you don't care, fine, but if you care about portability, you can define a π constant either as:

constexpr auto pi = 3.14159265358979f;

or, better, as:

template <typename T>
constexpr auto pi = T(3.14159265358979L); // add as many digits of precision as you please

in which case you can even make the conversion function a template:

template <typename T>
// possibly constrain T
constexpr auto to_radians(T const& v) noexcept((v * pi<T>) / T(180))
{
    return (v * pi<T>) / T(180);
}

Any way you do it, you should end up with:

const T1 azimuthAngle = to_radians(direction(i).first);
const T1 elevationAngle = to_radians(direction(i).second);

Now in the next few lines you need the sin and cos of the azimuth and elevation. You might as well precalculate them - especially since you reuse some of them:

auto const cos_azimuth = std::cos(azimuthAngle);
auto const sin_azimuth = std::sin(azimuthAngle);
auto const cos_elevation = std::cos(elevationAngle);
auto const sin_elevation = std::sin(elevationAngle);

transformedMatrix(0,i) = cos_azimuth * cos_elevation;
transformedMatrix(1,i) = sin_azimuth * cos_elevation;
transformedMatrix(2,i) = sin_elevation;

But transforming coordinate systems seems both like something you can reuse and - more importantly - something you can test in isolation. So this should be a function:

template <typename Pair>
auto spherical_to_cartesian(Pair const& spherical)
{
    const auto azimuthAngle = to_radians(std::get<0>(spherical));
    const auto elevationAngle = to_radians(std::get<1>(spherical));

    auto const cos_azimuth = std::cos(azimuthAngle);
    auto const sin_azimuth = std::sin(azimuthAngle);
    auto const cos_elevation = std::cos(elevationAngle);
    auto const sin_elevation = std::sin(elevationAngle);

    return std::tuple{cos_azimuth * cos_elevation, sin_azimuth * cos_elevation, sin_elevation};
}

Which makes your loop:

for (decltype(transformedMatrix.cols()) i = 0; i < transformedMatrix.cols(); ++i)
{
    std::tie(
            transformedMatrix(0, i),
            transformedMatrix(1, i),
            transformedMatrix(2, i))
        = spherical_to_cartesian(direction(i));
}

Now, if you're asking about how to do geometric transforms with Eigen, that seems more like a Stack Overflow question. As for vectorization, that depends on what SSPL::MatrixX is. But the loop above can be very easily parallelized, because each transform is independent. As for how, the standard way would require that Eigen::Array and SSPL::MatrixX could be used with standard algorithms, in which case the answer would simply be:

// Hypothetical code.
std::transform(std::par_unseq, begin(direction), end(direction), begin(transformedMatrix), spherical_to_cartesian);

Or you could look into OpenMP and parallel for.