1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // grid analysis //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // Grid_IMCORR.cpp //
17 // //
18 // Copyright (C) 2012 by //
19 // Magnus Bremer //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. SAGA is free software; you //
25 // can redistribute it and/or modify it under the terms //
26 // of the GNU General Public License as published by the //
27 // Free Software Foundation, either version 2 of the //
28 // License, or (at your option) any later version. //
29 // //
30 // SAGA is distributed in the hope that it will be //
31 // useful, but WITHOUT ANY WARRANTY; without even the //
32 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
33 // PARTICULAR PURPOSE. See the GNU General Public //
34 // License for more details. //
35 // //
36 // You should have received a copy of the GNU General //
37 // Public License along with this program; if not, see //
38 // <http://www.gnu.org/licenses/>. //
39 // //
40 //-------------------------------------------------------//
41 // //
42 // e-mail: magnus.bremer@uibk.ac.at //
43 // //
44 // contact: Magnus Bremer //
45 // Institute of Geography //
46 // University of Innsbruck //
47 // Innrain 52 //
48 // 6020 Innsbruck //
49 // Austria //
50 // //
51 ///////////////////////////////////////////////////////////
52
53 //---------------------------------------------------------
54 /* This implementation is based on the original FORTRAN code
55 available from <http://nsidc.org/data/velmap/imcorr.html>.
56 This implementation extends the functionality by introducing
57 the possibility to derive 3D displacement vectors once DTM
58 datasets are provided.
59 */
60 //---------------------------------------------------------
61
62
63 ///////////////////////////////////////////////////////////
64 // //
65 // //
66 // //
67 ///////////////////////////////////////////////////////////
68
69 //---------------------------------------------------------
70 #include "Grid_IMCORR.h"
71
72
73 ///////////////////////////////////////////////////////////
74 // //
75 // //
76 // //
77 ///////////////////////////////////////////////////////////
78
79 //---------------------------------------------------------
CGrid_IMCORR(void)80 CGrid_IMCORR::CGrid_IMCORR(void)
81 {
82 Set_Name (_TL("IMCORR - Feature Tracking"));
83
84 Set_Author ("Magnus Bremer (c) 2012");
85
86 Set_Description (_TW(
87 "The tool performs an image correlation based on two raster data sets. "
88 "Additionally, two DTMs can be given and used to optain 3D displacement vectors.\n"
89 "This is a SAGA implementation of the standalone IMCORR software provided by the "
90 "National Snow and Ice Data Center in Boulder, Colorado / US."
91 ));
92
93 Add_Reference(
94 "Scambos, T. A., Dutkiewicz, M. J., Wilson, J. C., and R. A. Bindschadler", "1992",
95 "Application of image cross-correlation to the measurement of glacier velocity using satellite image data",
96 "Remote Sensing Environ., 42(3), 177-186."
97 );
98
99 Add_Reference(
100 "Fahnestock, M. A., Scambos, T.A., and R. A. Bindschadler", "1992",
101 "Semi-automated ice velocity determination from satellite imagery",
102 "Eos, 73, 493."
103 );
104
105 Add_Reference(
106 "http://nsidc.org/data/velmap/imcorr.html", SG_T("standalone software and documentation")
107 );
108
109 Parameters.Add_Grid(
110 "", "GRID_1", _TL("Grid 1"),
111 _TL("The first grid to correlate"),
112 PARAMETER_INPUT
113 );
114
115 Parameters.Add_Grid(
116 "", "GRID_2" , _TL("Grid 2"),
117 _TL("The second grid to correlate"),
118 PARAMETER_INPUT
119 );
120
121 Parameters.Add_Grid(
122 "", "DTM_1", _TL("DTM 1"),
123 _TL("The first DTM used to assign height information to grid 1"),
124 PARAMETER_INPUT_OPTIONAL
125 );
126
127 Parameters.Add_Grid(
128 "", "DTM_2" , _TL("DTM 2"),
129 _TL("The second DTM used to assign height information to grid 2"),
130 PARAMETER_INPUT_OPTIONAL
131 );
132
133 Parameters.Add_Shapes(
134 "", "CORRPOINTS", _TL("Correlated Points"),
135 _TL("Correlated points with displacement and correlation information"),
136 PARAMETER_OUTPUT, SHAPE_TYPE_Point
137 );
138
139 Parameters.Add_Shapes(
140 "", "CORRLINES" , _TL("Displacement Vector"),
141 _TL("Displacement vectors between correlated points"),
142 PARAMETER_OUTPUT, SHAPE_TYPE_Line
143 );
144
145 Parameters.Add_Choice(
146 "", "SEARCH_CHIPSIZE", _TL("Search Chip Size (Cells)"),
147 _TL("Chip size of search chip, used to find correlating reference chip"),
148 CSG_String::Format("%s|%s|%s|%s|%s",
149 SG_T("16x16"),
150 SG_T("32x32"),
151 SG_T("64x64"),
152 SG_T("128x128"),
153 SG_T("256x256")
154 ), 2
155 );
156
157 Parameters.Add_Choice(
158 "", "REF_CHIPSIZE" , _TL("Reference Chip Size (Cells)"),
159 _TL("Chip size of reference chip to be found in search chip"),
160 CSG_String::Format("%s|%s|%s|%s",
161 SG_T("16x16"),
162 SG_T("32x32"),
163 SG_T("64x64"),
164 SG_T("128x128")
165 ), 1
166 );
167
168 Parameters.Add_Double(
169 "", "GRID_SPACING" , _TL("Grid Spacing (Map Units)"),
170 _TL("Grid spacing used for the construction of correlated points [map units]"),
171 10.0, 0.1, true, 256.0, true
172 );
173 }
174
175
176 ///////////////////////////////////////////////////////////
177 // //
178 ///////////////////////////////////////////////////////////
179
180 //---------------------------------------------------------
On_Execute(void)181 bool CGrid_IMCORR::On_Execute(void)
182 {
183 CSG_Grid *pGrid1, *pGrid2, *pDTM1, *pDTM2;
184 CSG_Shapes *pCorrPts, *pCorrLines;
185 CSG_Shape *pCorrPt, *pCorrLine;
186 int Search_Chipsize, Ref_Chipsize, Grid_Spacing;
187 double SpacingMetrics;
188
189 pGrid1 = Parameters("GRID_1")->asGrid();
190 pGrid2 = Parameters("GRID_2")->asGrid();
191 pDTM1 = Parameters("DTM_1")->asGrid();
192 pDTM2 = Parameters("DTM_2")->asGrid();
193 pCorrPts = Parameters("CORRPOINTS")->asShapes();
194 pCorrLines = Parameters("CORRLINES")->asShapes();
195 Search_Chipsize = Parameters("SEARCH_CHIPSIZE")->asInt();
196 Ref_Chipsize = Parameters("REF_CHIPSIZE")->asInt();
197 SpacingMetrics = Parameters("GRID_SPACING")->asDouble();
198
199 Search_Chipsize = (int)(pow(2.0,4+Search_Chipsize));
200 Ref_Chipsize = (int)(pow(2.0,4+Ref_Chipsize));
201 if (Search_Chipsize < Ref_Chipsize)
202 Search_Chipsize = Ref_Chipsize;
203
204 CSG_String Message = CSG_String::Format("%s: %d", _TL("Search chip size"), Search_Chipsize);
205 SG_UI_Msg_Add(Message,true);
206 Message = CSG_String::Format("%s: %d", _TL("Reference chip size"), Ref_Chipsize);
207 SG_UI_Msg_Add(Message,true);
208
209
210 if (pDTM1 == NULL || pDTM2 == NULL)
211 {
212 CSG_String name = CSG_String::Format("%s_CORRPOINTS", pGrid1->Get_Name());
213 pCorrPts->Create(SHAPE_TYPE_Point, name.c_str(), pCorrPts);
214 pCorrPts->Add_Field("ID" , SG_DATATYPE_Int);
215 pCorrPts->Add_Field("GX" , SG_DATATYPE_Double);
216 pCorrPts->Add_Field("GY" , SG_DATATYPE_Double);
217 pCorrPts->Add_Field("REALX" , SG_DATATYPE_Double);
218 pCorrPts->Add_Field("REALY" , SG_DATATYPE_Double);
219 pCorrPts->Add_Field("DISP" , SG_DATATYPE_Double);
220 pCorrPts->Add_Field("STRENGTH" , SG_DATATYPE_Double);
221 pCorrPts->Add_Field("FLAG" , SG_DATATYPE_Int);
222 pCorrPts->Add_Field("XDISP" , SG_DATATYPE_Double);
223 pCorrPts->Add_Field("YDISP" , SG_DATATYPE_Double);
224 pCorrPts->Add_Field("XDISP_UNIT" , SG_DATATYPE_Double);
225 pCorrPts->Add_Field("YDISP_UNIT" , SG_DATATYPE_Double);
226 pCorrPts->Add_Field("XTARG" , SG_DATATYPE_Double);
227 pCorrPts->Add_Field("YTARG" , SG_DATATYPE_Double);
228 pCorrPts->Add_Field("XERR" , SG_DATATYPE_Double);
229 pCorrPts->Add_Field("YERR" , SG_DATATYPE_Double);
230 pCorrPts->Add_Field("ASPECT" , SG_DATATYPE_Double);
231
232 CSG_String name2 = CSG_String::Format("%s_DISP_VEC", pGrid1->Get_Name());
233 pCorrLines->Create(SHAPE_TYPE_Line, name2.c_str(), pCorrLines);
234 pCorrLines->Add_Field("ID" , SG_DATATYPE_Int);
235 pCorrLines->Add_Field("GX" , SG_DATATYPE_Double);
236 pCorrLines->Add_Field("GY" , SG_DATATYPE_Double);
237 pCorrLines->Add_Field("REALX" , SG_DATATYPE_Double);
238 pCorrLines->Add_Field("REALY" , SG_DATATYPE_Double);
239 pCorrLines->Add_Field("DISP" , SG_DATATYPE_Double);
240 pCorrLines->Add_Field("STRENGTH" , SG_DATATYPE_Double);
241 pCorrLines->Add_Field("FLAG" , SG_DATATYPE_Int);
242 pCorrLines->Add_Field("XDISP" , SG_DATATYPE_Double);
243 pCorrLines->Add_Field("YDISP" , SG_DATATYPE_Double);
244 pCorrLines->Add_Field("XDISP_UNIT", SG_DATATYPE_Double);
245 pCorrLines->Add_Field("YDISP_UNIT", SG_DATATYPE_Double);
246 pCorrLines->Add_Field("XTARG" , SG_DATATYPE_Double);
247 pCorrLines->Add_Field("YTARG" , SG_DATATYPE_Double);
248 pCorrLines->Add_Field("XERR" , SG_DATATYPE_Double);
249 pCorrLines->Add_Field("YERR" , SG_DATATYPE_Double);
250 pCorrLines->Add_Field("ASPECT" , SG_DATATYPE_Double);
251 }
252 else // If DTM is given the Output gets 3D
253 {
254 CSG_String name = CSG_String::Format("%s_CORRPOINTS", pGrid1->Get_Name());
255 pCorrPts->Create(SHAPE_TYPE_Point, name.c_str(), pCorrPts, SG_VERTEX_TYPE_XYZ);
256 pCorrPts->Add_Field("ID" , SG_DATATYPE_Int);
257 pCorrPts->Add_Field("GX" , SG_DATATYPE_Double);
258 pCorrPts->Add_Field("GY" , SG_DATATYPE_Double);
259 pCorrPts->Add_Field("REALX" , SG_DATATYPE_Double);
260 pCorrPts->Add_Field("REALY" , SG_DATATYPE_Double);
261 pCorrPts->Add_Field("REALZ" , SG_DATATYPE_Double);
262 pCorrPts->Add_Field("DISP" , SG_DATATYPE_Double);
263 pCorrPts->Add_Field("DISP_REAL" , SG_DATATYPE_Double);
264 pCorrPts->Add_Field("STRENGTH" , SG_DATATYPE_Double);
265 pCorrPts->Add_Field("FLAG" , SG_DATATYPE_Int);
266 pCorrPts->Add_Field("XDISP" , SG_DATATYPE_Double);
267 pCorrPts->Add_Field("YDISP" , SG_DATATYPE_Double);
268 pCorrPts->Add_Field("ZDISP" , SG_DATATYPE_Double);
269 pCorrPts->Add_Field("XDISP_UNIT" , SG_DATATYPE_Double);
270 pCorrPts->Add_Field("YDISP_UNIT" , SG_DATATYPE_Double);
271 pCorrPts->Add_Field("XTARG" , SG_DATATYPE_Double);
272 pCorrPts->Add_Field("YTARG" , SG_DATATYPE_Double);
273 pCorrPts->Add_Field("ZTARG" , SG_DATATYPE_Double);
274 pCorrPts->Add_Field("XERR" , SG_DATATYPE_Double);
275 pCorrPts->Add_Field("YERR" , SG_DATATYPE_Double);
276 pCorrPts->Add_Field("ASPECT" , SG_DATATYPE_Double);
277 pCorrPts->Add_Field("SLOPE" , SG_DATATYPE_Double);
278
279 CSG_String name2 = CSG_String::Format("%s_DISP_VEC", pGrid1->Get_Name());
280 pCorrLines->Create(SHAPE_TYPE_Line, name2.c_str(), pCorrLines, SG_VERTEX_TYPE_XYZ);
281 pCorrLines->Add_Field("ID" , SG_DATATYPE_Int);
282 pCorrLines->Add_Field("GX" , SG_DATATYPE_Double);
283 pCorrLines->Add_Field("GY" , SG_DATATYPE_Double);
284 pCorrLines->Add_Field("REALX" , SG_DATATYPE_Double);
285 pCorrLines->Add_Field("REALY" , SG_DATATYPE_Double);
286 pCorrLines->Add_Field("REALZ" , SG_DATATYPE_Double);
287 pCorrLines->Add_Field("DISP" , SG_DATATYPE_Double);
288 pCorrLines->Add_Field("DISP_REAL" , SG_DATATYPE_Double);
289 pCorrLines->Add_Field("STRENGTH" , SG_DATATYPE_Double);
290 pCorrLines->Add_Field("FLAG" , SG_DATATYPE_Int);
291 pCorrLines->Add_Field("XDISP" , SG_DATATYPE_Double);
292 pCorrLines->Add_Field("YDISP" , SG_DATATYPE_Double);
293 pCorrLines->Add_Field("ZDISP" , SG_DATATYPE_Double);
294 pCorrLines->Add_Field("XDISP_UNIT", SG_DATATYPE_Double);
295 pCorrLines->Add_Field("YDISP_UNIT", SG_DATATYPE_Double);
296 pCorrLines->Add_Field("XTARG" , SG_DATATYPE_Double);
297 pCorrLines->Add_Field("YTARG" , SG_DATATYPE_Double);
298 pCorrLines->Add_Field("ZTARG" , SG_DATATYPE_Double);
299 pCorrLines->Add_Field("XERR" , SG_DATATYPE_Double);
300 pCorrLines->Add_Field("YERR" , SG_DATATYPE_Double);
301 pCorrLines->Add_Field("ASPECT" , SG_DATATYPE_Double);
302 pCorrLines->Add_Field("SLOPE" , SG_DATATYPE_Double);
303 }
304
305
306 // create containers
307 std::vector<std::vector<double> > SearchChip;
308 SearchChip.resize(Search_Chipsize);
309 for(int i = 0; i < Search_Chipsize; i++)
310 SearchChip[i].resize(Search_Chipsize);
311
312 std::vector<std::vector<double> > RefChip;
313 RefChip.resize(Ref_Chipsize);
314 for(int i = 0; i < Ref_Chipsize; i++)
315 RefChip[i].resize(Ref_Chipsize);
316
317 // defaults for values;
318 double mincorr = 2.0;
319 int fitmeth = 1, okparam = 1;
320 double maxdis = -1.0;
321
322 std::vector<double>ioffrq,nomoff;
323 nomoff.push_back(0.0);
324 nomoff.push_back((Search_Chipsize-Ref_Chipsize)*0.5);
325 nomoff.push_back((Search_Chipsize-Ref_Chipsize)*0.5);
326
327 ioffrq.push_back(0.0);
328 ioffrq.push_back(Search_Chipsize/2);
329 ioffrq.push_back(Search_Chipsize/2);
330
331
332 double disp=0.0;
333 double strength=0.0;
334 std::vector<double>best_fit, est_err;
335
336 Grid_Spacing = (int)(SpacingMetrics / (pGrid1->Get_Cellsize()));
337
338
339 // enshures that chips are always in grid
340 int ID = 0;
341 for(int gx1 = (Search_Chipsize/2-1); gx1 < pGrid1->Get_NX()-(Search_Chipsize/2) && Set_Progress(gx1, pGrid1->Get_NX()-(Search_Chipsize/2)); gx1 += Grid_Spacing)
342 {
343 for(int gy1 = (Search_Chipsize/2-1); gy1 < pGrid1->Get_NY()-(Search_Chipsize/2); gy1 += Grid_Spacing)
344 {
345 // get ref_chip
346 Get_This_Chip(RefChip, pGrid1, gx1, gy1, Ref_Chipsize);
347
348 // get search chip
349 Get_This_Chip(SearchChip, pGrid2, gx1, gy1, Search_Chipsize);
350 gcorr(SearchChip, RefChip, mincorr, fitmeth, maxdis, ioffrq, nomoff, okparam, strength, best_fit, est_err, disp);
351
352 if (okparam ==1)
353 {
354 disp = sqrt(best_fit[1]*best_fit[1] + best_fit[2]*best_fit[2]) * Get_Cellsize();
355 double DirNormX = (best_fit[2] * Get_Cellsize()) / disp;
356 double DirNormY = (best_fit[1] * Get_Cellsize()) / disp;
357 double Aspect;
358
359
360 if((DirNormY > 0.0 && DirNormX > 0))
361 Aspect = (atan(fabs(DirNormX)/fabs(DirNormY)))*M_RAD_TO_DEG;
362 else if ((DirNormY <= 0.0 && DirNormX > 0))
363 Aspect = 90+((atan(fabs(DirNormY)/fabs(DirNormX)))*M_RAD_TO_DEG);
364 else if ((DirNormY <= 0.0 && DirNormX <= 0))
365 Aspect = 180+((atan(fabs(DirNormX)/fabs(DirNormY)))*M_RAD_TO_DEG);
366 else
367 Aspect = 270+((atan(fabs(DirNormY)/fabs(DirNormX)))*M_RAD_TO_DEG);
368
369 double xReal = pGrid1->Get_System().Get_xGrid_to_World(gx1);
370 double yReal = pGrid1->Get_System().Get_yGrid_to_World(gy1);
371
372 double xReal2 = xReal + best_fit[2] * Get_Cellsize();
373 double yReal2 = yReal + best_fit[1] * Get_Cellsize();
374
375 if (pDTM1 == NULL || pDTM2 == NULL)
376 {
377 pCorrPt = pCorrPts->Add_Shape();
378 pCorrPt->Add_Point(xReal, yReal);
379 pCorrPt->Set_Value(0, ID);
380 pCorrPt->Set_Value(1, gx1);
381 pCorrPt->Set_Value(2, gy1);
382 pCorrPt->Set_Value(3, xReal);
383 pCorrPt->Set_Value(4, yReal);
384 pCorrPt->Set_Value(5, disp);
385 pCorrPt->Set_Value(6, strength);
386 pCorrPt->Set_Value(7, okparam);
387 pCorrPt->Set_Value(8, best_fit[2] * Get_Cellsize());
388 pCorrPt->Set_Value(9, best_fit[1] * Get_Cellsize());
389 pCorrPt->Set_Value(10, DirNormX);
390 pCorrPt->Set_Value(11, DirNormY);
391 pCorrPt->Set_Value(12, xReal2);
392 pCorrPt->Set_Value(13, yReal2);
393 pCorrPt->Set_Value(14, est_err[2]);
394 pCorrPt->Set_Value(15, est_err[1]);
395 pCorrPt->Set_Value(16, Aspect);
396
397 pCorrLine = pCorrLines->Add_Shape();
398 pCorrLine->Add_Point(xReal, yReal);
399 pCorrLine->Add_Point(xReal2, yReal2);
400 pCorrLine->Set_Value(0, ID);
401 pCorrLine->Set_Value(1, gx1);
402 pCorrLine->Set_Value(2, gy1);
403 pCorrLine->Set_Value(3, xReal);
404 pCorrLine->Set_Value(4, yReal);
405 pCorrLine->Set_Value(5, disp);
406 pCorrLine->Set_Value(6, strength);
407 pCorrLine->Set_Value(7, okparam);
408 pCorrLine->Set_Value(8, best_fit[2] * Get_Cellsize());
409 pCorrLine->Set_Value(9, best_fit[1] * Get_Cellsize());
410 pCorrLine->Set_Value(10, DirNormX);
411 pCorrLine->Set_Value(11, DirNormY);
412 pCorrLine->Set_Value(12, xReal2);
413 pCorrLine->Set_Value(13, yReal2);
414 pCorrLine->Set_Value(14, est_err[2]);
415 pCorrLine->Set_Value(15, est_err[1]);
416 pCorrLine->Set_Value(16, Aspect);
417 }
418 else
419 {
420 double zReal2 = pDTM2->Get_Value(xReal2, yReal2);
421 double zReal = pDTM1->asDouble(gx1, gy1);
422 double Slope = (atan((zReal2-zReal)/fabs(disp)))*M_RAD_TO_DEG;
423 double dispReal = sqrt(pow(zReal2-zReal,2) + disp*disp);
424 pCorrPt = pCorrPts->Add_Shape();
425 pCorrPt->Add_Point(xReal,yReal);
426 pCorrPt->Set_Z(zReal,ID,0);
427 pCorrPt->Set_Value(0, ID);
428 pCorrPt->Set_Value(1, gx1);
429 pCorrPt->Set_Value(2, gy1);
430 pCorrPt->Set_Value(3, xReal);
431 pCorrPt->Set_Value(4, yReal);
432 pCorrPt->Set_Value(5, zReal);
433 pCorrPt->Set_Value(6, disp);
434 pCorrPt->Set_Value(7, dispReal);
435 pCorrPt->Set_Value(8, strength);
436 pCorrPt->Set_Value(9, okparam);
437 pCorrPt->Set_Value(10, best_fit[2] * Get_Cellsize());
438 pCorrPt->Set_Value(11, best_fit[1] * Get_Cellsize());
439 pCorrPt->Set_Value(12, zReal2-zReal);
440 pCorrPt->Set_Value(13, DirNormX);
441 pCorrPt->Set_Value(14, DirNormY);
442 pCorrPt->Set_Value(15, xReal2);
443 pCorrPt->Set_Value(16, yReal2);
444 pCorrPt->Set_Value(17, zReal2);
445 pCorrPt->Set_Value(18, est_err[2]);
446 pCorrPt->Set_Value(19, est_err[1]);
447 pCorrPt->Set_Value(20, Aspect);
448 pCorrPt->Set_Value(21, Slope);
449
450 pCorrLine = pCorrLines->Add_Shape();
451 pCorrLine->Add_Point(xReal, yReal);
452 pCorrLine->Set_Z(zReal,0,0);
453 pCorrLine->Add_Point(xReal2, yReal2);
454 pCorrLine->Set_Z(zReal2,1,0);
455 pCorrLine->Set_Value(0, ID);
456 pCorrLine->Set_Value(1, gx1);
457 pCorrLine->Set_Value(2, gy1);
458 pCorrLine->Set_Value(3, xReal);
459 pCorrLine->Set_Value(4, yReal);
460 pCorrLine->Set_Value(5, zReal);
461 pCorrLine->Set_Value(6, disp);
462 pCorrLine->Set_Value(7, dispReal);
463 pCorrLine->Set_Value(8, strength);
464 pCorrLine->Set_Value(9, okparam);
465 pCorrLine->Set_Value(10, best_fit[2] * Get_Cellsize());
466 pCorrLine->Set_Value(11, best_fit[1] * Get_Cellsize());
467 pCorrLine->Set_Value(12, zReal2-zReal);
468 pCorrLine->Set_Value(13, DirNormX);
469 pCorrLine->Set_Value(14, DirNormY);
470 pCorrLine->Set_Value(15, xReal2);
471 pCorrLine->Set_Value(16, yReal2);
472 pCorrLine->Set_Value(17, zReal2);
473 pCorrLine->Set_Value(18, est_err[2]);
474 pCorrLine->Set_Value(19, est_err[1]);
475 pCorrLine->Set_Value(20, Aspect);
476 pCorrLine->Set_Value(21, Slope);
477 }
478 ID++;
479 }
480 }
481 }
482
483 return( true );
484 }
485
486
487
488 ///////////////////////////////////////////////////////////
489 // //
490 // Methods //
491 // //
492 ///////////////////////////////////////////////////////////
493
494 //---------------------------------------------------------
gcorr(std::vector<std::vector<double>> ChipSearch,std::vector<std::vector<double>> ChipRef,double csmin,int mfit,double ddmx,std::vector<double> ioffrq,std::vector<double> nomoff,int & iacrej,double & streng,std::vector<double> & bfoffs,std::vector<double> & tlerrs,double ddact)495 void CGrid_IMCORR::gcorr(std::vector<std::vector<double> >ChipSearch, std::vector<std::vector<double> >ChipRef, double csmin, int mfit, double ddmx, std::vector<double> ioffrq, std::vector<double> nomoff, int& iacrej, double& streng, std::vector<double>& bfoffs, std::vector<double>& tlerrs, double ddact)
496 {
497 bfoffs.resize(3);
498 std::vector<double>unormc;
499 // compute raw cross-product sums
500 cross(unormc,ChipSearch,ChipRef);
501
502 // compute normalized cross-correlation values and compile statistics
503 std::vector<double>ccnorm;
504 std::vector<double>pkval;
505 std::vector<int>ipkcol;
506 std::vector<int>ipkrow;
507 std::vector<double>sums;
508 gnorm(ccnorm,pkval, ipkcol, ipkrow,sums,ChipSearch, ChipRef,unormc);
509
510 int ncol = (int)(ChipSearch[0].size()-ChipRef[0].size()+1);
511 int nrow = (int)(ChipSearch.size()-ChipRef.size()+1);
512 // Evaluate Strength of correlation peak
513 std::vector<double>cpval;
514 eval(ncol, nrow, ccnorm,pkval, ipkcol, ipkrow,sums,csmin, streng, iacrej,cpval);
515
516 std::vector<double>pkoffs;
517 // determine offsets of peak relative to nominal location
518 if (iacrej==1)
519 {
520 if (mfit != 4)
521 {
522 fitreg(cpval, mfit, pkoffs, tlerrs);
523 bfoffs[1] = (ipkcol[1] - 1) - nomoff[1] + pkoffs[1];
524 bfoffs[2] = (ipkrow[1] - 1) - nomoff[2] + pkoffs[2];
525 }
526 else
527 {
528 bfoffs[1] = (ipkcol[1] - 1) - nomoff[1];
529 bfoffs[2] = (ipkrow[1] - 1) - nomoff[2];
530 tlerrs[1] = 0.5;
531 tlerrs[2] = 0.5;
532 }
533 //
534 // Determine diagonal displacement from nominal and check against maximum
535 // acceptable value
536 ddact = sqrt(bfoffs[1]*bfoffs[1] + bfoffs[2]*bfoffs[2]);
537 if (ddmx > 0.0)
538 {
539 if (ddact > ddmx)
540 iacrej = 5;
541 }
542 else
543 {
544 if ( (bfoffs[1]*bfoffs[1] > ioffrq[1]*ioffrq[1])|| (bfoffs[2]*bfoffs[2] > ioffrq[2]*ioffrq[2]) )
545 iacrej = 5;
546 }
547 }
548 return;
549 }
550
551
552 //---------------------------------------------------------
fitreg(std::vector<double> cpval,int mfit,std::vector<double> & pkoffs,std::vector<double> & tlerrs)553 void CGrid_IMCORR::fitreg(std::vector<double> cpval, int mfit, std::vector<double>& pkoffs, std::vector<double>& tlerrs)
554 {
555 pkoffs.resize(3);
556 tlerrs.resize(3);
557
558 int i, j;
559 std::vector<std::vector<float> >b;
560 std::vector<double> coeffs, vector, wghts, z;
561 coeffs.resize(7);
562 double denom;
563
564 // compute elements of matrix and column vector;
565 sums(cpval, mfit, z, wghts, b, vector);
566
567 // invert matrix to get best-fit peak offsets
568 kvert(b);
569 i = 1;
570 while (i <= 6)
571 {
572 coeffs[i] = 0.0;
573 j = 1;
574 while(j <= 6)
575 {
576 coeffs[i] += (b[i-1][j-1]*vector[j]);
577 j++;
578 }
579 i++;
580 }
581
582 denom = 4.0*coeffs[4]*coeffs[6] - coeffs[5]*coeffs[5];
583 pkoffs[1] = (coeffs[3] * coeffs[5] - 2.0 * coeffs[2] * coeffs[6])/denom;
584 pkoffs[2] = (coeffs[2]*coeffs[5] - 2.0 * coeffs[3] * coeffs[4])/denom;
585
586 // estimate errors in best fit offsets;
587 esterr(z, wghts, b, coeffs, pkoffs, tlerrs);
588
589 return;
590 }
591
592
593 //---------------------------------------------------------
esterr(std::vector<double> z,std::vector<double> wghts,std::vector<std::vector<float>> bnvrs,std::vector<double> coeffs,std::vector<double> & pkoffs,std::vector<double> & tlerrs)594 void CGrid_IMCORR::esterr(std::vector<double> z, std::vector<double> wghts,std::vector<std::vector<float> > bnvrs, std::vector<double> coeffs, std::vector<double>& pkoffs, std::vector<double>& tlerrs)
595 {
596 pkoffs.resize(3);
597 tlerrs.resize(4);
598
599 int i, ivalpt, j;
600 std::vector<double>du, dv;
601 du.resize(7);
602 dv.resize(7);
603 double c, denom, f, usum, var, vsum, x, xsum, y;
604
605 //Compute residual variance for elements of peak array;
606 ivalpt = 1;
607 var = 0.0;
608 y = -2.0;
609 while (y <= 2.0)
610 {
611 x= -2.0;
612 while(x <= 2.0)
613 {
614 f = coeffs[1] + coeffs[2]*x + coeffs[3]*y+ coeffs[4]*x+x + coeffs[5]*x*y + coeffs[6]*y*y;
615 var+= (wghts[ivalpt]*pow((f-z[ivalpt]),2));
616 ivalpt+=1;
617 x+=1.0;
618 }
619 y+=1.0;
620 }
621
622 // compute constants to use the weights
623 c = var/19.0;
624
625 // compute partial derivatives of peak offsets with respect to polynomal coefficients
626 denom = 4.0 * coeffs[4] * coeffs[6] - (coeffs[5] * coeffs[5]);
627 du[1] = 0.0;
628 du[2] = -2.0 * coeffs[6]/denom;
629 du[3] = coeffs[5]/denom;
630 du[4] = -4.0 * coeffs[6] * pkoffs[1]/denom;
631 du[5] = (coeffs[3]+2.0*coeffs[5]*pkoffs[1])/denom;
632 du[6] = (-2.0 * coeffs[2]-4.0*coeffs[4]*pkoffs[1])/denom;
633
634 dv[1] = 0.0;
635 dv[2] = du[3];
636 dv[3] = -2.0 * coeffs[4]/denom;
637 dv[4] = (-2.0*coeffs[3] - 4.0*coeffs[6]*pkoffs[2])/denom;
638 dv[5] = (coeffs[2] + 2.0*coeffs[5] * pkoffs[2])/denom;
639 dv[6] = -4.0 * coeffs[4] * pkoffs[2]/denom;
640
641 // compute estimated errors in offsets
642 usum = 0.0;
643 vsum = 0.0;
644 xsum = 0.0;
645 i=1;
646 while (i <= 6)
647 {
648 j=1;
649 while(j <= 6)
650 {
651 usum += (du[i] * du[j] * bnvrs[i-1][j-1]);
652 vsum += (dv[i] * dv[j] * bnvrs[i-1][j-1]);
653 xsum += (du[i] * dv[j] * bnvrs[i-1][j-1]);
654 j++;
655 }
656 i++;
657 }
658 tlerrs[1] = sqrt(abs(c*usum));
659 tlerrs[2] = sqrt(abs(c*vsum));
660 tlerrs[3] = c*xsum;
661
662 return;
663 }
664
665
666 //---------------------------------------------------------
sums(std::vector<double> cpval,int mfit,std::vector<double> & z,std::vector<double> & wghts,std::vector<std::vector<float>> & b,std::vector<double> & vector)667 void CGrid_IMCORR::sums(std::vector<double> cpval, int mfit, std::vector<double>& z, std::vector<double>& wghts,std::vector<std::vector<float> > & b, std::vector<double>& vector)
668 {
669 b.resize(6);
670 for(int i = 0; i <b.size(); i++)
671 b[i].resize(6);
672
673 vector.resize(26);
674 wghts.resize(26);
675 z.resize(26);
676
677 int i, ic, ir, ivalpt, j;
678
679 std::vector<double>term;
680 term.resize(7);
681
682 double val;
683
684 //initialize arrays and constants
685 i = 1;
686 while (i <= 6)
687 {
688 j=1;
689 while(j <= 6)
690 {
691 b[i-1][j-1] = 0.0;
692 j++;
693 }
694 vector[i] = 0.0;
695 i++;
696 }
697 term[1] = 1.0;
698 ivalpt = 0;
699
700 //compute function of correlation coefficient and assign weight for
701 // each location in neighbourhood of peak
702
703 ir = 1;
704 while(ir <= 5)
705 {
706 ic = 1;
707 while (ic <= 5)
708 {
709 ivalpt++;
710
711 if (cpval[ivalpt] > 1.0)
712 val = cpval[ivalpt];
713 else
714 val = 1.0;
715
716 if (mfit == 1) // eliptical paraboloid
717 {
718 z[ivalpt] = val;
719 wghts[ivalpt] = 1.0;
720 }
721 else if (mfit == 2) // eliptical gaussian
722 {
723 z[ivalpt] = log(val);
724 wghts[ivalpt] = pow(val,2);
725 }
726 else // reciprocal paraboloid
727 {
728 z[ivalpt] = 1.0/val;
729 wghts[ivalpt] = pow(val,4);
730 }
731
732 // compute matrix and vector elements in linear equations for best-fit coefficients
733 term[2] = ic -3;
734 term[3] = ir -3;
735 term[4] = term[2] * term[2];
736 term[5] = term[2] * term[3];
737 term[6] = term[3] * term[3];
738 i = 1;
739 while (i <= 6)
740 {
741 vector[i] += (wghts[ivalpt]*term[i]*z[ivalpt]);
742 j=1;
743 while (j <= 6)
744 {
745 b[i-1][j-1] += (float)(wghts[ivalpt]*term[i]*term[j]);
746 j++;
747 }
748 i++;
749 }
750 ic++;
751 }
752 ir++;
753 }
754
755 return;
756 }
757
758
759 //---------------------------------------------------------
cofact(float num[25][25],float f,std::vector<std::vector<float>> & INV)760 void CGrid_IMCORR::cofact( float num[ 25 ][ 25 ], float f, std::vector<std::vector<float> >& INV )
761 {
762 float b[ 25 ][ 25 ], fac[ 25 ][ 25 ];
763 int p, q, m, n, i, j;
764
765 for ( q = 0;q < f;q++ )
766 {
767 for ( p = 0;p < f;p++ )
768 {
769 m = 0;
770 n = 0;
771
772 for ( i = 0;i < f;i++ )
773 {
774 for ( j = 0;j < f;j++ )
775 {
776 b[ i ][ j ] = 0;
777
778 if ( i != q && j != p )
779 {
780 b[ m ][ n ] = num[ i ][ j ];
781
782 if ( n < ( f-2.0 ) )
783 n++;
784 else
785 {
786 n = 0;
787 m++;
788 }
789 }
790 }
791 }
792
793 fac[ q ][ p ] = (float)(pow( -1.0, (q + p) ) * detrm( b, (f-1) ));
794 }
795 }
796
797 trans( num, fac, f, INV );
798
799 return;
800 }
801
802
803 //---------------------------------------------------------
detrm(float a[25][25],float k)804 float CGrid_IMCORR::detrm( float a[ 25 ][ 25 ], float k )
805 {
806 float s = 1, det = 0, b[ 25 ][ 25 ];
807 int i, j, m, n, c;
808
809 if ( k == 1 )
810 {
811 return ( a[ 0 ][ 0 ] );
812 }
813 else
814 {
815 det = 0;
816
817 for ( c = 0;c < k;c++ )
818 {
819 m = 0;
820 n = 0;
821
822 for ( i = 0;i < k;i++ )
823 {
824 for ( j = 0;j < k;j++ )
825 {
826 b[ i ][ j ] = 0;
827
828 if ( i != 0 && j != c )
829 {
830 b[ m ][ n ] = a[ i ][ j ];
831
832 if ( n < ( k-2 ) )
833 n++;
834 else
835 {
836 n = 0;
837 m++;
838 }
839 }
840 }
841 }
842
843 det = det + s * ( a[ 0 ][ c ] * detrm( b, k-1 ) );
844 s = -1 * s;
845 }
846 }
847
848 return ( det );
849 }
850
851
852 //---------------------------------------------------------
trans(float num[25][25],float fac[25][25],float r,std::vector<std::vector<float>> & INV)853 void CGrid_IMCORR::trans( float num[ 25 ][ 25 ], float fac[ 25 ][ 25 ], float r , std::vector<std::vector<float> >& INV)
854 {
855 int i, j;
856 float b[ 25 ][ 25 ], inv[ 25 ][ 25 ], d;
857
858 for ( i = 0;i < r;i++ )
859 {
860 for ( j = 0;j < r;j++ )
861 {
862 b[ i ][ j ] = fac[ j ][ i ];
863 }
864 }
865
866 d = detrm( num, r );
867 inv[ i ][ j ] = 0;
868
869 for ( i = 0;i < r;i++ )
870 {
871 for ( j = 0;j < r;j++ )
872 {
873 inv[ i ][ j ] = b[ i ][ j ] / d;
874 }
875 }
876
877 INV.resize((int)r);
878 for (int ii = 0; ii < r; ii++)
879 {
880 INV[ii].resize((int)r);
881 }
882 for ( i = 0;i < r;i++ )
883 {
884 for ( j = 0;j < r;j++ )
885 {
886 INV[ i ][ j ] = inv[ i ][ j ];
887 }
888 }
889
890 return;
891 }
892
893
894 //---------------------------------------------------------
kvert(std::vector<std::vector<float>> & V)895 void CGrid_IMCORR::kvert(std::vector<std::vector<float> >& V)
896 {
897 const int k = (int)V[0].size();
898 float b[25][25];
899 for ( int i = 0;i < k;i++ )
900 {
901 for (int j = 0;j < k;j++ )
902 {
903 b[ i ][ j ] = V[ i ][ j ];
904 }
905 }
906 float d = detrm( b, (float)k );
907 if ( d != 0 )
908 cofact( b, (float)k, V);
909
910 return;
911 }
912
913
914 //---------------------------------------------------------
915 // code from original FORTRAN version
916 //void CGrid_IMCORR::kvert(std::vector<std::vector<double> >& V, std::vector<double>W)
917 //{
918 // int LV = V[0].size()-1;
919 // int N = V.size()-1;
920 // int H,I,J,K,L,M,N,O,P,Q;
921 // double S, T;
922 //
923 //
924 // if (N != 1)
925 // {
926 // O = N+1;
927 // L = 0;
928 // M = 1;
929 // if (L != N)
930 // {
931 // K = L;
932 // L = M;
933 // M = M + 1;
934 // // find pivot and start row swap
935 // P = L;
936 // Q = L;
937 // S = abs(V[L][L]);
938 // H = L;
939 // bool GoOut = false;
940 // while (H <= N)
941 // {
942 // I = L;
943 // while(I <= N)
944 // {
945 // T = abs(V[I][H]);
946 // if (T <= S)
947 // break;
948 // P = I;
949 // Q = H;
950 // S = T;
951 // if ( I>N)
952 // break;
953 // I++;
954 // }
955 // H++;
956 // }
957 // W[N+L] = P;
958 // W[O-L] = Q;
959 // I = 1;
960 // while(I <= N)
961 // {
962 // T = V[I][L];
963 // V[I][L] = V[I][Q];
964 // V[I][Q] = T;
965 // I++;
966 // }
967 // S = V[P][L];
968 // V[P][L] = V[L][L];
969 // if (S != 0.0)
970 // {
971 // // compute multipliers
972 // V[L][L] = -1;
973 // S = 1.0/S;
974 // I=1;
975 // while(I <= N)
976 // {
977 // V[I][L] = -S * V[I][L];
978 // I++;
979 // }
980 // J=L;
981 // J++;
982 // if(J > N)
983 // J = 1;
984 //
985 //
986 // }
987 //
988 // }
989 // }
990 // return;
991 //}
992
993
994 //---------------------------------------------------------
eval(int ncol,int nrow,std::vector<double> CCNORM,std::vector<double> pkval,std::vector<int> ipkcol,std::vector<int> ipkrow,std::vector<double> sums,double & csmin,double & streng,int & iacrej,std::vector<double> & cpval)995 void CGrid_IMCORR::eval(int ncol, int nrow, std::vector<double> CCNORM , std::vector<double> pkval, std::vector<int> ipkcol, std::vector<int> ipkrow, std::vector<double> sums, double& csmin, double& streng, int& iacrej, std::vector<double>& cpval)
996 {
997 cpval.resize(26,0.0);// shift to rebuild fortran;
998 std::vector<int>ipt5(3);
999 int i, icol, idist, iptr, ipt7, irow, krbase, lcol, lrow, npts, n5, n7, j;
1000 double bmean, bsigma;
1001 ipt5[0] = 0;
1002 ipt5[1] = 32;
1003 ipt5[2] = 32;
1004 iacrej = 1;
1005 streng = 0.0;
1006 ipt7 =1;
1007
1008 // check for peak value within two rows or columns of edge
1009 if ((ipkcol[1] <= 2) || (ipkcol[1] >= ncol-1) || (ipkrow[1] <= 2) || (ipkrow[1] >= nrow-1))
1010 {
1011 iacrej = 0;
1012 return;
1013 }
1014
1015 // Find largest values outside 5 by 5 and 7 by 7 neoghbourhoods of peak
1016 n5 = 0;
1017 n7 = 0;
1018 i=2;
1019 while ((n5 < 2) && (i <= 32))
1020 {
1021 if (abs(ipkcol[1]-ipkcol[i]) > abs(ipkrow[1]-ipkrow[i]))
1022 idist = abs(ipkcol[1]-ipkcol[i]);
1023 else
1024 idist = abs(ipkrow[1]-ipkrow[i]);
1025 if (idist > 2)
1026 {
1027 n5 = n5 +1;
1028 ipt5[n5] = i;
1029 if ((n7 == 0) && (idist < 3))
1030 {
1031 ipt7 = i;
1032 n7 = 1;
1033 }
1034 }
1035 i++;
1036 }
1037
1038 if((ipt5[1] <= 3) || ipt5[2] <=5)
1039 {
1040 iacrej = 3;
1041 return;
1042 }
1043
1044 //Find edges of 9 by 9 array centred on peak
1045 if (1 > ipkcol[1] -4)
1046 icol = 1;
1047 else
1048 icol = ipkcol[1] -4;
1049
1050 if (1 > ipkcol[1] -4)
1051 irow = 1;
1052 else
1053 irow = ipkcol[1] -4;
1054
1055 if (ncol < ipkcol[1] -4)
1056 lcol = ncol;
1057 else
1058 lcol = ipkcol[1] -4;
1059
1060 if (nrow < ipkcol[1] -4)
1061 lrow = nrow;
1062 else
1063 lrow = ipkcol[1] -4;
1064
1065 // eliminate points within 9 by 9 array from background statistics
1066 krbase = ncol*(nrow-1);
1067 i = irow;
1068 while (i <= lrow)
1069 {
1070 j=icol;
1071 while (j <= lcol)
1072 {
1073 sums[0] -= CCNORM[krbase+j];
1074 sums[1] -= (CCNORM[krbase+j]*CCNORM[krbase+j]);
1075 j++;
1076 }
1077 krbase+= ncol;
1078 i++;
1079 }
1080
1081 //Compute background mean and standard deviation
1082 npts = ncol*nrow - (lcol - icol +1) * (lrow - irow +1);
1083 bmean = sums[0]/npts;
1084 bsigma = sqrt(sums[1]/npts - bmean*bmean);
1085
1086 // Compute correlation strength and check against minimum. Note that
1087 // the LAS 4.x Version of the code contains an error... it often tries to
1088 // access pkval array element zero when there are no significant subsidiary peaks, resulting in
1089 // unpredictable results
1090
1091 if (n7 == 0)
1092 streng = 2 * ((pkval[1] -bmean)/bsigma) -0.2;
1093 else
1094 {
1095 streng = (pkval[1] -bmean)/bsigma + (pkval[0] - pkval[ipt7])/bsigma + 0.2*(n7-1.0);
1096 }
1097
1098 if (streng < csmin)
1099 {
1100 iacrej = 4;
1101 return;
1102 }
1103
1104 //convert 5 by 5 neighbourhood of peak to standard deviations above mean
1105 krbase = ncol*(ipkrow[1]-3);
1106 iptr = 1;
1107 int I = 1;
1108 while (I <= 5)
1109 {
1110 j = ipkcol[1]-2;
1111 while (j <= ipkcol[1]+2)
1112 {
1113 cpval[iptr] = ((CCNORM[krbase+j] - bmean)/bsigma);
1114 iptr++;
1115 j++;
1116 }
1117 krbase+= ncol;
1118 I++;
1119 }
1120
1121 return;
1122 }
1123
1124
1125 //---------------------------------------------------------
gnorm(std::vector<double> & CCNORM,std::vector<double> & pkval,std::vector<int> & ipkcol,std::vector<int> & ipkrow,std::vector<double> & sums,std::vector<std::vector<double>> ChipSearch,std::vector<std::vector<double>> ChipRef,std::vector<double> UNORMC)1126 void CGrid_IMCORR::gnorm(std::vector<double>& CCNORM , std::vector<double>& pkval, std::vector<int>& ipkcol, std::vector<int>& ipkrow, std::vector<double>& sums, std::vector<std::vector<double> >ChipSearch, std::vector<std::vector<double> >ChipRef, std::vector<double> UNORMC )
1127 {
1128 std::vector<double>ser, ref;
1129 //create arrays
1130 ser.push_back(0.0); // dummy index to allow fortran conform indexing
1131 for(int i = 0; i<ChipSearch.size(); i++)
1132 {
1133 for(int ii = 0; ii<ChipSearch[0].size(); ii++)
1134 {
1135 ser.push_back(ChipSearch[i][ii]);
1136 }
1137 }
1138
1139 ref.push_back(0.0); // dummy index to allow fortran conform indexing
1140 for(int i = 0; i<ChipRef.size(); i++)
1141 {
1142 for(int ii = 0; ii<ChipRef[0].size(); ii++)
1143 {
1144 ref.push_back(ChipRef[i][ii]);
1145 }
1146 }
1147
1148
1149 CCNORM.push_back(0.0);
1150
1151 // double vectors
1152 std::vector<double >xval, colsum, colsqr;
1153 xval.resize(33);
1154 colsum.resize(257,0.0);
1155 colsqr.resize(257,0.0);
1156
1157 // integer vectors
1158 std::vector<int>iptr, ixcol, ixrow;
1159 iptr.resize(33);
1160 ixcol.resize(33);
1161 ixrow.resize(33);
1162
1163 // in most cases search is bigger than ref (default = 64x64)
1164 std::vector<int>nsnew;
1165 nsnew.push_back((int)ChipSearch[0].size());
1166 nsnew.push_back((int)ChipSearch.size());
1167
1168 // in most cases ref is smaller than search (default = 32x32)
1169 std::vector<int>nrnew;
1170 nrnew.push_back((int)ChipRef[0].size());
1171 nrnew.push_back((int)ChipRef.size());
1172
1173 int ipfree, ipt, jstart, jstop, k, kol, koladd, kolsub, line, lnadd, lnsub, ncol, nrow, nrtot;
1174 double refsum, refsqr, rho, rmean, rsigma, rtotal, sigma1, sigmas, srchsm, srchsq, temp, tempmn, igl, iglnew, iglold;
1175 refsum = 0;
1176 refsqr = 0;
1177 nrtot = nrnew[0]*nrnew[1];
1178
1179 ipt =1;
1180 while (ipt <= nrtot)
1181 {
1182 igl = ref[ipt];
1183 refsum = refsum + igl;
1184 refsqr = refsqr +igl*igl;
1185 ipt++;
1186 }
1187 rtotal = nrtot;
1188 tempmn = 0.01 * (rtotal*rtotal);
1189 rmean = refsum/rtotal;
1190 if ((refsqr/rtotal - rmean*rmean) > 0.01)
1191 rsigma = sqrt(refsqr/rtotal - rmean*rmean);
1192 else
1193 rsigma = sqrt(0.01);
1194
1195 // clear sums and sums of squares of normalized correlation values
1196 sums.push_back(0.0);
1197 sums.push_back(0.0);
1198 k =1;
1199 while (k <= 32)
1200 {
1201 xval[k] = -1.0;
1202 iptr[k] = k;
1203 k++;
1204 }
1205 ncol = nsnew[0] - nrnew[0] +1;
1206 nrow = nsnew[1] - nrnew[1] +1;
1207
1208 // compute normalized cross-corr values for one row of alignments at a time
1209 jstart = 1;
1210 jstop = ncol;
1211 int i = 1;
1212 while (i <= nrow)
1213 {
1214 // Get column sums and sums of squares for portion of search
1215 // subimage overlain by reference in current row of alignments
1216 if (i == 1)
1217 {
1218 kol = 1;
1219 while(kol <= nsnew[0])
1220 {
1221 colsum[kol] = 0.0;
1222 colsqr[kol] = 0.0;
1223 kol++;
1224 }
1225 ipt = 1;
1226 line = 1;
1227 while (line <= nrnew[1])
1228 {
1229 kol = 1;
1230 while (kol <= nsnew[0])
1231 {
1232 igl = ser[ipt];
1233 colsum[kol]+=igl;
1234 colsqr[kol]+=(igl*igl);
1235 ipt++;
1236 kol++;
1237 }
1238 line++;
1239 }
1240 }
1241 else // if i != 1
1242 {
1243 lnsub = (i-2)*nsnew[0];
1244 lnadd = lnsub + nsnew[0] * nrnew[1];
1245 kol = 1;
1246 while (kol <= nsnew[0])
1247 {
1248 iglnew = ser[lnadd+kol];
1249 iglold = ser[lnsub+kol];
1250 colsum[kol]+= (iglnew-iglold);
1251 colsqr[kol]+= ((iglnew*iglnew)-(iglold*iglold));
1252 kol++;
1253 }
1254
1255 } // end if
1256
1257 // complete computation of search-subarea pixel statistics
1258 int j = jstart;
1259 while(j <= jstop)
1260 {
1261
1262 if(j == jstart)
1263 {
1264 srchsm = 0.0;
1265 srchsq = 0.0;
1266 kol = 1;
1267 while (kol <= nrnew[0])
1268 {
1269 srchsm+= colsum[kol];
1270 srchsq+= colsqr[kol];
1271 kol++;
1272 }
1273
1274 if (tempmn > (rtotal*srchsq - srchsm*srchsm))
1275 temp = tempmn;
1276 else
1277 temp = (rtotal*srchsq - srchsm*srchsm);
1278 sigmas = sqrt(temp);
1279 }
1280 else
1281 {
1282 kolsub = j -jstart;
1283 koladd = kolsub + nrnew[0];
1284 srchsm+= (colsum[koladd] - colsum[kolsub]);
1285 srchsq+= (colsqr[koladd] - colsqr[kolsub]);
1286
1287 if (tempmn > (rtotal*srchsq - srchsm*srchsm))
1288 temp = tempmn;
1289 else
1290 temp = (rtotal*srchsq - srchsm*srchsm);
1291
1292 sigma1 = 0.5 * (sigmas + temp/sigmas);
1293 sigmas = 0.5 *(sigma1 + temp/sigma1);
1294 }
1295
1296 // compute normalized cross-correlation value
1297 rho = (UNORMC[j] - rmean*srchsm)/(rsigma * sigmas);
1298 CCNORM.push_back(rho);
1299 sums[0]+= rho;
1300 sums[1]+= (rho*rho);
1301
1302 // check whether value is among top 32
1303 if (rho > xval[iptr[32]])
1304 {
1305 k = 32;
1306 ipfree = iptr[32];
1307 while ( k > 1 && rho > xval[iptr[k-1]])
1308 {
1309 iptr[k] = iptr[k-1];
1310 k = k-1;
1311 }
1312 iptr[k] = ipfree;
1313 xval[ipfree] = rho;
1314 ixcol[ipfree] = j - jstart + 1;
1315 ixrow[ipfree] = i;
1316 }
1317 j++;
1318 }
1319 jstart+= ncol;
1320 jstop+= ncol;
1321 i++;
1322 }
1323
1324 // copy peak values and coordinates in correct sequence
1325 k= 1;
1326 pkval.push_back(0.0);
1327 ipkcol.push_back(0);
1328 ipkrow.push_back(0);
1329 while (k <= 32)
1330 {
1331 pkval.push_back(xval[iptr[k]]);
1332 ipkcol.push_back(ixcol[iptr[k]]);
1333 ipkrow.push_back(ixrow[iptr[k]]);
1334 k++;
1335 }
1336
1337 return;
1338 }
1339
1340
1341 //---------------------------------------------------------
cross(std::vector<double> & UNORMC,std::vector<std::vector<double>> ChipSearch,std::vector<std::vector<double>> ChipRef)1342 void CGrid_IMCORR::cross(std::vector<double>& UNORMC , std::vector<std::vector<double> >ChipSearch, std::vector<std::vector<double> >ChipRef)
1343 {
1344 // in most cases search is bigger than ref (default = 64x64)
1345 std::vector<int>nsnew;
1346 nsnew.push_back((int)ChipSearch[0].size());
1347 nsnew.push_back((int)ChipSearch.size());
1348
1349 // in most cases ref is smaller than search (default = 32x32)
1350 std::vector<int>nrnew;
1351 nrnew.push_back((int)ChipRef[0].size());
1352 nrnew.push_back((int)ChipRef.size());
1353
1354 std::vector<std::vector<double> > ChipRef2;
1355 ChipRef2.resize(ChipSearch[0].size());
1356 for(int i = 0; i < ChipSearch[0].size(); i++)
1357 ChipRef2[i].resize(ChipSearch.size(),0.0);
1358
1359 // zero extent chipref to search chip size
1360 for(int i = 0; i<ChipRef.size(); i++)
1361 {
1362 for(int ii = 0; ii<ChipRef[0].size(); ii++)
1363 {
1364 ChipRef2[i][ii] = ChipRef[i][ii];
1365 }
1366 }
1367
1368
1369 std::vector<double>ser, ref;
1370 //create arrays
1371 ser.push_back(0.0); // dummy index to allow fortran conform indexing
1372 for(int i = 0; i<ChipSearch.size(); i++)
1373 {
1374 for(int ii = 0; ii<ChipSearch[0].size(); ii++)
1375 {
1376 ser.push_back(ChipSearch[i][ii]);
1377 }
1378 }
1379
1380 ref.push_back(0.0); // dummy index to allow fortran conform indexing
1381 for(int i = 0; i<ChipRef2.size(); i++)
1382 {
1383 for(int ii = 0; ii<ChipRef2[0].size(); ii++)
1384 {
1385 ref.push_back(ChipRef2[i][ii]);
1386 }
1387 }
1388
1389 int lnstrt = 1;
1390 int imgptr = 1;
1391 int line = 1;
1392 int lncont = 0;
1393
1394 // pseudo complex arrays
1395 std::vector<double>cser, cref;
1396 cser.push_back(0.0);
1397 cref.push_back(0.0);
1398 for (int i = 1; i < ser.size(); i++)
1399 {
1400 cser.push_back(ser[i]);
1401 cser.push_back(0.0);
1402 cref.push_back(ref[i]);
1403 cref.push_back(0.0);
1404 }
1405
1406 // take fast fourier transform of search and reference data
1407 fft2(cser,nsnew,1); // from "house" format to frequency format
1408 fft2(cref,nsnew,1);
1409
1410 // make point by multiplication of ft of search ft with conjugate of reference image
1411 for(int i = 1; i<cser.size(); i+=2)
1412 {
1413 // cser[i] = cser[i] * conjg(cref[i])
1414 double temp_cser_real = cser[i];
1415 cser[i] = cser[i]*cref[i] - cser[i+1]*(-cref[i+1]);
1416 cser[i+1] = temp_cser_real*(-cref[i+1]) + cref[i]*(cser[i+1]);
1417 }
1418
1419 // take inverse fft of cser
1420 fft2(cser,nsnew,-1); // real signal format (house) again
1421
1422 //Extract useful valid correlation
1423 int ncol = nsnew[0] - nrnew[0] +1;
1424 int nrow = nsnew[1] - nrnew[1] +1;
1425 int denom = (int)(ChipSearch[0].size() * ChipSearch.size());
1426 int ndxout = 1;
1427 int i =1;
1428
1429 UNORMC.push_back(0.0); // shift for fortran compatibility;
1430 std::vector<int>WhichValues;
1431 while(i <= nrow)
1432 {
1433 int j = 1;
1434 while(j <= ncol*2)
1435 {
1436
1437 //UNORMC.push_back((double)(cser[(j-1)*nsnew[1]+i])/(double) (denom));
1438 UNORMC.push_back((double)(cser[(i-1)*(nsnew[1]*2)+j])/(double) (denom));
1439 WhichValues.push_back((i-1)*(nsnew[1]*2)+j);
1440 ndxout++;
1441 j+=2; // because of pseudo complex configuration
1442 }
1443 i+=1; // I think rows are not influenced by complex padding
1444 }
1445
1446 return;
1447 }
1448
1449
1450 //---------------------------------------------------------
fft2(std::vector<double> & data,std::vector<int> nel,int isign)1451 void CGrid_IMCORR::fft2(std::vector<double>& data, std::vector<int>nel, int isign)
1452 {
1453 double wr,wi,wpr,wpi,wtemp,theta;
1454 int ntot = nel[0]*nel[1];
1455
1456 int nprev = 1;
1457 int idim = 1;
1458 while(idim <= 2)
1459 {
1460 int n = nel[idim-1];
1461 int nrem = ntot/(n*nprev);
1462 int ip1 = 2*nprev;
1463 int ip2 = ip1*n;
1464 int ip3 = ip2 * nrem;
1465 int i2rev = 1;
1466
1467 // perform bit reversal
1468 int i2 = 1;
1469 while(i2 <= ip2)
1470 {
1471 if(i2 < i2rev)
1472 {
1473 int i1= i2;
1474 while (i1 <= (i2+ip1-2))
1475 {
1476 int i3=i1;
1477 while (i3 <= ip3)
1478 {
1479 int i3rev = i2rev+i3-i2;
1480 double tempr = data[i3];
1481 double tempi = data[i3+1];
1482
1483 data[i3] = data[i3rev];
1484 data[i3+1] = data[i3rev+1];
1485 data[i3rev] = tempr;
1486 data[i3rev+1] = tempi;
1487 i3+=ip2;
1488 }
1489 i1+=2;
1490 }
1491 }
1492 int ibit = ip2/2;
1493 while (ibit >= ip1 && i2rev > ibit)
1494 {
1495 i2rev = i2rev-ibit;
1496 ibit = ibit/2;
1497 }
1498 i2rev = i2rev+ibit;
1499 i2+=ip1;
1500 }
1501 // end of bit reversal
1502
1503
1504 // Danielson-Lanczos Section
1505 int ifp1 = ip1;
1506 while(ifp1 < ip2)
1507 {
1508 int ifp2 = 2* ifp1;
1509 // Initialize for the trigonometric recurrence
1510 theta = isign *(M_PI*2) / (ifp2/ip1);
1511 wtemp = sin(0.5*theta);
1512 wpr = -2.0 * (wtemp*wtemp);
1513 wpi = sin(theta);
1514 wr = 1.0;
1515 wi = 0.0;
1516 int i3=1;
1517 while (i3 <= ifp1)
1518 {
1519 int i1 = i3;
1520 while (i1 <= (i3+ip1-2))
1521 {
1522 int i2 = i1;
1523 while (i2 <= ip3)
1524 {
1525 int k1 = i2;
1526 int k2 = k1+ifp1;
1527 double tempr = wr * data[k2] - wi * data[k2+1];
1528 double tempi = wr * data[k2+1] + wi * data[k2];
1529 data[k2] = data[k1]-tempr;
1530 data[k2+1] = data[k1+1]-tempi;
1531 data[k1] = data[k1]+tempr;
1532 data[k1+1] = data[k1+1]+tempi;
1533 i2+=ifp2;
1534 }
1535 i1+=2;
1536 }
1537 // the trigonometric recurrence
1538 wtemp = wr;
1539 wr = wr*wpr-wi*wpi+wr;
1540 wi=wi*wpr+wtemp*wpi+wi;
1541 i3+=ip1;
1542 }
1543 ifp1=ifp2;
1544 }
1545 nprev = n*nprev;
1546 idim++;
1547 }
1548
1549 return;
1550 }
1551
1552
1553 //---------------------------------------------------------
decimal(std::vector<int> Bin)1554 int CGrid_IMCORR::decimal(std::vector<int> Bin)
1555 {
1556 int Sum = 0;
1557 for (int i = 0; i < Bin.size(); i++)
1558 {
1559 int Exponent = (int)(Bin.size()-1)-i;
1560 Sum+= Bin[i] * (int)pow(2.0, Exponent);
1561 }
1562 return Sum;
1563 }
1564
1565
1566 //---------------------------------------------------------
binary(std::vector<int> & Bin,int number)1567 void CGrid_IMCORR::binary(std::vector<int>& Bin, int number)
1568 {
1569 int remainder;
1570
1571 if(number <= 1)
1572 {
1573 Bin.push_back(number);
1574 return;
1575 }
1576 remainder = number%2;
1577 binary(Bin, number >> 1);
1578 Bin.push_back(remainder);
1579 }
1580
1581
1582 //---------------------------------------------------------
Get_This_Chip(std::vector<std::vector<double>> & Chip,CSG_Grid * pGrid,int gx,int gy,int Ref_Chipsize)1583 void CGrid_IMCORR::Get_This_Chip(std::vector<std::vector<double> >& Chip, CSG_Grid *pGrid, int gx, int gy, int Ref_Chipsize)
1584 {
1585 int ref_chipX = 0;
1586 int ref_chipY = 0;
1587 for (int refx = gx - (Ref_Chipsize/2-1); refx < gx - (Ref_Chipsize/2-1) + Ref_Chipsize; refx++)
1588 {
1589 ref_chipY = 0;
1590 for (int refy = gy - (Ref_Chipsize/2-1); refy < gy - (Ref_Chipsize/2-1) + Ref_Chipsize; refy++)
1591 {
1592 Chip[ref_chipX][ref_chipY] = pGrid->asDouble(refx,refy);
1593 ref_chipY++;
1594 }
1595 ref_chipX++;
1596 }
1597
1598 return;
1599 }
1600
1601 //---------------------------------------------------------
1602