Search code examples
matrixeigen

Multiply a (n x 1) matrix by (n x m) matrix coefficient wise

I am trying to implement the following Bézier curve degree elevation equation using Eigen library:

enter image description here

The following code snippet is working to calculate the new control points. In this code degree is the variable n from the equation.

const size_t dimension = 3; // 2d or 3d points
const size_t degree = 3;
const size_t order = degree + 1;

// Create and fill Eigen::Matrix with original control points
Eigen::Matrix<double, order, dimension> P;

    // Fill matrix with original control points. Should be degree + 1 points.

// Calculate the new control points
Eigen::Matrix<double, degree, 1> M1 = FillElevationMatrix<double, degree>();
Eigen::Matrix<double, degree, 1> M2;
M2.setOnes();
M2 -= M1;
Eigen::Matrix<double, degree, dimension> Q;
for (size_t i = 0; i < degree; ++i) {
    Q.block(i, 0, 1, dimension) = (M1.row(i) * P.row(i)) + (M2.row(i) * P.row(i + 1));
}

Is there is a way to eliminate the loop and do the calculation in one shot? Or, more generically, How do I multiply a column of scalars (n x 1 matrix) by a n x m matrix so that just the corresponding row of the first matrix is multiplied by each element of the corresponding row in the second matrix in one operation?. The loop is doing it one row at a time. What I would like is something like this:

Q = (M1 * P.block(0, 0, degree, dimension)) + (M2 * P.block(1, 0, degree, dimension));

Solution

  • Your loop corresponds to multiplying with diagonal martrices, i.e:

    Q = M1.asDiagonal() * P.topRows<degree>() + M2.asDiagonal() * P.bottomRows<degree>();

    For initializing M1 and M2, also have a look at LinSpaced and reverse.