154 {
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){
182 double p=static_cast<double>(rand())/(RAND_MAX);
183 p*=100.0;
184 if(p>m_threshold)
185 continue;
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
194 double value1=rowBuffer1[icol1];
195 double value2=rowBuffer2[icol2];
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
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}