pktools  2.6.7
Processing Kernel for geospatial data
Vector2d.h
1 /**********************************************************************
2 Vector2d.h: 2-dimensional vector class (inherits from stl vector class)
3 Copyright (C) 2008-2012 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #ifndef _VECTOR2D_H_
21 #define _VECTOR2D_H_
22 
23 #include <vector>
24 #include <list>
25 #include <algorithm>
26 #include <numeric>
27 #include <gsl/gsl_matrix.h>
28 #include "IndexValue.h"
29 #include "algorithms/StatFactory.h"
30 
31 template<class T> class Vector2d: public std::vector<std::vector <T> >
32 {
33 public:
34  Vector2d();
35  Vector2d(const Vector2d<T>& v1);//copy constructor
36  ~Vector2d();
37  Vector2d(int nrow);
38  Vector2d(int nrow, int ncol);
39  Vector2d(int nrow, int ncol, const T& value);
40  Vector2d(const gsl_matrix* gsl_m);
41  void resize(int nrow)
42  {
43  std::vector< std::vector<T> >::resize(nrow);
44  };
45  void resize(int nrow, int ncol);
46  int nRows() const {return this->size();};
47  int nCols() const {if(this->size()) return this->begin()->size(); else return 0;};
48  void selectCol(int col, std::vector<T> &output) const;
49  void selectCol(int col, T* output) const;
50  std::vector<T> selectCol(int col);
51  void selectCols(const std::list<int> &cols, Vector2d<T> &output) const;
52  void setMask(const Vector2d<T> &mask, T msknodata, T nodata=0);
53  void transpose(Vector2d<T> &output) const{
54  output.resize(nCols(),nRows());
55  for(int irow=0;irow<nRows();++irow){
56  for(int icol=0;icol<nCols();++icol){
57  output[icol][irow]=(*this)[irow][icol];
58  }
59  }
60  };
61  void selectCols(const std::list<int> &cols);
62  void sort(Vector2d<T>& output);
63  void scale(const std::vector<double> &scaleVector, const std::vector<double> &offsetVector, Vector2d<T>& scaledOutput);
64  void scale(const T lbound, const T ubound, std::vector<double> &scaleVector, std::vector<double> &offsetVector, Vector2d<T>& scaledOutput);
65  Vector2d<T> operator=(const Vector2d<T>& v1);
66  Vector2d<T> operator+=(const Vector2d<T>& v1);
67 // std::ostream& operator<<(std::ostream& os, const Vector2d<T>& v);
68 // template<class T> std::ostream& operator<<(std::ostream& os, const Vector2d<T>& v);
69  template<class T1> friend std::ostream& operator<<(std::ostream & os, const Vector2d<T1>& v);
70  Vector2d<T> sum(const Vector2d<T>& v1, const Vector2d<T>& v2) const;
71  T mymax(int& x, int& y, double maxValue) const;
72 
73  T sum() const;
74 };
75 
76 template<class T> Vector2d<T>::Vector2d()
77  : std::vector< std::vector<T> >()
78 {
79 }
80 
81 template<class T> Vector2d<T>::~Vector2d()
82 {
83 }
84 
85 //copy constructor
86 template<class T> Vector2d<T>::Vector2d(const Vector2d<T>& v1){
87  this->resize(v1.size());
88  for(int irow=0;irow<v1.size();++irow)
89  this->at(irow)=v1[irow];
90 }
91 
92 template<class T> Vector2d<T> Vector2d<T>::operator=(const Vector2d<T>& v1){
93  //check for assignment to self (of the form v=v)
94  if(this==&v1)
95  return *this;
96  else{
97  this->resize(v1.size());
98  for(int irow=0;irow<v1.size();++irow)
99  this->at(irow)=v1[irow];
100  return *this;
101  }
102 }
103 
104 template<class T> Vector2d<T> Vector2d<T>::operator+=(const Vector2d<T>& v1){
105  assert(v1.nRows()==nRows());
106  assert(v1.nCols()==nCols());
107  for(int irow=0;irow<nRows();++irow)
108  for(int icol=0;icol<nCols();++icol)
109  (*this)[irow][icol]+=v1[irow][icol];
110  return *this;
111 }
112 
113 template<class T> Vector2d<T>::Vector2d(int nrow)
114  : std::vector< std::vector<T> >(nrow)
115 {
116 }
117 
118 template<class T> Vector2d<T>::Vector2d(int nrow, int ncol)
119 // : std::vector< std::vector<T> >(nrow)
120 {
121  this->resize(nrow);
122  for(int irow=0;irow<nrow;++irow){
123  (this->operator[](irow)).resize(ncol);
124 // (*this)[irow].resize(ncol);
125  }
126 }
127 
128 template<class T> Vector2d<T>::Vector2d(int nrow, int ncol, const T& value)
129 {
130  this->resize(nrow);
131  for(int irow=0;irow<nrow;++irow){
132  (this->operator[](irow)).resize(ncol);
133  for(int icol=0;icol<ncol;++icol)
134  (this->operator[](irow))[icol]=value;
135  }
136 }
137 
138 template<class T> Vector2d<T>::Vector2d(const gsl_matrix* gsl_m)
139 {
140  this->resize(gsl_m->size1);
141  for(int irow=0;irow<this->size();++irow){
142  (this->operator[](irow)).resize(gsl_m->size2);
143  for(int icol=0;icol<this->operator[](irow).size();++icol)
144  (this->operator[](irow))[icol]=gsl_matrix_get(gsl_m,irow,icol);
145  }
146 }
147 
148 
149 template<class T> void Vector2d<T>::resize(int nrow, int ncol)
150 {
151  this->std::vector< std::vector<T> >::resize(nrow);
152  for(int irow=0;irow<nrow;++irow){
153  (this->operator[](irow)).resize(ncol);
154  }
155 }
156 
157 template<class T> void Vector2d<T>::selectCols(const std::list<int> &cols, Vector2d<T> &output) const
158 {
159  output.resize(this->size());
160  std::list<int>::const_iterator it;
161  for(int irow=0;irow<this->size();++irow){
162  output[irow].resize(cols.size());
163  it=cols.begin();
164  for(int icol=0;icol<cols.size();++icol)
165  output[irow][icol]=(*this)[irow][*(it++)];
166  }
167 }
168 
169 template<class T> void Vector2d<T>::selectCol(int col, std::vector<T> &output) const
170 {
171  assert(col>=0);
172  assert(col<(*this)[0].size());
173  output.resize(this->size());
174  for(int irow=0;irow<this->size();++irow){
175  output[irow]=(*this)[irow][col];
176  }
177 }
178 
179 template<class T> std::vector<T> Vector2d<T>::selectCol(int col)
180 {
181  assert(col>=0);
182  assert(col<(*this)[0].size());
183  std::vector<T> output(this->size());
184  for(int irow=0;irow<this->size();++irow)
185  output[irow]=(*this)[irow][col];
186  return(output);
187 }
188 
189 template<class T> void Vector2d<T>::selectCol(int col, T* output) const
190 {
191  assert(col>=0);
192  assert(col<(*this)[0].size());
193  for(int irow=0;irow<this->size();++irow){
194  output[irow]=(*this)[irow][col];
195  }
196 }
197 
198 template<class T> void Vector2d<T>::selectCols(const std::list<int> &cols)
199 {
200  for(int irow=0;irow<this->size();++irow)
201  for(int icol=((*this)[irow]).size()-1;icol>=0;--icol)
202  if(find(cols.begin(),cols.end(),icol)==cols.end())
203  (*this)[irow].erase(((*this)[irow]).begin()+icol);
204 }
205 
206 template<class T> void Vector2d<T>::setMask(const Vector2d<T> &mask, T msknodata, T nodata)
207 {
208  assert(mask.nRows()==nRows());
209  assert(mask.nCols()==nCols());
210  for(int irow=0;irow<this->size();++irow)
211  for(int icol=0;icol<((*this)[irow]).size()-1;++icol)
212  if(mask[irow][icol]==msknodata)
213  (*this)[irow][icol]=nodata;
214 }
215 
216 template<class T1> std::ostream& operator<<(std::ostream& os, const Vector2d<T1>& v)
217 {
218  for(int irow=0;irow<v.size();++irow){
219  for(int icol=0;icol<v[irow].size();++icol){
220  os << v[irow][icol] << "\t";
221  }
222  os << std::endl;
223  }
224  return os;
225  // os << theOption.getLongName() << ": ";
226  // for(int index=0;index<theOption.size();++index)
227  // os << type2string<T>(theOption[index]) << " ";
228  // os << std::endl;
229  // return os;
230 }
231 
232 template<class T> void Vector2d<T>::sort(Vector2d<T>& output)
233 {
234  //sort according to first sample (ex. wavelength)
235  int nsample=this->size();//including first sample (ex. wavelength)
236  int nband=(*this)[0].size();
237  std::vector<IndexValue> sortW(nband);
238  for(int ilevel=0;ilevel<nband;++ilevel){
239  IndexValue pv;
240  pv.position=ilevel;
241  pv.value=(*this)[0][ilevel];
242  sortW[ilevel]=pv;
243  }
244  std::sort(sortW.begin(),sortW.end(),Increase_IndexValue());
245  output.resize(nsample);
246  for(int isample=0;isample<nsample;++isample){
247  output[isample].resize(nband);
248  for(int iband=0;iband<nband;++iband)
249  output[isample][iband]=(*this)[isample][sortW[iband].position];
250  }
251 }
252 
253 template<class T> void Vector2d<T>::scale(const std::vector<double> &scaleVector,const std::vector<double> &offsetVector, Vector2d<T>& scaledOutput)
254 {
255  int nsample=this->size();//including first sample (ex. wavelength)
256  int nband=(*this)[0].size();
257  assert(scaleVector.size()==nband);
258  assert(offsetVector.size()==nband);
259  std::vector<T> pixel(nband);
260  scaledOutput.resize(nsample,nband);
261  for(int isample=0;isample<nsample;++isample)
262  for(int iband=0;iband<nband;++iband)
263  scaledOutput[isample][iband]=((*this)[isample][iband])*scaleVector[iband]+offsetVector[iband];
264 }
265 
266 template<class T> void Vector2d<T>::scale(const T lbound, const T ubound, std::vector<double> &scaleVector, std::vector<double> &offsetVector, Vector2d<T>& scaledOutput)
267 {
268  //scale to lbound and ubound
269  int nsample=this->size();//including first sample (ex. wavelength)
270  int nband=(*this)[0].size();
271  scaleVector.resize(nband);
272  offsetVector.resize(nband);
273  std::vector<T> pixel(nsample);
274  T theMin;
275  T theMax;
277  scaledOutput.resize(nsample,nband);
278  for(int iband=0;iband<nband;++iband){
279  pixel=selectCol(iband);
280  stat.minmax(pixel, pixel.begin(), pixel.end(), theMin, theMax);
281  scaleVector[iband]=static_cast<double>(ubound-lbound)/(theMax-theMin);
282  offsetVector[iband]=static_cast<double>(-theMin*scaleVector[iband])-lbound;
283  for(int isample=0;isample<pixel.size();++isample)
284  scaledOutput[isample][iband]=((*this)[isample][iband])*scaleVector[iband]+offsetVector[iband];
285  }
286 }
287 
288 template<class T> Vector2d<T> Vector2d<T>::sum(const Vector2d<T>& v1, const Vector2d<T>& v2) const{
289  Vector2d<T> vsum(v1.size());
290  assert(v1.size()==v2.size());
291  for(int irow=0;irow<v1.size();++irow){
292  assert(v1[irow].size()==v2[irow].size());
293  vsum[irow].resize(v1[irow].size());
294  for(int icol=0;icol<v1.size();++icol)
295  vsum[irow][icol]=v1[irow][icol]+v2[irow][icol];
296  }
297  return vsum;
298 }
299 
300 template<class T> T Vector2d<T>::sum() const{
301  double theSum=0;
302  for(int irow=0;irow<this->size();++irow){
303  for(int icol=0;icol<this->operator[](irow).size();++icol)
304  theSum+=(this->operator[](irow))[icol];
305  }
306  return theSum;
307 }
308 
309 template<class T> T Vector2d<T>::mymax(int& x, int& y, double maxValue) const{
310  //todo: what if this->operator[](0)[0] >=maxValue?
311  // double theMax=(this->operator[](0))[0];
312  double theMax=0;
313  for(int irow=0;irow<this->size();++irow){
314  for(int icol=0;icol<(this->operator[](irow)).size();++icol){
315  double currentValue=(this->operator[](irow))[icol];
316  if(currentValue<maxValue&&currentValue>theMax){
317  assert(theMax<maxValue);
318  y=irow;
319  x=icol;
320  theMax=currentValue;
321  }
322  }
323  }
324  assert(theMax<maxValue);
325  return theMax;
326 }
327 
328 #endif /* _VECTOR2D_H_ */