Snooping clean (#4565) (ab2249a2) · Commits · aflab / astrogeology / ISIS3

isis/src/control/objs/BundleAdjust/BundleAdjust.cpp

+235 −7

Original line number	Diff line number	Diff line
		@@ -24,6 +24,7 @@ find files of those names at the top level of this repository. **/
		#include <boost/numeric/ublas/io.hpp>
		#include <boost/numeric/ublas/matrix_sparse.hpp>
		#include <boost/numeric/ublas/vector_proxy.hpp>
		#include <boost/math/distributions/fisher_f.hpp>

		// Isis lib
		#include "Application.h"
		@@ -850,12 +851,22 @@ namespace Isis {
		// compute residuals
		vtpv = computeResiduals();

		// flag outliers
		// variance of unit weight (also reference variance, variance factor, etc.)
		m_bundleResults.computeSigma0(vtpv, m_bundleSettings->convergenceCriteria());

		// flag outliers based on median absolute deviation
		if ( m_bundleSettings->outlierRejection() ) {
		computeRejectionLimit();
		flagOutliers();
		}

		// flag outliers based on normalized residuals
		if ( m_bundleSettings->grossOutlierRejection() ) {
		computeGrossRejectionLimit();
		computeGrossOutlierTestStatistics();
		flagGrossOutliers();
		}

		// testing
		if (m_abort) {
		m_bundleResults.setConverged(false);
		@@ -865,9 +876,6 @@ namespace Isis {
		}
		// testing

		// variance of unit weight (also reference variance, variance factor, etc.)
		m_bundleResults.computeSigma0(vtpv, m_bundleSettings->convergenceCriteria());

		// Set up formatting for status updates with doubles (e.g. Sigma0, Elapsed Time)
		int fieldWidth = 20;
		char format = 'f';
		@@ -1895,6 +1903,7 @@ namespace Isis {

		if (m_bundleSettings->solveTargetBody()) {
		observation->computeTargetPartials(coeffTarget, measure, m_bundleSettings, m_bundleTargetBody);
		// measure->setTargetPartial(&coeffTarget);
		}

		observation->computeImagePartials(coeffImage, measure);
		@@ -2365,6 +2374,225 @@ namespace Isis {
		return true;
		}

		/**
		* @brief Compute rejection limit for normalized residual outlier rejection.
		*
		* Computes a critical value from the f-distribution based on the systems degrees of freedom
		* and probability level. Then, sets the rejection limit in m_bundleResults.
		*
		* @return @b bool If the rejection limit was successfully computed and set.
		*
		* @TODO should this be in BundleResults?
		*
		*/
		bool BundleAdjust::computeGrossRejectionLimit() {
		int dof = m_bundleResults.degreesOfFreedom();
		double a0 = m_bundleSettings->grossOutlierProbabilityLevel();

		// define the cdf distribution as a function of x
		auto fDistributionProbability = [&] (double x) {
		// define the f distribution for a single point against entire observation set
		boost::math::fisher_f_distribution<double> fDist(1, dof-1);
		return boost::math::cdf(fDist, x) - a0;
		};

		// define termination condition
		auto terminationCondition = [] (double min, double max) {
		return abs(min - max) <= 0.0001;
		};

		// feed cdf distirbution function-a0 into bisection algorithm
		using boost::math::tools::bisect;
		std::pair<double, double> result = bisect(fDistributionProbability, 0.001, 100.0, terminationCondition);
		double c = (result.first + result.second) / 2;

		QString status = "\nProbability level: ";
		status.append(QString("%1").arg(a0));
		status.append("\n");
		outputBundleStatus(status);

		m_bundleResults.setGrossRejectionLimit(c);

		status = "\nGross Rejection Limit: ";
		status.append(QString("%1").arg(m_bundleResults.grossRejectionLimit()));
		status.append("\n");
		outputBundleStatus(status);

		return true;
		}


		bool BundleAdjust::computeGrossOutlierTestStatistics() {
		int numObjectPoints = m_bundleControlPoints.size();
		double sigma0 = m_bundleResults.sigma0();

		for (int i = 0; i < numObjectPoints; i++) {
		BundleControlPointQsp point = m_bundleControlPoints.at(i);

		int numMeasures = point->numberOfMeasures();
		for (int j = 0; j < numMeasures; j++) {
		BundleMeasureQsp measure = point->at(j);

		// compute ai*xi
		static LinearAlgebra::Vector aixi(2);

		// compute qvv=(1-ak*xi)/pi
		static LinearAlgebra::Vector qvv(2);

		// calculate test statistics
		double vx = measure->sampleResidual();
		double vy = measure->lineResidual();
		//double wx = fabs(vx)/(sigma0*sqrt(qvv[0]));
		//double wy = fabs(vy)/(sigma0*sqrt(qvv[1]));
		// attach wi to BundleMeasure -> grossOutlierTestStatistic
		//measure.setGrossOutlierTestStatistic(wx,wy);
		}
		}

		return true;
		}


		/**
		* Flags outlier measures and control points by comparing
		* the observations normalized residual to a critical
		* value determined by F-distribution and a user-defined
		* probability level.
		*
		* @return @b bool If the flagging was successful.
		*
		* @TODO How should we handle points with few measures.
		*/
		bool BundleAdjust::flagGrossOutliers() {
		int numRejected;
		int totalNumRejected = 0;

		int maxTestStatisticIndex;
		double maxTestStatistic;

		double wx, wy;
		double usedRejectionLimit = m_bundleResults.grossRejectionLimit();


		// TODO What to do if usedRejectionLimit is too low?

		int numComingBack = 0;

		int numObjectPoints = m_bundleControlPoints.size();

		outputBundleStatus("\n");
		for (int i = 0; i < numObjectPoints; i++) {
		BundleControlPointQsp point = m_bundleControlPoints.at(i);

		point->zeroNumberOfRejectedMeasures();

		numRejected = 0;
		maxTestStatisticIndex = -1;
		maxTestStatistic = -1.0;

		int numMeasures = point->numberOfMeasures();
		for (int j = 0; j < numMeasures; j++) {

		BundleMeasureQsp measure = point->at(j);

		// grab test statistics
		//wx = measure.grossOutlierLineTestStatistic();
		//wy = measure.grossOutlierSampleTestStatistic();
		wx = 1.34;
		wy = 2.58;

		// measure is good
		if ( wx <= usedRejectionLimit &&
		wy <= usedRejectionLimit) {

		// was it previously rejected?
		if ( measure->isRejected() ) {
		QString status = "Coming back in: ";
		status.append(QString("%1").arg(point->id().toLatin1().data()));
		status.append("\r");
		outputBundleStatus(status);
		numComingBack++;
		m_controlNet->DecrementNumberOfRejectedMeasuresInImage(measure->cubeSerialNumber());
		}

		measure->setRejected(false);
		continue;
		}

		// if it's still rejected, skip it
		if ( measure->isRejected() ) {
		numRejected++;
		totalNumRejected++;
		continue;
		}

		if ( wx > maxTestStatistic ) {
		maxTestStatistic = wx;
		maxTestStatisticIndex = j;
		}
		if ( wy > maxTestStatistic ) {
		maxTestStatistic = wy;
		maxTestStatisticIndex = j;
		}

		}

		// no observations above the current rejection limit for this 3D point
		if ( maxTestStatistic == -1.0 \|\| maxTestStatistic <= usedRejectionLimit ) {
		point->setNumberOfRejectedMeasures(numRejected);
		continue;
		}

		// this is another kluge - if we only have two observations
		// we won't reject (for now)
		if ((numMeasures - (numRejected + 1)) < 2) {
		point->setNumberOfRejectedMeasures(numRejected);
		continue;
		}

		// otherwise, we have at least one observation
		// for this point whose residual is above the
		// current rejection limit - we'll flag the
		// worst of these as rejected
		BundleMeasureQsp rejected = point->at(maxTestStatisticIndex);
		rejected->setRejected(true);
		numRejected++;
		point->setNumberOfRejectedMeasures(numRejected);
		m_controlNet->IncrementNumberOfRejectedMeasuresInImage(rejected->cubeSerialNumber());
		totalNumRejected++;

		// do we still have sufficient remaining observations for this 3D point?
		if ( ( numMeasures-numRejected ) < 2 ) {
		point->setRejected(true);
		QString status = "Rejecting Entire Point: ";
		status.append(QString("%1").arg(point->id().toLatin1().data()));
		status.append("\r");
		outputBundleStatus(status);
		}
		else
		point->setRejected(false);

		}

		int numberRejectedObservations = 2*totalNumRejected;

		QString status = "\nRejected Observations:";
		status.append(QString("%1").arg(numberRejectedObservations));
		status.append(" (Rejection Limit:");
		status.append(QString("%1").arg(usedRejectionLimit));
		status.append(")\n");
		outputBundleStatus(status);

		m_bundleResults.setNumberRejectedObservations(numberRejectedObservations);

		status = "\nMeasures that came back: ";
		status.append(QString("%1").arg(numComingBack));
		status.append("\n");
		outputBundleStatus(status);

		return true;
		}


		/**
		* This method returns the image list used in the bundle adjust. If a QList<ImageList *> was passed

isis/src/control/objs/BundleAdjust/BundleAdjust.h

+3 −0

Original line number	Diff line number	Diff line
		@@ -389,7 +389,10 @@ namespace Isis {
		bool errorPropagation();
		double computeResiduals();
		bool computeRejectionLimit();
		bool computeGrossRejectionLimit();
		bool computeGrossOutlierTestStatistics();
		bool flagOutliers();
		bool flagGrossOutliers();

		// normal equation matrices methods

isis/src/control/objs/BundleUtilities/BundleMeasure.cpp

+127 −10

Original line number	Diff line number	Diff line
		@@ -36,6 +36,18 @@ namespace Isis {
		* Destructor
		*/
		BundleMeasure::~BundleMeasure() {
		if (m_imagePartial) {
		delete m_imagePartial;
		m_imagePartial = nullptr;
		}
		if (m_point3DPartial) {
		delete m_point3DPartial;
		m_point3DPartial = nullptr;
		}
		if (m_targetPartial) {
		delete m_targetPartial;
		m_targetPartial = nullptr;
		}
		}


		@@ -108,6 +120,61 @@ namespace Isis {
		}


		/**
		* Sets the BundleMeasure's image partial
		*
		* @param imagePartial BundleMeasure's imagePartial matrix
		*/
		void BundleMeasure::setImagePartial(const LinearAlgebra::Matrix &imagePartial) {
		if (m_imagePartial) {
		delete m_imagePartial;
		m_imagePartial = nullptr;
		}

		m_imagePartial = new LinearAlgebra::Matrix(imagePartial);
		}


		/**
		* Sets the BundleMeasure's 3D point partial
		*
		* @param point3DPartial BundleMeasure's point3DPartial matrix
		*/
		void BundleMeasure::setPoint3DPartial(const LinearAlgebra::Matrix &point3DPartial) {
		if (m_point3DPartial) {
		delete m_point3DPartial;
		m_point3DPartial = nullptr;
		}

		m_point3DPartial = new LinearAlgebra::Matrix(point3DPartial);
		}


		/**
		* Sets the BundleMeasure's target partial
		*
		* @param targetPartial BundleMeasure's targetPartial matrix
		*/
		void BundleMeasure::setTargetPartial(const LinearAlgebra::Matrix &targetPartial) {
		if (m_targetPartial) {
		delete m_targetPartial;
		m_targetPartial = nullptr;
		}

		m_targetPartial = new LinearAlgebra::Matrix(targetPartial);
		}

		/**
		* Sets the test statisitc used for data snooping outlier rejection method
		*
		* @param w vector of sample and line test statistic
		*/
		void BundleMeasure::setGrossOutlierTestStatistic(double sampleOutlierTestStatistic, double lineOutlierTestStatistic) {
		m_sampleOutlierTestStatistic = sampleOutlierTestStatistic;
		m_lineOutlierTestStatistic = lineOutlierTestStatistic;
		}


		/**
		* Determines whether or not this BundleMeasure is rejected
		*
		@@ -194,26 +261,55 @@ namespace Isis {


		/**
		* Accesses the sample residual for this control measure
		* Accesses the current line measurement for this control measure
		*
		* @see ControlMeasure::GetSampleResidual()
		* @see ControlMeasure::GetLine()
		*
		* @return @b double Returns the sample residual
		* @return @b double Returns the line measurement for this control measure
		*/
		double BundleMeasure::sampleResidual() const {
		return m_controlMeasure->GetSampleResidual();
		double BundleMeasure::line() const {
		return m_controlMeasure->GetLine();
		}


		/**
		* Accesses the current line measurement for this control measure
		* Acessess the pointer to the BundleMeasure's image partial
		*
		* @see ControlMeasure::GetLine()
		* @return @b image partial pointer
		*/
		LinearAlgebra::Matrix *BundleMeasure::imagePartial() {
		return m_imagePartial;
		}


		/**
		* Acessess the pointer to the BundleMeasure's point3D partial
		*
		* @return @b double Returns the line measurement for this control measure
		* @return @b 3D point partial pointer
		*/
		double BundleMeasure::line() const {
		return m_controlMeasure->GetLine();
		LinearAlgebra::Matrix *BundleMeasure::point3DPartial() {
		return m_point3DPartial;
		}


		/**
		* Acessess the pointer to the BundleMeasure's target partial
		*
		* @return @b target partial pointer
		*/
		LinearAlgebra::Matrix *BundleMeasure::targetPartial() {
		return m_targetPartial;
		}

		/**
		* Accesses the sample residual for this control measure
		*
		* @see ControlMeasure::GetSampleResidual()
		*
		* @return @b double Returns the sample residual
		*/
		double BundleMeasure::sampleResidual() const {
		return m_controlMeasure->GetSampleResidual();
		}


		@@ -241,6 +337,27 @@ namespace Isis {
		}


		/**
		* Accesses the outlier test statistic used for the data snooping
		* outlier rejection method associate with the BundleMeasure's sample.
		*
		* @return @b BundleMeasure sample's outlier test statistic
		*/
		double BundleMeasure::sampleOutlierTestStatistic() const {
		return m_sampleOutlierTestStatistic;
		}

		/**
		* Accesses the outlier test statistic used for the data snooping
		* outlier rejection method associate with the BundleMeasure's line.
		*
		* @return @b BundleMeasure line's outlier test statistic
		*/
		double BundleMeasure::lineOutlierTestStatistic() const {
		return m_lineOutlierTestStatistic;
		}


		/**
		* Accesses the serial number of the cube containing this control measure
		*

isis/src/control/objs/BundleUtilities/BundleMeasure.h

+19 −1

Original line number	Diff line number	Diff line
		@@ -10,6 +10,7 @@ find files of those names at the top level of this repository. **/
		/* SPDX-License-Identifier: CC0-1.0 */

		#include <QSharedPointer>
		#include "LinearAlgebra.h"

		namespace Isis {
		class BundleControlPoint;
		@@ -68,6 +69,10 @@ namespace Isis {
		void setParentObservation(QSharedPointer<BundleObservation> observation);
		void setParentImage(QSharedPointer<BundleImage> image);
		void setRejected(bool reject);
		void setImagePartial(const LinearAlgebra::Matrix &imagePartials);
		void setPoint3DPartial(const LinearAlgebra::Matrix &point3DPartial);
		void setTargetPartial(const LinearAlgebra::Matrix &targetPartial);
		void setGrossOutlierTestStatistic(double sampleOutlierTestStatistic, double lineOutlierTestStatistic);

		bool isRejected() const;
		Camera *camera() const;
		@@ -77,10 +82,15 @@ namespace Isis {
		const QSharedPointer<BundleObservationSolveSettings> observationSolveSettings();

		double sample() const;
		double sampleResidual() const;
		double line() const;
		LinearAlgebra::Matrix *imagePartial();
		LinearAlgebra::Matrix *point3DPartial();
		LinearAlgebra::Matrix *targetPartial();
		double sampleResidual() const;
		double lineResidual() const;
		double residualMagnitude() const;
		double sampleOutlierTestStatistic() const;
		double lineOutlierTestStatistic() const;
		QString cubeSerialNumber() const;
		double focalPlaneComputedX() const;
		double focalPlaneComputedY() const;
		@@ -94,6 +104,14 @@ namespace Isis {
		bundle control measure **/
		QSharedPointer<BundleImage> m_parentBundleImage; /< Parent image of this bundle control measure /
		QSharedPointer<BundleObservation> m_parentObservation; /< Parent bundle observation /

		LinearAlgebra::Matrix m_imagePartial = nullptr; /< Image partials associated with the control measure */
		LinearAlgebra::Matrix m_point3DPartial = nullptr; /< Point partials associated with the control measure */
		LinearAlgebra::Matrix m_targetPartial = nullptr; /< Target partials associated with the control measure */

		// Test statistic used for determining gross outliers
		double m_sampleOutlierTestStatistic;
		double m_lineOutlierTestStatistic;
		};
		//! Definition for BundleMeasureQsp, a shared pointer to a BundleMeasure.
		typedef QSharedPointer<BundleMeasure> BundleMeasureQsp;