# Transforming from spherical coordinates to Cartesian coordinates using Eigen

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.

• Don't post the exact same question on multiple Stack Exchange sites. – Mike Borkland Jul 27 '18 at 21:39
• @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. – Soo Jul 28 '18 at 17:32

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);


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)
{

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};
}


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.

• 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. – Soo Jul 29 '18 at 18:58
• Is it advisable to copy the azimuth and elevation into Eigen::VectorX so that I can get vectorization advantages? – Soo Jul 29 '18 at 19:08
• @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. – indi Jul 31 '18 at 5:34