pktools  2.6.7
Processing Kernel for geospatial data
ImgRegression.cc
1 /**********************************************************************
2 ImgRegression.cc: class to calculate regression between two raster datasets
3 Copyright (C) 2008-2014 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 #include "ImgRegression.h"
21 #include <iostream>
22 
23 using namespace imgregression;
24 
25 ImgRegression::ImgRegression(void)
26 : m_threshold(0), m_down(1)
27 {}
28 
29 ImgRegression::~ImgRegression(void)
30 {}
31 
32 double ImgRegression::getRMSE(ImgReaderGdal& imgReader1, ImgReaderGdal& imgReader2, double& c0, double& c1, unsigned short band1, unsigned short band2, short verbose) const{
33  c0=0;
34  c1=1;
35  int icol1=0,irow1=0;
36  std::vector<double> rowBuffer1(imgReader1.nrOfCol());
37  std::vector<double> rowBuffer2(imgReader2.nrOfCol());
38  std::vector<double> buffer1;
39  std::vector<double> buffer2;
40 
41  srand(time(NULL));
42  for(irow1=0;irow1<imgReader1.nrOfRow();++irow1){
43  if(irow1%m_down)
44  continue;
45  icol1=0;
46  double icol2=0,irow2=0;
47  double geox=0,geoy=0;
48  imgReader1.readData(rowBuffer1,GDT_Float64,irow1,band1);
49  imgReader1.image2geo(icol1,irow1,geox,geoy);
50  imgReader2.geo2image(geox,geoy,icol2,irow2);
51  icol2=static_cast<int>(icol2);
52  irow2=static_cast<int>(irow2);
53  if(irow2<0||irow2>=imgReader2.nrOfRow())
54  continue;
55  imgReader2.readData(rowBuffer2,GDT_Float64,irow2,band2);
56  for(icol1=0;icol1<imgReader1.nrOfCol();++icol1){
57  if(icol1%m_down)
58  continue;
59  if(m_threshold>0){//percentual value
60  double p=static_cast<double>(rand())/(RAND_MAX);
61  p*=100.0;
62  if(p>m_threshold)
63  continue;//do not select for now, go to next column
64  }
65  imgReader1.image2geo(icol1,irow1,geox,geoy);
66  imgReader2.geo2image(geox,geoy,icol2,irow2);
67  if(icol2<0||icol2>=imgReader2.nrOfCol())
68  continue;
69  icol2=static_cast<int>(icol2);
70  irow2=static_cast<int>(irow2);
71  //check for nodata
72  double value1=rowBuffer1[icol1];
73  double value2=rowBuffer2[icol2];
74  if(imgReader1.isNoData(value1)||imgReader2.isNoData(value2))
75  continue;
76 
77  buffer1.push_back(value1);
78  buffer2.push_back(value2);
79  if(verbose>1)
80  std::cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << std::endl;
81  }
82  }
83  double err=0;
84  if(buffer1.size()&&buffer2.size()){
85  statfactory::StatFactory stat;
86  err=stat.linear_regression_err(buffer1,buffer2,c0,c1);
87  }
88  if(verbose)
89  std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with rmse: " << err << std::endl;
90  return err;
91 }
92 
93 double ImgRegression::getR2(ImgReaderGdal& imgReader1, ImgReaderGdal& imgReader2, double& c0, double& c1, unsigned short band1, unsigned short band2, short verbose) const{
94  c0=0;
95  c1=1;
96  int icol1=0,irow1=0;
97  std::vector<double> rowBuffer1(imgReader1.nrOfCol());
98  std::vector<double> rowBuffer2(imgReader2.nrOfCol());
99  std::vector<double> buffer1;
100  std::vector<double> buffer2;
101 
102  srand(time(NULL));
103  for(irow1=0;irow1<imgReader1.nrOfRow();++irow1){
104  if(irow1%m_down)
105  continue;
106  icol1=0;
107  double icol2=0,irow2=0;
108  double geox=0,geoy=0;
109  imgReader1.readData(rowBuffer1,GDT_Float64,irow1,band1);
110  imgReader1.image2geo(icol1,irow1,geox,geoy);
111  imgReader2.geo2image(geox,geoy,icol2,irow2);
112  icol2=static_cast<int>(icol2);
113  irow2=static_cast<int>(irow2);
114  if(irow2<0||irow2>=imgReader2.nrOfRow())
115  continue;
116  imgReader2.readData(rowBuffer2,GDT_Float64,irow2,band2);
117  for(icol1=0;icol1<imgReader1.nrOfCol();++icol1){
118  if(icol1%m_down)
119  continue;
120  if(m_threshold>0){//percentual value
121  double p=static_cast<double>(rand())/(RAND_MAX);
122  p*=100.0;
123  if(p>m_threshold)
124  continue;//do not select for now, go to next column
125  }
126  imgReader1.image2geo(icol1,irow1,geox,geoy);
127  imgReader2.geo2image(geox,geoy,icol2,irow2);
128  if(icol2<0||icol2>=imgReader2.nrOfCol())
129  continue;
130  icol2=static_cast<int>(icol2);
131  irow2=static_cast<int>(irow2);
132  //check for nodata
133  double value1=rowBuffer1[icol1];
134  double value2=rowBuffer2[icol2];
135  if(imgReader1.isNoData(value1)||imgReader2.isNoData(value2))
136  continue;
137 
138  buffer1.push_back(value1);
139  buffer2.push_back(value2);
140  if(verbose>1)
141  std::cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << std::endl;
142  }
143  }
144  double r2=0;
145  if(buffer1.size()&&buffer2.size()){
146  statfactory::StatFactory stat;
147  r2=stat.linear_regression(buffer1,buffer2,c0,c1);
148  }
149  if(verbose)
150  std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r2 << std::endl;
151  return r2;
152 }
153 
154 double ImgRegression::pgetR2(ImgReaderGdal& imgReader1, ImgReaderGdal& imgReader2, double& c0, double& c1, unsigned short band1, unsigned short band2, short verbose) const{
155  c0=0;
156  c1=1;
157  int icol1=0,irow1=0;
158  std::vector<double> rowBuffer1(imgReader1.nrOfCol());
159  std::vector<double> rowBuffer2(imgReader2.nrOfCol());
160  std::vector<double> buffer1;
161  std::vector<double> buffer2;
162 
163  srand(time(NULL));
164  for(irow1=0;irow1<imgReader1.nrOfRow();++irow1){
165  if(irow1%m_down)
166  continue;
167  icol1=0;
168  double icol2=0,irow2=0;
169  double geox=0,geoy=0;
170  imgReader1.readData(rowBuffer1,GDT_Float64,irow1,band1);
171  imgReader1.image2geo(icol1,irow1,geox,geoy);
172  imgReader2.geo2image(geox,geoy,icol2,irow2);
173  icol2=static_cast<int>(icol2);
174  irow2=static_cast<int>(irow2);
175  if(irow2<0||irow2>=imgReader2.nrOfRow())
176  continue;
177  imgReader2.readData(rowBuffer2,GDT_Float64,irow2,band2);
178  for(icol1=0;icol1<imgReader1.nrOfCol();++icol1){
179  if(icol1%m_down)
180  continue;
181  if(m_threshold>0){//percentual value
182  double p=static_cast<double>(rand())/(RAND_MAX);
183  p*=100.0;
184  if(p>m_threshold)
185  continue;//do not select for now, go to next column
186  }
187  imgReader1.image2geo(icol1,irow1,geox,geoy);
188  imgReader2.geo2image(geox,geoy,icol2,irow2);
189  if(icol2<0||icol2>=imgReader2.nrOfCol())
190  continue;
191  icol2=static_cast<int>(icol2);
192  irow2=static_cast<int>(irow2);
193  //check for nodata
194  double value1=rowBuffer1[icol1];
195  double value2=rowBuffer2[icol2];
196  if(imgReader1.isNoData(value1)||imgReader2.isNoData(value2))
197  continue;
198 
199  buffer1.push_back(value1);
200  buffer2.push_back(value2);
201  if(verbose>1)
202  std::cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << std::endl;
203  }
204  }
205  double r=0;
206  if(buffer1.size()&&buffer2.size()){
207  statfactory::StatFactory stat;
208  r=stat.correlation(buffer1,buffer2);
209  // r=stat.gsl_correlation(buffer1,buffer2);
210  double m1=0;
211  double v1=0;
212  double m2=0;
213  double v2=0;
214  stat.meanVar(buffer1,m1,v1);
215  stat.meanVar(buffer2,m2,v2);
216  if(v1>0){
217  if(r>=0)
218  c1=v2/v1;
219  else
220  c1=-v2/v1;
221  }
222  c0=m2-c1*m1;
223  }
224  if(verbose)
225  std::cout << "orthogonal regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r*r << std::endl;
226  return r*r;
227 }
228 
229 double ImgRegression::getRMSE(ImgReaderGdal& imgReader, unsigned short band1, unsigned short band2, double& c0, double& c1, short verbose) const{
230  c0=0;
231  c1=1;
232  int icol=0,irow=0;
233  std::vector<double> rowBuffer1(imgReader.nrOfCol());
234  std::vector<double> rowBuffer2(imgReader.nrOfCol());
235  std::vector<double> buffer1;
236  std::vector<double> buffer2;
237 
238  srand(time(NULL));
239  assert(band1>=0);
240  assert(band1<imgReader.nrOfBand());
241  assert(band2>=0);
242  assert(band2<imgReader.nrOfBand());
243  for(irow=0;irow<imgReader.nrOfRow();++irow){
244  if(irow%m_down)
245  continue;
246  icol=0;
247  imgReader.readData(rowBuffer1,GDT_Float64,irow,band1);
248  imgReader.readData(rowBuffer2,GDT_Float64,irow,band2);
249  for(icol=0;icol<imgReader.nrOfCol();++icol){
250  if(icol%m_down)
251  continue;
252  if(m_threshold>0){//percentual value
253  double p=static_cast<double>(rand())/(RAND_MAX);
254  p*=100.0;
255  if(p>m_threshold)
256  continue;//do not select for now, go to next column
257  }
258  //check for nodata
259  double value1=rowBuffer1[icol];
260  double value2=rowBuffer2[icol];
261  if(imgReader.isNoData(value1)||imgReader.isNoData(value2))
262  continue;
263 
264  buffer1.push_back(value1);
265  buffer2.push_back(value2);
266  if(verbose>1)
267  std::cout << icol << " " << irow << " " << buffer1.back() << " " << buffer2.back() << std::endl;
268  }
269  }
270  double err=0;
271  if(buffer1.size()&&buffer2.size()){
272  statfactory::StatFactory stat;
273  err=stat.linear_regression_err(buffer1,buffer2,c0,c1);
274  }
275  if(verbose)
276  std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with rmse: " << err << std::endl;
277  return err;
278 }
279 
280 double ImgRegression::getR2(ImgReaderGdal& imgReader, unsigned short band1, unsigned short band2, double& c0, double& c1, short verbose) const{
281  c0=0;
282  c1=1;
283  int icol=0,irow=0;
284  std::vector<double> rowBuffer1(imgReader.nrOfCol());
285  std::vector<double> rowBuffer2(imgReader.nrOfCol());
286  std::vector<double> buffer1;
287  std::vector<double> buffer2;
288 
289  srand(time(NULL));
290  assert(band1>=0);
291  assert(band1<imgReader.nrOfBand());
292  assert(band2>=0);
293  assert(band2<imgReader.nrOfBand());
294  for(irow=0;irow<imgReader.nrOfRow();++irow){
295  if(irow%m_down)
296  continue;
297  icol=0;
298  imgReader.readData(rowBuffer1,GDT_Float64,irow,band1);
299  imgReader.readData(rowBuffer2,GDT_Float64,irow,band2);
300  for(icol=0;icol<imgReader.nrOfCol();++icol){
301  if(icol%m_down)
302  continue;
303  if(m_threshold>0){//percentual value
304  double p=static_cast<double>(rand())/(RAND_MAX);
305  p*=100.0;
306  if(p>m_threshold)
307  continue;//do not select for now, go to next column
308  }
309  //check for nodata
310  double value1=rowBuffer1[icol];
311  double value2=rowBuffer2[icol];
312  if(imgReader.isNoData(value1)||imgReader.isNoData(value2))
313  continue;
314 
315  buffer1.push_back(value1);
316  buffer2.push_back(value2);
317  if(verbose>1)
318  std::cout << icol << " " << irow << " " << buffer1.back() << " " << buffer2.back() << std::endl;
319  }
320  }
321  double r2=0;
322  if(buffer1.size()&&buffer2.size()){
323  statfactory::StatFactory stat;
324  r2=stat.linear_regression(buffer1,buffer2,c0,c1);
325  }
326  if(verbose)
327  std::cout << "linear regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r2 << std::endl;
328  return r2;
329 }
330 
331 double ImgRegression::pgetR2(ImgReaderGdal& imgReader, unsigned short band1, unsigned short band2, double& c0, double& c1, short verbose) const{
332  c0=0;
333  c1=1;
334  int icol=0,irow=0;
335  std::vector<double> rowBuffer1(imgReader.nrOfCol());
336  std::vector<double> rowBuffer2(imgReader.nrOfCol());
337  std::vector<double> buffer1;
338  std::vector<double> buffer2;
339 
340  srand(time(NULL));
341  assert(band1>=0);
342  assert(band1<imgReader.nrOfBand());
343  assert(band2>=0);
344  assert(band2<imgReader.nrOfBand());
345  for(irow=0;irow<imgReader.nrOfRow();++irow){
346  if(irow%m_down)
347  continue;
348  icol=0;
349  imgReader.readData(rowBuffer1,GDT_Float64,irow,band1);
350  imgReader.readData(rowBuffer2,GDT_Float64,irow,band2);
351  for(icol=0;icol<imgReader.nrOfCol();++icol){
352  if(icol%m_down)
353  continue;
354  if(m_threshold>0){//percentual value
355  double p=static_cast<double>(rand())/(RAND_MAX);
356  p*=100.0;
357  if(p>m_threshold)
358  continue;//do not select for now, go to next column
359  }
360  //check for nodata
361  double value1=rowBuffer1[icol];
362  double value2=rowBuffer2[icol];
363  if(imgReader.isNoData(value1)||imgReader.isNoData(value2))
364  continue;
365 
366  buffer1.push_back(value1);
367  buffer2.push_back(value2);
368  if(verbose>1)
369  std::cout << icol << " " << irow << " " << buffer1.back() << " " << buffer2.back() << std::endl;
370  }
371  }
372  double r=0;
373  if(buffer1.size()&&buffer2.size()){
374  statfactory::StatFactory stat;
375  r=stat.correlation(buffer1,buffer2);
376  // r=stat.gsl_correlation(buffer1,buffer2);
377  double m1=0;
378  double v1=0;
379  double m2=0;
380  double v2=0;
381  stat.meanVar(buffer1,m1,v1);
382  stat.meanVar(buffer2,m2,v2);
383  if(v1>0){
384  if(r>=0)
385  c1=v2/v1;
386  else
387  c1=-v2/v1;
388  }
389  c0=m2-c1*m1;
390  }
391  if(verbose)
392  std::cout << "orthogonal regression based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with r^2: " << r*r << std::endl;
393  return r*r;
394 }
ImgRasterGdal::geo2image
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
Definition: ImgRasterGdal.cc:387
ImgRasterGdal::isNoData
bool isNoData(double value) const
Check if value is nodata in this dataset.
Definition: ImgRasterGdal.h:136
ImgRasterGdal::nrOfCol
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
ImgReaderGdal::readData
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::image2geo
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y) ...
Definition: ImgRasterGdal.cc:416
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
Generated on Tue Feb 23 2021 08:50:20 for pktools by   doxygen 1.8.13