Let's say I have a function templated according to Eigen's documentation, in order to use it both from C++ and Python using pybind11.
The main goal of that function if to perform a cartesian -> polar coordinate conversion. Here is the code to that function:
template <typename VertexType>
typename Eigen::MatrixBase<VertexType>::PlainObject cartesian_to_polar(
const Eigen::MatrixBase<VertexType>& cartesian
) {
using MatrixType = typename Eigen::MatrixBase<VertexType>::PlainObject;
using ScalarType = typename MatrixType::Scalar;
using ReturnMatrixType = Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic, MatrixType::Options>;
ReturnMatrixType polar;
polar.resize(cartesian.rows(), 3);
/* (1) */ const auto z = cartesian.col(2).eval().array();
/* (2) */ const auto yx = (cartesian.col(1) / cartesian.col(0)).eval().array();
polar.col(0) = Eigen::acos(z);
polar.col(1) = Eigen::atan(yx);
/* snip... compute radius and store in polar.col(2) */
return polar;
}
According to Eigen's rules on lazy evaluation, the line (1) takes the data from the third cartesian matrix's column, and eval()-uates it into a new matrix, probably of type MatrixBase<Scalar, Eigen::Dynamic, 1>. However, I get compile errors on line (2), where the error message points out I cannot apply operator/ between two ConstColXpr. Why? Shouldn't the type of an expression cast away the const-ness of it, since expressions will only read from their source?
I understand that I cannot modify the original (cartesian) matrix, but wasn't the goal of expressions to apply symbolic operations on matrices ? In that case, the const-ness of the original block's matrix should not prevent me from evaluating a series of expressions into a temporary matrix and modifying it.
And, as an aside: I also get compile errors on the Eigen::{acos,atan}(...) lines. How can I explicitly cast MatrixBase expressions to ArrayBase? The .array() method does not seem to change anything when using that function either from a compile-time-known Matrix<...> or a Python-dependent pybind11::EigenDRef<...> variable.
As requested, here's a minimal example using that function:
#include "./cartesian_to_polar.h" /* Where the function above is defined. */
#include <Eigen/Eigen>
int main() {
Eigen::Matrix3f vertices;
vertices << 1.0f, .0f, .0f, .0f, 1.0f, .0f, .0f, .0f, 1.0f;
auto polar_vertices = cartesian_to_polar(vertices);
}
First of all, Eigen and C++11/14's auto don't mix very well. This is because Eigen in nearly all cases returns special expression templates when invoking methods, and storing those expression templates can lead to undefined behavior (for example when lifetimes of temporaries isn't extended). If you know what you are doing, it will work, but you should be careful. In your specific case this isn't the source of your problems, and your specific use of auto here is safe, but I would strongly recommend not using auto in conjunction with Eigen's APIs but always specifying the exact types you want to use. (The auto in your example main() would be OK because your API guarantees to return a PlainObject.)
Second, your issue has nothing to do with the fact that you've templated the code. The line
const auto yx = (cartesian.col(1) / cartesian.col(0)).eval().array();
fails because cartesian.col(1) and cartesian.col(0) are column vector expressions, and not 1d array expressions, so you can't divide them. (The same way you can also not do VectorXd{...} / VectorXd{...}.) What you want to do is use .array() before the division itself:
const auto yx = (cartesian.array().col(1) / cartesian.array().col(0)).eval();
And finally, when calculating an angle, you shouldn't use atan(); instead there is the atan2() function that always generates the correct values regardless of the signs of the input variables. (atan(y/x) will yield the wrong angle for both y and x negative, for example.) Unfortunately Eigen doesn't provide a predefined wrapper for that function, so you're going to have to use the standard library:
polar.col(1) = cartesian.array().col(1).binaryExpr(cartesian.array().col(0), [] (ScalarType y, ScalarType x) { return std::atan2(y, x); });
And also unfortunately, atan2() is typically slower than division and plain atan(), but unless you are absolutely sure that your angles are always within a predetermined quadrant, you'll still need to use atan2() to get the correct result.
Btw. why are you calculating acos(z)? Are you attempting to use cylindrical or spherical coordinates? In either case acos(z) is not correct. (In cylindrical you just keep the z as-is, and in spherical you calculate acos(z/r) with r being the total radius in 3D.) And the order of your coordinates is weird, too - I personally would sort them canonically as (rho, phi, z) instead of (z, phi, rho) (cylindrical case), or (r, theta, phi) instead of (phi, theta, r) (spherical case).