Loading [MathJax]/jax/output/NativeMML/config.js
Isis 3 Programmer Reference
Gruen.cpp
1
5
6/* SPDX-License-Identifier: CC0-1.0 */
7#include "Chip.h"
8#include "GruenTypes.h"
9#include "Gruen.h"
10#include "DbProfile.h"
11#include "Pvl.h"
12
13#include "tnt/tnt_array1d_utils.h"
14
15#include <cmath>
16#include <iostream>
17#include <sstream>
18#include <iomanip>
19
20#include <QTime>
21
22
23using namespace std;
24
25namespace Isis {
26
35
47 Gruen::Gruen(Pvl &pvl) : AutoReg(pvl) {
48 init(pvl);
49 }
50
92 void Gruen::WriteSubsearchChips(const QString &pattern) {
93 m_filePattern = pattern;
94 return;
95 }
96
111 void Gruen::setAffineRadio(const AffineRadio &affrad) {
112 m_affine = affrad;
113 return;
114 }
115
128 m_affine = getDefaultAffineRadio();
129 }
130
156 int Gruen::algorithm(Chip &pattern, Chip &subsearch,
157 const Radiometric &radio, BigInt &ptsUsed,
158 double &resid, GMatrix &atai, AffineRadio &affrad) {
159
160 m_totalIterations++; // Bump iteration counter
161
162 // Initialize loop variables
163 int tackSamp = pattern.TackSample();
164 int tackLine = pattern.TackLine();
165
166 // Internal variables
167 double rshift = radio.Shift();
168 double rgain = radio.Gain();
169
170 int maxPnts((pattern.Samples()-2) * (pattern.Lines()-2));
171 GMatrix a(maxPnts, 8);
172 GVector lvec(maxPnts);
173
174 // pattern chip is rh image , subsearch chip is lh image
175 resid = 0.0;
176 int npts = 0;
177 for(int line = 2 ; line < pattern.Lines() ; line++) {
178 for(int samp = 2; samp < pattern.Samples() ; samp++) {
179
180 if(!pattern.IsValid(samp, line)) continue;
181 if(!subsearch.IsValid(samp, line)) continue;
182 if(!subsearch.IsValid(samp + 1, line)) continue;
183 if(!subsearch.IsValid(samp - 1, line)) continue;
184 if(!subsearch.IsValid(samp, line - 1)) continue;
185 if(!subsearch.IsValid(samp, line + 1)) continue;
186
187 // Sample/Line numbers
188 double x0 = (double)(samp - tackSamp);
189 double y0 = (double)(line - tackLine);
190
191 // Discrete derivatives (delta sample #)
192 double gxtemp = subsearch.GetValue(samp + 1, line) -
193 subsearch.GetValue(samp - 1, line);
194 double gytemp = subsearch.GetValue(samp, line + 1) -
195 subsearch.GetValue(samp, line - 1);
196
197 a[npts][0] = gxtemp;
198 a[npts][1] = gxtemp * x0;
199 a[npts][2] = gxtemp * y0;
200 a[npts][3] = gytemp;
201 a[npts][4] = gytemp * x0;
202 a[npts][5] = gytemp * y0;
203 a[npts][6] = 1.0;
204 a[npts][7] = subsearch.GetValue(samp, line);
205
206 double ell = pattern.GetValue(samp, line) -
207 (((1.0 + rgain) * subsearch.GetValue(samp, line)) +
208 rshift);
209
210 lvec[npts] = ell;
211 resid += (ell * ell);
212 npts++;
213 }
214 }
215
216 // Check for enough points
217 ptsUsed = npts;
218 if(!ValidPoints(maxPnts, npts)) {
219 std::ostringstream mess;
220 mess << "Minimum points (" << MinValidPoints(maxPnts)
221 << ") criteria not met (" << npts << ")";
222 return (logError(NotEnoughPoints, mess.str().c_str()));
223 }
224
225 // Create the ATA array
226 GMatrix ata(8,8);
227 for(int i = 0 ; i < 8 ; i++) {
228 for(int j = 0 ; j < 8 ; j++) {
229 double temp(0.0);
230 for (int k = 0 ; k < npts ; k++) {
231 temp += (a[k][i] * a[k][j]);
232 }
233 ata[i][j] = temp;
234 }
235 }
236
237 // Solve the ATAI with Cholesky
238 try {
239 GVector p(8);
240 atai = Choldc(ata, p);
241 GMatrix b = Identity(8);
242 GMatrix x = b;
243 atai = Cholsl(atai, p, b, x);
244 }
245 catch(IException &ie) {
246 QString mess = "Cholesky Failed:: " + ie.toString();
247 return (logError(CholeskyFailed, mess));
248 }
249
250 // Compute the affine update if result are requested by caller.
251 GVector atl(8);
252 for (int i = 0 ; i < 8 ; i++) {
253 double temp(0.0);
254 for (int k = 0 ; k < npts ; k++) {
255 temp += a[k][i] * lvec[k];
256 }
257 atl[i] = temp;
258 }
259
260 GVector alpha(8, 0.0);
261 for (int i = 0 ; i < 8 ; i++) {
262 for (int k = 0 ; k < 8 ; k++) {
263 alpha[i] += atai[i][k] * atl[k];
264 }
265 }
266
267 try {
268 affrad = AffineRadio(alpha);
269 }
270 catch(IException &ie) {
271 QString mess = "Affine failed: " + ie.toString();
272 return (logError(AffineNotInvertable, mess));
273 }
274
275 return (0);
276 }
277
289 Analysis Gruen::errorAnalysis(const BigInt &npts, const double &resid,
290 const GMatrix &atai) {
291
292 Analysis results;
293 results.m_npts = npts;
294
295 // Converged, compute error anaylsis
296 double variance = resid / DegreesOfFreedom(npts);
297 GMatrix kmat(8,8);
298 for(int r = 0 ; r < 8 ; r++) {
299 for(int c = 0 ; c < 8 ; c++) {
300 kmat[r][c] = variance * atai[r][c];
301 }
302 }
303 results.m_variance = variance;
304
305
306 // Set up submatrix
307 GMatrix skmat(2,2);
308 skmat[0][0] = kmat[0][0];
309 skmat[0][1] = kmat[0][3];
310 skmat[1][0] = kmat[3][0];
311 skmat[1][1] = kmat[3][3];
312 try {
313 GVector eigen(2);
314 Jacobi(skmat, eigen, skmat);
315 EigenSort(eigen, skmat);
316 for (int i = 0 ; i < 2 ; i++) {
317 results.m_sevals[i] = eigen[i];
318 results.m_kmat[i] = kmat[i*3][i*3];
319 }
320 }
321 catch(IException &ie) {
322 QString errmsg = "Eigen Solution Failed:: " + ie.toString();
323 results.m_status = logError(EigenSolutionFailed, errmsg);
324 return (results);
325 }
326
327 results.m_status = 0;
328 return (results);
329 }
330
341 GMatrix Gruen::Choldc(const GMatrix &a, GVector &p) const {
342 int nrows(a.dim1());
343 int ncols(a.dim2());
344
345 GMatrix aa = a.copy();
346 p = GVector(ncols);
347
348 for(int i = 0 ; i < nrows ; i++) {
349 for(int j = i ; j < ncols ; j++) {
350 double sum = aa[i][j];
351 for(int k = i-1 ; k >= 0 ; k--) {
352 sum -= (aa[i][k] * aa[j][k]);
353 }
354 // Handle diagonal special
355 if (i == j) {
356 if (sum <= 0.0) {
358 "Choldc failed - matrix not postive definite",
359 _FILEINFO_);
360 }
361 p[i] = sqrt(sum);
362 }
363 else {
364 aa[j][i] = sum / p[i];
365 }
366 }
367 }
368 return (aa);
369 }
370
383 GMatrix Gruen::Cholsl(const GMatrix &a,const GVector &p, const GMatrix &b,
384 const GMatrix &x) const {
385 assert(b.dim1() == x.dim1());
386 assert(b.dim2() == x.dim2());
387 assert(p.dim1() == b.dim2());
388
389 int nrows(a.dim1());
390 int ncols(a.dim2());
391
392 GMatrix xout = x.copy();
393 for(int j = 0 ; j < nrows ; j++) {
394
395 for(int i = 0 ; i < ncols ; i++) {
396 double sum = b[j][i];
397 for(int k = i-1 ; k >= 0 ; k--) {
398 sum -= (a[i][k] * xout[j][k]);
399 }
400 xout[j][i] = sum / p[i];
401 }
402
403 for (int i = ncols-1 ; i >= 0 ; i--) {
404 double sum = xout[j][i];
405 for(int k = i+1 ; k < ncols ; k++) {
406 sum -= (a[k][i] * xout[j][k]);
407 }
408 xout[j][i] = sum / p[i];
409 }
410
411 }
412 return (xout);
413 }
414
427 int Gruen::Jacobi(const GMatrix &a, GVector &evals,
428 GMatrix &evecs, const int &MaxIters) const {
429
430 int nrows(a.dim1());
431 int ncols(a.dim2());
432 GMatrix v = Identity(nrows);
433 GVector d(nrows),b(nrows), z(nrows);
434
435 for(int ip = 0 ; ip < nrows ; ip++) {
436 b[ip] = a[ip][ip];
437 d[ip] = b[ip];
438 z[ip] = 0.0;
439 }
440
441 double n2(double(nrows) * double(nrows));
442 GMatrix aa = a.copy();
443 int nrot(0);
444 for ( ; nrot < MaxIters ; nrot++) {
445 double sm(0.0);
446 for(int ip = 0 ; ip < nrows-1 ; ip++) {
447 for(int iq = ip+1 ; iq < nrows ; iq++) {
448 sm += fabs(aa[ip][iq]);
449 }
450 }
451
452 // Test for termination condition
453 if (sm == 0.0) {
454 evals = d;
455 evecs = v;
456 return (nrot);
457 }
458
459 double thresh = (nrot < 3) ? (0.2 * sm/n2 ): 0.0;
460 for (int ip = 0 ; ip < nrows-1 ; ip++) {
461 for (int iq = ip+1 ; iq < nrows ; iq++) {
462 double g = 100.0 * fabs(aa[ip][iq]);
463 if ( (nrot > 3) &&
464 ( (fabs(d[ip]+g) == fabs(d[ip])) ) &&
465 ( (fabs(d[iq]+g) == fabs(d[iq])) ) ) {
466 aa[ip][iq] = 0.0;
467 }
468 else if ( fabs(aa[ip][iq]) > thresh ) {
469 double h = d[iq] - d[ip];
470 double t;
471 if ( (fabs(h)+g) == fabs(h) ) {
472 t = aa[ip][iq]/h;
473 }
474 else {
475 double theta = 0.5 * h / aa[ip][iq];
476 t = 1.0 / (fabs(theta) + sqrt(1.0 + theta * theta));
477 if (theta < 0.0) t = -1.0 * t;
478 }
479 double c = 1./sqrt(1.0 + t * t);
480 double s = t * c;
481 double tau = s / (1.0 + c);
482
483 h = t * aa[ip][iq];
484 z[ip] = z[ip] - h;
485 z[iq] = z[iq] + h;
486 d[ip] = d[ip] - h;
487 d[iq] = d[iq] + h;
488 aa[ip][iq] = 0.0;
489
490 double g;
491 for (int j = 0 ; j < ip-1 ; j++) {
492 g = aa[j][ip];
493 h = aa[j][iq];
494 aa[j][ip] = g - s * (h + g * tau);
495 aa[j][iq] = h + s * (g - h * tau);
496 }
497
498 for (int j = ip+1 ; j < iq-1 ; j++) {
499 g = aa[ip][j];
500 h = aa[j][iq];
501 aa[ip][j] = g - s * (h + g * tau);
502 aa[j][iq] = h + s * (g - h * tau);
503 }
504
505 for (int j = iq+1 ; j < ncols ; j++) {
506 g = aa[ip][j];
507 h = aa[j][iq];
508 aa[ip][j] = g - s * (h + g * tau);
509 aa[j][iq] = h + s * (g - h * tau);
510 }
511
512 for (int j = 0 ; j < ncols ; j++) {
513 g = v[j][ip];
514 h = v[j][iq];
515 v[j][ip] = g - s * (h + g * tau);
516 v[j][iq] = h + s * (g - h * tau);
517 }
518 nrot++;
519 }
520 }
521 }
522
523 for (int ip = 0 ; ip < nrows ; ip++) {
524 b[ip] = b[ip] + z[ip];
525 d[ip] = b[ip];
526 z[ip] = 0.0;
527 }
528 }
529
530 // Reach here and we have too many iterations
531 evals = d;
532 evecs = v;
534 "Too many iterations in Jacobi",
535 _FILEINFO_);
536 return (nrot);
537 }
538
539
540 GMatrix Gruen::Identity(const int &ndiag) const {
541 GMatrix ident(ndiag, ndiag, 0.0);
542 for (int i = 0 ; i < ndiag ; i++) {
543 ident[i][i] = 1.0;
544 }
545 return (ident);
546 }
547
556 void Gruen::EigenSort(GVector &evals, GMatrix &evecs) const {
557 assert(evals.dim1() == evecs.dim1());
558 int n(evals.dim1());
559 for (int i = 0 ; i < n-1 ; i++ ) {
560 int k = i;
561 double p = evals[i];
562 for (int j = i+1 ; j < n ; j++) {
563 if (evals[j] >= p) {
564 k = j;
565 p = evals[j];
566 }
567 }
568 if (k != i) {
569 evals[k] = evals[i];
570 evals[i] = p;
571 for (int j = 0 ; j < n ; j++) {
572 p = evecs[j][i];
573 evecs[j][i] = evecs[j][k];
574 evecs[j][k] = p;
575 }
576 }
577 }
578 return;
579 }
580
581
608 double Gruen::MatchAlgorithm(Chip &pattern, Chip &subsearch) {
609 // reset();
610
611 BigInt npts;
612 double resid;
613 GMatrix atai;
614 AffineRadio affrad;
615 int status = algorithm(pattern, subsearch, getDefaultRadio(),
616 npts, resid, atai, affrad);
617 if (status == 0) {
618 // Compute fit quality
619 Analysis analysis = errorAnalysis(npts, resid, atai);
620 if (analysis.isValid()) {
621 return (analysis.getEigen());
622 }
623 }
624
625 // Error conditions return failure
626 return (Null);
627 }
628
636 bool Gruen::CompareFits(double fit1, double fit2) {
637 return (fit1 <= fit2);
638 }
639
681 Chip &fChip, int startSamp,
682 int startLine, int endSamp,
683 int endLine, int bestSamp,
684 int bestLine) {
685 // Set conditions for writing subsearch chip states per call. This code
686 // implementation ensures caller must request it per call to this routine.
687 // See the WriteSubsearchChips() method.
688 m_callCount++;
689 QString chipOut = m_filePattern;
690 m_filePattern = "";
691
692 // Initialize match point. Ensure points are centered to get real cube
693 // line/sample positions
694 pChip.SetChipPosition(pChip.TackSample(), pChip.TackLine());
695 sChip.SetChipPosition(sChip.TackSample(), sChip.TackLine());
696 MatchPoint matchpt = MatchPoint(PointPair(Coordinate(pChip),
697 Coordinate(sChip)));
698
699 // Create the fit chip whose size is the same as the pattern chip.
700 // This chip will contain the final image at the last iteration. This
701 // usage differs from the non-adaptive purpose.
702 // It is critical to use the original search chip to create the subsearch.
703 // Copying the original search chip and then resizing preserves established
704 // minimum/maximum value ranges. Then, establish chip convergence
705 // condition for Gruen's affine.
706 fChip = sChip;
707 fChip.SetSize(pChip.Samples(), pChip.Lines());
708 Threshold thresh(fChip, getAffineTolerance());
709
710 // Set up Affine transform by establishing search tack point
711 AffineRadio affine = m_affine;
712 m_affine = AffineRadio(); // Reset initial condition for next call
713
714 // Set up bestLine/bestSample position. Do this using the local affine
715 // and not the search chip.
716 Coordinate best(bestLine-sChip.TackLine(), bestSamp-sChip.TackSample());
717 affine.Translate(best);
718
719
720 // Algorithm parameters
721 bool done = false;
722 m_nIters = 0;
723 do {
724
725 // Extract the sub search chip. The algorithm method handles the
726 // determination of good data volumes
727 Affine extractor(affine.m_affine);
728 sChip.Extract(fChip, extractor);
729
730 // If requested for this run, write the current subsearch chip state
731 if (!chipOut.isEmpty()) {
732 std::ostringstream ss;
733 ss << "C" << std::setw(6) << std::setfill('0') << m_callCount << "I"
734 << std::setw(3) << std::setfill('0') << m_nIters;
735 QString sfname = chipOut + ss.str().c_str() + ".cub";
736 fChip.Write(sfname);
737 }
738
739 // Try to match the two subchips
740 AffineRadio alpha;
741 BigInt npts(0);
742 GMatrix atai;
743 double resid;
744 int status = algorithm(pChip, fChip, affine.m_radio, npts, resid,
745 atai, alpha);
746 if (status != 0) {
747 // Set failed return condition and give up!
748 return (Status(matchpt.setStatus(status)));
749 }
750
751 // Test for termination conditions - errors or convergence
752 matchpt.m_nIters = ++m_nIters;
753 if (m_nIters > m_maxIters) {
754 QString errmsg = "Maximum Iterations exceeded";
755 matchpt.setStatus(logError(MaxIterationsExceeded, errmsg));
756 return (Status(matchpt)); // Error condition
757 }
758
759 // Check for convergence after the first pass
760 if (m_nIters > 1) {
761 if(thresh.hasConverged(alpha)) {
762 // Compute error analysis
763 matchpt.m_affine = affine;
764 matchpt.m_analysis = errorAnalysis(npts, resid, atai);
765 matchpt.setStatus(matchpt.m_analysis.m_status);
766 if (matchpt.isValid()) {
767 // Update the point even if constraints don't pass
768 Coordinate uCoord = getChipUpdate(sChip, matchpt);
769 SetChipSample(uCoord.getSample());
770 SetChipLine(uCoord.getLine());
771 SetGoodnessOfFit(matchpt.getEigen());
772
773 // Check constraints
774 matchpt.setStatus(CheckConstraints(matchpt));
775 }
776
777 // Set output point
778 m_point = matchpt;
779 return (Status(matchpt)); // with AutoReg status
780 }
781 }
782 // Not done yet - update the affine transform for next loop
783 try {
784 affine += alpha;
785 }
786 catch (IException &ie) {
787 QString mess = "Affine invalid/not invertable";
788 matchpt.setStatus(logError(AffineNotInvertable, mess));
789 return (Status(matchpt)); // Another error condition to return
790 }
791 } while (!done);
792
793 return (Status(matchpt));
794 }
795
804 static Pvl regdef;
805 regdef = Pvl("$ISISROOT/appdata/templates/autoreg/coreg.adaptgruen.p1515s3030.def");
806 return (regdef);
807 }
808
825 PvlGroup algo("GruenFailures");
826 algo += PvlKeyword("Name", AlgorithmName());
827 algo += PvlKeyword("Mode", "Adaptive");
828
829 // Log errors
830 for (int e = 0 ; e < m_errors.size() ; e++) {
831 algo += m_errors.getNth(e).LogIt();
832 }
833
834 if (m_unclassified > 0) {
835 algo += PvlKeyword("UnclassifiedErrors", toString(m_unclassified));
836 }
837 pvl.addGroup(algo);
838 pvl.addGroup(StatsLog());
839 pvl.addGroup(ParameterLog());
840 return (pvl);
841 }
842
854 PvlGroup stats("GruenStatistics");
855
856 stats += PvlKeyword("TotalIterations", toString(m_totalIterations));
857 stats += ValidateKey("IterationMinimum", m_iterStat.Minimum());
858 stats += ValidateKey("IterationAverage", m_iterStat.Average());
859 stats += ValidateKey("IterationMaximum", m_iterStat.Maximum());
860 stats += ValidateKey("IterationStandardDeviation", m_iterStat.StandardDeviation());
861
862 stats += ValidateKey("EigenMinimum", m_eigenStat.Minimum());
863 stats += ValidateKey("EigenAverage", m_eigenStat.Average());
864 stats += ValidateKey("EigenMaximum", m_eigenStat.Maximum());
865 stats += ValidateKey("EigenStandardDeviation", m_eigenStat.StandardDeviation());
866
867 stats += ValidateKey("RadioShiftMinimum", m_shiftStat.Minimum());
868 stats += ValidateKey("RadioShiftAverage", m_shiftStat.Average());
869 stats += ValidateKey("RadioShiftMaximum", m_shiftStat.Maximum());
870 stats += ValidateKey("RadioShiftStandardDeviation", m_shiftStat.StandardDeviation());
871
872 stats += ValidateKey("RadioGainMinimum", m_gainStat.Minimum());
873 stats += ValidateKey("RadioGainAverage", m_gainStat.Average());
874 stats += ValidateKey("RadioGainMaximum", m_gainStat.Maximum());
875 stats += ValidateKey("RadioGainStandardDeviation", m_gainStat.StandardDeviation());
876
877 return (stats);
878
879 }
880
891 PvlGroup parms("GruenParameters");
892
893 parms += PvlKeyword("MaximumIterations", toString(m_maxIters));
894 parms += ValidateKey("AffineScaleTolerance", m_scaleTol);
895 parms += ValidateKey("AffineShearTolerance", m_shearTol);
896 parms += ValidateKey("AffineTranslationTolerance", m_transTol);
897
898 parms += ParameterKey("AffineTolerance", m_affineTol);
899 parms += ParameterKey("SpiceTolerance", m_spiceTol);
900
901 parms += ParameterKey("RadioShiftTolerance", m_shiftTol);
902
903 parms += ParameterKey("RadioGainMinTolerance", m_rgainMinTol);
904 parms += ParameterKey("RadioGainMaxTolerance", m_rgainMaxTol);
905
906 parms += ValidateKey("DefaultRadioGain", m_defGain);
907 parms += ValidateKey("DefaultRadioShift", m_defShift);
908
909 return (parms);
910 }
911
912
923 ErrorList elist;
924 elist.add(1, ErrorCounter(1, "NotEnoughPoints"));
925 elist.add(2, ErrorCounter(2, "CholeskyFailed"));
926 elist.add(3, ErrorCounter(3, "EigenSolutionFailed"));
927 elist.add(4, ErrorCounter(4, "AffineNotInvertable"));
928 elist.add(5, ErrorCounter(5, "MaxIterationsExceeded"));
929 elist.add(6, ErrorCounter(6, "RadShiftExceeded"));
930 elist.add(7, ErrorCounter(7, "RadGainExceeded"));
931 elist.add(8, ErrorCounter(8, "MaxEigenExceeded"));
932 elist.add(9, ErrorCounter(9, "AffineDistExceeded"));
933 return (elist);
934 }
935
936
948 int Gruen::logError(int gerrno, const QString &gerrmsg) {
949 if (!m_errors.exists(gerrno)) {
950 m_unclassified++;
951 }
952 else {
953 m_errors.get(gerrno).BumpIt();
954 }
955 return (gerrno);
956 }
957
967 void Gruen::init(Pvl &pvl) {
968 // Establish the parameters
969 if (pvl.hasObject("AutoRegistration")) {
970 m_prof = DbProfile(pvl.findGroup("Algorithm", Pvl::Traverse));
971 }
972 else {
973 m_prof = DbProfile(pvl);
974 }
975
976 if (m_prof.Name().isEmpty()) m_prof.setName("Gruen");
977
978 // Define internal parameters
979 m_maxIters = toInt(ConfKey(m_prof, "MaximumIterations", toString(30)));
980
981 m_transTol = toDouble(ConfKey(m_prof, "AffineTranslationTolerance", toString(0.1)));
982 m_scaleTol = toDouble(ConfKey(m_prof, "AffineScaleTolerance", toString(0.3)));
983 m_shearTol = toDouble(ConfKey(m_prof, "AffineShearTolerance", toString(m_scaleTol)));
984 m_affineTol = toDouble(ConfKey(m_prof, "AffineTolerance", toString(DBL_MAX)));
985
986 m_spiceTol = toDouble(ConfKey(m_prof, "SpiceTolerance", toString(DBL_MAX)));
987
988 m_shiftTol = toDouble(ConfKey(m_prof, "RadioShiftTolerance", toString(DBL_MAX)));
989 m_rgainMinTol = toDouble(ConfKey(m_prof, "RadioGainMinTolerance", toString(-DBL_MAX)));
990 m_rgainMaxTol = toDouble(ConfKey(m_prof, "RadioGainMaxTolerance", toString(DBL_MAX)));
991
992 // Set radiometric defaults
993 m_defGain = toDouble(ConfKey(m_prof, "DefaultRadioGain", toString(0.0)));
994 m_defShift = toDouble(ConfKey(m_prof, "DefaultRadioShift", toString(0.0)));
995
996 m_callCount = 0;
997 m_filePattern = "";
998
999 m_nIters = 0;
1000 m_totalIterations = 0;
1001
1002 m_errors = initErrorList();
1003 m_unclassified = 0;
1004
1005 m_defAffine = AffineRadio(getDefaultRadio());
1006 m_affine = getAffineRadio();
1007 m_point = MatchPoint(m_affine);
1008
1009 //reset();
1010 return;
1011 }
1012
1013
1019 m_eigenStat.Reset();
1020 m_iterStat.Reset();
1021 m_shiftStat.Reset();
1022 m_gainStat.Reset();
1023 return;
1024 }
1025
1039 double pts = (double) totalPoints * (PatternValidPercent() / 100.0);
1040 return ((BigInt) pts);
1041 }
1042
1055 bool Gruen::ValidPoints(BigInt totalPoints, BigInt nPoints) const {
1056 return (nPoints > MinValidPoints(totalPoints));
1057 }
1058
1067 return (AffineTolerance(m_transTol, m_scaleTol, m_shearTol));
1068 }
1069
1091
1092 // Point must be good for check to occur
1093 if (point.isValid()) {
1094 if (point.m_nIters > m_maxIters) {
1095 QString errmsg = "Maximum Iterations exceeded";
1096 return (logError(MaxIterationsExceeded, errmsg));
1097 }
1098 m_iterStat.AddData(point.m_nIters);
1099
1100 if (point.getEigen() > Tolerance()) {
1101 QString errmsg = "Maximum Eigenvalue exceeded";
1102 return (logError(MaxEigenExceeded, errmsg));
1103 }
1104 m_eigenStat.AddData(point.getEigen());
1105
1106 double shift = point.m_affine.m_radio.Shift();
1107 if ( shift > m_shiftTol) {
1108 QString errmsg = "Radiometric shift exceeds tolerance";
1109 return (logError(RadShiftExceeded, errmsg));
1110 }
1111 m_shiftStat.AddData(shift);
1112
1113 double gain = point.m_affine.m_radio.Gain();
1114 if (((1.0 + gain) > m_rgainMaxTol) ||
1115 ((1.0 + gain) < m_rgainMinTol)) {
1116 QString errmsg = "Radiometric gain exceeds tolerances";
1117 return (logError(RadGainExceeded, errmsg));
1118 }
1119 m_gainStat.AddData(gain);
1120
1121
1122 double dist = point.getAffinePoint(Coordinate(0.0, 0.0)).getDistance();
1123 if (dist > getAffineConstraint()) {
1124 QString errmsg = "Affine distance exceeded";
1125 return (logError(AffineDistExceeded, errmsg));
1126 }
1127 }
1128 return (point.getStatus());
1129 }
1130
1142 Coordinate chippt = point.getAffinePoint(Coordinate(0.0, 0.0));
1143 chip.SetChipPosition(chip.TackSample(), chip.TackLine());
1144 chip.TackCube(chip.CubeSample()+chippt.getSample(),
1145 chip.CubeLine()+chippt.getLine());
1146 return (Coordinate(chip.ChipLine(), chip.ChipSample()));
1147 }
1148
1149
1162 if ( mpt.isValid() ) { return (AutoReg::SuccessSubPixel); }
1164 }
1165
1166} // namespace Isis
Affine basis function.
Definition Affine.h:65
Container for affine and radiometric parameters.
Definition GruenTypes.h:242
void Translate(const Coordinate &offset)
Apply a translation to the given offset.
Definition GruenTypes.h:263
RegisterStatus
Enumeration of the Register() method's return status.
Definition AutoReg.h:179
@ SuccessSubPixel
Success registering to sub-pixel accuracy.
Definition AutoReg.h:181
@ AdaptiveAlgorithmFailed
Error occured in Adaptive algorithm.
Definition AutoReg.h:189
void SetChipLine(double line)
Sets the search chip subpixel line that matches the pattern tack line.
Definition AutoReg.h:404
AutoReg(Pvl &pvl)
Create AutoReg object.
Definition AutoReg.cpp:66
double PatternValidPercent() const
Return pattern chip valid percent. The default value is.
Definition AutoReg.h:280
void SetChipSample(double sample)
Sets the search chip subpixel sample that matches the pattern tack sample.
Definition AutoReg.h:393
void SetGoodnessOfFit(double fit)
Sets the goodness of fit for adaptive algorithms.
Definition AutoReg.h:413
double Tolerance() const
Return match algorithm tolerance.
Definition AutoReg.h:290
A small chip of data used for pattern matching.
Definition Chip.h:86
double CubeLine() const
Definition Chip.h:210
void SetSize(const int samples, const int lines)
Change the size of the Chip.
Definition Chip.cpp:134
void Write(const QString &filename)
Writes the contents of the Chip to a cube.
Definition Chip.cpp:985
double GetValue(int sample, int line)
Loads a Chip with a value.
Definition Chip.h:145
int Samples() const
Definition Chip.h:99
int TackSample() const
This method returns a chip's fixed tack sample; the middle of the chip.
Definition Chip.h:176
int Lines() const
Definition Chip.h:106
void SetChipPosition(const double sample, const double line)
Compute the position of the cube given a chip coordinate.
Definition Chip.cpp:643
double ChipSample() const
Definition Chip.h:219
double ChipLine() const
Definition Chip.h:226
void TackCube(const double cubeSample, const double cubeLine)
This sets which cube position will be located at the chip tack position.
Definition Chip.cpp:182
double CubeSample() const
Definition Chip.h:203
bool IsValid(int sample, int line)
Definition Chip.h:240
Chip Extract(int samples, int lines, int samp, int line)
Extract a sub-chip from a chip.
Definition Chip.cpp:727
int TackLine() const
This method returns a chip's fixed tack line; the middle of the chip.
Definition Chip.h:187
void add(const K &key, const T &value)
Adds the element to the list.
Define a generic Y/X container.
Definition GruenTypes.h:53
double getDistance(const Coordinate &pntA=Coordinate(0.0, 0.0)) const
Computes the distance from this point and the point provided.
Definition GruenTypes.h:100
A DbProfile is a container for access parameters to a database.
Definition DbProfile.h:51
void EigenSort(GVector &evals, GMatrix &evecs) const
Sort eigenvectors from highest to lowest.
Definition Gruen.cpp:556
double getAffineConstraint() const
Returns the Affine tolerance constraint as read from config file.
Definition Gruen.h:95
void WriteSubsearchChips(const QString &pattern="SubChip")
Set up for writing subsearch for a a given registration call.
Definition Gruen.cpp:92
void resetStats()
Reset Gruen statistics as needed.
Definition Gruen.cpp:1018
QString ConfKey(const DbProfile &conf, const QString &keyname, const QString &defval, int index=0) const
Helper method to initialize parameters.
Definition Gruen.h:241
Gruen()
Default constructor.
Definition Gruen.cpp:32
int Jacobi(const GMatrix &a, GVector &evals, GMatrix &evecs, const int &MaxIters=50) const
Compute the Jacobian of a covariance matrix.
Definition Gruen.cpp:427
const AffineRadio & getDefaultAffineRadio() const
Returns default settings for Affine/Radiometric parameters.
Definition Gruen.h:101
CollectorMap< int, ErrorCounter > ErrorList
Declaration of error count list.
Definition Gruen.h:182
void init(Pvl &pvl)
Initialize the object.
Definition Gruen.cpp:967
Coordinate getChipUpdate(Chip &sChip, MatchPoint &point) const
Compute the chip coordinate of the registered pixel.
Definition Gruen.cpp:1141
virtual double MatchAlgorithm(Chip &pattern, Chip &subsearch)
Minimization of data set using Gruen algorithm.
Definition Gruen.cpp:608
virtual QString AlgorithmName() const
Returns the default name of the algorithm as Gruen.
Definition Gruen.h:124
Analysis errorAnalysis(const BigInt &npts, const double &resid, const GMatrix &atai)
Compute the error analysis of convergent Gruen matrix.
Definition Gruen.cpp:289
static Pvl & getDefaultParameters()
Load default Gruen parameter file in $ISISROOT/appdata/templates.
Definition Gruen.cpp:803
int logError(int gerrno, const QString &gerrmsg)
Logs a Gruen error.
Definition Gruen.cpp:948
bool ValidPoints(BigInt totalPoints, BigInt nPoints) const
Determines if number of points is valid percentage of all points.
Definition Gruen.cpp:1055
void setAffineRadio()
Set affine parameters to defaults.
Definition Gruen.cpp:127
virtual bool CompareFits(double fit1, double fit2)
This virtual method must return if the 1st fit is equal to or better than the second fit.
Definition Gruen.cpp:636
GMatrix Cholsl(const GMatrix &a, const GVector &p, const GMatrix &b, const GMatrix &x) const
Compute Cholesky solution matrix from correlation.
Definition Gruen.cpp:383
AutoReg::RegisterStatus Status(const MatchPoint &result)
Returns the proper status given a Gruen result container.
Definition Gruen.cpp:1161
PvlGroup ParameterLog() const
Create a PvlGroup with the Gruen specific parameters.
Definition Gruen.cpp:890
PvlKeyword ValidateKey(const QString keyname, const double &value, const QString &unit="") const
Checks value of key, produces appropriate value.
Definition Gruen.h:323
virtual AutoReg::RegisterStatus Registration(Chip &sChip, Chip &pChip, Chip &fChip, int startSamp, int startLine, int endSamp, int endLine, int bestSamp, int bestLine)
Applies the adaptive Gruen algorithm to pattern and search chips.
Definition Gruen.cpp:680
PvlGroup StatsLog() const
Create a PvlGroup with the Gruen specific statistics.
Definition Gruen.cpp:853
Radiometric getDefaultRadio() const
Returns the default radiometric gain value.
Definition Gruen.h:342
int algorithm(Chip &pattern, Chip &subsearch, const Radiometric &radio, BigInt &ptsUsed, double &resid, GMatrix &atai, AffineRadio &affrad)
Real workhorse of the computational Gruen algorithm.
Definition Gruen.cpp:156
int CheckConstraints(MatchPoint &point)
Test user limits/contraints after the algorithm has converged.
Definition Gruen.cpp:1090
ErrorList initErrorList() const
Creates an error list from know Gruen errors.
Definition Gruen.cpp:922
GMatrix Choldc(const GMatrix &a, GVector &p) const
Compute Cholesky solution.
Definition Gruen.cpp:341
PvlKeyword ParameterKey(const QString &keyname, const T &value, const QString &unit="") const
Keyword formatter for Gruen parameters.
Definition Gruen.h:301
double DegreesOfFreedom(const int npts) const
Returns number of degrees of freedom of points.
Definition Gruen.h:347
virtual Pvl AlgorithmStatistics(Pvl &pvl)
Create Gruen error and processing statistics Pvl output.
Definition Gruen.cpp:824
BigInt MinValidPoints(BigInt totalPoints) const
Computes the number of minimum valid points.
Definition Gruen.cpp:1038
const AffineRadio & getAffineRadio() const
Return current state of Affine/Radio state.
Definition Gruen.h:104
AffineTolerance getAffineTolerance() const
Return set of tolerances for affine convergence.
Definition Gruen.cpp:1066
Isis exception class.
Definition IException.h:91
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition IException.h:146
QString toString() const
Returns a string representation of this exception.
Structure containing comprehensive registration info/results.
Definition GruenTypes.h:433
Coordinate getAffinePoint(const Coordinate &coord=Coordinate(0.0, 0.0)) const
Return registration offset of a given chip coordinate from center.
Definition GruenTypes.h:452
Define a point set of left, right and geometry at that location.
Definition GruenTypes.h:171
Contains multiple PvlContainers.
Definition PvlGroup.h:41
Container for cube-like labels.
Definition Pvl.h:119
A single keyword-value pair.
Definition PvlKeyword.h:87
bool hasObject(const QString &name) const
Returns a boolean value based on whether the object exists in the current PvlObject or not.
Definition PvlObject.h:323
@ Traverse
Search child objects.
Definition PvlObject.h:158
void addGroup(const Isis::PvlGroup &group)
Add a group to the object.
Definition PvlObject.h:186
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition PvlObject.h:129
Store for radiometric gain and shift parameters.
Definition GruenTypes.h:206
Compute/test the Affine convergence from given parameters/chip.
Definition GruenTypes.h:348
bool hasConverged(const AffineRadio &affine) const
Determines convergence from an affine/radiometric fit.
Definition GruenTypes.h:363
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition IString.cpp:211
int toInt(const QString &string)
Global function to convert from a string to an integer.
Definition IString.cpp:93
const double Null
Value for an Isis Null pixel.
long long int BigInt
Big int.
Definition Constants.h:49
double toDouble(const QString &string)
Global function to convert from a string to a double.
Definition IString.cpp:149
Namespace for the standard library.
Container for Affine limits parameters.
Definition GruenTypes.h:319
Error analysis of Gruen match point solution.
Definition GruenTypes.h:383
double getEigen() const
Returns the square of the of sum of the squares of eigenvalues.
Definition GruenTypes.h:396
Struct that maintains error counts.
Definition Gruen.h:166