Topic: Matrix math problem - Align coordinate system

[pkohut]

I'm trying to figure out how the arx api for AcGeMatrix2d::alignCoordSystem works and just can't figure it out. Can someone help me out or provide a website with instructions to do this calculation? TIA.

The docs say:

Returns matrix that maps the coordinate system with the origin fromOrigin and the axes fromE0 and fromE1 to the coordinate system with the origin toOrigin and axes toE0 and toE1.

Code: [Select]

     (2d input vectors)                                  (matrix produced from origin e0 e1, row major)
fromOrigin = 10.0, 10.0   --->   mat1 (from coorinate system) = 1, 5, 10
    fromE0 = 1.0, 3.0                                           3, 7, 10
    fromE1 = 5.0, 7.0                                           0, 0,  1

  toOrigin = 20.0,  5.0   --->  mat2 (to coorinate system)   = 11, 17, 20
      toE0 = 11.0, 13.0                                        13, 19,  5
      toE1 = 17.0, 19.0                                     0,  0,  1


Matrix2d result = SetToAlignCoordSys(fromOrigin, fromE0, fromE1, toOrigin, toE0, toE1)
result = -3.25, 4.75,   5.0
         -4.25, 5.75, -10.0
          0.0 , 0.0 ,   1.0

What are the steps required to get the result?

[Lee Mac]

Linear Algebra isn't my strong point, so I welcome others to chime in, but this is how I see it:

You need to find a change of basis matrix, say P, from Coordinate system E0 to E1, given the basis vectors that span E0 and E1 (the coordinate axes), lets call them e0₁, e0₂, e0₃ (similarly for E1)

The basis vectors are linearly independent, so each of the basis vectors of E0 can be written as a linear combination of the basis vectors of E1, hence:

Code: [Select]


e01 = a11 e11 + a12 e12 + a13 e13
e02 = a21 e12 + a22 e12 + a23 e13
.
.
.

Where the a11, a12, etc are unique constants. So here we are saying, what do the basis vectors of E0 look like in terms of the new basis of E1?

Solving the above equations, we reach a change of basis matrix P, with columns:

Code: [Select]

( a11 a21 a31 )
( a12 a22 a13 )
( a13 a23 a33 )

Now:  E1 = P E0

[pkohut]

Linear Algebra isn't my strong point, so I welcome others to chime in, but this is how I see it:

Thanks Lee. Well be diving in further when I get back in a few hours.

[gile]

Hi,

I'm not sure to undersand the task and certainly not able to explain the method due to my English more poor than may math and coding knowlege.
But it looks like you can use the Gauss Jordan method (search in an English Wiki)

Here's a LISP implementation of the method I used to Inverse a matrix

(GaussJordan '((1.0 3.0 0.0) (5.0 7.0 0.0) (10.0 10.0 1.0)) '((11.0 13.0 0.0) (17.0 19.0 0.0) (20.0 5.0 1.0)))
returns
((-3.25 -4.25 0.0) (4.75 5.75 0.0) (5.0 -10.0 1.0))

Code: [Select]

;; GaussJordan (gile)
;; Apply the Gauss-Jordan method to two matrices
;;
;; Arguments : 2 matrices

