Собственный: оценить операции от постоянной матрицы до новой матрицы/массива
Допустим, у меня есть функция , шаблонированная в соответствии с документацией Эйгена , чтобы использовать ее как из C++, так и из Python с помощью pybind11.
Основная цель этой функции - выполнить преобразование декартовой -> полярной координаты. Вот код этой функции:
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;
}
Согласно правилам Эйгена о ленивых вычислениях, строка (1) берет данные из третьего столбца cartesian матрицы и eval()-u преобразует их в новую матрицу, вероятно, типа MatrixBase<Scalar, Eigen::Dynamic, 1>. Однако я получаю ошибки компиляции в строке (2), где сообщение об ошибке указывает, что я не могу применить operator/ между двумя ConstColXpr. Почему? Разве тип выражения не должен отбрасывать его константность, поскольку выражения будут считываться только из своего источника?
Я понимаю, что не могу изменить исходную (cartesian) матрицу, но разве цель выражений не заключалась в применении символических операций к матрицам? В этом случае константность матрицы исходного блока не должна мешать мне вычислить серию выражений во временной матрице и изменить ее.
И кстати: я также получаю ошибки компиляции в строках Eigen::{acos,atan}(...). Как я могу явно привести выражения MatrixBase к ArrayBase? Кажется, что метод .array() ничего не меняет при использовании этой функции ни из известной во время компиляции Matrix<...>, ни из зависящей от Python переменной pybind11::EigenDRef<...>.
В соответствии с просьбой, вот минимальный пример использования этой функции:
#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);
}
Ах я вижу. Моя операция копирования и вставки просто удалила ключевое слово template :)
Ответы 1
Прежде всего, Eigen и C++11/14 auto не очень хорошо сочетаются. Это связано с тем, что Eigen почти во всех случаях возвращает специальные шаблоны выражений при вызове методов, и сохранение этих шаблонов выражений может привести к неопределенному поведению (например, когда время жизни временных объектов не увеличивается). Если вы знаете, что делаете, это сработает, но будьте осторожны. В вашем конкретном случае это не является источником ваших проблем, и конкретное использование auto здесь безопасно, но я настоятельно рекомендую не использовать auto в сочетании с API-интерфейсами Eigen, а всегда указывать точные типы, которые вы хотите использовать. (В вашем примере auto подойдет, потому что ваш API гарантирует возврат main().)
Во-вторых, ваша проблема не имеет ничего общего с тем, что вы создали шаблон кода. Линия
const auto yx = (cartesian.col(1) / cartesian.col(0)).eval().array();
не удается, потому что PlainObject и cartesian.col(1) являются выражениями вектора-столбца, а не выражениями массива 1d, поэтому вы не можете их разделить. (Точно так же вы можете не делать cartesian.col(0).) Вам нужно использовать VectorXd{...} / VectorXd{...} перед самим разделением:
const auto yx = (cartesian.array().col(1) / cartesian.array().col(0)).eval();
И, наконец, при расчете угла не следует использовать .array(); вместо этого есть функция atan(), которая всегда генерирует правильные значения независимо от знаков входных переменных. (Например, atan2() даст неверный угол как для отрицательного значения y, так и для x.) К сожалению, Eigen не предоставляет предопределенную оболочку для этой функции, поэтому вам придется использовать стандартную библиотеку:
polar.col(1) = cartesian.array().col(1).binaryExpr(cartesian.array().col(0), [] (ScalarType y, ScalarType x) { return std::atan2(y, x); });
А также, к сожалению, atan(y/x) обычно медленнее, чем деление и простое atan2(), но если вы не абсолютно уверены, что ваши углы всегда находятся в пределах заранее определенного квадранта, вам все равно придется использовать atan(), чтобы получить правильный результат.
Кстати. зачем ты рассчитываешь atan2()? Вы пытаетесь использовать цилиндрические или сферические координаты? В любом случае acos(z) неверно. (В цилиндрическом случае вы просто сохраняете acos(z) как есть, а в сферическом вы вычисляете z, где acos(z/r) — это общий радиус в 3D.) И порядок ваших координат тоже странный — лично я бы отсортировал их канонически как r вместо (rho, phi, z) (цилиндрический случай) или (z, phi, rho) вместо (r, theta, phi) (сферический случай).
Спасибо за подсказку про atan2(), я забыл об этом. Да, я пытаюсь вычислить сферические координаты для набора точек, и я взял преобразования из PBRTv4, здесь - Я не осознавал, что на данный момент они уже нормализованы. Если вы выполняете оценку cartesian.array().col(x), копирует ли она всю матрицу во временном режиме?
Нет, большинство собственных методов по умолчанию являются выражениями. cartesian.array().col(x) — это выражение, которое ссылается на хранилище, стоящее за cartesian, и не создает копию. (По сути, это неявная ссылка.) Вот почему auto — такая плохая идея в сочетании с Эйгеном. Есть некоторые собственные методы, которые по умолчанию оценивают новое хранилище, но они встречаются редко. Даже такие вещи, как продукты матрицы, являются выражениями, но в отличие от .array() или .col() они будут оцениваться временно, если над ними выполняются другие операции.
Присвоение (expression).eval().array(); переменной auto безопасно только в том случае, если .eval() ничего не делает или expression уже является выражением массива — в противном случае вы получите ArrayWrapper вокруг временного объекта.
@TedLyngmo Эта функция создана по шаблону для использования с любой собственной матрицей. Для минимального воспроизводимого примера заполните любой
Eigen::Matrixи вызовите функцию, используя его в качестве единственного аргумента. Я обновил свой вопрос, чтобы добавить его.