BundleAdjust.h

#ifndef BundleAdjust_h
#define BundleAdjust_h
 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include <QObject> // parent class
 
// std lib
#include <vector>
#include <fstream>
 
// cholmod lib
#include <cholmod.h>
 
// Isis lib
#include "BundleControlPoint.h"
#include "BundleLidarPointVector.h"
#include "BundleObservationSolveSettings.h"
#include "BundleObservationVector.h"
#include "BundleResults.h"
#include "BundleSettings.h"
#include "BundleSolutionInfo.h"
#include "BundleTargetBody.h"
#include "Camera.h"
#include "CameraGroundMap.h"
#include "ControlMeasure.h"
#include "ControlNet.h"
#include "LidarData.h"
#include "LinearAlgebra.h"
#include "MaximumLikelihoodWFunctions.h"
#include "ObservationNumberList.h"
#include "Progress.h"
#include "SerialNumberList.h"
#include "SparseBlockMatrix.h"
#include "Statistics.h"
 
template< typename T > class QList;
template< typename A, typename B > class QMap;
 
namespace Isis {
  class Control;
  class ImageList;
 
  class BundleAdjust : public QObject {
      Q_OBJECT
    public:
      BundleAdjust(BundleSettingsQsp bundleSettings,
                   const QString &cnetFile,
                   const QString &cubeList,
                   bool printSummary = true);
      BundleAdjust(BundleSettingsQsp bundleSettings,
                   const QString &cnetFile,
                   const QString &cubeList,
                   const QString &lidarDataFile,
                   bool printSummary = true);
      BundleAdjust(BundleSettingsQsp bundleSettings,
                   QString &cnet,
                   SerialNumberList &snlist,
                   bool printSummary = true);
      BundleAdjust(BundleSettingsQsp bundleSettings,
                   Control &cnet,
                   SerialNumberList &snlist,
                   bool bPrintSummary);
      BundleAdjust(BundleSettingsQsp bundleSettings,
                   ControlNet &cnet,
                   SerialNumberList &snlist,
                   bool printSummary = true);
      BundleAdjust(BundleSettingsQsp bundleSettings,
                   ControlNetQsp cnet,
                   const QString &cubeList,
                   bool printSummary = true);
      BundleAdjust(BundleSettingsQsp bundleSettings,
                   Control &control,
                   QList<ImageList *> imgList,
                   bool printSummary);
      ~BundleAdjust();
      BundleSolutionInfo*    solveCholeskyBR();
 
      QList<ImageList *> imageLists();
      bool isAborted();
 
    public slots:
      bool solveCholesky();
      void abortBundle();
      void outputBundleStatus(QString status);
 
      // accessors
 
      ControlNetQsp    controlNet();
      LidarDataQsp     lidarData();
      SerialNumberList *serialNumberList();
      QString          fileName(int index);
      QString          iterationSummaryGroup() const;
      bool             isConverged();
      Table            cMatrix(int index);
      Table            spVector(int index);
      QString          modelState(int index);
      int              numberOfImages() const;
      double           iteration() const;
 
    signals:
      void statusUpdate(QString);
      void error(QString);
      void iterationUpdate(int);
      void pointUpdate(int);
      void statusBarUpdate(QString);
      void resultsReady(BundleSolutionInfo *bundleSolveInformation);
      void finished();
 
    private:
      //TODO Should there be a resetBundle(BundleSettings bundleSettings) method
      //     that allows for rerunning with new settings? JWB
      void init(Progress *progress = 0);
      bool initializeNormalEquationsMatrix();
      bool validateNetwork();
      bool solveSystem();
      void iterationSummary();
      BundleSolutionInfo* bundleSolveInformation();
      bool computeBundleStatistics();
      void applyParameterCorrections();
      bool errorPropagation();
      void computeResiduals();
      double computeVtpv();
      bool computeRejectionLimit();
      bool flagOutliers();
 