(defun GaussJordan (m1 m2 / len mat todo cnt row col piv new)
  (setq len (length m1))
  (if (= len (length m2))
    (progn
      (setq mat (mapcar '(lambda (x1 x2) (append x1 x2)) m1 m2)
    todo mat
    cnt 0
      )
      (while todo
(setq row (nth cnt mat)
      col (mapcar '(lambda (x) (abs (car x))) todo)
)
(repeat (vl-position 0 (vl-sort-i col '>))
  (setq mat (append (vl-remove row mat) (list row))
row (nth cnt mat)
  )
)
(if (equal (setq piv (car row)) 0.0 1e-14)
  (setq mat  nil
todo nil
  )
  (setq piv  (/ 1.0 piv)
new  (mapcar '(lambda (x) (* x piv)) row)
mat  (mapcar
       '(lambda (r / e)
  (setq e (car r))
  (if (equal r row)
    (cdr new)
    (cdr (mapcar '(lambda (x n) (- x (* n e))) r new)
    )
  )
)
       mat
     )
todo mat
cnt  (1+ cnt)
  )
)
(and todo (repeat cnt (setq todo (cdr todo))))
      )
      mat
    )
  )
)

[LE3]

hi Paul,
don't think i have used AcGeMatrix2d::alignCoordSystem, but AcGeMatrix3d::setToAlignCoordSys a lot... maybe by looking on the implementation or usage, you could figure out if could works for you, i can post some partial code, if you want...

[MickD]

here's some helper methods in C# that may give some clues
basically the class does the grunt work, you feed it a point for the origin and ortho vectors for the x,y vectors for the coordinate system you want to xform from and feed the new origin and vectors for the coordinate system you want to set, basically it's the same as setting a ucs if you like.

Don't get hung up on how they work internally as there's a lot to soak up, experiment with the out put to see how you need to use them.
Just remember that your vectors need to be orthographic (ie, at perp to each other) or it will fail.
hth

Code: [Select]

namespace DCSGeo
{
/// <summary>
/// Xform - transformation/translation and geometry helper classes
/// </summary>
public class Xform
{
public static Point3d UcsToWcs(Point3d point, Database db)
{
Matrix3d mat = new Matrix3d();
mat = Matrix3d.AlignCoordinateSystem(db.Ucsorg,db.Ucsxdir,db.Ucsydir,db.Ucsxdir.CrossProduct(db.Ucsydir),
new Point3d(),Vector3d.XAxis,Vector3d.YAxis,Vector3d.ZAxis);
Point3d retPoint = point.TransformBy(mat.Inverse());
return retPoint;
}

public static Vector3d UcsToWcs(Vector3d vec, Database db)
{
Matrix3d mat = new Matrix3d();
mat = Matrix3d.AlignCoordinateSystem(db.Ucsorg,db.Ucsxdir,db.Ucsydir,db.Ucsxdir.CrossProduct(db.Ucsydir),
new Point3d(),Vector3d.XAxis,Vector3d.YAxis,Vector3d.ZAxis);
Vector3d retVec = vec.TransformBy(mat.Inverse());
return retVec;
}

public static Point3d WcsToUcs(Point3d point, Database db)
{
Matrix3d mat = new Matrix3d();
mat = Matrix3d.AlignCoordinateSystem(db.Ucsorg,db.Ucsxdir,db.Ucsydir,db.Ucsxdir.CrossProduct(db.Ucsydir),
new Point3d(),Vector3d.XAxis,Vector3d.YAxis,Vector3d.ZAxis);
Point3d retPoint = point.TransformBy(mat);
return retPoint;
}

public static void WcsToUcs(Entity ent, Database db)
{
Matrix3d mat = new Matrix3d();
mat = Matrix3d.AlignCoordinateSystem(db.Ucsorg,db.Ucsxdir,db.Ucsydir,db.Ucsxdir.CrossProduct(db.Ucsydir),
new Point3d(),Vector3d.XAxis,Vector3d.YAxis,Vector3d.ZAxis);
ent.TransformBy(mat);
}
}
}

[pkohut]

Thanks everyone. Definitely should be able to figure this out now.  I'm using Eigen as a matrix and vector library for a stand alone application (was looking at Blitz++ and TvMet), and have wrapped most of the AcGeVector2d and AcGeMatrix2d functionality using it.  Initially I was doing this because the Blitz++ and TvMet libaries are so cryptic to use, and I wanted to wrap them into something familiar. Anyways, found Eigen on Wednesday, and have re-wrapped most everything.  So far alignCoordSys (setToAlignCoordSys) has been the only one I've not been able to wrap my head around.

http://eigen.tuxfamily.org/
http://www.oonumerics.org/blitz/
http://tvmet.sourceforge.net/

[pkohut]

I've posted a question to this problem in the Eigen forum, hopefully someone will reply and say "Oh ya, just do this and that, bang."

In the meantime I found a working method using SysMat (http://pagesperso-orange.fr/jean-pierre.moreau/c_matrices.html), written in C.  I've got it converted and using the Eigen matrix library, just need to clean it up a bit more and switch it to use templates.

Thanks for heads up on what to look for guys, this has been a pretty fun research project.

[gile]

Paul,

As far as I understand the method shown in the link you posted uses the Gauss Jordan elimination method as the LISP I posted up there.

[pkohut]

Paul,

As far as I understand the method shown in the link you posted uses the Gauss Jordan elimination method as the LISP I posted up there.

Ya, they produce the same results. 

[pkohut]

Thought it might interest some comparing Giles lisp version to a C++ version.  This version still needs some code refactoring done, for one it calcs column major instead of row major, hence the 2 transpose calls at the end.  Also, as written, the GaussJordanElimination function requires matrices 1 x 1 bigger than the input matrices, requiring extra matrix copying.  It's based on the SysMat function linked to in a earlier post.
Attached is a compile executable of the code below.

Code: [Select]

#include "stdafx.h"
#include <eigen\Eigen>
#include <iostream>

using namespace std;

template<class T, int NSize>
int GaussJordanElimination(Eigen::Matrix<T, NSize, NSize> & mm1, Eigen::Matrix<T, NSize, NSize> & mm2, T & det);

int _tmain(int argc, _TCHAR* argv[])
{
    double det;
    // Create 3 x 3 matrix and populate
    Eigen::Matrix<double, 3, 3> m1, m2;
    m1 << 1.0, 3.0, 0.0,
        5.0, 7.0, 0.0,
        10.0, 10.0, 1.0;

    m2 << 11.0, 13.0, 0.0,
        17.0, 19.0, 0.0,
        20.0, 5.0, 1.0;

    cout << "Calculating 3 x 3 Gauss Jordan Elimination" << endl;
    cout << "Input matrix m1" << endl << m1 << endl;
    cout << "Input matrix m2" << endl << m2 << endl;
    if(GaussJordanElimination<double, 3>(m1, m2, det)) {
        cout << "Gauss Jordan Elimination Calculated" << endl;
        cout << "Output matrix: Inverse m1" << endl << m1 << endl;
        cout << "Output matrix: Solution" << endl << m2 << endl;
        cout << "Output: Determinant = " << det << endl;
    }
    else
        cout << "Gauss Jordan Elimination not computed" << endl;
    cout << endl;

    Eigen::Matrix<double, 4, 4> m3, m4;
    m3 << 1.0, 3.0, 5.0, 7.0,
        9.0, 11.0, 13.0, 17.0,
        19.0, 23.0, 29.0, 31.0,
        32.0, 34.0, 35.0, 36.0;

    m4 << 33.0, 37.0, 39.0, 41.0,
        43.0, 47.0, 49.0, 51.0,
        53.0, 57.0, 59.0, 71.0,
        72.0, 73.0, 74.0, 75.0;

    cout << "Calculating 4 x 4 Gauss Jordan Elimination" << endl;
    cout << "Input matrix m3" << endl << m3 << endl;
    cout << "Input matrix m4" << endl << m4 << endl;
    if(GaussJordanElimination<double, 4>(m3, m4, det)) {
        cout << "Gauss Jordan Elimination Calculated" << endl;
        cout << "Output matrix: Inverse m3" << endl << m3 << endl;
        cout << "Output matrix: Solution" << endl << m4 << endl;
        cout << "Output: Determinant = " << det << endl;
    }
    else
        cout << "Gauss Jordan Elimination not computed" << endl;
    cout << endl;


    system("PAUSE");

    return 0;
}


Code: [Select]

template<int NSize>
inline void SetNaN(Eigen::Matrix<double, NSize, NSize> & m)
{
    m.fill(_Nan._Double);
}

template<int NSize>
inline void SetNaN(Eigen::Matrix<float, NSize, NSize> & m)
{
    m.fill(_Nan._Float);
}

template<class T, int NSize>
int GaussJordanElimination(Eigen::Matrix<T, NSize, NSize> & mm1, Eigen::Matrix<T, NSize, NSize> & mm2, T & det)
{
    enum { nCols = NSize, mCols = NSize, };

    Eigen::Matrix<T, NSize + 1, 1> pc, pl, cs;
    pc.fill(0.0); pl.fill(0.0); cs.fill(0.0);

    // Create matrix +1 bigger than input matrix
    Eigen::Matrix<T, NSize + 1, NSize + 1> m1, m2;
    // Copy mm1 matrix to m1
    m1.block<NSize, NSize>(0, 0) = mm1;
    // fill last row and column with 0.0
    m1.row(NSize).fill(0.0);
    m1.col(NSize).fill(0.0);

    m2.block<NSize, NSize>(0, 0) = mm2;
    m2.row(NSize).fill(0.0);
    m2.col(NSize).fill(0.0);

    det = 1.0;

    for(int k = 0; k < nCols; ++k) {
        T pv = m1(k, k);
        int ik = k;
        int jk = k;
        T pav = fabs(pv);
        for(int j, i = k; i < nCols; ++i) {
            for(j = k; j < nCols; ++j) {
                if(fabs(m1(i, j)) > pav) {
                    pv = m1(i, j);
                    pav = fabs(pv);
                    ik = i;
                    jk = j;
                }
            }

            // Search terminated, the pivot is in location i=ik, j=jk
            // Memorizing pivot location;
            pc(k) = (T) jk;
            pl(k) = (T) ik;

            // Determinant det is actualised
            // If det = 0, error message and stop
            // Machine dependant EPSMACH equals here is 1e-20

            if(ik != k) det = -det;
            if(jk != k) det = -det;
            det = det * pv;
            if(fabs(det) < 1e-20) {
                SetNaN<nCols>(mm1);
                SetNaN<nCols>(mm2);
                return 0;
            }

            // Positioning pivot in k,k
            if(ik != k)
                for(int i = 0; i < nCols; ++i) {
                    // exchange lines ik and k of matrix m1
                    T dTemp = m1(ik, i);
                    m1(ik, i) = m1(k,i);
                    m1(k,i) = dTemp;
                }

                if(mCols) {
                    for(i = 0; i < mCols; ++i) {
                        T dTemp = m2(ik, i);
                        m2(ik, i) = m2(k, i);
                        m2(k, i) = dTemp;
                    }

                    // Pivot is a correct line
                    if(jk != k)
                        for(i = 0; i < nCols; ++i) {
                            // exchange columns jk and k of matrix m1
                            T dTemp = m1(i, jk);
                            m1(i, jk) = m1(i, k);
                            m1(i, k) = dTemp;
                        }

                        // Pivot is at correct column and is located in k,k

                        // Column k of matrix m1 is stored in cs vector
                        // then column k is set to zero.
                        for(i = 0; i < nCols; ++i) {
                            cs(i) = m1(i, k);
                            m1(i, k) = 0.0;
                        }

                        cs(k) = 0.0;
                        m1(k,k) = 1.0;
                        // Line k of matrix m1 is modified
                        if(fabs(pv) < 1e-20) {
                            SetNaN<nCols>(mm1);
                            SetNaN<nCols>(mm2);
                            return 0;
                        }
                        for(i = 0; i < nCols; ++i)
                            m1(k,i) = m1(k,i) / pv;

                        if(mCols)
                            for(i = 0; i < mCols; ++i)
                                m2(k, i) = m2(k, i) / pv;

                        // Other lines of matrix m1 are modified
                        for(j = 0; j < nCols; ++j) {
                            if(j == k) j++;
                            for(i = 0; i < nCols; ++i)
                                // Line j of matrix m1 is modified
                                m1(j,i)= m1(j,i) - cs(j) * m1(k,i);
                            if(mCols)
                                for(i = 0; i < mCols; ++i)
                                    m2(j, i) = m2(j, i) - cs(j) * m2(k,i);
                        }
                        // Line k is ready

                } // mCols
        }
    } // end of k loop

    // exchange lines
    for(int i = nCols - 1; i >= 0; --i) {
        int ik = (int) pc(i);
        if(ik == i)
            continue;
        // exchange lines i and pc(i) of matrix m1
        for(int j = 0; j < nCols; ++j) {
            T dTemp = m1(i,j);
            m1(i,j) = m1(ik,j);
            m1(ik,j) = dTemp;
        }
        if(mCols) {
            for(int j = 0; j < mCols; ++j) {
                T dTemp = m2(i,j);
                m2(i,j) = m2(ik,j);
                m2(ik,j) = dTemp;
            }
        }
        // no more exchange is necessary
        // goto next line
    } // for i

    // exchange columns
    for(int j = nCols - 1; j >= 0; --j) {
        int jk = (int) pl(j);
        if(jk == j)
            continue;
        // exchange columns j et pl(j) of matrix m1
        for(int i = 0; i < nCols; ++i) {
            T dTemp = m1(i,j);
            m1(i,j) = m1(i, jk);
            m1(i, jk) = dTemp;
        }
        // No more exchange is necessary
        // go to next column
    }
    // rearrangement terminated
    // for(int i = 0; i < )

    mm1 = m1.block<NSize,NSize>(0,0).transpose();
    mm2 = m2.block<NSize,NSize>(0,0).transpose();

    return 1;
}


[pkohut]

Code: [Select]

;; GaussJordan (gile)
;; Apply the Gauss-Jordan method to two matrices
.....
(repeat (vl-position 0 (vl-sort-i col '>))
  (setq mat (append (vl-remove row mat) (list row))
row (nth cnt mat)
  )
)
.....
)

Gile super nice code!  I'm trying to figure out its magic, but am unclear what this sort does (I see what it does, just don't know why).  Can you comment on that? TIA.

[LE3]

want to play... but this is as far I am now... let's see if the next part can be resolved to find the rotation matrix that will map one matrix onto another
hope am i right [added the resultant see if is right]

Code: [Select]

;;;;;;;;;;;;;;;;;;;;;
;;; matrix original
(setq a  1.0 b  3.0 c 0.0)
(setq d  5.0 e  7.0 f 0.0)
(setq g 10.0 h 10.0 i 1.0)

(setq det (- (* a (- (* e i) (* h f))) (+ (* b (- (* d i) (* g f))) (* c (- (* d h) (* g e))))))

;;; inverse from the original
(setq ia (/ (- (* e i) (* f h)) det))
(setq ib (/ (* (- (* b i) (* h c)) -1) det))
(setq ic (/ (- (* b f) (* e c)) det))

(setq id (/ (* (- (* d i) (* f g)) -1) det))
(setq ie (/ (- (* a i) (* g c)) det))
(setq _if (/ (* (- (* a f) (* d c)) -1) det))

(setq ig (/ (- (* d h) (* g e)) det))
(setq ih (/ (* (- (* a h) (* g b)) -1) det))
(setq ii (/ (- (* a e) (* b d)) det))

(setq invmat (list (list ia ib ic) (list id ie _if) (list ig ih ii)))
return: ((-0.875 0.375 0.0) (0.625 -0.125 0.0) (2.5 -2.5 1.0))

;; and to obtain the resultant rotation matrix by multiplying both matrices (the inverse by the final one):
;;; final matrix
(setq p 11.0 q  13.0 r 0.0)
(setq s 17.0 tt 19.0 u 0.0)
(setq v 20.0 w  5.0  x 1.0)

(setq rotMat
(list
  (list
(+ (* ia p) (* ib s) (* ic v))
(+ (* ia q) (* ib tt) (* ic w))
(+ (* ia r) (* ib u) (* ic x)))

  (list
(+ (* id p) (* ie s) (* _if v))
(+ (* id q) (* ie tt) (* _if w))
(+ (* id r) (* ie u) (* _if x)))

  (list
(+ (* ig p) (* ih s) (* ii v))
(+ (* ig q) (* ih tt) (* ii w))
(+ (* ig r) (* ih u) (* ii x)))))

return=((-3.25 -4.25 0.0) (4.75 5.75 0.0) (5.0 -10.0 1.0))
;;;;;;;;;;;;;;;;;;;;;


[pkohut]

want to play... but this is as far I am now...

 

I've started converting Gile's to C++.  Just got done implementing the vl-sort-i routine. So that's about the line I'm on in his code.

[LE3]

he he... after giving a try (last night and a few hours today) ... it is not hard at all - end up being easy 
did my first part in C++ too, but as noticed for a 3x3 matrix - but it can easy to add another column and row...

want to play... but this is as far I am now...

 

I've started converting Gile's to C++.  Just got done implementing the vl-sort-i routine. So that's about the line I'm on in his code.

[pkohut]

This is what I've got so far this morning, maybe 1 hour coding/API lookup and 2 hours studing Gile's code.  Any time to work on it this weekend will have to be sneaked in, as the kids are over.  One of them has been over my shoulder watching me for the last 10 minutes, and keeps asking when I'll be done.  Right Now!

Code: [Select]

struct SortedIndex
{
    SortedIndex(int nPos, const double * pVal) : m_nPos(nPos), m_pVal(pVal) {}
    int Position(void) const { return m_nPos; }

    int m_nPos;
    const double * m_pVal;
};

inline bool sort(const SortedIndex & si1, const SortedIndex & si2)
{
    return (*si1.m_pVal > *si2.m_pVal);
}

inline void Sort_I(const double * pData, std::vector<SortedIndex> & index)
{
    std::sort(index.begin(), index.end(), sort);
}

inline int Find(std::vector<SortedIndex> & index, int nVal)
{
    for(std::vector<SortedIndex>::iterator it = index.begin(); it != index.end(); ++it) {
        if(it->Position() == nVal) {
            return it->Position();
        }
    }
    return -1;
}

template<class T, int NSize>
int Giles(const Eigen::Matrix<T, NSize, NSize> & m1, const Eigen::Matrix<T, NSize, NSize> & m2, Eigen::Matrix<T, NSize, NSize> & m3)
{
    enum { nCols = NSize, mCols = NSize, };
    Eigen::Matrix<T, NSize, NSize * 2> mm;
    mm.block<NSize, NSize>(0,0) = m1;//.transpose();
    mm.block<NSize, NSize>(0, NSize) = m2;//.transpose();

    for(int i = 0; i < nCols; ++i) {

        Eigen::Matrix<T, NSize * 2, 1> row;
        Eigen::Matrix<T, NSize, 1> col;
        row = mm.row(i);
        col = mm.col(i);

        std::vector<SortedIndex> index;
        for(int i = 0; i < NSize; ++i)
            index.push_back(SortedIndex(i, &col(i)));
        Sort_I(row.data(), index);

        int nPos = Find(index, 0);

//        for(int i = 0; i < index.)
    }
    return 0;
}


Code: [Select]

int _tmain(int argc, _TCHAR * argv[])
{
    Eigen::Matrix<double, 3, 3> m1, m2, m3;

    m1 << 107, 105, 103, 99, 115, 87, 123, 125, 127;
    m2 << 253, 10, 257, 263, 50, 267, 1, 275, 277;

    int nVal = Giles<double, 3>(m1, m2, m3);
    return 0;
}


[LE3]

give a few minutes, and will post what I have in C++.... here it goes (does not verify if the determinant is a non-zero and zero error checking):

Code: [Select]

typedef float M3[9];
static void test1 (void)
{
M3 origMat;
origMat[0] = 1.0;
origMat[1] = 3.0;
origMat[2] = 0.0;
origMat[3] = 5.0;
origMat[4] = 7.0;
origMat[5] = 0.0;
origMat[6] = 10.0;
origMat[7] = 10.0;
origMat[8] = 1.0;

M3 finalMat;
finalMat[0] = 11.0;
finalMat[1] = 13.0;
finalMat[2] = 0.0;
finalMat[3] = 17.0;
finalMat[4] = 19.0;
finalMat[5] = 0.0;
finalMat[6] = 20.0;
finalMat[7] = 5.0;
finalMat[8] = 1.0;

float det;
det = origMat[0] * (origMat[4] * origMat[8] - origMat[7] * origMat[5])
- origMat[1] * (origMat[3] * origMat[8] - origMat[6] * origMat[5])
+ origMat[2] * (origMat[3] * origMat[7] - origMat[6] * origMat[4]);

M3 invMat;
invMat[0]=(((origMat[4] * origMat[8]) - (origMat[5] * origMat[7])) / det);
invMat[1]=(-(((origMat[1] * origMat[8]) - (origMat[7] * origMat[2]))) / det);
invMat[2]=(((origMat[1] * origMat[5]) - (origMat[4] * origMat[2])) / det);

invMat[3]=(-(((origMat[3] * origMat[8]) - (origMat[5] * origMat[6]))) / det);
invMat[4]=(((origMat[0] * origMat[8]) - (origMat[6] * origMat[2])) / det);
invMat[5]=(-(((origMat[0] * origMat[5]) - (origMat[3] * origMat[2]))) / det);

invMat[6]=(((origMat[3] * origMat[7]) - (origMat[6] * origMat[4])) / det);
invMat[7]=(-(((origMat[0] * origMat[7]) - (origMat[6] * origMat[1]))) / det);
invMat[8]=(((origMat[0] * origMat[4]) - (origMat[1] * origMat[3])) / det);

M3 m;
m[0]=invMat[0] * finalMat[0] + invMat[1] * finalMat[3] + invMat[2] * finalMat[6];
m[1]=invMat[0] * finalMat[1] + invMat[1] * finalMat[4] + invMat[2] * finalMat[7];
m[2]=invMat[0] * finalMat[2] + invMat[1] * finalMat[5] + invMat[2] * finalMat[8];

m[3]=invMat[3] * finalMat[0] + invMat[4] * finalMat[3] + invMat[5] * finalMat[6];
m[4]=invMat[3] * finalMat[1] + invMat[4] * finalMat[4] + invMat[5] * finalMat[7];
m[5]=invMat[3] * finalMat[2] + invMat[4] * finalMat[5] + invMat[5] * finalMat[8];

m[6]=invMat[6] * finalMat[0] + invMat[7] * finalMat[3] + invMat[8] * finalMat[6];
m[7]=invMat[6] * finalMat[1] + invMat[7] * finalMat[4] + invMat[8] * finalMat[7];
m[8]=invMat[6] * finalMat[2] + invMat[7] * finalMat[5] + invMat[8] * finalMat[8];

acutPrintf(_T("\n (%g %g %g) (%g %g %g) (%g %g %g) \n"), m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8]);
}

[gile]

Code: [Select]

;; GaussJordan (gile)
;; Apply the Gauss-Jordan method to two matrices
.....
(repeat (vl-position 0 (vl-sort-i col '>))
  (setq mat (append (vl-remove row mat) (list row))
row (nth cnt mat)
  )
)
.....
)

Gile super nice code!  I'm trying to figure out its magic, but am unclear what this sort does (I see what it does, just don't know why).  Can you comment on that? TIA.

Thanks Paul.

As you certainly know, while treating a row in the process, the 'pivot' can't be 0 (and if all rows first items are 0 the equations can't be solved).
So, rather than just testing if the 'pivot' is 0 and permute rows until the the 'pivot' of this index row is different from 0 (as I did in the first releases of this routine), I permute rows until the one with the greater 'pivot' absolute value becomes this index row.
This is a little more time expensive but allows a better accuracy and avoid some errors if the 'pivot' isn't 0 but very closed to (thanks to Tim Willey for reporting this bug).

[gile]

I tried to write a little C# Matrix class there with including a GaussJordan method.

[pkohut]

I've posted a question to this problem in the Eigen forum, hopefully someone will reply and say "Oh ya, just do this and that, bang."

Response from the Eigen forum:  

you just have to solve a linear system of equation m1 * result = m2. For instance you can do:
result = m1.inverse() * m2;
or:
result = m1.lu().solve(m2);

For 3x3 matrices I don't know which option is the fastest.
In both cases you need to #include <Eigen/LU>

Code: [Select]

#include <eigen\Eigen>
#include <iostream>

int _tmain(int argc, _TCHAR * argv[])
{
    Eigen::Matrix<double, 3, 3> m1, m2, m3;

    m1 << 1, 3, 0, 5, 7, 0, 10, 10, 1;
    m2 << 11, 13, 0, 17, 19, 0, 20, 5, 1;

    m3 = m1.inverse() * m2;
    cout << m3 << endl << endl;

    system("PAUSE");
    return 0;
}


Code: [Select]

Results -
-3.25 -4.25 0
4.75 5.75 0
5 -10 1

Press any key to continue . . .

Now, just need to be careful when visualizing how the data is stored in memory, and remember that Eigen defaults to column major, OpenGL prefers column major, and DirectX prefers row major.

edit: Now I feel kind of foolish for not asking yesterday if this was a Ax=B problem. I assume because I'd tried result = m1.inverse() * m2; and bVal = m1.lu().solve(m2, &m3); yesterday and didn't get the correct results that it wasn't, though my research and Lee's hint suggested it was it was. And the reason for the wrong results is I'd stuffed the arrays in the wrong order.  Thanks again everyone.

[pkohut]

I tried to write a little C# Matrix class there with including a GaussJordan method.

If we could hold pointers to the rows and columns and shuffle those around instead of the actual rows and columns, then that would be a pretty good optimization.  I had been thinking about that all day for the port of your lisp code.

[LE3]

Do not know if you tried what I posted (just need the check of the determinant value ie: if (fabs(det) < 0.0005) return;) but it is really way to easy (after seeing the results, and some coffees)

Matrix rotation = Matrix final x (Matrix original) -1 = Then you get the inverse from the original and multiplying by the final one and making sure the determinant from the original is a non-zero
Mr = Mf * (Mo)-1

Did not tested with any other values... hth

edit: Now I feel kind of foolish for not asking yesterday if this was a Ax=B problem. I assume because I'd tried result = m1.inverse() * m2; and bVal = m1.lu().solve(m2, &m3); yesterday and didn't get the correct results that it wasn't, though my research and Lee's hint suggested it was it was. And the reason for the wrong results is I'd stuffed the arrays in the wrong order. 

[pkohut]

Do not know if you tried what I posted (just need the check of the determinant value ie: if (fabs(det) < 0.0005) return;) but it is really way to easy (after seeing the results, and some coffees)

Matrix rotation = Matrix final x (Matrix original) -1 = Then you get the inverse from the original and multiplying by the final one and making sure the determinant from the original is a non-zero
Mr = Mf * (Mo)-1

I'd looked at it but thought it looked way to easy, especially since I'd tried something similar and didn't get the expected results, that being the case of "garbage in - garbage out"

[LE3]

I'd looked at it but thought it looked way to easy, especially since I'd tried something similar and didn't get the expected results, that being the case of "garbage in - garbage out"

yep... and does return the same output as the ones posted here... here is again maybe with some better format, just to have a different approach

Code: [Select]

typedef float M3[9];

static BOOL mapCoordSysM3 (M3 Mo, M3 Mf, M3 &m)
{
float det;
det = Mo[0] * (Mo[4] * Mo[8] - Mo[7] * Mo[5])
- Mo[1] * (Mo[3] * Mo[8] - Mo[6] * Mo[5])
+ Mo[2] * (Mo[3] * Mo[7] - Mo[6] * Mo[4]);

if (fabs(det) < 0.0005) return FALSE;

M3 iM;
iM[0] = (((Mo[4] * Mo[8]) - (Mo[5] * Mo[7]))    / det);
iM[1] = (-(((Mo[1] * Mo[8]) - (Mo[7] * Mo[2]))) / det);
iM[2] = (((Mo[1] * Mo[5]) - (Mo[4] * Mo[2]))    / det);

iM[3] = (-(((Mo[3] * Mo[8]) - (Mo[5] * Mo[6]))) / det);
iM[4] = (((Mo[0] * Mo[8]) - (Mo[6] * Mo[2]))    / det);
iM[5] = (-(((Mo[0] * Mo[5]) - (Mo[3] * Mo[2]))) / det);

iM[6] = (((Mo[3] * Mo[7]) - (Mo[6] * Mo[4]))    / det);
iM[7] = (-(((Mo[0] * Mo[7]) - (Mo[6] * Mo[1]))) / det);
iM[8] = (((Mo[0] * Mo[4]) - (Mo[1] * Mo[3]))    / det);

m[0] = iM[0] * Mf[0] + iM[1] * Mf[3] + iM[2] * Mf[6];
m[1] = iM[0] * Mf[1] + iM[1] * Mf[4] + iM[2] * Mf[7];
m[2] = iM[0] * Mf[2] + iM[1] * Mf[5] + iM[2] * Mf[8];

m[3] = iM[3] * Mf[0] + iM[4] * Mf[3] + iM[5] * Mf[6];
m[4] = iM[3] * Mf[1] + iM[4] * Mf[4] + iM[5] * Mf[7];
m[5] = iM[3] * Mf[2] + iM[4] * Mf[5] + iM[5] * Mf[8];

m[6] = iM[6] * Mf[0] + iM[7] * Mf[3] + iM[8] * Mf[6];
m[7] = iM[6] * Mf[1] + iM[7] * Mf[4] + iM[8] * Mf[7];
m[8] = iM[6] * Mf[2] + iM[7] * Mf[5] + iM[8] * Mf[8];
return TRUE;
}

static void usageMapCoordSysM3 (void)
{
M3 Mo;
Mo[0] = 1.0;
Mo[1] = 3.0;
Mo[2] = 0.0;
Mo[3] = 5.0;
Mo[4] = 7.0;
Mo[5] = 0.0;
Mo[6] = 10.0;
Mo[7] = 10.0;
Mo[8] = 1.0;

M3 Mf;
Mf[0] = 11.0;
Mf[1] = 13.0;
Mf[2] = 0.0;
Mf[3] = 17.0;
Mf[4] = 19.0;
Mf[5] = 0.0;
Mf[6] = 20.0;
Mf[7] = 5.0;
Mf[8] = 1.0;

M3 m;
if (mapCoordSysM3(Mo, Mf, m) == TRUE)
{
acutPrintf(_T("\n (%g %g %g) (%g %g %g) (%g %g %g) \n"),
m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8]);
}
}

static void SDSArxCmdsTester_TEST(void)
{
usageMapCoordSysM3();
}

Command: TEST
(-3.250000 -4.250000 0.000000) (4.750000 5.750000 0.000000) (5.000000 -10.000000 1.000000)
 

[pkohut]

Finally finished the drop in replacements for AcGeMatrix2d, AcGeVector2d, and validation tests. 

As mentioned before, needed something familiar for my standalone app.  Thanks everyone for your help.