Сокращение строк расширенной матрицы - 3D сплайн расчеты - PullRequest
Новая
       40

Сокращение строк расширенной матрицы - 3D сплайн расчеты

0 голосов
/ 01 ноября 2019

Я пытаюсь реализовать Интерполяцию по расслабленным кубическим сплайнам , которую можно найти в 5-й главе этой статьи (стр. 9): https://www.math.ucla.edu/~baker/149.1.02w/handouts/dd_splines.pdf

Пока что яесть следующее: 

auto GetControlPoints = [](const std::vector<Vector3d>& S) {
    int n = S.size();
    float var = n - 1.0f;
    MatrixXd M(n - 1, n - 1);
    VectorXd C[3] = {
        VectorXd(n - 1),
        VectorXd(n - 1),
        VectorXd(n - 1)
    };

    for (int i = 0; i < n - 1; ++i) {
        auto r = RowVectorXd(n - 1);

        for (int j = 0; j < n - 1; ++j) {
            if (j == i)
                r[j] = var;
            else if (j == i - 1 || j == i + 1)
                r[j] = 1.f;
            else
                r[j] = 0.f;
        }

        M.row(i) = r;

        if (i == 0) {
            for (int j = 0; j < 3; ++j) {
                C[j] << (n + 1) * S[1][j] - S[0][j];
            }
        }
        else if (i == n - 1) {
            for (int j = 0; j < 3; ++j) {
                C[j] << (n + 1) * S[n - 1][j] - S[n][j];
            }
        }
        else {
            for (int j = 0; j < 3; ++j) {
                C[j] << (n + 1) * S[i][j];
            }
        }
    }

    MatrixXd augMC[3] = {
        MatrixXd(n - 1, n),
        MatrixXd(n - 1, n),
        MatrixXd(n - 1, n)
    };

    for (int i = 0; i < 3; ++i) {
        augMC[i].block(0, 0, n - 1, n - 1) = M;
        augMC[i].block(n - 1, n - 1, n - 1, 1) = C[i].transpose();
    }
};

Я дошел до того, что я создал расширенную матрицу, используя M и C, но я понятия не имею, как ее уменьшить. Есть мысли?

1 Ответ

0 голосов
/ 01 ноября 2019

Вы можете использовать вариант размещения PartialPivLU - но похоже, что вы действительно хотите решить систему M*B = C, для которой вам нужно просто разложить M (так как она симметрична)вы можете использовать декомпозицию LLt или LDLt), а затем использовать метод декомпозиции solve.

Для настройки M вам также следует использовать метод diagonal (не проверено): 

MatrixXd M(n - 1, n - 1);
M.setZero();
M.diagonal().setConstant(n - 1.0);
M.diagonal<1>().setOnes();
M.diagonal<-1>().setOnes();

LLT<MatrixXd> lltOfM(M);
for (int i = 0; i < 3; ++i) { B[i] = lltOfM.solve(C[i]); }

Для больших n это неоптимально, поскольку в нем не используется трехдиагональная структура M. Вы можете попробовать модуль разреженности для этого, но на самом деле должен быть прямой алгоритм (хотя у Eigen явно нет трехдиагонального типа матрицы).

Для C вы, вероятно, также можете использовать MatrixX3d (Я не совсем понимаю, как вы заполняете свои C векторы - я думаю, что ваш текущий код должен утверждаться во время выполнения, если только n==2 или вы не отключили утверждения).

...