3
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Don't post the exact same question on multiple Stack Exchange sites. \$\endgroup\$ – Mike Borkland Jul 27 '18 at 21:39
  • \$\begingroup\$ @MikeBorkland I agree with you. But only after posting it on StackOverflow, I realised that the question may be more relevant here! Sorry for the trouble. \$\endgroup\$ – Soo Jul 28 '18 at 17:32
2
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Thanks for the answer. As you guessed correctly it was Eigen::Matrix3Xf. Please also comment whether I can use some vectorization for calculating the "sin" and "cos". I forgot to mention that I need to use C++ 11 only. Hence I will not be able to use std::par_unseq. \$\endgroup\$ – Soo Jul 29 '18 at 18:58
  • \$\begingroup\$ Is it advisable to copy the azimuth and elevation into Eigen::VectorX so that I can get vectorization advantages? \$\endgroup\$ – Soo Jul 29 '18 at 19:08
  • \$\begingroup\$ @Soo It has been a long time since I dealt with SSE, but I don't think it will help you at all because (last time I checked) there are no sin/cos operations. Your biggest gain will probably be with parallelization, not vectorization. \$\endgroup\$ – indi Jul 31 '18 at 5:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.