Code Monkey home page Code Monkey logo

px4-matrix's Introduction

matrix

ARCHIVED - PX4-Matrix is now maintained within PX4-Autopilot. https://github.com/PX4/PX4-Autopilot/tree/4a3d64f1d76856d22323d1061ac6e560efda0a05/src/lib/matrix

A simple and efficient template based matrix library.

License

  • (BSD) The Matrix library is licensed under a permissive 3-clause BSD license. Contributions must be made under the same license.

Features

  • Compile time size checking.
  • Euler angle (321), DCM, Quaternion conversion through constructors.
  • High testing coverage.

Limitations

  • No dynamically sized matrices.

Library Overview

  • matrix/math.hpp : Provides matrix math routines.

    • Matrix (MxN)
    • Square Matrix (MxM, has inverse)
    • Vector (Mx1)
    • Scalar (1x1)
    • Quaternion
    • Dcm
    • Euler (Body 321)
    • Axis Angle
  • matrix/filter.hpp : Provides filtering routines.

    • kalman_correct
  • matrix/integrate.hpp : Provides integration routines.

    • integrate_rk4 (Runge-Kutta 4th order)

Testing

The tests can be executed as following:

mkdir build
cd build
cmake  -DTESTING=ON ..
make check

Formatting

The format can be checked as following:

mkdir build
cd build
cmake  -DFORMAT=ON ..
make check_format

Example

See the test directory for detailed examples. Some simple examples are included below:

    // define an euler angle (Body 3(yaw)-2(pitch)-1(roll) rotation)
    float roll = 0.1f;
    float pitch = 0.2f;
    float yaw = 0.3f;
    Eulerf euler(roll, pitch, yaw);

    // convert to quaternion from euler
    Quatf q_nb(euler);

    // convert to DCM from quaternion
    Dcmf dcm(q_nb);

    // you can assign a rotation instance that already exist to another rotation instance, e.g.
    dcm = euler;

    // you can also directly create a DCM instance from euler angles like this
    dcm = Eulerf(roll, pitch, yaw);

    // create an axis angle representation from euler angles
    AxisAngle<float> axis_angle = euler;

    // use axis angle to initialize a DCM
    Dcmf dcm2(AxisAngle(1, 2, 3));
    
    // use axis angle with axis/angle separated to init DCM
    Dcmf dcm3(AxisAngle(Vector3f(1, 0, 0), 0.2));

    // do some kalman filtering
    const size_t n_x = 5;
    const size_t n_y = 3;

    // define matrix sizes
    SquareMatrix<float, n_x> P;
    Vector<float, n_x> x;
    Vector<float, n_y> y;
    Matrix<float, n_y, n_x> C;
    SquareMatrix<float, n_y> R;
    SquareMatrix<float, n_y> S;
    Matrix<float, n_x, n_y> K;

    // define measurement matrix
    C = zero<float, n_y, n_x>(); // or C.setZero()
    C(0,0) = 1;
    C(1,1) = 2;
    C(2,2) = 3;

    // set x to zero
    x = zero<float, n_x, 1>(); // or x.setZero()

    // set P to identity * 0.01
    P = eye<float, n_x>()*0.01;

    // set R to identity * 0.1
    R = eye<float, n_y>()*0.1;

    // measurement
    y(0) = 1;
    y(1) = 2;
    y(2) = 3;

    // innovation
    r = y - C*x;

    // innovations variance
    S = C*P*C.T() + R;

    // Kalman gain matrix
    K = P*C.T()*S.I();
    // S.I() is the inverse, defined for SquareMatrix

    // correction
    x += K*r;

    // slicing
    float data[9] = {0, 2, 3,
                     4, 5, 6,
                     7, 8, 10
                    };
    SquareMatrix<float, 3> A(data);

    // Slice a 3,3 matrix starting at row 1, col 0,
    // with size 2 x 3, warning, no size checking
    Matrix<float, 2, 3> B(A.slice<2, 3>(1, 0));

    // this results in:
    // 4, 5, 6
    // 7, 8, 10

px4-matrix's People

Contributors

bartslinger avatar baumanta avatar birchera avatar bkueng avatar bresch avatar bugobliterator avatar dagar avatar jacobcrabill avatar jaredkw avatar jgoppert avatar jkflying avatar jlecoeur avatar julianoes avatar kamilritz avatar lorenzmeier avatar madcowswe avatar maetugr avatar mcharleb avatar mitchell-lee-93 avatar mortenfyhn avatar mrivi avatar natergator avatar nicolasm0 avatar pavel-kirienko avatar romain-chiap avatar romanbapst avatar thomasgubler avatar tsc21 avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

px4-matrix's Issues

Fix Unit test properly

