Go to the documentation of this file.00001
00002 #include "mbl_mod_gram_schmidt.h"
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <vcl_vector.h>
00023 #include <vnl/vnl_vector.h>
00024
00025
00026
00027
00028 void mbl_mod_gram_schmidt(const vnl_matrix<double>& v,
00029 vnl_matrix<double>& e)
00030 {
00031 unsigned N=v.cols();
00032
00033
00034
00035 vcl_vector<vnl_vector<double > > vbasis(N);
00036 vcl_vector<vnl_vector<double > > evecs(N);
00037
00038
00039
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
00049
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
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
00068
00069
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
00077
00078 vcl_vector<vnl_vector<double > > vbasis(N);
00079 vcl_vector<vnl_vector<double > > evecs(N);
00080
00081
00082
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
00093
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
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 }