pktools 2.6.7
Processing Kernel for geospatial data
ImgRegression.cc
1/**********************************************************************
2ImgRegression.cc: class to calculate regression between two raster datasets
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#include "ImgRegression.h"
21#include <iostream>
22
23using namespace imgregression;
24
25ImgRegression::ImgRegression(void)
26: m_threshold(0), m_down(1)
27{}
28
29ImgRegression::~ImgRegression(void)
30{}
31
32double 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()){
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
93double 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()){
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
154double 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()){
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
229double 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()){
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
280double 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()){
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
331double 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()){
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}
int nrOfRow(void) const
Get the number of rows of this dataset.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
bool isNoData(double value) const
Check if value is nodata in this dataset.
int nrOfBand(void) const
Get the number of bands of this dataset.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
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