lsqnonlin (Matlab) или эквивалент SciPy наименьших квадратов в C ++ - PullRequest
0 голосов
/ 28 февраля 2020

Я уже несколько дней пытаюсь понять, как перенести код с MatLab или Python на C ++ ...

Моя проблема с переносом функции оптимизации нелинейных наименьших квадратов и Я в основном преодолел блокпост, поскольку до сих пор все работает отлично.

Реализация MathLab: https://github.com/SuwoongHeo/Deformation-Transfer-Matlab

Python Реализация: https://github.com/ziyeshanwai/python-deformation-transfer

Точный бит кода I `у меня проблемы с переносом ...

Код MatLab:

[b, res, resi] = lsqnonlin(@(b) resSimXform( b,double(A),double(B) ),(double(b0)),[],[],options);

function result = resSimXform( b,A,B )

r = b(1:4);
t = b(5:7);
s = b(8);
% s = b(8:10);
n = size(A,2);

if ~isreal(r)
    a = 1;
end
R = vrrotvec2mat(r);

test =repmat(t', 1, n);
rot_A = diag(s) * R * A + repmat(t', 1, n);

result = sum(sum((B-rot_A).^2,2));
end

и эквивалентный python код:

b = least_squares(fun=resSimXform, x0=b0, jac='3-point', method='lm', args=(Points_A, Points_B),
                      ftol=1e-12, xtol=1e-12, gtol=1e-12, max_nfev=100000)

def resSimXform(b, A, B):
    print("resSimXform function")
    t = b[4:7]
    R = np.zeros((3, 3))
    R = R_axis_angle(R, b[0:3], b[3])
    rot_A = b[7]*R.dot(A) + t[:, np.newaxis]
    result = np.sqrt(np.sum((B-rot_A)**2, axis=0))
    return result

Я пробовал различные оптимизация библиотек на C ++ без удачи ... Видимо, они работают иначе, чем в Python или Matlab.

Я пробовал Eigen LevenbergMarquardt и DLibs solve_least_squares_lm, но не смог ... Я не мог понять, как их настроить, поскольку, очевидно, они отличаются от своих аналогов.

Я уже реализовал функцию «Остаток сходства» (resSimXform) в C ++, и она отлично работает. Я постоянно отлаживаю результаты и проверяю на соответствие выводам из MathLab и Python.

Моя версия на c ++ выглядит примерно так: (Я попытаюсь очистить ее больше, когда выясню, как используйте его с функцией наименьшего квадрата ... но в данный момент он возвращает точное значение как Matlab)

double residual(MatrixXd x, MatrixXd A, MatrixXd B) {
    MatrixXd r(1, 4);
    r(0, 0) = x(0, 0);
    r(0, 1) = x(0, 1);
    r(0, 2) = x(0, 2);
    r(0, 3) = x(0, 3);

    MatrixXd t(1, 3);
    t(0, 0) = x(0, 4);
    t(0, 1) = x(0, 5);
    t(0, 2) = x(0, 6);

    double s = x(0, 7);


    MatrixXd R(3, 3);
    R = R_axis_angle(r);


    MatrixXd rot_A(A.rows(), A.cols());

    t.transposeInPlace();




    t = t.col(0).replicate(1, A.cols());


    rot_A = s * R * A + t;


    MatrixXd fvecM = B - rot_A;
    for (int i = 0; i < fvecM.rows(); i++) {
        for (int j = 0; j < fvecM.cols(); j++) {
            fvecM(i, j) = pow(fvecM(i, j), 2);
        }
    }

    Eigen::VectorXd sums(3);
    sums(0) = fvecM.row(0).sum();
    sums(1) = fvecM.row(1).sum();
    sums(2) = fvecM.row(2).sum();
    return fvecM.sum();
}

Я был бы очень признателен, если кто-нибудь может помочь мне выяснить, какая библиотека оптимизации для c ++ будет работать в этом сценарии и как его реализовать ... если есть огромные отличия от таковых в Matlab или Python, так как мой опыт работы с такого рода оптимизаторами ограничен.

Заранее спасибо ... Я уже пробовал гуглить пару дней и не смог найти решение, поэтому я здесь:)

1 Ответ

0 голосов
/ 29 февраля 2020

Разобрался с Eigen`s не поддерживается LevenbergMarquardt. Вы создаете функтор ... my_functor (void): Functor (8, 8) {} .. Если вы хотите решить 8 значений ... это должно быть 8x8 ...

Операторная функция ... это, в основном, остаток сходства ...

также, если ваша функция функтора возвращает только одно значение ... вам нужно передать его 8 раз, иначе оно не сойдет, если вы просто передадите его, как fve c (0) = результат

for (int i = 0; i < x.rows(); i++) {

    fvec(i) = sum;
}

И затем вы используете режим NumericDiff ...

Eigen::NumericalDiff<my_functor> numDiff(functor);
Eigen::LevenbergMarquardt<Eigen::NumericalDiff<my_functor>, double> lm(numDiff);

Надеюсь, это кому-нибудь поможет ... Ура

...