contrib/mul/mbl/mbl_mod_gram_schmidt.cxx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_mod_gram_schmidt.cxx
00002 #include "mbl_mod_gram_schmidt.h"
00003 //:
00004 // \file
00005 // \brief Orthogonalise a basis using modified Gram-Schmidt (and normalise)
00006 // \author Martin Roberts
00007 //
00008 // Transform basis {vk} to orthonormal basis {ek} with k in range 1..N
00009 // for j = 1 to N
00010 //     ej = vj
00011 //     for k = 1 to j-1
00012 //         ej = ej - <ej,ek>ek   //NB Classical GS has vj in inner product
00013 //     end
00014 //     ej = ej/|ej|
00015 //  end
00016 //
00017 // \verbatim
00018 //  Modifications
00019 //   Mar. 2011 - Patrick Sauer - added variant that returns normalisation multipliers
00020 // \endverbatim
00021 
00022 #include <vcl_vector.h>
00023 #include <vnl/vnl_vector.h>
00024 
00025 //=======================================================================
00026 //: Convert input basis {v} to orthonormal basis {e}
00027 // Each basis vector is a column of v, and likewise the orthonormal bases are returned as columns of e
00028 void mbl_mod_gram_schmidt(const vnl_matrix<double>& v,
00029                           vnl_matrix<double>& e)
00030 {
00031     unsigned N=v.cols();
00032 
00033     //Note internally it is easier to deal with the basis as a vector of vnl_vectors
00034     //As matrices are stored row-wise
00035     vcl_vector<vnl_vector<double > > vbasis(N);
00036     vcl_vector<vnl_vector<double > > evecs(N);
00037 
00038     //Copy into more convenient holding storage as vector of vectors
00039     //And also initialise output basis to input
00040     for (unsigned jcol=0;jcol<N;++jcol)
00041     {
00042         evecs[jcol] = vbasis[jcol] = v.get_column(jcol);
00043     }
00044     evecs[0].normalize();
00045 
00046     for (unsigned j=1;j<N;++j)
00047     {
00048         //Loop over previously created bases and subtract off partial projections
00049         //Thus producing orthogonality
00050 
00051         unsigned n2 = j-1;
00052         for (unsigned k=0;k<=n2;++k)
00053         {
00054             evecs[j] -= dot_product(evecs[j],evecs[k]) * evecs[k];
00055         }
00056         evecs[j].normalize();
00057     }
00058 
00059     //And copy into column-wise matrix (kth base is the kth column)
00060     e.set_size(v.rows(),N);
00061     for (unsigned jcol=0;jcol<N;++jcol)
00062     {
00063         e.set_column(jcol,evecs[jcol]);
00064     }
00065 }
00066 
00067 //: Convert input basis {v} to orthonormal basis {e}
00068 // Each basis vector is a column of v, and likewise the orthonormal bases are returned as columns of e
00069 // The multipliers used to normalise each vector in {e} are returned in np
00070 void mbl_mod_gram_schmidt(const vnl_matrix<double>& v,
00071                           vnl_matrix<double>& e, vnl_vector<double>& np)
00072 {
00073     unsigned N=v.cols();
00074     np.set_size(N);
00075 
00076     //Note internally it is easier to deal with the basis as a vector of vnl_vectors
00077     //As matrices are stored row-wise
00078     vcl_vector<vnl_vector<double > > vbasis(N);
00079     vcl_vector<vnl_vector<double > > evecs(N);
00080 
00081     //Copy into more convenient holding storage as vector of vectors
00082     //And also initialise output basis to input
00083     for (unsigned jcol=0;jcol<N;++jcol)
00084     {
00085         evecs[jcol] = vbasis[jcol] = v.get_column(jcol);
00086     }
00087     np[0]     = evecs[0].magnitude();
00088     evecs[0] /= np[0];
00089 
00090     for (unsigned j=1;j<N;++j)
00091     {
00092         //Loop over previously created bases and subtract off partial projections
00093         //Thus producing orthogonality
00094 
00095         unsigned n2 = j-1;
00096         for (unsigned k=0;k<=n2;++k)
00097         {
00098             evecs[j] -= dot_product(evecs[j],evecs[k]) * evecs[k];
00099         }
00100         np[j]     = evecs[j].magnitude();
00101         evecs[j] /= np[j];
00102     }
00103 
00104     //And copy into column-wise matrix (kth base is the kth column)
00105     e.set_size(v.rows(),N);
00106     for (unsigned jcol=0;jcol<N;++jcol)
00107     {
00108         e.set_column(jcol,evecs[jcol]);
00109     }
00110 }