      // normal equation matrices methods
 
      bool formNormalEquations();
      bool computePartials(LinearAlgebra::Matrix  &coeffTarget,
                           LinearAlgebra::Matrix  &coeffImage,
                           LinearAlgebra::Matrix  &coeffPoint3D,
                           LinearAlgebra::Vector  &coeffRHS,
                           BundleMeasure          &measure,
                           BundleControlPoint     &point);
      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);
      int formPointNormals(LinearAlgebra::MatrixUpperTriangular &N22,
                           SparseBlockColumnMatrix              &N12,
                           LinearAlgebra::Vector                &n2,
                           LinearAlgebra::Vector                &nj,
                           BundleControlPointQsp                &point);
      int formLidarPointNormals(LinearAlgebra::MatrixUpperTriangular &N22,
                                SparseBlockColumnMatrix              &N12,
                                LinearAlgebra::Vector                &n2,
                                LinearAlgebra::Vector                &nj,
                                BundleLidarControlPointQsp           &point);
      bool formWeightedNormals(LinearAlgebra::VectorCompressed  &n1,
                               LinearAlgebra::Vector            &nj);
 
      // dedicated matrix functions
 
      void productAB(SparseBlockColumnMatrix &A,
                     SparseBlockRowMatrix    &B);
      void accumProductAlphaAB(double                alpha,
                               SparseBlockRowMatrix  &A,
                               LinearAlgebra::Vector &B,
                               LinearAlgebra::Vector &C);
      bool invert3x3(LinearAlgebra::MatrixUpperTriangular &m);
      bool productATransB(LinearAlgebra::MatrixUpperTriangular &N22,
                          SparseBlockColumnMatrix              &N12,
                          SparseBlockRowMatrix                 &Q);
 
      // CHOLMOD library methods
 
      bool initializeCHOLMODLibraryVariables();
      bool freeCHOLMODLibraryVariables();
      bool loadCholmodTriplet();
 
      // member variables
 
      BundleSettingsQsp m_bundleSettings;                    
      BundleResults  m_bundleResults;                        
      ControlNetQsp m_controlNet;                            
      QString m_cnetFileName;                                
 
      QVector <BundleControlPointQsp> m_bundleControlPoints; 
      BundleLidarPointVector m_bundleLidarControlPoints;     
 
      QString m_lidarFileName;                               
      LidarDataQsp m_lidarDataSet;                           
      int m_numLidarConstraints;                             
 
      BundleObservationVector m_bundleObservations;          
      SerialNumberList *m_serialNumberList;                  
      BundleTargetBodyQsp m_bundleTargetBody;                
      bool m_abort;                                          
      QString m_iterationSummary;                            
      bool m_printSummary;                                   
      bool m_cleanUp;                                        
      int m_rank;                                            
      int m_iteration;                                       
      double m_iterationTime;                                
      int m_numberOfImagePartials;                           
      QList<ImageList *> m_imageLists;                        
      // ==========================================================================================
      // === BEYOND THIS PLACE (THERE BE DRAGONS) all refers to the folded bundle solution.     ===
      // === Currently, everything uses the CHOLMOD library,                                    ===
      // === there is no dependence on the least-squares class.                                 ===
      // ==========================================================================================
 
      boost::numeric::ublas::symmetric_matrix<
          double,
          boost::numeric::ublas::upper,
          boost::numeric::ublas::column_major > m_normalInverse;
      cholmod_common  m_cholmodCommon;                       
      LinearAlgebra::Vector m_RHS;                           
      SparseBlockMatrix m_sparseNormals;                     
      cholmod_triplet  *m_cholmodTriplet;                    
      cholmod_sparse *m_cholmodNormal;                       
      cholmod_factor *m_L;                                   
      LinearAlgebra::Vector m_imageSolution;                 
      int m_previousNumberImagePartials;                     
  };
}
 
#endif