[pktools] 99/375: removed dependency of gslwrap: replaced gslwrap with armadillo in SensorModel.h and pksensormodel, replaced gslwrap with stl::vector and Vector2d
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:03 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit 8e25ca92a1e4f2751965103729316bd7ae94ab9a
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Thu Apr 25 17:38:35 2013 +0200
removed dependency of gslwrap: replaced gslwrap with armadillo in SensorModel.h and pksensormodel, replaced gslwrap with stl::vector and Vector2d
---
src/algorithms/Filter.h | 34 ++----
src/algorithms/Makefile.am | 2 +-
src/apps/Makefile.am | 5 +-
src/apps/pksensormodel.cc | 182 +++++++++++++--------------
src/apps/pksensormodel.h | 25 ++--
src/models/SensorModel.h | 298 ++++++++++++++++++++++-----------------------
6 files changed, 265 insertions(+), 281 deletions(-)
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 9d2eb9d..3da2601 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -22,9 +22,6 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <vector>
#include <iostream>
-#include <gslwrap/vector_double.h>
-#include <gslwrap/matrix_double.h>
-#include <gslwrap/matrix_vector_operators.h>
#include "StatFactory.h"
#include "imageclasses/ImgReaderGdal.h"
#include "imageclasses/ImgWriterGdal.h"
@@ -120,9 +117,6 @@ private:
gsl_spline *spline;
stat.getSpline(interpolationType,nband,spline);
stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
- // gsl_interp_accel *acc=gsl_interp_accel_alloc();
- // gsl_spline *spline=gsl_spline_alloc(gsl_interp_linear,nband);
- // gsl_spline_init(spline,&(srf[0][0]),&(srf[1][0]),nband);
if(verbose)
std::cout << "calculating norm of srf" << std::endl << std::flush;
double norm=0;
@@ -201,9 +195,6 @@ private:
gsl_spline *spline;
stat.getSpline(interpolationType,nband,spline);
stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
- // gsl_interp_accel *acc=gsl_interp_accel_alloc();
- // gsl_spline *spline=gsl_spline_alloc(gsl_interp_linear,nband);
- // gsl_spline_init(spline,&(srf[0][0]),&(srf[1][0]),nband);
if(verbose)
std::cout << "calculating norm of srf" << std::endl << std::flush;
double norm=0;
@@ -292,22 +283,23 @@ template<class T> void Filter::applyFwhm(const vector<double> &wavelengthIn, con
int nbandOut=wavelengthOut.size();
output.resize(nbandOut);
- gsl::matrix tf(nbandIn,nbandOut);
+ Vector2d<double> tf(nbandIn,nbandOut);
for(int indexOut=0;indexOut<nbandOut;++indexOut){
double norm=0;
for(int indexIn=0;indexIn<nbandIn;++indexIn){
- tf(indexIn,indexOut)=
+ // tf(indexIn,indexOut)=
+ tf[indexIn][indexOut]=
exp((wavelengthOut[indexOut]-wavelength_fine[indexIn])
*(wavelength_fine[indexIn]-wavelengthOut[indexOut])
/2.0/stddev[indexOut]
/stddev[indexOut]);
- tf(indexIn,indexOut)/=sqrt(2.0*M_PI);
- tf(indexIn,indexOut)/=stddev[indexOut];
- norm+=tf(indexIn,indexOut);
+ tf[indexIn][indexOut]/=sqrt(2.0*M_PI);
+ tf[indexIn][indexOut]/=stddev[indexOut];
+ norm+=tf[indexIn][indexOut];
}
output[indexOut]=0;
for(int indexIn=0;indexIn<nbandIn;++indexIn)
- output[indexOut]+=input_fine[indexIn]*tf(indexIn,indexOut)/norm;
+ output[indexOut]+=input_fine[indexIn]*tf[indexIn][indexOut]/norm;
}
}
@@ -335,19 +327,19 @@ template<class T> void Filter::applyFwhm(const vector<double> &wavelengthIn, con
int nbandOut=wavelengthOut.size();
output.resize(nbandOut,(input[0].size()+down-1)/down);
- gsl::matrix tf(nbandIn,nbandOut);
+ Vector2d<double> tf(nbandIn,nbandOut);
vector<double> norm(nbandOut);
for(int indexOut=0;indexOut<nbandOut;++indexOut){
norm[indexOut]=0;
for(int indexIn=0;indexIn<nbandIn;++indexIn){
- tf(indexIn,indexOut)=
+ tf[indexIn][indexOut]=
exp((wavelengthOut[indexOut]-wavelength_fine[indexIn])
*(wavelength_fine[indexIn]-wavelengthOut[indexOut])
/2.0/stddev[indexOut]
/stddev[indexOut]);
- tf(indexIn,indexOut)/=sqrt(2.0*M_PI);
- tf(indexIn,indexOut)/=stddev[indexOut];
- norm[indexOut]+=tf(indexIn,indexOut);
+ tf[indexIn][indexOut]/=sqrt(2.0*M_PI);
+ tf[indexIn][indexOut]/=stddev[indexOut];
+ norm[indexOut]+=tf[indexIn][indexOut];
}
}
@@ -362,7 +354,7 @@ template<class T> void Filter::applyFwhm(const vector<double> &wavelengthIn, con
stat.interpolateUp(wavelengthIn,inputValues,wavelength_fine,interpolationType,input_fine,verbose);
output[indexOut][(isample+down-1)/down]=0;
for(int indexIn=0;indexIn<nbandIn;++indexIn){
- output[indexOut][(isample+down-1)/down]+=input_fine[indexIn]*tf(indexIn,indexOut)/norm[indexOut];
+ output[indexOut][(isample+down-1)/down]+=input_fine[indexIn]*tf[indexIn][indexOut]/norm[indexOut];
}
}
}
diff --git a/src/algorithms/Makefile.am b/src/algorithms/Makefile.am
index d9b35a9..b3f9028 100644
--- a/src/algorithms/Makefile.am
+++ b/src/algorithms/Makefile.am
@@ -37,4 +37,4 @@ libalgorithms_la_SOURCES = $(libalgorithms_la_HEADERS) Egcs.cc Filter2d.cc Filte
# list of sources for the binaries
#pktestStat_SOURCES = pktestStat.cc
-#pktestStat_LDADD = $(GSL_LIBS) -lgslwrap $(top_builddir)/src/algorithms/libalgorithms.la $(top_builddir)/src/imageclasses/libimageClasses.la -lgslwrap
+#pktestStat_LDADD = $(GSL_LIBS) $(top_builddir)/src/algorithms/libalgorithms.la $(top_builddir)/src/imageclasses/libimageClasses.la
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index b64c749..f80400e 100644
--- a/src/apps/Makefile.am
+++ b/src/apps/Makefile.am
@@ -29,8 +29,9 @@ pklas2img_LDADD = -llas $(AM_LDFLAGS)
endif
if USE_NLOPT
bin_PROGRAMS += pksensormodel pkopt_svm
+pksensormodel_LDFLAGS = -O1
pksensormodel_SOURCES = $(top_srcdir)/src/models/SensorModel.h pksensormodel.h pksensormodel.cc
-pksensormodel_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lnlopt -lm -lgslwrap
+pksensormodel_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lnlopt -lm -larmadillo #-lgslwrap
pkopt_svm_SOURCES = $(top_srcdir)/src/algorithms/OptFactory.h pkclassify_nn.h pkopt_svm.cc
pkopt_svm_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lnlopt
endif
@@ -53,7 +54,7 @@ pkegcs_LDADD = -lgdal $(AM_LDFLAGS) -lgdal
pkextract_SOURCES = pkextract.cc
pkfillnodata_SOURCES = pkfillnodata.cc
pkfilter_SOURCES = pkfilter.cc
-pkfilter_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgslwrap -lgdal
+pkfilter_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal #-lgslwrap
pkdsm2shadow_SOURCES = pkdsm2shadow.cc
pkmosaic_SOURCES = pkmosaic.cc
pkndvi_SOURCES = pkndvi.cc
diff --git a/src/apps/pksensormodel.cc b/src/apps/pksensormodel.cc
index a63a589..24594c7 100644
--- a/src/apps/pksensormodel.cc
+++ b/src/apps/pksensormodel.cc
@@ -32,10 +32,10 @@ double objFunction(const std::vector<double> &x, std::vector<double> &grad, void
assert(grad.empty());
double error=0;
DataModel *dm=reinterpret_cast<DataModel*> (my_func_data);
- gsl::vector bc_att(3);
- bc_att[0]=x[0];
- bc_att[1]=x[1];
- bc_att[2]=x[2];
+ arma::vec bc_att(3);
+ bc_att(0)=x[0];
+ bc_att(1)=x[1];
+ bc_att(2)=x[2];
dm->setBoresightAtt(bc_att);
for(unsigned int ipoint=0;ipoint<dm->getSize();++ipoint){
double e=dm->getDistGeo(ipoint);
@@ -190,8 +190,8 @@ int main(int argc, char *argv[])
theModel.setPPy(ppy_opt[0]);
theModel.setPolynome(polynome_opt);
theModel.setDatum(datum_opt[0]);
- // gsl::vector bc_pos(3);
- gsl::vector bc_att(3);
+ // arma::vec bc_pos(3);
+ arma::vec bc_att(3);
// while(bcpos_opt.size()<3)
// bcpos_opt.push_back(bcpos_opt[0]);
while(bcatt_opt.size()<3)
@@ -199,16 +199,16 @@ int main(int argc, char *argv[])
if(bcrad_opt[0]){
// bc_pos[0]=theModel.rad2deg(bcpos_opt[0]);
// bc_pos[1]=theModel.rad2deg(bcpos_opt[1]);
- bc_att[0]=theModel.rad2deg(bcatt_opt[0]);
- bc_att[1]=theModel.rad2deg(bcatt_opt[1]);
- bc_att[2]=theModel.rad2deg(bcatt_opt[2]);
+ bc_att(0)=theModel.rad2deg(bcatt_opt[0]);
+ bc_att(1)=theModel.rad2deg(bcatt_opt[1]);
+ bc_att(2)=theModel.rad2deg(bcatt_opt[2]);
}
else{
- // bc_pos[0]=bcpos_opt[0];
- // bc_pos[1]=bcpos_opt[1];
- bc_att[0]=bcatt_opt[0];
- bc_att[1]=bcatt_opt[1];
- bc_att[2]=bcatt_opt[2];
+ // bc_pos(0)=bcpos_opt(0);
+ // bc_pos(1)=bcpos_opt(1);
+ bc_att(0)=bcatt_opt[0];
+ bc_att(1)=bcatt_opt[1];
+ bc_att(2)=bcatt_opt[2];
}
// bc_pos[2]=bcpos_opt[2];
// theModel.setBoresightPos(bc_pos);
@@ -246,9 +246,9 @@ int main(int argc, char *argv[])
istringstream csvstream(csvRecord);
string item;
int icol=0;//first column is id
- gsl::vector thePosGCP(3);
- gsl::vector thePosPlatform(3);
- gsl::vector theAttPlatform(3);
+ arma::vec thePosGCP(3);
+ arma::vec thePosPlatform(3);
+ arma::vec theAttPlatform(3);
while(getline(csvstream,item,fs_opt[0])){//read a column
if(verbose_opt[0]>1)
std::cout << item << " ";
@@ -268,59 +268,59 @@ int main(int argc, char *argv[])
case(3)://X
dvalue=atof(item.c_str());
if(gcprad_opt[0])
- thePosGCP[0]=theModel.rad2deg(dvalue);
+ thePosGCP(0)=theModel.rad2deg(dvalue);
else
- thePosGCP[0]=dvalue;
+ thePosGCP(0)=dvalue;
break;
case(4)://Y
dvalue=atof(item.c_str());
if(gcprad_opt[0])
- thePosGCP[1]=theModel.rad2deg(dvalue);
+ thePosGCP(1)=theModel.rad2deg(dvalue);
else
- thePosGCP[1]=dvalue;
+ thePosGCP(1)=dvalue;
break;
case(5)://Z
dvalue=atof(item.c_str());
- thePosGCP[2]=dvalue;
+ thePosGCP(2)=dvalue;
break;
case(6)://XL
dvalue=atof(item.c_str());
if(pplrad_opt[0])
- thePosPlatform[0]=theModel.rad2deg(dvalue);
+ thePosPlatform(0)=theModel.rad2deg(dvalue);
else
- thePosPlatform[0]=dvalue;
+ thePosPlatform(0)=dvalue;
break;
case(7)://YL
dvalue=atof(item.c_str());
if(pplrad_opt[0])
- thePosPlatform[1]=theModel.rad2deg(dvalue);
+ thePosPlatform(1)=theModel.rad2deg(dvalue);
else
- thePosPlatform[1]=dvalue;
+ thePosPlatform(1)=dvalue;
break;
case(8)://ZL
dvalue=atof(item.c_str());
- thePosPlatform[2]=dvalue;
+ thePosPlatform(2)=dvalue;
break;
case(9)://Roll
dvalue=atof(item.c_str());
if(aplrad_opt[0])
- theAttPlatform[0]=theModel.rad2deg(dvalue);
+ theAttPlatform(0)=theModel.rad2deg(dvalue);
else
- theAttPlatform[0]=dvalue;
+ theAttPlatform(0)=dvalue;
break;
case(10)://Pitch
dvalue=atof(item.c_str());
if(aplrad_opt[0])
- theAttPlatform[1]=theModel.rad2deg(dvalue);
+ theAttPlatform(1)=theModel.rad2deg(dvalue);
else
- theAttPlatform[1]=dvalue;
+ theAttPlatform(1)=dvalue;
break;
case(11)://Yaw
dvalue=atof(item.c_str());
if(aplrad_opt[0])
- theAttPlatform[2]=theModel.rad2deg(dvalue);
+ theAttPlatform(2)=theModel.rad2deg(dvalue);
else
- theAttPlatform[2]=dvalue;
+ theAttPlatform(2)=dvalue;
break;
case(12)://track
break;
@@ -332,10 +332,10 @@ int main(int argc, char *argv[])
}
if(s_srs_opt[0]!=4326){
if(verbose_opt[0]>1)
- std::cout << "transforming: " << thePosGCP[0] << ", " << thePosGCP[1] << ", " << thePosGCP[2] << std::endl;
- poCT->Transform(1,&(thePosGCP[0]),&(thePosGCP[1]),&(thePosGCP[2]));
+ std::cout << "transforming: " << thePosGCP(0) << ", " << thePosGCP(1) << ", " << thePosGCP(2) << std::endl;
+ poCT->Transform(1,&(thePosGCP(0)),&(thePosGCP(1)),&(thePosGCP(2)));
if(verbose_opt[0]>1)
- std::cout << "into: " << thePosGCP[0] << ", " << thePosGCP[1] << ", " << thePosGCP[2] << std::endl;
+ std::cout << "into: " << thePosGCP(0) << ", " << thePosGCP(1) << ", " << thePosGCP(2) << std::endl;
// poCT->Transform(1,&(thePosPlatform[0]),&(thePosPlatform[1]),&(thePosPlatform[2]));
}
theDataModel.pushPosGCP(thePosGCP);
@@ -354,9 +354,9 @@ int main(int argc, char *argv[])
istringstream lineStream(spaceRecord);
string item;
int icol=0;
- gsl::vector thePosGCP(3);
- gsl::vector thePosPlatform(3);
- gsl::vector theAttPlatform(3);
+ arma::vec thePosGCP(3);
+ arma::vec thePosPlatform(3);
+ arma::vec theAttPlatform(3);
while(lineStream >> item){
if(verbose_opt[0]>1)
std::cout << item << " ";
@@ -377,59 +377,59 @@ int main(int argc, char *argv[])
case(3)://X
itemStream >> dvalue;
if(gcprad_opt[0])
- thePosGCP[0]=theModel.rad2deg(dvalue);
+ thePosGCP(0)=theModel.rad2deg(dvalue);
else
- thePosGCP[0]=dvalue;
+ thePosGCP(0)=dvalue;
break;
case(4)://Y
itemStream >> dvalue;
if(gcprad_opt[0])
- thePosGCP[1]=theModel.rad2deg(dvalue);
+ thePosGCP(1)=theModel.rad2deg(dvalue);
else
- thePosGCP[1]=dvalue;
+ thePosGCP(1)=dvalue;
break;
case(5)://Z
itemStream >> dvalue;
- thePosGCP[2]=dvalue;
+ thePosGCP(2)=dvalue;
break;
case(6)://XL
itemStream >> dvalue;
if(pplrad_opt[0])
- thePosPlatform[0]=theModel.rad2deg(dvalue);
+ thePosPlatform(0)=theModel.rad2deg(dvalue);
else
- thePosPlatform[0]=dvalue;
+ thePosPlatform(0)=dvalue;
break;
case(7)://YL
itemStream >> dvalue;
if(pplrad_opt[0])
- thePosPlatform[1]=theModel.rad2deg(dvalue);
+ thePosPlatform(1)=theModel.rad2deg(dvalue);
else
- thePosPlatform[1]=dvalue;
+ thePosPlatform(1)=dvalue;
break;
case(8)://ZL
itemStream >> dvalue;
- thePosPlatform[2]=dvalue;
+ thePosPlatform(2)=dvalue;
break;
case(9)://Roll
itemStream >> dvalue;
if(aplrad_opt[0])
- theAttPlatform[0]=theModel.rad2deg(dvalue);
+ theAttPlatform(0)=theModel.rad2deg(dvalue);
else
- theAttPlatform[0]=dvalue;
+ theAttPlatform(0)=dvalue;
break;
case(10)://Pitch
itemStream >> dvalue;
if(aplrad_opt[0])
- theAttPlatform[1]=theModel.rad2deg(dvalue);
+ theAttPlatform(1)=theModel.rad2deg(dvalue);
else
- theAttPlatform[1]=dvalue;
+ theAttPlatform(1)=dvalue;
break;
case(11)://Yaw
itemStream >> dvalue;
if(aplrad_opt[0])
- theAttPlatform[2]=theModel.rad2deg(dvalue);
+ theAttPlatform(2)=theModel.rad2deg(dvalue);
else
- theAttPlatform[2]=dvalue;
+ theAttPlatform(2)=dvalue;
break;
case(12)://track
break;
@@ -441,10 +441,10 @@ int main(int argc, char *argv[])
}
if(s_srs_opt[0]!=4326){
if(verbose_opt[0]>1)
- std::cout << "transforming: " << thePosGCP[0] << ", " << thePosGCP[1] << ", " << thePosGCP[2] << std::endl;
- poCT->Transform(1,&(thePosGCP[0]),&(thePosGCP[1]),&(thePosGCP[2]));
+ std::cout << "transforming: " << thePosGCP(0) << ", " << thePosGCP(1) << ", " << thePosGCP(2) << std::endl;
+ poCT->Transform(1,&(thePosGCP(0)),&(thePosGCP(1)),&(thePosGCP(2)));
if(verbose_opt[0]>1)
- std::cout << "into: " << thePosGCP[0] << ", " << thePosGCP[1] << ", " << thePosGCP[2] << std::endl;
+ std::cout << "into: " << thePosGCP(0) << ", " << thePosGCP(1) << ", " << thePosGCP(2) << std::endl;
// poCT->Transform(1,&(thePosPlatform[0]),&(thePosPlatform[1]),&(thePosPlatform[2]));
}
theDataModel.pushPosGCP(thePosGCP);
@@ -465,8 +465,8 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>1)
std::cout << "ipoint: " << ipoint << std::endl;
vID.push_back(ipoint);
- gsl::vector thePosPlatform(3);
- gsl::vector theAttPlatform(3);
+ arma::vec thePosPlatform(3);
+ arma::vec theAttPlatform(3);
assert(row_opt.size()>ipoint);
assert(col_opt.size()>ipoint);
theDataModel.pushRow(row_opt[ipoint]);
@@ -474,33 +474,33 @@ int main(int argc, char *argv[])
if(z_opt.size()>ipoint){
// assert(x_opt.size()>ipoint);
// assert(y_opt.size()>ipoint);
- gsl::vector thePosGCP(3);
+ arma::vec thePosGCP(3);
if(x_opt.size()>ipoint&&y_opt.size()>ipoint){
if(gcprad_opt[0]){
- thePosGCP[0]=theModel.rad2deg(x_opt[ipoint]);
- thePosGCP[1]=theModel.rad2deg(y_opt[ipoint]);
+ thePosGCP(0)=theModel.rad2deg(x_opt[ipoint]);
+ thePosGCP(1)=theModel.rad2deg(y_opt[ipoint]);
}
else{
- thePosGCP[0]=x_opt[ipoint];
- thePosGCP[1]=y_opt[ipoint];
+ thePosGCP(0)=x_opt[ipoint];
+ thePosGCP(1)=y_opt[ipoint];
}
}
- thePosGCP[2]=z_opt[ipoint];
+ thePosGCP(2)=z_opt[ipoint];
if(s_srs_opt[0]!=4326)
- poCT->Transform(1,&(thePosGCP[0]),&(thePosGCP[1]),&(thePosGCP[2]));
+ poCT->Transform(1,&(thePosGCP(0)),&(thePosGCP(1)),&(thePosGCP(2)));
theDataModel.pushPosGCP(thePosGCP);
}
assert(yl_opt.size()>ipoint);
assert(zl_opt.size()>ipoint);
if(pplrad_opt[0]){
- thePosPlatform[0]=theModel.rad2deg(xl_opt[ipoint]);
- thePosPlatform[1]=theModel.rad2deg(yl_opt[ipoint]);
+ thePosPlatform(0)=theModel.rad2deg(xl_opt[ipoint]);
+ thePosPlatform(1)=theModel.rad2deg(yl_opt[ipoint]);
}
else{
- thePosPlatform[0]=xl_opt[ipoint];
- thePosPlatform[1]=yl_opt[ipoint];
+ thePosPlatform(0)=xl_opt[ipoint];
+ thePosPlatform(1)=yl_opt[ipoint];
}
- thePosPlatform[2]=zl_opt[ipoint];
+ thePosPlatform(2)=zl_opt[ipoint];
// if(s_srs_opt[0]!=4326)
// poCT->Transform(1,&(thePosPlatform[0]),&(thePosPlatform[1]),&(thePosPlatform[2]));
theDataModel.pushPosPlatform(thePosPlatform);
@@ -509,14 +509,14 @@ int main(int argc, char *argv[])
assert(pitch_opt.size()>ipoint);
assert(yaw_opt.size()>ipoint);
if(aplrad_opt[0]){
- theAttPlatform[0]=theModel.rad2deg(roll_opt[ipoint]);
- theAttPlatform[1]=theModel.rad2deg(pitch_opt[ipoint]);
- theAttPlatform[2]=theModel.rad2deg(yaw_opt[ipoint]);
+ theAttPlatform(0)=theModel.rad2deg(roll_opt[ipoint]);
+ theAttPlatform(1)=theModel.rad2deg(pitch_opt[ipoint]);
+ theAttPlatform(2)=theModel.rad2deg(yaw_opt[ipoint]);
}
else{
- theAttPlatform[0]=roll_opt[ipoint];
- theAttPlatform[1]=pitch_opt[ipoint];
- theAttPlatform[2]=yaw_opt[ipoint];
+ theAttPlatform(0)=roll_opt[ipoint];
+ theAttPlatform(1)=pitch_opt[ipoint];
+ theAttPlatform(2)=yaw_opt[ipoint];
}
theDataModel.pushAttPlatform(theAttPlatform);
}
@@ -527,10 +527,10 @@ int main(int argc, char *argv[])
for(int ipoint=0;ipoint<vID.size();++ipoint){
if(verbose_opt[0]>1)
std::cout << "point: " << ipoint << std::endl;
- gsl::vector pos_platform=theDataModel.getPosPlatform(ipoint);
- gsl::vector att_platform=theDataModel.getAttPlatform(ipoint);
- gsl::vector pos_calc=theDataModel.getPos(ipoint);
- gsl::vector pos_gcp=theDataModel.getPosGCP(ipoint);
+ arma::vec pos_platform=theDataModel.getPosPlatform(ipoint);
+ arma::vec att_platform=theDataModel.getAttPlatform(ipoint);
+ arma::vec pos_calc=theDataModel.getPos(ipoint);
+ arma::vec pos_gcp=theDataModel.getPosGCP(ipoint);
ostringstream gcpss;
double e=theDataModel.getModel().getDistGeo(pos_gcp,pos_calc);
if(e>=theDataModel.getThreshold()){
@@ -608,11 +608,11 @@ int main(int argc, char *argv[])
for(int ipoint=0;ipoint<vID.size();++ipoint){
if(verbose_opt[0]>1)
std::cout << "point: " << ipoint << std::endl;
- gsl::vector pos_platform=theDataModel.getPosPlatform(ipoint);
- gsl::vector att_platform=theDataModel.getAttPlatform(ipoint);
+ arma::vec pos_platform=theDataModel.getPosPlatform(ipoint);
+ arma::vec att_platform=theDataModel.getAttPlatform(ipoint);
- gsl::vector pos_calc;
- gsl::vector pos_gcp;
+ arma::vec pos_calc;
+ arma::vec pos_gcp;
if(input_opt.size()||x_opt.size())
pos_gcp=theDataModel.getPosGCP(ipoint);
else
@@ -628,20 +628,20 @@ int main(int argc, char *argv[])
ostringstream gcpss;
gcpss.precision(12);
- poCT->Transform(1,&(pos_gcp[0]),&(pos_gcp[1]),&(pos_gcp[2]));
- gcpss << pos_gcp[0] << " " << pos_gcp[1] << " " << pos_gcp[2] << " ";
+ poCT->Transform(1,&(pos_gcp(0)),&(pos_gcp(1)),&(pos_gcp(2)));
+ gcpss << pos_gcp(0) << " " << pos_gcp(1) << " " << pos_gcp(2) << " ";
gcpss << errorv.back() << " ";
- poCT->Transform(1,&(pos_calc[0]),&(pos_calc[1]),&(pos_calc[2]));
+ poCT->Transform(1,&(pos_calc(0)),&(pos_calc(1)),&(pos_calc(2)));
if(output_opt.size()){
- outputStream << setprecision(12) << vID[ipoint] << " " << theDataModel.getRow(ipoint) << " " << theDataModel.getCol(ipoint) << " " << pos_calc[0] << " " << pos_calc[1] << " " << pos_calc[2] << " " << gcpss.str();
+ outputStream << setprecision(12) << vID[ipoint] << " " << theDataModel.getRow(ipoint) << " " << theDataModel.getCol(ipoint) << " " << pos_calc(0) << " " << pos_calc(1) << " " << pos_calc(2) << " " << gcpss.str();
if(getzenith_opt[0])
outputStream << " " << theDataModel.getModel().getZenith(att_platform,theDataModel.getRow(ipoint),theDataModel.getCol(ipoint));
outputStream << std::endl;
}
else{
- std::cout << setprecision(12) << vID[ipoint] << " " << theDataModel.getRow(ipoint) << " " << theDataModel.getCol(ipoint) << " " << pos_calc[0] << " " << pos_calc[1] << " " << pos_calc[2] << " " << gcpss.str();
+ std::cout << setprecision(12) << vID[ipoint] << " " << theDataModel.getRow(ipoint) << " " << theDataModel.getCol(ipoint) << " " << pos_calc(0) << " " << pos_calc(1) << " " << pos_calc(2) << " " << gcpss.str();
if(getzenith_opt[0])
std::cout << " " << theDataModel.getModel().getZenith(att_platform,theDataModel.getRow(ipoint),theDataModel.getCol(ipoint));
std::cout << std::endl;
diff --git a/src/apps/pksensormodel.h b/src/apps/pksensormodel.h
index fd13ecb..09a80f4 100644
--- a/src/apps/pksensormodel.h
+++ b/src/apps/pksensormodel.h
@@ -20,7 +20,8 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#ifndef _PKSENSORMODEL_H_
#define _PKSENSORMODEL_H_
#include <vector>
-#include <gslwrap/matrix_double.h>
+// #include <gslwrap/matrix_double.h>
+#include <armadillo>
#include "models/SensorModel.h"
double objFunction(const std::vector<double> &x, std::vector<double> &grad, void *my_func_data);
@@ -42,15 +43,15 @@ class DataModel{
m_row.erase(m_row.begin()+index);
m_col.erase(m_col.begin()+index);
};
- int pushAttPlatform(const gsl::vector& atp){m_attPlatform.push_back(atp); return m_attPlatform.size();};
- int pushPosPlatform(const gsl::vector& ppl){m_posPlatform.push_back(ppl); return m_posPlatform.size();};
- int pushPosGCP(const gsl::vector& pgcp){m_posGCP.push_back(pgcp); return m_posGCP.size();};
+ int pushAttPlatform(const arma::vec& atp){m_attPlatform.push_back(atp); return m_attPlatform.size();};
+ int pushPosPlatform(const arma::vec& ppl){m_posPlatform.push_back(ppl); return m_posPlatform.size();};
+ int pushPosGCP(const arma::vec& pgcp){m_posGCP.push_back(pgcp); return m_posGCP.size();};
int pushRow(int r){m_row.push_back(r); return m_row.size();};
int pushCol(int c){m_col.push_back(c); return m_col.size();};
- gsl::vector getPosPlatform(int index) const{assert(index>=0);assert(index<m_posPlatform.size());return(m_posPlatform[index]);};
- gsl::vector getAttPlatform(int index) const{assert(index>=0);assert(index<m_attPlatform.size());return(m_attPlatform[index]);};
- gsl::vector getPosGCP(int index) const{assert(index>=0);assert(index<m_posGCP.size());return(m_posGCP[index]);};
- gsl::vector getPos(int index) const{
+ arma::vec getPosPlatform(int index) const{assert(index>=0);assert(index<m_posPlatform.size());return(m_posPlatform[index]);};
+ arma::vec getAttPlatform(int index) const{assert(index>=0);assert(index<m_attPlatform.size());return(m_attPlatform[index]);};
+ arma::vec getPosGCP(int index) const{assert(index>=0);assert(index<m_posGCP.size());return(m_posGCP[index]);};
+ arma::vec getPos(int index) const{
assert(index>=0);
assert(index<m_posPlatform.size());
assert(index<m_attPlatform.size());
@@ -63,16 +64,16 @@ class DataModel{
int getRow(int index) const{assert(index>=0);assert(index<m_row.size());return(m_row[index]);};
int getCol(int index) const{assert(index>=0);assert(index<m_col.size());return(m_col[index]);};
double getHeight(int index) const{assert(index>=0);assert(index<m_posGCP.size());return(m_posGCP[index][2]);};
- void setBoresightAtt(const gsl::vector& bc_att){
+ void setBoresightAtt(const arma::vec& bc_att){
m_model.setBoresightAtt(bc_att);
// for(int index=0;index<m_attPlatform.size();++index)
// m_attPlatform[index]+=bc_att;
};
private:
SensorModel::SensorModel m_model;
- vector<gsl::vector> m_posPlatform;
- vector<gsl::vector> m_posGCP;
- vector<gsl::vector> m_attPlatform;
+ vector<arma::vec> m_posPlatform;
+ vector<arma::vec> m_posGCP;
+ vector<arma::vec> m_attPlatform;
vector<int> m_row;
vector<int> m_col;
double m_threshold;
diff --git a/src/models/SensorModel.h b/src/models/SensorModel.h
index 6bebf9a..77209e2 100644
--- a/src/models/SensorModel.h
+++ b/src/models/SensorModel.h
@@ -23,16 +23,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <math.h>
#include <vector>
#include <iostream>
-// #include <gsl/gsl_matrix.h>
-// #include <gsl/gsl_vector.h>
-// #include <gsl/gsl_permutation.h>
-// #include <gsl/gsl_linalg.h>
-// #include <gsl/gsl_blas.h>
-// #include <gsl/gsl_rng.h>
-// #include <gsl/gsl_randist.h>
-#include <gslwrap/vector_double.h>
-#include <gslwrap/matrix_double.h>
-#include <gslwrap/matrix_vector_operators.h>
+#include <armadillo>
#include "ogr_spatialref.h"
#ifndef PI
@@ -50,8 +41,8 @@ namespace SensorModel
class SensorModel{
public:
SensorModel(void){
- m_bcpos.resize(3);
- m_bcatt.resize(3);
+ m_bcpos.set_size(3);
+ m_bcatt.set_size(3);
m_bcpos[0]=0;
m_bcpos[1]=0;
m_bcpos[2]=0;
@@ -60,8 +51,8 @@ namespace SensorModel
m_bcatt[2]=0;
};
SensorModel(int theModel) : m_model(theModel){
- m_bcpos.resize(3);
- m_bcatt.resize(3);
+ m_bcpos.set_size(3);
+ m_bcatt.set_size(3);
m_bcpos[0]=0;
m_bcpos[1]=0;
m_bcpos[2]=0;
@@ -84,117 +75,117 @@ namespace SensorModel
double getF() const {return m_fc;};
void setPPx(double ppx){m_ppx=ppx;};
void setPPy(double ppy){m_ppy=ppy;};
- void setBoresightPos(const gsl::vector& bcpos){m_bcpos=bcpos;};
- void setBoresightAtt(const gsl::vector& bcatt){m_bcatt=bcatt;};
- gsl::vector getBoresightPos() const{return m_bcpos;};
- gsl::vector getBoresightAtt() const{return m_bcatt;};
+ void setBoresightPos(const arma::vec& bcpos){m_bcpos=bcpos;};
+ void setBoresightAtt(const arma::vec& bcatt){m_bcatt=bcatt;};
+ arma::vec getBoresightPos() const{return m_bcpos;};
+ arma::vec getBoresightAtt() const{return m_bcatt;};
void setPolynome(const std::vector<double>& polynome) {m_polynome=polynome;};
void setDatum(const std::string& theDatum="WGS84"){m_spatialref.SetWellKnownGeogCS(theDatum.c_str());};
- double getZenith(const gsl::vector& att_platform, int row, int column) const{
- gsl::vector normallevel(3);
- gsl::vector normalplatform(3);
+ double getZenith(const arma::vec& att_platform, int row, int column) const{
+ arma::vec normallevel(3);
+ arma::vec normalplatform(3);
- gsl::vector apl_deg=att_platform;
+ arma::vec apl_deg=att_platform;
apl_deg+=m_bcatt;
- gsl::vector apl_rad(3);
- apl_rad[0]=deg2rad(apl_deg[0]);//roll
- apl_rad[1]=deg2rad(apl_deg[1]);//pitch
- apl_rad[2]=deg2rad(apl_deg[2]);//yaw
+ arma::vec apl_rad(3);
+ apl_rad(0)=deg2rad(apl_deg(0));//roll
+ apl_rad(1)=deg2rad(apl_deg(1));//pitch
+ apl_rad(2)=deg2rad(apl_deg(2));//yaw
if(getModel()==PUSHBROOM){
- gsl::vector scanAngle(2);
+ arma::vec scanAngle(2);
scanAngle=scanAngle_PB(row,column);
- apl_rad[0]+=scanAngle[1];
- apl_rad[1]+=scanAngle[0];
+ apl_rad(0)+=scanAngle(1);
+ apl_rad(1)+=scanAngle(0);
}
else if(getModel()==WHISKBROOM){
- apl_rad[0]+=0;
- apl_rad[1]=tan(scanAngle_WB(column));
+ apl_rad(0)+=0;
+ apl_rad(1)=tan(scanAngle_WB(column));
}
else//not implemented?
assert(0);
- normallevel[0]=0;
- normallevel[1]=0;
- normallevel[2]=-1;
- gsl::matrix rotM(3,3);
- rotM=getRz(apl_rad[2])*getRy(apl_rad[1])*getRx(apl_rad[0]);
+ normallevel(0)=0;
+ normallevel(1)=0;
+ normallevel(2)=-1;
+ arma::Mat<double> rotM(3,3);
+ rotM=getRz(apl_rad(2))*getRy(apl_rad(1))*getRx(apl_rad(0));
normalplatform=rotM*normallevel;
- return rad2deg(acos(-normalplatform[2]));
+ return rad2deg(acos(-normalplatform(2)));
};
- gsl::vector getPos(const gsl::vector& pos_platform, const gsl::vector& att_platform, int row, int column, double elevation) const{
- gsl::vector thePosition(3);
- gsl::vector pos_ellips(3);
- gsl::vector ppl_deg=pos_platform;
- gsl::vector apl_deg=att_platform;
+ arma::vec getPos(const arma::vec& pos_platform, const arma::vec& att_platform, int row, int column, double elevation) const{
+ arma::vec thePosition(3);
+ arma::vec pos_ellips(3);
+ arma::vec ppl_deg=pos_platform;
+ arma::vec apl_deg=att_platform;
ppl_deg+=m_bcpos;
apl_deg+=m_bcatt;
- gsl::vector ppl_rad(3);
- gsl::vector apl_rad(3);
- ppl_rad[0]=deg2rad(ppl_deg[0]);
- ppl_rad[1]=deg2rad(ppl_deg[1]);
- ppl_rad[2]=ppl_deg[2];//add geoid elevation if necessary...
- apl_rad[0]=deg2rad(apl_deg[0]);//roll
- apl_rad[1]=deg2rad(apl_deg[1]);//pitch
- apl_rad[2]=deg2rad(apl_deg[2]);//yaw
- gsl::vector pos_ecef=pECEF(ppl_rad,apl_rad,row,column);
+ arma::vec ppl_rad(3);
+ arma::vec apl_rad(3);
+ ppl_rad(0)=deg2rad(ppl_deg(0));
+ ppl_rad(1)=deg2rad(ppl_deg(1));
+ ppl_rad(2)=ppl_deg(2);//add geoid elevation if necessary...
+ apl_rad(0)=deg2rad(apl_deg(0));//roll
+ apl_rad(1)=deg2rad(apl_deg(1));//pitch
+ apl_rad(2)=deg2rad(apl_deg(2));//yaw
+ arma::vec pos_ecef=pECEF(ppl_rad,apl_rad,row,column);
pos_ellips=ecef2geo(pos_ecef);
- pos_ellips[2]=0;
+ pos_ellips(2)=0;
thePosition=getXYatZ(elevation,ppl_deg,pos_ellips);
return thePosition;
};
- double getDistGeo(const gsl::vector& pos1, const gsl::vector& pos2) const{
- double lon1=pos1[0];
- double lat1=pos1[1];
- double lon2=pos2[0];
- double lat2=pos2[1];
+ double getDistGeo(const arma::vec& pos1, const arma::vec& pos2) const{
+ double lon1=pos1(0);
+ double lat1=pos1(1);
+ double lon2=pos2(0);
+ double lat2=pos2(1);
double result;
//simplified formula (spherical approximation)
// result=2*asin(sqrt(pow(sin((lat1-lat2)/2),2) + cos(lat1)*cos(lat2)*pow(sin((lon1-lon2)/2),2)));
//using loxodromes
- result=(pos1[1]==pos2[1]) ? lox2(deg2rad(pos1[1]),deg2rad(pos1[0]),deg2rad(pos2[0])) : lox1(deg2rad(pos1[1]),deg2rad(pos1[0]),deg2rad(pos2[1]),deg2rad(pos2[0]));
+ result=(pos1(1)==pos2(1)) ? lox2(deg2rad(pos1(1)),deg2rad(pos1(0)),deg2rad(pos2(0))) : lox1(deg2rad(pos1(1)),deg2rad(pos1(0)),deg2rad(pos2(1)),deg2rad(pos2(0)));
return result;
};
- gsl::vector geo2ecef(const gsl::vector& pos_geo) const{
- gsl::vector pos_ecef(3);
- double nu=getNu(pos_geo[1]);
- double f1=(nu+pos_geo[2])*cos(pos_geo[1]);
+ arma::vec geo2ecef(const arma::vec& pos_geo) const{
+ arma::vec pos_ecef(3);
+ double nu=getNu(pos_geo(1));
+ double f1=(nu+pos_geo(2))*cos(pos_geo(1));
double e1=getE1();
- pos_ecef[0]=f1*cos(pos_geo[0]);
- pos_ecef[1]=f1*sin(pos_geo[0]);
- pos_ecef[2]=(nu*(1-e1*e1)+pos_geo[2])*sin(pos_geo[1]);
+ pos_ecef(0)=f1*cos(pos_geo(0));
+ pos_ecef(1)=f1*sin(pos_geo(0));
+ pos_ecef(2)=(nu*(1-e1*e1)+pos_geo(2))*sin(pos_geo(1));
return pos_ecef;
};
- gsl::vector ecef2geo(const gsl::vector& pos_ecef) const{
- gsl::vector pos_geo(3);
- pos_geo[0]=atan2(pos_ecef[1],pos_ecef[0]);
- while(pos_geo[0]>2*PI)
- pos_geo[0]-=2*PI;
- double r=sqrt(pos_ecef[0]*pos_ecef[0]+pos_ecef[1]*pos_ecef[1]);
+ arma::vec ecef2geo(const arma::vec& pos_ecef) const{
+ arma::vec pos_geo(3);
+ pos_geo(0)=atan2(pos_ecef(1),pos_ecef(0));
+ while(pos_geo(0)>2*PI)
+ pos_geo(0)-=2*PI;
+ double r=sqrt(pos_ecef(0)*pos_ecef(0)+pos_ecef(1)*pos_ecef(1));
double f_earth=1.0/m_spatialref.GetInvFlattening();
double e2=f_earth*(2-f_earth);
double Ae=m_spatialref.GetSemiMajor();
- pos_geo[1]=atan(pos_ecef[2]/r);
+ pos_geo(1)=atan(pos_ecef(2)/r);
double c=1;
double lat=1E+30;
int iterations=0;
- while(sqrt((pos_geo[1]-lat)*(pos_geo[1]-lat))>1E-15){
+ while(sqrt((pos_geo(1)-lat)*(pos_geo(1)-lat))>1E-15){
++iterations;
- lat=pos_geo[1];
+ lat=pos_geo(1);
c=1.0/sqrt(1-e2*sin(lat)*sin(lat));
- pos_geo[1]=atan((pos_ecef[2]+Ae*c*e2*sin(lat))/r);
+ pos_geo(1)=atan((pos_ecef(2)+Ae*c*e2*sin(lat))/r);
}
- pos_geo[2]=r/cos(pos_geo[1])-Ae*c;
- pos_geo[0]=rad2deg(pos_geo[0]);
- pos_geo[1]=rad2deg(pos_geo[1]);
+ pos_geo(2)=r/cos(pos_geo(1))-Ae*c;
+ pos_geo(0)=rad2deg(pos_geo(0));
+ pos_geo(1)=rad2deg(pos_geo(1));
return pos_geo;
};
private:
//line function; In this function, point 0 is the platform position, point 1 is the ellipsoid position
- gsl::vector getXYatZ (double inZ, gsl::vector p0, gsl::vector p1) const{
- gsl::vector posatz(3);
- posatz[0]=p0[0]+(p0[0]-p1[0])/(p0[2]-p1[2])*(inZ-p0[2]);
- posatz[1]=p0[1]+(p0[1]-p1[1])/(p0[2]-p1[2])*(inZ-p0[2]);
- posatz[2]=inZ;
+ arma::vec getXYatZ (double inZ, arma::vec p0, arma::vec p1) const{
+ arma::vec posatz(3);
+ posatz(0)=p0(0)+(p0(0)-p1(0))/(p0(2)-p1(2))*(inZ-p0(2));
+ posatz(1)=p0(1)+(p0(1)-p1(1))/(p0(2)-p1(2))*(inZ-p0(2));
+ posatz(2)=inZ;
return posatz;
};
//get First eccentricity of the Earth's surface
@@ -228,68 +219,66 @@ namespace SensorModel
poCT->Transform(1,&lon,&lat);
};
- gsl::vector pECEF(const gsl::vector& pos, const gsl::vector& attitude, int row, int column) const{
- gsl::vector A(3);
- gsl::matrix B(3,3);
- B.set_element(0,0,-sin(pos[1])*cos(pos[0]));
- B.set_element(0,1,-sin(pos[0]));
- B.set_element(0,2,-cos(pos[1])*cos(pos[0]));
- B.set_element(1,0,-sin(pos[1])*sin(pos[0]));
- B.set_element(1,1,cos(pos[0]));
- B.set_element(1,2,-cos(pos[1])*sin(pos[0]));
- B.set_element(2,0,cos(pos[1]));
- B.set_element(2,1,0);
- B.set_element(2,2,-sin(pos[1]));
+ arma::vec pECEF(const arma::vec& pos, const arma::vec& attitude, int row, int column) const{
+ arma::vec A(3);
+ arma::Mat<double> B(3,3);
+ B(0,0)=-sin(pos(1))*cos(pos(0));
+ B(0,1)=-sin(pos(0));
+ B(0,2)=-cos(pos(1))*cos(pos(0));
+ B(1,0)=-sin(pos(1))*sin(pos(0));
+ B(1,1)=cos(pos(0));
+ B(1,2)=-cos(pos(1))*sin(pos(0));
+ B(2,0)=cos(pos(1));
+ B(2,1)=0;
+ B(2,2)=-sin(pos(1));
A=geo2ecef(pos);
- gsl::vector C(3);
- A+=B*SV(row,column,attitude,pos[2]);
+ arma::vec C(3);
+ A+=B*SV(row,column,attitude,pos(2));
return A;
};
- gsl::vector getIV(int row, int column) const{
- gsl::vector iv(3);
- iv[0]=0;
- gsl::vector scanAngle(2);
+ arma::vec getIV(int row, int column) const{
+ arma::vec iv(3);
+ iv(0)=0;
+ arma::vec scanAngle(2);
if(getModel()==PUSHBROOM){
scanAngle=scanAngle_PB(row,column);
- iv[0]=scanAngle[1];
- iv[1]=tan(scanAngle[0]);
+ iv(0)=scanAngle(1);
+ iv(1)=tan(scanAngle(0));
}
else if(getModel()==WHISKBROOM){
- iv[0]=0;
- iv[1]=tan(scanAngle_WB(column));
+ iv(0)=0;
+ iv(1)=tan(scanAngle_WB(column));
}
else//not implemented?
assert(0);
- iv[2]=-1.0;
+ iv(2)=-1.0;
return iv;
};
double scanAngle_WB(int column) const{
return (-m_fov/2.0+column*(m_fov/(m_ncol-1)));
};
- gsl::vector SV(int row, int column, const gsl::vector& attitude, double z) const{
- gsl::vector rv(3);
- gsl::matrix rotM(3,3);
- rotM=getRz(attitude[2])*getRy(attitude[1])*getRx(attitude[0]);
+ arma::vec SV(int row, int column, const arma::vec& attitude, double z) const{
+ arma::vec rv(3);
+ arma::Mat<double> rotM(3,3);
+ rotM=getRz(attitude(2))*getRy(attitude(1))*getRx(attitude(0));
rv=rotM*getIV(row,column);
- double height=rv[2];
- // gsl_blas_dscal(z/height,&rv);
- // return rv;
- rv*=z/height;
+ double height=rv(2);
+ rv=rv*z/height;
return rv;
};
- gsl::vector scanAngle_PB(int row, int column) const{
- gsl::vector alpha(2);
+ arma::vec scanAngle_PB(int row, int column) const{
+ arma::vec alpha(2);
double r=corr_along(column);
double theAx=(column<m_ppx) ? m_dx*0.000001*(m_ppx-(column+corr_across(column))) : m_dx*0.000001*(m_ppx-(column-corr_across(column)));
double theAy=m_dy*0.000001*(m_ppy-r);
- alpha[0]=(getModel()==FRAME) ? atan(theAx/getF()) : atan(-theAx/getF());
- alpha[1]=atan(-theAy/getF());
+ alpha(0)=(getModel()==FRAME) ? atan(theAx/getF()) : atan(-theAx/getF());
+ alpha(1)=atan(-theAy/getF());
return alpha;
};
//geocentric coordinate system is right-handed, orthogogal Cartesian system with its origin at the centre of the earth. The direct Helmert transformation from lat/lon and ellipsoid height to geocentric x,y,z is returned in posX, posY, posZ
//Euclidian distance between points on sphere
//Euclidian distance between two vectors
- double getDist(const gsl::vector& v1, const gsl::vector& v2) const{gsl::vector dv=v1;dv-=v2;return dv.norm2();};
+ double getDist(const arma::vec& v1, const arma::vec& v2) const{arma::vec dv=v1;dv-=v2;return arma::norm(dv,2);};
double lox1(double lat1, double lon1, double lat2, double lon2) const{
double result=(M(lat2)-M(lat1))/cos(Az(lat1,lon1,lat2,lon2));
return result;
@@ -348,43 +337,44 @@ namespace SensorModel
}
return result;
};
- gsl::matrix getRx(double theta) const{
- gsl::matrix Rx(3,3);
- Rx.set_element(0,0,1);
- Rx.set_element(0,1,0);
- Rx.set_element(0,2,0);
- Rx.set_element(1,0,0);
- Rx.set_element(1,1,cos(theta));
- Rx.set_element(1,2,-sin(theta));
- Rx.set_element(2,0,0);
- Rx.set_element(2,1,sin(theta));
- Rx.set_element(2,2,cos(theta));
+ arma::Mat<double> getRx(double theta) const{
+ arma::Mat<double> Rx(3,3);
+ Rx(0,0)=1;
+ Rx(0,1)=0;
+ Rx(0,2)=0;
+ Rx(1,0)=0;
+ Rx(1,1)=cos(theta);
+ Rx(1,2)=-sin(theta);
+ Rx(2,0)=0;
+ Rx(2,1)=sin(theta);
+ Rx(2,2)=cos(theta);
return Rx;
};
- gsl::matrix getRy(double theta) const{
- gsl::matrix Ry(3,3);
- Ry.set_element(0,0,cos(theta));
- Ry.set_element(0,1,0);
- Ry.set_element(0,2,sin(theta));
- Ry.set_element(1,0,0);
- Ry.set_element(1,1,1);
- Ry.set_element(1,2,0);
- Ry.set_element(2,0,-sin(theta));
- Ry.set_element(2,1,0);
- Ry.set_element(2,2,cos(theta));
+ arma::Mat<double> getRy(double theta) const{
+ arma::Mat<double> Ry(3,3);
+ Ry(0,0)=cos(theta);
+ Ry(0,1)=0;
+ Ry(0,2)=sin(theta);
+ Ry(1,0)=0;
+ Ry(1,1)=1;
+ Ry(1,2)=0;
+ Ry(2,0)=-sin(theta);
+ Ry(2,1)=0;
+ Ry(2,2)=cos(theta);
return Ry;
};
- gsl::matrix getRz(double theta) const{
- gsl::matrix Rz(3,3);
- Rz.set_element(0,0,cos(theta));
- Rz.set_element(0,1,-sin(theta));
- Rz.set_element(0,2,0);
- Rz.set_element(1,0,sin(theta));
- Rz.set_element(1,1,cos(theta));
- Rz.set_element(1,2,0);
- Rz.set_element(2,0,0);
- Rz.set_element(2,1,0);
- Rz.set_element(2,2,1);
+
+ arma::Mat<double> getRz(double theta) const{
+ arma::Mat<double> Rz(3,3);
+ Rz(0,0)=cos(theta);
+ Rz(0,1)=-sin(theta);
+ Rz(0,2)=0;
+ Rz(1,0)=sin(theta);
+ Rz(1,1)=cos(theta);
+ Rz(1,2)=0;
+ Rz(2,0)=0;
+ Rz(2,1)=0;
+ Rz(2,2)=1;
return Rz;
};
int m_model;
@@ -401,8 +391,8 @@ namespace SensorModel
int m_nrow;
double m_dx;
double m_dy;
- gsl::vector m_bcpos;
- gsl::vector m_bcatt;
+ arma::vec m_bcpos;
+ arma::vec m_bcatt;
std::vector<double> m_polynome;
};
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list