I changed the hard-coded cast to float to a cast to the type, which now explodes in the unit test:
https://travis-ci.org/PX4/Matrix#L139

However, if I don't do this (and it looks right) it explodes in the PX4 unit test. I think we need to fix it for good here without hard-forcing a specific type in a templated function.

Assert Function Coverage 100%

We need a way to tell coveralls to assert that function coverage is 100%. For template functions that are not used they are thrown away in the source and it reports 100% coverage when it doesn't even call the functions.

platform testing

I was thinking about actually running the matrix tests on a pixhawk, snapdragon, etc. This ensures we don't hit compiler or other platform surprises.

Could we modify the existing matrix test so it continues to work standalone on travis-ci, but also could be built as a px4 module and in px4fmu-v2_test (PX4/PX4-Autopilot#4532). Alternatively we could just take a copy and hack them together into a test that lives in Firmware.

Problem with calculating floating point numbers with PseudoInverse library

Hello,
I had an error while using the pseudo inv library so I want to share the bug.
There were some cases that pseudo inv becomes a null matrix when I use a floating-point number matrix as an input. It was because of the tolerance in fullRankCholesky function. In case of using float, even though L(k,r) is near zero value in floating-point number (ex:5.96e-8), the tolerance might has even smaller value(ex:8.59e-12). If it happens, it makes an error and outputs the null matrix(which makes the output a zero vector).

스크린샷, 2020-08-26 13-44-21

To fix this, In my case which only using the floating-number matrix, I made the following change in void fullRankCholeskyTolerance(Type &tol) function

'tol /= 1000000000;' to 'tol = 0.000000001f;'

and added the following code right before 'if (L(k, r) > tol) {' in fullRankCholesky function

'if (L(k,r)<0.0000001f) L(k,r)=0.0f;'

the result has been changed
from

image

to

image

I don't think this is a good solution for the general case but in my case, I guess it works.

Best regards,

Mitchell Lee

Need to be 100% clear about what quaternion convention is used

There are many different combinations of quaternion conventions. All have their own peculiarities. Hamilton, STS, JPL, etc. to name a few. We need to be 100% clear about what version of operations is used.

With this said, hamilton multiplication formula would look like this:

    r(0) = p(0) * q(0) - p(1) * q(1) - p(2) * q(2) - p(3) * q(3); 
    r(1) = p(0) * q(1) + p(1) * q(0) + p(2) * q(3) - p(3) * q(2); 
    r(2) = p(0) * q(2) - p(1) * q(3) + p(2) * q(0) + p(3) * q(1); 
    r(3) = p(0) * q(3) + p(1) * q(2) - p(2) * q(1) + p(3) * q(0);

Your formula looked like this:
r(0) = p(0) * q(0) - p(1) * q(1) - p(2) * q(2) - p(3) * q(3);
r(1) = p(0) * q(1) + p(1) * q(0) - p(2) * q(3) + p(3) * q(2);
r(2) = p(0) * q(2) + p(1) * q(3) + p(2) * q(0) - p(3) * q(1);
r(3) = p(0) * q(3) - p(1) * q(2) + p(2) * q(1) + p(3) * q(0);

I'm looking into whether this confusion was causing issues for me. It definitely caused me to mix wrong formulas together. If we agree to use hamilton then following rules would apply:

  • component order: wxyz
  • algebra ij = k
  • handedness: right
  • right to left product means local to global transformation
  • default operation: vglobal = q * q(0, vlocal) * q.inversed()

This is important because if someone starts using JPL for example and mix the two in the same code, things will get really really weird in very nontrivial ways.

cmake header only library

Matrix as a header only library should still use cmake and create an INTERFACE library.

This would enable the use of find_library and handle include paths, any particular compiler flags, and versioning (if desired).

Dcm Type Error

If the data type is not float, line 533 of the file Quaternion.hpp will appear a compile error.
If Change ‘return Quaternion(Dcmf(dcm));’ to ‘return Quaternion(Dcm(dcm));’
It can adapt to more data types.

May be a bug in matrix inverse

It may be a bug in matirx inverse function.
In SquareMatrix.inv() :

// if diagonal is zero, swap with row below
if (fabsf(U(n, n)) < 1e-8f) {
    //printf("trying pivot for row %d\n",n);
    for (size_t i = 0; i < M; i++) {
        if (i == n) {
            continue;
        }

        //printf("\ttrying row %d\n",i);
        if (fabsf(U(i, n)) > 1e-8f) {
            //printf("swapped %d\n",i);
            U.swapRows(i, n);
            P.swapRows(i, n);
        }
    }
}   

may be change to----〉

// if diagonal is zero, swap with row below
if (fabsf(U(n, n)) < 1e-8f) {
    //printf("trying pivot for row %d\n",n);
    for (size_t i = n + 1; i < M; i++) {
    
        //printf("\ttrying row %d\n",i);
        if (fabsf(U(i, n)) > 1e-8f) {
            //printf("swapped %d\n",i);
            U.swapRows(i, n);
            P.swapRows(i, n);
            L.swapRows(i, n);
            L.swapCols(i, n);   
            break;                     
        }
    }
}  

The matrix my test is:

  -2.000000   1.000000   1.000000  -1.000000  -5.000000   1.000000   2.000000  -1.000000   0.000000
  -3.000000   2.000000  -1.000000   0.000000   2.000000   2.000000  -1.000000  -5.000000   3.000000
   0.000000   0.000000   0.000000   1.000000   4.000000  -3.000000   3.000000   0.000000  -2.000000
   2.000000   2.000000  -1.000000  -2.000000  -1.000000   0.000000   3.000000   0.000000   1.000000
  -1.000000   2.000000  -1.000000  -1.000000  -3.000000   3.000000   0.000000  -2.000000   3.000000
   0.000000   1.000000   1.000000  -3.000000   3.000000  -2.000000   0.000000  -4.000000   0.000000
   1.000000   0.000000   0.000000   0.000000   0.000000   0.000000  -2.000000   4.000000  -3.000000
   1.000000  -1.000000   0.000000  -1.000000  -1.000000   1.000000  -1.000000  -3.000000   4.000000
   0.000000   3.000000  -1.000000  -2.000000   2.000000   1.000000  -2.000000   0.000000  -1.000000

Quaternion multiplication

I was reviewing px4 ecl library.

I found there is difference between matlab script and source code; Quaternion multiplication order..

and I was realized that this from Matrix library.

I think in this library quaternion multiplication is inversed..

/** * Quaternion multiplication operator * * @param q quaternion to multiply with * @return product */ Quaternion operator*(const Quaternion &q) const { const Quaternion &p = *this; Quaternion r; r(0) = p(0) * q(0) - p(1) * q(1) - p(2) * q(2) - p(3) * q(3); r(1) = p(0) * q(1) + p(1) * q(0) - p(2) * q(3) + p(3) * q(2); r(2) = p(0) * q(2) + p(1) * q(3) + p(2) * q(0) - p(3) * q(1); r(3) = p(0) * q(3) - p(1) * q(2) + p(2) * q(1) + p(3) * q(0); return r; }

I don't know why. I think it could errupt many mistakes..
Is there special reason for this inversion of multiplication order?

Add class for axis angle representation

This code breaks the style of the rest of the library:

    // get rotation axis from quaternion (nonzero rotation)
    q = Quatf(cosf(1.0f / 2), 0.0f, sinf(1.0f / 2), 0.0f);
    rot = q.to_axis_angle();
    TEST(fabsf(rot(0)) < eps);
    TEST(fabsf(rot(1) -1.0f) < eps);
    TEST(fabsf(rot(2)) < eps);

    // get rotation axis from quaternion (zero rotation)
    q = Quatf(1.0f, 0.0f, 0.0f, 0.0f);
    rot = q.to_axis_angle();
    TEST(fabsf(rot(0)) < eps);
    TEST(fabsf(rot(1)) < eps);
    TEST(fabsf(rot(2)) < eps);

    // from axis angle (zero rotation)
    rot(0) = rot(1) = rot(2) = 0.0f;
    q.from_axis_angle(rot, 0.0f);
    q_true = Quatf(1.0f, 0.0f, 0.0f, 0.0f);
    TEST(fabsf(q(0) - q_true(0)) < eps);
    TEST(fabsf(q(1) - q_true(1)) < eps);
    TEST(fabsf(q(2) - q_true(2)) < eps);
    TEST(fabsf(q(3) - q_true(3)) < eps);

How to deduce the maximum force generated in X, Y and Z directions for the fully-actuated UAV?

Thank you for your reply.But How to isolate fmax in the equation?There are six variables, but there are only five equations.Therefore, I cannot infer the relationship between Fmax and rpm。
the M matrix is shown below:
image
In the above case, what is the relation of Fmax(fxmax/fymax/fzmax) with respect to the tilt Angle alpha and Beta?

I would be very grateful if you could help meThank you very much!

Implementation of quaternion derivative

Can somebody explain the reason behind the implementation of the derivative function here: https://github.com/PX4/Matrix/blob/master/matrix/Quaternion.hpp#L238-L253

Why does it have to be implemented that way, as opposed to just:
matrix::Quaternion w(0, _w(0), _w(1), _w(2));
matrix::Quaternion dq = (0.5f * q * w);

(The above code is "textbook" example but it gives wrong results (seems like it results in angles being applied in the wrong order..). Yet it is what I was able to find in the literature so far so I want to know why it does not work).

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.