9#include "BundleAdjust.h"
17#include <QCoreApplication>
23#include <boost/lexical_cast.hpp>
24#include <boost/numeric/ublas/io.hpp>
25#include <boost/numeric/ublas/matrix_proxy.hpp>
26#include <boost/numeric/ublas/matrix_sparse.hpp>
27#include <boost/numeric/ublas/vector_proxy.hpp>
30#include "Application.h"
31#include "BundleLidarControlPoint.h"
32#include "BundleMeasure.h"
33#include "BundleObservation.h"
34#include "IsisBundleObservation.h"
35#include "BundleObservationSolveSettings.h"
36#include "BundleResults.h"
37#include "BundleSettings.h"
38#include "BundleSolutionInfo.h"
39#include "BundleTargetBody.h"
41#include "CameraDetectorMap.h"
42#include "CameraDistortionMap.h"
43#include "CameraFocalPlaneMap.h"
44#include "CameraGroundMap.h"
47#include "ControlPoint.h"
48#include "CorrelationMatrix.h"
55#include "MaximumLikelihoodWFunctions.h"
56#include "SpecialPixel.h"
57#include "StatCumProbDistDynCalc.h"
58#include "SurfacePoint.h"
61using namespace boost::numeric::ublas;
78 const char* message) {
92 errlog +=
". (See print.prt for details)";
108 const QString &cnetFile,
109 const QString &cubeList,
121 bundleSettings->controlPointCoordTypeReports()) );
129 m_bundleTargetBody = bundleSettings->bundleTargetBody();
146 const QString &cnetFile,
147 const QString &cubeList,
148 const QString &lidarDataFile,
179 m_bundleTargetBody = bundleSettings->bundleTargetBody();
216 m_bundleTargetBody = bundleSettings->bundleTargetBody();
253 m_bundleTargetBody = bundleSettings->bundleTargetBody();
289 m_bundleTargetBody = bundleSettings->bundleTargetBody();
305 const QString &cubeList,
320 m_bundleTargetBody = bundleSettings->bundleTargetBody();
337 QList<ImageList *> imgLists,
350 m_imageLists = imgLists;
358 foreach (
Image *image, *imgList) {
364 m_bundleTargetBody = bundleSettings->bundleTargetBody();
426 emit(statusUpdate(
"Initialization"));
435 m_iterationSummary =
"";
470 for (
int i = 0; i < numImages; i++) {
482 QString msg =
"In BundleAdjust::init(): image " +
fileName +
"is null." +
"\n";
490 QString msg =
"In BundleAdjust::init(): observation "
491 + observationNumber +
"is null." +
"\n";
498 for (
int i = 0; i < numControlPoints; i++) {
500 if (point->IsIgnored()) {
509 int numMeasures = bundleControlPoint->size();
510 for (
int j = 0; j < numMeasures; j++) {
512 QString cubeSerialNumber = measure->cubeSerialNumber();
516 BundleImageQsp image = observation->imageByCubeSerialNumber(cubeSerialNumber);
518 measure->setParentObservation(observation);
519 measure->setParentImage(image);
520 measure->setSigma(1.4);
527 int numLidarPoints = 0;
531 for (
int i = 0; i < numLidarPoints; i++) {
533 if (lidarPoint->IsIgnored()) {
542 int numMeasures = bundleLidarPoint->size();
543 for (
int j = 0; j < numMeasures; j++) {
545 QString cubeSerialNumber = measure->cubeSerialNumber();
549 BundleImageQsp image = observation->imageByCubeSerialNumber(cubeSerialNumber);
551 measure->setParentObservation(observation);
552 measure->setParentImage(image);
553 measure->setSigma(30.0);
565 lidarPoint->ComputeApriori();
568 bundleLidarPoint->initializeRangeConstraints();
605 if (m_bundleTargetBody->solveMeanRadius() || m_bundleTargetBody->solveTriaxialRadii()) {
607 "is NOT possible and likely increases error in the solve.\n");
610 if (m_bundleTargetBody->solveMeanRadius()){
612 bool isMeanRadiusValid =
true;
613 double localRadius, aprioriRadius;
617 SurfacePoint surfacepoint = point->adjustedSurfacePoint();
619 localRadius = surfacepoint.GetLocalRadius().meters();
620 aprioriRadius = m_bundleTargetBody->meanRadius().meters();
624 if (aprioriRadius >= 2 * localRadius || aprioriRadius <= localRadius / 2) {
625 isMeanRadiusValid =
false;
629 if (!isMeanRadiusValid) {
631 "This can cause a bundle failure.\n");
674 int imagesWithInsufficientMeasures = 0;
675 QString msg =
"Images with one or less measures:\n";
676 int numObservations = m_bundleObservations.size();
677 for (
int i = 0; i < numObservations; i++) {
678 int numImages = m_bundleObservations.at(i)->size();
679 for (
int j = 0; j < numImages; j++) {
682 GetNumberOfValidMeasuresInImage(bundleImage->serialNumber());
684 if (numMeasures > 1) {
688 imagesWithInsufficientMeasures++;
689 msg += bundleImage->fileName() +
": " +
toString(numMeasures) +
"\n";
693 if ( imagesWithInsufficientMeasures > 0 ) {
716 cholmod_l_start(&m_cholmodCommon);
722 m_cholmodCommon.nmethods = 1;
723 m_cholmodCommon.method[0].ordering = CHOLMOD_AMD;
741 cholmod_l_free_factor(&
m_L, &m_cholmodCommon);
743 cholmod_l_finish(&m_cholmodCommon);
761 int nBlockColumns = m_bundleObservations.size();
776 for (
int i = 2; i < nBlockColumns; i++) {
783 for (
int i = 0; i < nBlockColumns; i++) {
832 emit(statusBarUpdate(
"Solving"));
851 if (m_bundleTargetBody && m_bundleTargetBody->solveMeanRadius()) {
853 for (
int i = 0; i < numControlPoints; i++) {
855 SurfacePoint surfacepoint = point->adjustedSurfacePoint();
857 surfacepoint.ResetLocalRadius(m_bundleTargetBody->meanRadius());
859 point->setAdjustedSurfacePoint(surfacepoint);
864 if (m_bundleTargetBody && m_bundleTargetBody->solveTriaxialRadii()
867 for (
int i = 0; i < numControlPoints; i++) {
869 SurfacePoint surfacepoint = point->adjustedSurfacePoint();
871 Distance localRadius = m_bundleTargetBody->localRadius(surfacepoint.GetLatitude(),
872 surfacepoint.GetLongitude());
873 surfacepoint.ResetLocalRadius(localRadius);
875 point->setAdjustedSurfacePoint(surfacepoint);
882 double previousSigma0 = 0.0;
885 clock_t solveStartClock = clock();
894 emit statusUpdate(
"\n aborting...");
900 emit statusUpdate( QString(
"\nstarting iteration %1\n").arg(
m_iteration) );
902 clock_t iterationStartClock = clock();
918 emit statusUpdate(
"\n aborting...");
934 emit statusUpdate(
"\n aborting...");
946 emit statusUpdate(
"\n aborting...");
953 emit(statusBarUpdate(
"Computing Residuals"));
968 emit statusUpdate(
"\n aborting...");
984 emit statusUpdate(QString(
"Iteration: %1 \n")
986 emit statusUpdate(QString(
"Sigma0: %1 \n")
991 emit statusUpdate(QString(
"Observations: %1 \n")
993 emit statusUpdate(QString(
"Constrained Parameters:%1 \n")
995 emit statusUpdate(QString(
"Unknowns: %1 \n")
997 emit statusUpdate(QString(
"Degrees of Freedom: %1 \n")
1024 emit statusUpdate(
"Bundle has converged\n");
1025 emit statusBarUpdate(
"Converged");
1027 clock_t iterationStopClock = clock();
1029 / (
double)CLOCKS_PER_SEC;
1036 int numConvergedParams = 0;
1038 for (
int ij = 0; ij < numImgParams; ij++) {
1043 numConvergedParams++;
1046 if ( numConvergedParams == numImgParams ) {
1048 emit statusUpdate(
"Bundle has converged\n");
1049 emit statusBarUpdate(
"Converged");
1055 clock_t iterationStopClock = clock();
1057 / (
double)CLOCKS_PER_SEC;
1058 emit statusUpdate( QString(
"End of Iteration %1 \n").arg(
m_iteration) );
1059 emit statusUpdate( QString(
"Elapsed Time: %1 \n").arg(
m_iterationTime,
1066 emit(statusBarUpdate(
"Max Iterations Reached"));
1079 cholmod_l_free_factor(&
m_L, &m_cholmodCommon);
1090 clock_t errorPropStartClock = clock();
1095 emit statusUpdate(
"\n\nError Propagation Complete\n");
1096 clock_t errorPropStopClock = clock();
1098 / (
double)CLOCKS_PER_SEC);
1101 clock_t solveStopClock = clock();
1103 / (
double)CLOCKS_PER_SEC);
1117 emit statusUpdate(
"\nBundle Complete\n");
1123 emit statusUpdate(
"\n aborting...");
1124 emit statusBarUpdate(
"Failed to Converge");
1126 QString msg =
"Could not solve bundle adjust.";
1127 throw IException(e, e.errorType(), msg, _FILEINFO_);
1141 emit(statusBarUpdate(
"Computing Measure Residuals"));
1148 emit(statusBarUpdate(
"Computing Lidar Measure Residuals"));
1171 return bundleSolutionInfo;
1187 emit(statusBarUpdate(
"Forming Normal Equations"));
1188 bool status =
false;
1205 int numTargetBodyParameters =
m_bundleSettings->numberTargetBodyParameters();
1207 coeffTarget.resize(2,numTargetBodyParameters);
1216 coeffPoint3D.clear();
1222 int numObservations = 0;
1223 int numGood3DPoints = 0;
1224 int numRejected3DPoints = 0;
1225 int numConstrainedCoordinates = 0;
1230 for (
int i = 0; i < num3DPoints; i++) {
1231 emit(pointUpdate(i+1));
1234 if (point->isRejected()) {
1235 numRejected3DPoints++;
1246 int numMeasures = point->size();
1247 for (
int j = 0; j < numMeasures; j++) {
1252 if (measure->isRejected()) {
1256 status =
computePartials(coeffTarget, coeffImage, coeffPoint3D, coeffRHS, *measure,
1266 numObservations += 2;
1268 formMeasureNormals(N22, N12, n1, n2, coeffTarget, coeffImage, coeffPoint3D, coeffRHS,
1269 measure->observationIndex());
1281 int numRejectedLidarPoints = 0.0;
1282 int numGoodLidarPoints = 0.0;
1283 numObservations = 0;
1284 numConstrainedCoordinates = 0;
1291 for (
int i = 0; i < numLidarPoints; i++) {
1292 emit(pointUpdate(i+1));
1295 if (point->isRejected()) {
1296 numRejectedLidarPoints++;
1305 int numMeasures = point->size();
1306 for (
int j = 0; j < numMeasures; j++) {
1309 if (measure->isRejected()) {
1313 status =
computePartials(coeffTarget, coeffImage, coeffPoint3D, coeffRHS, *measure,
1322 numObservations += 2;
1324 formMeasureNormals(N22, N12, n1, n2, coeffTarget, coeffImage, coeffPoint3D, coeffRHS,
1325 measure->observationIndex());
1333 numGoodLidarPoints++;
1377 int observationIndex) {
1379 int blockIndex = observationIndex;
1383 int numTargetPartials = coeffTarget.size2();
1393 numTargetPartials, coeffImage.size2());
1394 (*(*
m_sparseNormals[blockIndex])[0]) += prod(trans(coeffTarget),coeffImage);
1397 N12.insertMatrixBlock(0, numTargetPartials, 3);
1398 *N12[0] += prod(trans(coeffTarget), coeffPoint3D);
1401 vector_range<LinearAlgebra::VectorCompressed> n1_range(n1, range(0, numTargetPartials));
1403 n1_range += prod(trans(coeffTarget), coeffRHS);
1407 int numImagePartials = coeffImage.size2();
1411 numImagePartials, numImagePartials);
1413 (*(*
m_sparseNormals[blockIndex])[blockIndex]) += prod(trans(coeffImage), coeffImage);
1416 N12.insertMatrixBlock(blockIndex, numImagePartials, 3);
1417 *N12[blockIndex] += prod(trans(coeffImage), coeffPoint3D);
1420 vector_range<LinearAlgebra::VectorCompressed> vr(
1426 vr += prod(trans(coeffImage), coeffRHS);
1429 N22 += prod(trans(coeffPoint3D), coeffPoint3D);
1432 n2 += prod(trans(coeffPoint3D), coeffRHS);
1461 boost::numeric::ublas::bounded_vector<double, 3> &NIC = bundleControlPoint->nicVector();
1467 int numConstrainedCoordinates = 0;
1471 boost::numeric::ublas::bounded_vector<double, 3> &weights
1472 = bundleControlPoint->weights();
1473 boost::numeric::ublas::bounded_vector<double, 3> &corrections
1474 = bundleControlPoint->corrections();
1476 if (weights(0) > 0.0) {
1477 N22(0,0) += weights(0);
1478 n2(0) += (-weights(0) * corrections(0));
1479 numConstrainedCoordinates++;
1482 if (weights(1) > 0.0) {
1483 N22(1,1) += weights(1);
1484 n2(1) += (-weights(1) * corrections(1));
1485 numConstrainedCoordinates++;
1488 if (weights(2) > 0.0) {
1489 N22(2,2) += weights(2);
1490 n2(2) += (-weights(2) * corrections(2));
1491 numConstrainedCoordinates++;
1500 bundleControlPoint->setAdjustedSurfacePoint(
SurfacePoint);
1506 NIC = prod(N22, n2);
1514 return numConstrainedCoordinates;
1541 boost::numeric::ublas::bounded_vector<double, 3> &NIC = bundleLidarControlPoint->nicVector();
1547 int numConstrainedCoordinates = 0;
1551 boost::numeric::ublas::bounded_vector<double, 3> &weights
1552 = bundleLidarControlPoint->weights();
1553 boost::numeric::ublas::bounded_vector<double, 3> &corrections
1554 = bundleLidarControlPoint->corrections();
1556 if (weights(0) > 0.0) {
1557 N22(0,0) += weights(0);
1558 n2(0) += (-weights(0) * corrections(0));
1559 numConstrainedCoordinates++;
1562 if (weights(1) > 0.0) {
1563 N22(1,1) += weights(1);
1564 n2(1) += (-weights(1) * corrections(1));
1565 numConstrainedCoordinates++;
1568 if (weights(2) > 0.0) {
1569 N22(2,2) += weights(2);
1570 n2(2) += (-weights(2) * corrections(2));
1571 numConstrainedCoordinates++;
1580 bundleLidarControlPoint->setAdjustedSurfacePoint(
SurfacePoint);
1586 NIC = prod(N22, n2);
1594 return numConstrainedCoordinates;
1617 vector<double> &nj) {
1625 if ( !diagonalBlock )
1632 vector<double> weights = m_bundleTargetBody->parameterWeights();
1633 vector<double> corrections = m_bundleTargetBody->parameterCorrections();
1635 int blockSize = diagonalBlock->size1();
1636 for (
int j = 0; j < blockSize; j++) {
1637 if (weights[j] > 0.0) {
1638 (*diagonalBlock)(j,j) += weights[j];
1639 nj[n] -= weights[j] * corrections(j);
1649 observation = m_bundleObservations.at(i-1);
1652 observation = m_bundleObservations.at(i);
1659 int blockSize = diagonalBlock->size1();
1660 for (
int j = 0; j < blockSize; j++) {
1661 if (weights(j) > 0.0) {
1662 (*diagonalBlock)(j,j) += weights(j);
1663 nj[n] -= weights(j) * corrections(j);
1691 QMapIterator< int, LinearAlgebra::Matrix * > N12it(N12);
1693 while ( N12it.hasNext() ) {
1696 int rowIndex = N12it.key();
1699 Q.insertMatrixBlock(rowIndex, 3, N12it.value()->size1());
1701 *(Q[rowIndex]) = prod(N22,trans(*(N12it.value())));
1720 QMapIterator<int, LinearAlgebra::Matrix*> N12it(N12);
1721 QMapIterator<int, LinearAlgebra::Matrix*> Qit(Q);
1724 while ( N12it.hasNext() ) {
1727 int rowIndex = N12it.key();
1730 while ( Qit.hasNext() ) {
1733 int columnIndex = Qit.key();
1735 if ( rowIndex > columnIndex ) {
1743 N12block->size1(), Qblock->size2());
1745 (*(*
m_sparseNormals[columnIndex])[rowIndex]) -= prod(*N12block,*Qblock);
1765 vector<double> &nj) {
1773 QMapIterator<int, LinearAlgebra::Matrix*> Qit(Q);
1775 while ( Qit.hasNext() ) {
1778 int columnIndex = Qit.key();
1785 for (
unsigned i = 0; i < blockProduct.size(); i++) {
1786 nj(numParams+i) += alpha*blockProduct(i);
1805 QString msg =
"CHOLMOD: Failed to load Triplet matrix";
1823 if (m_cholmodCommon.status == CHOLMOD_NOT_POSDEF) {
1824 QString msg =
"Matrix NOT positive-definite: failure at column " +
toString((
int)
m_L->minor);
1832 cholmod_dense *x, *b;
1838 double *px = (
double*)b->x;
1839 for (
int i = 0; i <
m_rank; i++) {
1844 x = cholmod_l_solve(CHOLMOD_A,
m_L, b, &m_cholmodCommon);
1847 double *sx = (
double*)x->x;
1848 for (
int i = 0; i <
m_rank; i++) {
1854 cholmod_l_free_dense(&b, &m_cholmodCommon);
1855 cholmod_l_free_dense(&x, &m_cholmodCommon);
1877 -1, CHOLMOD_REAL, &m_cholmodCommon);
1892 long numEntries = 0;
1895 for (
int columnIndex = 0; columnIndex < numBlockcolumns; columnIndex++) {
1899 if ( !normalsColumn ) {
1900 QString status =
"\nSparseBlockColumnMatrix retrieval failure at column " +
1901 QString::number(columnIndex);
1906 int numLeadingColumns = normalsColumn->startColumn();
1908 QMapIterator< int, LinearAlgebra::Matrix * > it(*normalsColumn);
1910 while ( it.hasNext() ) {
1913 int rowIndex = it.key();
1920 if ( !normalsBlock ) {
1921 QString status =
"\nmatrix block retrieval failure at column ";
1922 status.append(columnIndex);
1923 status.append(
", row ");
1924 status.append(rowIndex);
1926 status =
"Total # of block columns: " + QString::number(numBlockcolumns);
1933 if ( columnIndex == rowIndex ) {
1934 for (
unsigned ii = 0; ii < normalsBlock->size1(); ii++) {
1935 for (
unsigned jj = ii; jj < normalsBlock->size2(); jj++) {
1936 entryValue = normalsBlock->at_element(ii,jj);
1937 int entryColumnIndex = jj + numLeadingColumns;
1938 int entryRowIndex = ii + numLeadingRows;
1941 tripletColumns[numEntries] = entryColumnIndex;
1942 tripletRows[numEntries] = entryRowIndex;
1946 tripletValues[numEntries] = entryValue;
1953 for (
unsigned ii = 0; ii < normalsBlock->size1(); ii++) {
1954 for (
unsigned jj = 0; jj < normalsBlock->size2(); jj++) {
1955 entryValue = normalsBlock->at_element(ii,jj);
1956 int entryColumnIndex = jj + numLeadingColumns;
1957 int entryRowIndex = ii + numLeadingRows;
1960 tripletColumns[numEntries] = entryRowIndex;
1961 tripletRows[numEntries] = entryColumnIndex;
1965 tripletValues[numEntries] = entryValue;
1994 boost::numeric::ublas::symmetric_matrix< double, upper > c = m;
1996 den = m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))
1997 - m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0))
1998 + m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));
2001 if (fabs(den) < 1.0e-100) {
2007 m(0, 0) = (c(1, 1) * c(2, 2) - c(1, 2) * c(2, 1)) * det;
2008 m(0, 1) = (c(0, 2) * c(2, 1) - c(0, 1) * c(2, 2)) * det;
2009 m(0, 2) = (c(0, 1) * c(1, 2) - c(0, 2) * c(1, 1)) * det;
2010 m(1, 1) = (c(0, 0) * c(2, 2) - c(0, 2) * c(2, 0)) * det;
2011 m(1, 2) = (c(0, 2) * c(1, 0) - c(0, 0) * c(1, 2)) * det;
2012 m(2, 2) = (c(0, 0) * c(1, 1) - c(0, 1) * c(1, 0)) * det;
2038 matrix<double> &coeffImage,
2039 matrix<double> &coeffPoint3D,
2040 vector<double> &coeffRHS,
2044 Camera *measureCamera = measure.camera();
2047 int numImagePartials = observation->numberParameters();
2053 coeffImage.resize(2,numImagePartials);
2063 measureCamera->SetImage(measure.sample(), measure.line());
2067 if (measureCamera->GetCameraType() !=
Camera::Csm) {
2071 double computedX, computedY;
2073 &computedX, &computedY,
false))) {
2074 QString msg =
"Unable to map apriori surface point for measure ";
2075 msg += measure.cubeSerialNumber() +
" on point " + point.
id() +
" into focal plane";
2081 observation->computeTargetPartials(coeffTarget, measure,
m_bundleSettings, m_bundleTargetBody);
2084 observation->computeImagePartials(coeffImage, measure);
2090 observation->computePoint3DPartials(coeffPoint3D, measure, coordType);
2093 observation->computeRHSPartials(coeffRHS, measure);
2095 double deltaX = coeffRHS(0);
2096 double deltaY = coeffRHS(1);
2104 double residualR2ZScore = sqrt(deltaX * deltaX + deltaY * deltaY) / sqrt(2.0);
2112 coeffImage *= observationWeight;
2113 coeffPoint3D *= observationWeight;
2114 coeffRHS *= observationWeight;
2117 coeffTarget *= observationWeight;
2129 emit(statusBarUpdate(
"Updating Parameters"));
2135 int numTargetBodyParameters = m_bundleTargetBody->numberParameters();
2137 m_bundleTargetBody->applyParameterCorrections(subrange(
m_imageSolution,0,
2138 numTargetBodyParameters));
2140 t += numTargetBodyParameters;
2144 int numObservations = m_bundleObservations.size();
2145 for (
int i = 0; i < numObservations; i++) {
2150 observation->applyParameterCorrections(subrange(
m_imageSolution,t,t+numParameters));
2155 QSharedPointer<IsisBundleObservation> isisObservation = qSharedPointerDynamicCast<IsisBundleObservation>(observation);
2156 isisObservation->updateBodyRotation();
2166 if (point->isRejected()) {
2171 m_bundleTargetBody);
2179 if (point->isRejected()) {
2184 m_bundleTargetBody);
2197 emit(statusBarUpdate(
"Computing Residuals"));
2219 for (
int i = 0; i < m_bundleObservations.size(); i++) {
2220 vtpv += m_bundleObservations.at(i)->vtpv();
2224 if ( m_bundleTargetBody) {
2225 vtpv += m_bundleTargetBody->vtpv();
2256 std::vector<double> residuals;
2258 residuals.resize(numResiduals);
2261 int residualIndex = 0;
2263 for (
int i = 0; i < numObjectPoints; i++) {
2267 if ( point->isRejected() ) {
2271 int numMeasures = point->numberOfMeasures();
2272 for (
int j = 0; j < numMeasures; j++) {
2276 if ( measure->isRejected() ) {
2280 vx = measure->sampleResidual();
2281 vy = measure->lineResidual();
2283 residuals[residualIndex] = sqrt(vx*vx + vy*vy);
2290 std::sort(residuals.begin(), residuals.end());
2296 int midpointIndex = numResiduals / 2;
2298 if ( numResiduals % 2 == 0 ) {
2299 median = (residuals[midpointIndex - 1] + residuals[midpointIndex]) / 2;
2302 median = residuals[midpointIndex];
2306 for (
int i = 0; i < numResiduals; i++) {
2307 residuals[i] = fabs(residuals[i] - median);
2310 std::sort(residuals.begin(), residuals.end());
2312 if ( numResiduals % 2 == 0 ) {
2313 medianDev = (residuals[midpointIndex - 1] + residuals[midpointIndex]) / 2;
2316 medianDev = residuals[midpointIndex];
2319 QString status =
"\nmedian deviation: ";
2320 status.append(QString(
"%1").arg(medianDev));
2321 status.append(
"\n");
2324 mad = 1.4826 * medianDev;
2327 status.append(QString(
"%1").arg(mad));
2328 status.append(
"\n");
2334 status =
"\nRejection Limit: ";
2336 status.append(
"\n");
2353 int totalNumRejected = 0;
2355 int maxResidualIndex;
2362 int numComingBack = 0;
2367 for (
int i = 0; i < numObjectPoints; i++) {
2370 point->zeroNumberOfRejectedMeasures();
2373 maxResidualIndex = -1;
2376 int numMeasures = point->numberOfMeasures();
2377 for (
int j = 0; j < numMeasures; j++) {
2381 vx = measure->sampleResidual();
2382 vy = measure->lineResidual();
2384 sumSquares = sqrt(vx*vx + vy*vy);
2387 if ( sumSquares <= usedRejectionLimit ) {
2390 if ( measure->isRejected() ) {
2391 QString status =
"Coming back in: ";
2392 status.append(QString(
"%1").arg(point->id().toLatin1().data()));
2393 status.append(
"\r");
2396 m_controlNet->DecrementNumberOfRejectedMeasuresInImage(measure->cubeSerialNumber());
2399 measure->setRejected(
false);
2404 if ( measure->isRejected() ) {
2410 if ( sumSquares > maxResidual ) {
2411 maxResidual = sumSquares;
2412 maxResidualIndex = j;
2417 if ( maxResidual == -1.0 || maxResidual <= usedRejectionLimit ) {
2418 point->setNumberOfRejectedMeasures(numRejected);
2424 if ((numMeasures - (numRejected + 1)) < 2) {
2425 point->setNumberOfRejectedMeasures(numRejected);
2434 rejected->setRejected(
true);
2436 point->setNumberOfRejectedMeasures(numRejected);
2437 m_controlNet->IncrementNumberOfRejectedMeasuresInImage(rejected->cubeSerialNumber());
2441 if ( ( numMeasures-numRejected ) < 2 ) {
2442 point->setRejected(
true);
2443 QString status =
"Rejecting Entire Point: ";
2444 status.append(QString(
"%1").arg(point->id().toLatin1().data()));
2445 status.append(
"\r");
2449 point->setRejected(
false);
2453 int numberRejectedObservations = 2*totalNumRejected;
2455 QString status =
"\nRejected Observations:";
2456 status.append(QString(
"%1").arg(numberRejectedObservations));
2457 status.append(
" (Rejection Limit:");
2458 status.append(QString(
"%1").arg(usedRejectionLimit));
2459 status.append(
")\n");
2464 status =
"\nMeasures that came back: ";
2465 status.append(QString(
"%1").arg(numComingBack));
2466 status.append(
"\n");
2482 if (m_imageLists.count() > 0) {
2483 return m_imageLists;
2490 imgList->append(image);
2493 m_imageLists.append(imgList);
2496 QString msg =
"Invalid image in serial number list\n";
2501 QString msg =
"No images used in bundle adjust\n";
2505 return m_imageLists;
2530 emit(statusBarUpdate(
"Error Propagation"));
2545 QString status =
" Time: ";
2546 status.append(currentTime.c_str());
2547 status.append(
"\n\n");
2551 std::vector< symmetric_matrix<double> > pointCovariances(numObjectPoints,
2552 symmetric_matrix<double>(3));
2553 for (
int d = 0; d < numObjectPoints; d++) {
2554 pointCovariances[d].clear();
2560 b = cholmod_l_zeros (
m_rank, 1, CHOLMOD_REAL, &m_cholmodCommon );
2561 double *pb = (
double*)b->x;
2572 QFile matrixOutput(matrixFile.expanded());
2577 matrixOutput.open(QIODevice::WriteOnly);
2579 QDataStream outStream(&matrixOutput);
2582 int columnIndex = 0;
2585 for (i = 0; i < numBlockColumns; i++) {
2590 numColumns = normalsColumn->numberOfColumns();
2591 int numRows = normalsColumn->numberOfRows();
2593 inverseMatrix.zeroBlocks();
2596 if (normalsColumn->numberOfColumns() == numColumns) {
2597 int numRows = normalsColumn->numberOfRows();
2598 inverseMatrix.insertMatrixBlock(i, numRows, numColumns);
2599 inverseMatrix.zeroBlocks();
2602 numColumns = normalsColumn->numberOfColumns();
2605 inverseMatrix.wipe();
2608 for (j = 0; j < (i+1); j++) {
2610 int numRows = normalsRow->numberOfRows();
2620 for (j = 0; j < numColumns; j++) {
2621 if ( columnIndex > 0 ) {
2622 pb[columnIndex - 1] = 0.0;
2624 pb[columnIndex] = 1.0;
2626 x = cholmod_l_solve ( CHOLMOD_A,
m_L, b, &m_cholmodCommon );
2631 for (k = 0; k < inverseMatrix.size(); k++) {
2634 int sz1 = matrix->size1();
2636 for (
int ii = 0; ii < sz1; ii++) {
2637 (*matrix)(ii,localCol) = px[ii + rp];
2639 rp += matrix->size1();
2645 cholmod_l_free_dense(&x,&m_cholmodCommon);
2650 vector< double > &adjustedSigmas = m_bundleTargetBody->adjustedSigmas();
2651 matrix< double > *targetCovMatrix = inverseMatrix.value(i);
2653 for (
int z = 0; z < numColumns; z++)
2660 observation = m_bundleObservations.at(i-1);
2663 observation = m_bundleObservations.at(i);
2665 vector< double > &adjustedSigmas = observation->adjustedSigmas();
2666 matrix< double > *imageCovMatrix = inverseMatrix.value(i);
2667 for (
int z = 0; z < numColumns; z++) {
2674 outStream << inverseMatrix;
2679 for (j = 0; j < numObjectPoints; j++) {
2680 emit(pointUpdate(j+1));
2682 if ( point->isRejected() ) {
2688 QString status =
"\rError Propagation: Inverse Block ";
2689 status.append(QString::number(i+1));
2690 status.append(
" of ");
2691 status.append(QString::number(numBlockColumns));
2692 status.append(
"; Point ");
2693 status.append(QString::number(j+1));
2694 status.append(
" of ");
2695 status.append(QString::number(numObjectPoints));
2707 boost::numeric::ublas::symmetric_matrix<double> &covariance = pointCovariances[pointIndex];
2718 QMapIterator< int, LinearAlgebra::Matrix * > it(Q);
2719 while ( it.hasNext() ) {
2722 int nKey = it.key();
2730 if ( !secondQBlock ) {
2736 if ( !inverseBlock ) {
2740 T = prod(*inverseBlock, trans(*firstQBlock));
2741 T = prod(*secondQBlock,T);
2751 catch (std::exception &e) {
2753 QString msg =
"Input data and settings are not sufficiently stable "
2754 "for error propagation.";
2764 matrixOutput.close();
2773 cholmod_l_free_dense(&b,&m_cholmodCommon);
2779 status =
"\rFilling point covariance matrices: Time ";
2780 status.append(currentTime.c_str());
2787 for (j = 0; j < numObjectPoints; j++) {
2791 if ( point->isRejected() ) {
2796 status =
"\rError Propagation: Filling point covariance matrices ";
2797 status.append(QString(
"%1").arg(j+1));
2798 status.append(
" of ");
2799 status.append(QString(
"%1").arg(numObjectPoints));
2800 status.append(
"\r");
2805 boost::numeric::ublas::symmetric_matrix<double> &covariance = pointCovariances[pointIndex];
2813 boost::numeric::ublas::symmetric_matrix <double,boost::numeric::ublas::upper> pCovar;
2816 pCovar =
SurfacePoint.GetSphericalMatrix(SurfacePoint::Kilometers);
2820 pCovar =
SurfacePoint.GetRectangularMatrix(SurfacePoint::Kilometers);
2822 pCovar += covariance;
2823 pCovar *= sigma0Squared;
2947 return m_controlNet->Camera(i)->instrumentRotation()->Cache(
"InstrumentPointing");
2962 return m_controlNet->Camera(i)->instrumentPosition()->Cache(
"InstrumentPosition");
2977 QString msg =
"Cannot get model state for image [" +
toString(i) +
2978 "] because it is not a CSM camera model.";
2996 QString iterationNumber;
3005 PvlGroup summaryGroup(iterationNumber);
3011 summaryGroup +=
PvlKeyword(
"Constrained_Point_Parameters",
3013 summaryGroup +=
PvlKeyword(
"Constrained_Image_Parameters",
3016 summaryGroup +=
PvlKeyword(
"Constrained_Target_Parameters",
3019 summaryGroup +=
PvlKeyword(
"Unknown_Parameters",
3021 summaryGroup +=
PvlKeyword(
"Degrees_of_Freedom",
3023 summaryGroup +=
PvlKeyword(
"Rejected_Measures",
3030 summaryGroup +=
PvlKeyword(
"Maximum_Likelihood_Tier: ",
3032 summaryGroup +=
PvlKeyword(
"Median_of_R^2_residuals: ",
3037 summaryGroup +=
PvlKeyword(
"Converged",
"TRUE");
3041 summaryGroup +=
PvlKeyword(
"ErrorPropagationElapsedTime",
3050 std::ostringstream ostr;
3051 ostr << summaryGroup << std::endl;
3052 m_iterationSummary += QString::fromStdString( ostr.str() );
3058 std::cout << summaryGroup << std::endl;
3091 return m_iterationSummary;
3106 printf(
"%s", status.toStdString().c_str());
3108 else if (QCoreApplication::applicationName() !=
"ipce") {
3109 printf(
"%s", status.toStdString().c_str());
3156 QVector<Statistics> rmsImageSampleResiduals(numberImages);
3157 QVector<Statistics> rmsImageLineResiduals(numberImages);
3158 QVector<Statistics> rmsImageResiduals(numberImages);
3161 for (
int i = 0; i < numObjectPoints; i++) {
3165 if (point->isRejected()) {
3169 int numMeasures = point->numberOfMeasures();
3170 for (
int j = 0; j < numMeasures; j++) {
3174 if (measure->isRejected()) {
3178 double sampleResidual = fabs(measure->sampleResidual());
3179 double lineResidual = fabs(measure->lineResidual());
3185 rmsImageSampleResiduals[imageIndex].AddData(sampleResidual);
3186 rmsImageLineResiduals[imageIndex].AddData(lineResidual);
3187 rmsImageResiduals[imageIndex].AddData(lineResidual);
3188 rmsImageResiduals[imageIndex].AddData(sampleResidual);
3192 QVector<Statistics> rmsLidarImageSampleResiduals(numberImages);
3193 QVector<Statistics> rmsLidarImageLineResiduals(numberImages);
3194 QVector<Statistics> rmsLidarImageResiduals(numberImages);
3198 for (
int i = 0; i < numLidarPoints; i++) {
3202 if (point->isRejected()) {
3206 int numMeasures = point->numberOfMeasures();
3207 for (
int j = 0; j < numMeasures; j++) {
3211 if (measure->isRejected()) {
3215 double sampleResidual = fabs(measure->sampleResidual());
3216 double lineResidual = fabs(measure->lineResidual());
3222 rmsLidarImageSampleResiduals[imageIndex].AddData(sampleResidual);
3223 rmsLidarImageLineResiduals[imageIndex].AddData(lineResidual);
3224 rmsLidarImageResiduals[imageIndex].AddData(lineResidual);
3225 rmsLidarImageResiduals[imageIndex].AddData(sampleResidual);
3235 QString minSigmaCoord1PointId =
"";
3238 QString maxSigmaCoord1PointId =
"";
3242 QString minSigmaCoord2PointId =
"";
3245 QString maxSigmaCoord2PointId =
"";
3249 QString minSigmaCoord3PointId =
"";
3252 QString maxSigmaCoord3PointId =
"";
3259 Distance sigmaCoord1Dist, sigmaCoord2Dist, sigmaCoord3Dist;
3264 for (
int i = 0; i < numPoints; i++) {
3268 maxSigmaCoord1Dist = point->adjustedSurfacePoint().GetSigmaDistance(type,
3270 minSigmaCoord1Dist = maxSigmaCoord1Dist;
3272 maxSigmaCoord2Dist = point->adjustedSurfacePoint().GetSigmaDistance(type,
3274 minSigmaCoord2Dist = maxSigmaCoord2Dist;
3276 maxSigmaCoord1PointId = point->id();
3277 maxSigmaCoord2PointId = maxSigmaCoord1PointId;
3278 minSigmaCoord1PointId = maxSigmaCoord1PointId;
3279 minSigmaCoord2PointId = maxSigmaCoord1PointId;
3283 maxSigmaCoord3Dist = point->adjustedSurfacePoint().GetSigmaDistance(type,
3284 SurfacePoint::Three);
3285 minSigmaCoord3Dist = maxSigmaCoord3Dist;
3287 maxSigmaCoord3PointId = maxSigmaCoord1PointId;
3288 minSigmaCoord3PointId = maxSigmaCoord1PointId;
3293 for (
int i = 0; i < numPoints; i++) {
3297 sigmaCoord1Dist = point->adjustedSurfacePoint().GetSigmaDistance(type,
3299 sigmaCoord2Dist = point->adjustedSurfacePoint().GetSigmaDistance(type,
3301 sigmaCoord3Dist = point->adjustedSurfacePoint().GetSigmaDistance(type,
3302 SurfacePoint::Three);
3304 sigmaCoord1Stats.AddData(sigmaCoord1Dist.meters());
3305 sigmaCoord2Stats.AddData(sigmaCoord2Dist.meters());
3306 sigmaCoord3Stats.AddData(sigmaCoord3Dist.meters());
3308 if (sigmaCoord1Dist > maxSigmaCoord1Dist) {
3309 maxSigmaCoord1Dist = sigmaCoord1Dist;
3310 maxSigmaCoord1PointId = point->id();
3312 if (sigmaCoord2Dist > maxSigmaCoord2Dist) {
3313 maxSigmaCoord2Dist = sigmaCoord2Dist;
3314 maxSigmaCoord2PointId = point->id();
3317 if (sigmaCoord3Dist > maxSigmaCoord3Dist) {
3318 maxSigmaCoord3Dist = sigmaCoord3Dist;
3319 maxSigmaCoord3PointId = point->id();
3322 if (sigmaCoord1Dist < minSigmaCoord1Dist) {
3323 minSigmaCoord1Dist = sigmaCoord1Dist;
3324 minSigmaCoord1PointId = point->id();
3326 if (sigmaCoord2Dist < minSigmaCoord2Dist) {
3327 minSigmaCoord2Dist = sigmaCoord2Dist;
3328 minSigmaCoord2PointId = point->id();
3331 if (sigmaCoord3Dist < minSigmaCoord3Dist) {
3332 minSigmaCoord3Dist = sigmaCoord3Dist;
3333 minSigmaCoord3PointId = point->id();
3342 minSigmaCoord1PointId, maxSigmaCoord1PointId);
3345 minSigmaCoord2PointId, maxSigmaCoord2PointId);
3348 minSigmaCoord3PointId, maxSigmaCoord3PointId);
3351 sigmaCoord2Stats.Rms(),
3352 sigmaCoord3Stats.Rms());
3355 rmsImageSampleResiduals.toList(),
3356 rmsImageResiduals.toList());
3359 rmsLidarImageSampleResiduals.toList(),
3360 rmsLidarImageResiduals.toList());
static void Log(PvlGroup &results)
Writes Pvl results to sessionlog and printfile.
QString modelState(int index)
Return the updated model state for the ith cube in the cube list given to the constructor.
bool m_cleanUp
!< If the iteration summaries should be output to the log file.
int m_rank
!< If the serial number lists need to be deleted by the destructor.
SparseBlockMatrix m_sparseNormals
!< The right hand side of the normal equations.
int numberOfImages() const
Returns the number of images.
void productAB(SparseBlockColumnMatrix &A, SparseBlockRowMatrix &B)
Perform the matrix multiplication C = N12 x Q.
bool solveCholesky()
Compute the least squares bundle adjustment solution using Cholesky decomposition.
bool computeBundleStatistics()
Compute Bundle statistics and store them in m_bundleResults.
double m_iterationTime
Time for last iteration.
int formLidarPointNormals(LinearAlgebra::MatrixUpperTriangular &N22, SparseBlockColumnMatrix &N12, LinearAlgebra::Vector &n2, LinearAlgebra::Vector &nj, BundleLidarControlPointQsp &point)
Compute the Q matrix and NIC vector for a control point.
bool formNormalEquations()
Form the least-squares normal equations matrix via cholmod.
SerialNumberList * serialNumberList()
Returns a pointer to the serial number list.
void init(Progress *progress=0)
Initialize all solution parameters.
bool computePartials(LinearAlgebra::Matrix &coeffTarget, LinearAlgebra::Matrix &coeffImage, LinearAlgebra::Matrix &coeffPoint3D, LinearAlgebra::Vector &coeffRHS, BundleMeasure &measure, BundleControlPoint &point)
Compute partial derivatives and weighted residuals for a measure.
bool isConverged()
Returns if the BundleAdjust converged.
QString fileName(int index)
Return the ith filename in the cube list file given to constructor.
void abortBundle()
Flag to abort when bundle is threaded.
QString m_lidarFileName
Input lidar point filename.
int m_previousNumberImagePartials
!< The image parameter solution vector.
double iteration() const
Returns what iteration the BundleAdjust is currently on.
void outputBundleStatus(QString status)
Slot for deltack and jigsaw to output the bundle status.
bool flagOutliers()
Flags outlier measures and control points.
QVector< BundleControlPointQsp > m_bundleControlPoints
Vector of control points.
LidarDataQsp lidarData()
Returns a pointer to the output lidar data file.
bool isAborted()
Returns if the BundleAdjust has been aborted.
BundleSettingsQsp m_bundleSettings
Contains the solve settings.
bool invert3x3(LinearAlgebra::MatrixUpperTriangular &m)
Dedicated quick inverse of 3x3 matrix.
int m_numLidarConstraints
TODO: temp.
Table cMatrix(int index)
Return the updated instrument pointing table for the ith cube in the cube list given to the construct...
double computeVtpv()
Computes vtpv, the weighted sum of squares of residuals.
QList< ImageList * > imageLists()
This method returns the image list used in the bundle adjust.
void applyParameterCorrections()
Apply parameter corrections for current iteration.
bool initializeCHOLMODLibraryVariables()
Initializations for CHOLMOD sparse matrix package.
boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper, boost::numeric::ublas::column_major > m_normalInverse
!< The lists of images used in the bundle.
LinearAlgebra::Vector m_imageSolution
!< The lower triangular L matrix from Cholesky decomposition.
cholmod_triplet * m_cholmodTriplet
!< The sparse block normal equations matrix.
ControlNetQsp m_controlNet
Output control net.
bool validateNetwork()
control network validation - on the very real chance that the net has not been checked before running...
SerialNumberList * m_serialNumberList
!< Vector of observations.
bool loadCholmodTriplet()
Load sparse normal equations matrix into CHOLMOD triplet.
LidarDataQsp m_lidarDataSet
Output lidar data.
cholmod_factor * m_L
!< The CHOLMOD sparse normal equations matrix used by cholmod_factorize to solve the system.
bool solveSystem()
Compute the solution to the normal equations using the CHOLMOD library.
BundleAdjust(BundleSettingsQsp bundleSettings, const QString &cnetFile, const QString &cubeList, bool printSummary=true)
Construct a BundleAdjust object from the given settings, control network file, and cube list.
void accumProductAlphaAB(double alpha, SparseBlockRowMatrix &A, LinearAlgebra::Vector &B, LinearAlgebra::Vector &C)
Performs the matrix multiplication nj = nj + alpha (Q x n2).
bool initializeNormalEquationsMatrix()
Initialize Normal Equations matrix (m_sparseNormals).
QString iterationSummaryGroup() const
Returns the iteration summary string.
LinearAlgebra::Vector m_RHS
!< Contains object parameters, statistics, and workspace used by the CHOLMOD library.
bool errorPropagation()
Error propagation for solution.
BundleResults m_bundleResults
Stores the results of the bundle adjust.
QString m_cnetFileName
The control net filename.
bool formMeasureNormals(LinearAlgebra::MatrixUpperTriangular &N22, SparseBlockColumnMatrix &N12, LinearAlgebra::VectorCompressed &n1, LinearAlgebra::Vector &n2, LinearAlgebra::Matrix &coeffTarget, LinearAlgebra::Matrix &coeffImage, LinearAlgebra::Matrix &coeffPoint3D, LinearAlgebra::Vector &coeffRHS, int observationIndex)
Form the auxilary normal equation matrices for a measure.
BundleSolutionInfo * bundleSolveInformation()
Create a BundleSolutionInfo containing the settings and results from the bundle adjustment.
int formPointNormals(LinearAlgebra::MatrixUpperTriangular &N22, SparseBlockColumnMatrix &N12, LinearAlgebra::Vector &n2, LinearAlgebra::Vector &nj, BundleControlPointQsp &point)
Compute the Q matrix and NIC vector for a control point.
bool m_abort
!< Contains information about the target body.
BundleLidarPointVector m_bundleLidarControlPoints
Vector of lidar points.
void computeResiduals()
Compute image measure residuals.
bool m_printSummary
!< Summary of the most recently completed iteration.
bool computeRejectionLimit()
Compute rejection limit.
bool productATransB(LinearAlgebra::MatrixUpperTriangular &N22, SparseBlockColumnMatrix &N12, SparseBlockRowMatrix &Q)
Perform the matrix multiplication Q = N22 x N12(transpose)
void iterationSummary()
Creates an iteration summary and an iteration group for the solution summary.
int m_iteration
The current iteration.
bool freeCHOLMODLibraryVariables()
Free CHOLMOD library variables.
bool formWeightedNormals(LinearAlgebra::VectorCompressed &n1, LinearAlgebra::Vector &nj)
Apply weighting for spacecraft position, velocity, acceleration and camera angles,...
BundleSolutionInfo * solveCholeskyBR()
Compute the least squares bundle adjustment solution using Cholesky decomposition.
ControlNetQsp controlNet()
Returns a pointer to the output control network.
cholmod_sparse * m_cholmodNormal
!< The CHOLMOD triplet representation of the sparse normal equations matrix.
~BundleAdjust()
Destroys BundleAdjust object, deallocates pointers (if we have ownership), and frees variables from c...
Table spVector(int index)
Return the updated instrument position table for the ith cube in the cube list given to the construct...
This class holds information about a control point that BundleAdjust needs to run correctly.
SurfacePoint adjustedSurfacePoint() const
Accesses the adjusted SurfacePoint associated with this BundleControlPoint.
QString id() const
Accesses the Point ID associated with this BundleControlPoint.
This class hold image information that BundleAdjust needs to run correctly.Definition for a BundleIma...
This class holds information about a lidar control point that BundleAdjust requires.
double vtpvRangeContribution()
Compute vtpv of lidar range constraints.
A container class for a ControlMeasure.
BundleObservationQsp addNew(BundleImageQsp image, QString observationNumber, QString instrumentId, BundleSettingsQsp bundleSettings)
Adds a new BundleObservation to this vector or fetches an existing BundleObservation if this vector a...
int numberParameters()
Returns the sum of the position parameters and pointing parameters for the contained BundleObservatio...
BundleObservationQsp observationByCubeSerialNumber(QString cubeSerialNumber)
Accesses a BundleObservation associated with the passed serial number.
void addProbabilityDistributionObservation(double obsValue)
Adds an observation to the cumulative probability distribution of |R^2 residuals|.
void incrementNumberConstrainedImageParameters(int incrementAmount)
Increase the number of constrained image parameters.
double elapsedTimeErrorProp() const
Returns the elapsed time for error propagation.
int numberRejectedObservations() const
Returns the number of observation that were rejected.
double rejectionLimit() const
Returns the rejection limit.
int numberObservations() const
Returns the number of observations.
void setElapsedTime(double time)
Sets the elapsed time for the bundle adjustment.
MaximumLikelihoodWFunctions maximumLikelihoodModelWFunc(int modelIndex) const
Returns the maximum likelihood model at the given index.
void setBundleControlPoints(QVector< BundleControlPointQsp > controlPoints)
Sets the bundle control point vector.
void setNumberConstrainedPointParameters(int numberParameters)
Set number of contrained point parameters.
void setBundleLidarPoints(QVector< BundleLidarControlPointQsp > lidarPoints)
Sets the bundle lidar point vector.
double sigma0() const
Returns the Sigma0 of the bundle adjustment.
void setIterations(int iterations)
Sets the number of iterations taken by the BundleAdjust.
bool converged() const
Returns whether or not the bundle adjustment converged.
void setSigmaCoord1Range(Distance minCoord1Dist, Distance maxCoord1Dist, QString minCoord1PointId, QString maxCoord1PointId)
Sets the min and max sigma distances and point ids for coordinate 1.
void incrementMaximumLikelihoodModelIndex()
Increases the value that indicates which stage the maximum likelihood adjustment is currently on.
void setOutputControlNet(ControlNetQsp outNet)
Sets the output ControlNet.
void setSigmaCoord2Range(Distance minCoord2Dist, Distance maxCoord2Dist, QString minCoord2PointId, QString maxCoord2PointId)
Sets the min and max sigma distances and point ids for coordinate 2.
int numberConstrainedPointParameters() const
Returns the number of constrained point parameters.
void resetNumberConstrainedImageParameters()
Resets the number of constrained image parameters to 0.
void setNumberRejectedObservations(int numberObservations)
Sets the number of rejected observations.
void maximumLikelihoodSetUp(QList< QPair< MaximumLikelihoodWFunctions::Model, double > > modelsWithQuantiles)
This method steps up the maximum likelihood estimation solution.
int numberUnknownParameters() const
Returns the number of unknown parameters.
void addResidualsProbabilityDistributionObservation(double obsValue)
Adds an observation to the cumulative probability distribution of residuals used for reporting.
int numberMaximumLikelihoodModels() const
Returns how many maximum likelihood models were used in the bundle adjustment.
void setNumberLidarRangeConstraints(int numberLidarRangeConstraints)
Sets the total number of lidar range constraints.
int maximumLikelihoodModelIndex() const
Returns which step the bundle adjustment is on.
void setRejectionLimit(double rejectionLimit)
Sets the rejection limit.
void incrementNumberConstrainedTargetParameters(int incrementAmount)
Increases the number of constrained target parameters.
void setOutputLidarData(LidarDataQsp outLidarData)
Sets the output LidarData object.
void setNumberImageObservations(int numberObservations)
Sets the number of photogrammetric image observations.
void setRmsImageResidualLists(QList< Statistics > rmsImageLineResiduals, QList< Statistics > rmsImageSampleResiduals, QList< Statistics > rmsImageResiduals)
Sets the root mean square image residual Statistics lists.
double maximumLikelihoodMedianR2Residuals() const
Returns the median of the |R^2 residuals|.
void setNumberLidarImageObservations(int numberLidarObservations)
Sets the number of lidar observations.
int numberConstrainedTargetParameters() const
Return the number of constrained target parameters.
double elapsedTime() const
Returns the elapsed time for the bundle adjustment.
void setNumberUnknownParameters(int numberParameters)
Sets the total number of parameters to solve for.
void setSigmaCoord3Range(Distance minCoord3Dist, Distance maxCoord3Dist, QString minCoord3PointId, QString maxCoord3PointId)
Sets the min and max sigma distances and point ids for coordinate 3.
void setNumberConstrainedLidarPointParameters(int numberParameters)
Set number of contrained point parameters.
void setObservations(BundleObservationVector observations)
Sets the vector of BundleObservations.
void initializeResidualsProbabilityDistribution(unsigned int nodes=20)
Initializes or resets the cumulative probability distribution of residuals used for reporting.
void setRmsFromSigmaStatistics(double rmsFromSigmaCoord1Stats, double rmsFromSigmaCoord2Stats, double rmsFromSigmaCoord3Stats)
Sets the root mean square values of the adjusted sigmas for all three coordinates.
void setRmsXYResiduals(double rx, double ry, double rxy)
Sets the root mean square of the x and y residuals.
void printMaximumLikelihoodTierInformation()
Prints out information about which tier the solution is in and the status of the residuals.
void resizeSigmaStatisticsVectors(int numberImages)
Resizes all image sigma vectors.
void setCorrMatCovFileName(FileName name)
Set the covariance file name for the matrix used to calculate the correlation matrix.
void setRmsLidarImageResidualLists(QList< Statistics > rmsLidarImageLineResiduals, QList< Statistics > rmsLidarImageSampleResiduals, QList< Statistics > rmsLidarImageResiduals)
Sets the root mean square lidar image residual Statistics lists.
void computeSigma0(double dvtpv, BundleSettings::ConvergenceCriteria criteria)
Computes the sigma0 and stores it internally.
int numberConstrainedImageParameters() const
Returns the number of constrained image parameters.
void setElapsedTimeErrorProp(double time)
Sets the elapsed time for error propegation.
void setConverged(bool converged)
Sets if the bundle adjustment converged.
void resetNumberConstrainedTargetParameters()
Resets the number of constrained target parameters to 0.
int degreesOfFreedom() const
Returns the degrees of freedom.
@ Sigma0
The value of sigma0 will be used to determine that the bundle adjustment has converged.
Container class for BundleAdjustment results.
void setRunTime(QString runTime)
Sets the run time, and the name if a name is not already set.
QString getModelState() const
Get the CSM Model state string to re-create the CSM Model.
@ Csm
Community Sensor Model Camera.
This represents an ISIS control net in a project-based GUI interface.
QString fileName() const
Access the name of the control network file associated with this Control.
Status ComputeApriori()
Computes a priori lat/lon/radius point coordinates by determining the average lat/lon/radius of all m...
Distance measurement, usually in meters.
File name manipulation and expansion.
@ User
A type of error that could only have occurred due to a mistake on the user's part (e....
@ Programmer
This error is for when a programmer made an API call that was illegal.
This represents a cube in a project-based GUI interface.
void closeCube()
Cleans up the Cube pointer.
QString fileName() const
Get the file name of the cube that this image represents.
Internalizes a list of images and allows for operations on the entire list.
boost::numeric::ublas::compressed_vector< double > VectorCompressed
Definition for an Isis::LinearAlgebra::VectorCompressed of doubles.
boost::numeric::ublas::vector< double > Vector
Definition for an Isis::LinearAlgebra::Vector of doubles.
boost::numeric::ublas::matrix< double > Matrix
Definition for an Isis::LinearAlgebra::Matrix of doubles.
boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper > MatrixUpperTriangular
Definition for an Isis::LinearAlgebra::MatrixUpperTriangular of doubles with an upper configuration.
double sqrtWeightScaler(double residualZScore)
This provides the scaler to the sqrt of the weight, which is very useful for building normal equation...
Program progress reporter.
Contains multiple PvlContainers.
A single keyword-value pair.
Serial Number list generator.
QString observationNumber(int index)
Return a observation number given an index.
QString serialNumber(const QString &filename)
Return a serial number given a filename.
void add(const QString &filename, bool def2filename=false)
Adds a new filename / serial number pair to the SerialNumberList.
QString spacecraftInstrumentId(int index)
Return the spacecraftname/instrumentid at the given index.
int size() const
How many serial number / filename combos are in the list.
int serialNumberIndex(const QString &sn)
Return a list index given a serial number.
QString fileName(const QString &sn)
Return a filename given a serial number.
void zeroBlocks()
Sets all elements of all matrix blocks to zero.
int numberOfBlocks()
Returns total number of blocks in matrix.
void wipe()
Deletes all pointer elements and removes them from the list.
bool insertMatrixBlock(int nColumnBlock, int nRowBlock, int nRows, int nCols)
Inserts a "newed" boost LinearAlgebra::Matrix pointer of size (nRows, nCols) into the matrix at nColu...
LinearAlgebra::Matrix * getBlock(int column, int row)
Returns pointer to boost matrix at position (column, row).
bool setNumberOfColumns(int n)
Initializes number of columns (SparseBlockColumnMatrix).
int numberOfElements()
Returns number of matrix elements in matrix.
void zeroBlocks()
Sets all elements of all matrix blocks to zero.
This class is used to accumulate statistics on double arrays.
This class defines a body-fixed surface point.
CoordinateType
Defines the coordinate typ, units, and coordinate index for some of the output methods.
@ Latitudinal
Planetocentric latitudinal (lat/lon/rad) coordinates.
@ Rectangular
Body-fixed rectangular x/y/z coordinates.
void SetMatrix(CoordinateType type, const boost::numeric::ublas::symmetric_matrix< double, boost::numeric::ublas::upper > &covar)
Set the covariance matrix.
Class for storing Table blobs information.
static QString CurrentLocalTime()
Returns the current local time This time is taken directly from the system clock, so if the system cl...
This is free and unencumbered software released into the public domain.
QSharedPointer< LidarData > LidarDataQsp
Definition for a shared pointer to a LidarData object.
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
static void cholmodErrorHandler(int nStatus, const char *file, int nLineNo, const char *message)
Custom error handler for CHOLMOD.
QSharedPointer< ControlNet > ControlNetQsp
Typedef for QSharedPointer to control network. This typedef is for future implementation of target bo...