Added the ability to refit a piecewise polynomial over and existing polynomial. Fixes #5011. (14758c9f) · Commits · aflab / astrogeology / ISIS3

isis/src/base/objs/PiecewisePolynomial/PiecewisePolynomial.cpp

+144 −52

Original line number	Diff line number	Diff line
		@@ -105,7 +105,6 @@ namespace Isis {
		* @return @b std::vector<double> The output values.
		*/
		std::vector<double> PiecewisePolynomial::evaluate(double value) {
		try {
		int knotIndex = segmentIndex(value);

		std::vector<double> output( dimensions() );
		@@ -114,12 +113,6 @@ namespace Isis {
		}
		return output;
		}
		catch (IException &e) {
		QString msg = "Failed evaluating piecewise polynomial at value ["
		+ toString(value) + "].";
		throw IException(e, IException::Programmer, msg, _FILEINFO_);
		}
		}


		/**
		@@ -143,40 +136,75 @@ namespace Isis {

		/**
		* Returns the index of the segment that a given value belongs to.
		* If the value is less than the minimum value or greater than
		* the maximum value, then an error is thrown.
		* If the value is less than the lowest knot, then the first index is returned.
		* If the value is greater than the highest knot, then the last index is returned.
		*
		* @param value The value to find the segment for.
		*
		* @return @b int The zero-based index of the segment.
		*
		* @throws IException IException::Programmer "Value is not within valid range"
		*
		* @TODO Do we really want to throw an error? Maybe we should just return
		* the first/last index? JAM.
		*/
		int PiecewisePolynomial::segmentIndex(double value) {
		if ( value < m_knots.front() - 1e-10 \|\| value > m_knots.back() + 1e-10 ) {
		QString msg = "Value [" + toString(value) + "] is not within valid range [" +
		toString( m_knots.front() ) + ", " + toString( m_knots.back() ) +
		"].";
		throw IException(IException::Programmer, msg, _FILEINFO_);
		std::vector<double>::const_iterator knotIt;
		knotIt = std::upper_bound(m_knots.begin(), m_knots.end(), value);
		if (knotIt == m_knots.begin()) {
		return 0;
		}

		if ( value == m_knots.back() ) {
		else if (knotIt == m_knots.end()) {
		return m_knots.size() - 2;
		}
		else if ( value == m_knots.front() ) {
		return 0;
		}

		std::vector<double>::const_iterator knotIt;
		knotIt = std::upper_bound(m_knots.begin(), m_knots.end(), value);
		int knotIndex = std::distance<std::vector<double>::const_iterator>(m_knots.begin(), knotIt);
		return knotIndex - 1;
		}


		/**
		* Fit a new set of polynomials over the input number of segments. The
		* existing polynomials are sampled to generate a new data set that the new
		* polynomials are fit over.
		*
		* @param segments The number of segments for the new polynomials.
		*
		* @todo Is there a better way to refit the polynomials? Potentially
		* something to do with projecting the existing polynomials onto
		* the set of polynomials over the new knots. This would take some
		* serious calculations.
		*/
		void PiecewisePolynomial::refitPolynomials(int segments, int samples) {

		// Generate sample data.
		std::vector<double> sampleValues;
		std::vector< std::vector<double> > sampleData;
		double firstSample = m_knots.front();
		double lastSample = m_knots.back();

		// If fitting over the zero polynomial, just take two samples
		if ( isZero() ) {
		sampleValues.push_back(firstSample);
		sampleData.push_back(evaluate(firstSample));
		sampleValues.push_back(lastSample);
		sampleData.push_back(evaluate(lastSample));
		}

		// Otherwise sample as normal
		double sampleRate = (lastSample - firstSample) / samples;
		for (int i = 0; i < samples; i++) {
		double currentValue = firstSample + sampleRate * i;
		sampleValues.push_back(currentValue);
		sampleData.push_back(evaluate(currentValue));
		}

		// Fit over the sample data.
		try {
		fitPolynomials(sampleValues, sampleData, segments);
		}
		catch(IException &e) {
		QString msg = "Failed to refit polynomials to [" + toString(segments)
		+ "] segments.";
		throw IException(e, IException::Programmer, msg, _FILEINFO_);
		}
		}


		/**
		* @brief Compute knots and fit polynomials to a data set.
		*
		@@ -211,6 +239,7 @@ namespace Isis {
		*
		* Checks the following:
		* <ul>
		* <li>At least one segment is being used</li>
		* <li>The input values are sorted</li>
		* <li>There are the same number of input values and data points</li>
		* <li>The data points have the correct dimensions</li>
		@@ -234,12 +263,21 @@ namespace Isis {
		void PiecewisePolynomial::validateData(const std::vector<double> &values,
		const std::vector< std::vector<double> > &data,
		int segmentCount) {
		// Check that there is a valid segment count
		if (segmentCount < 1) {
		QString msg = "Cannot fit polynomials over [" + toString(segmentCount)
		+ "] segments. There must be at least one segment.";
		}

		// Check that each value has an associated data point
		if ( values.size() != data.size() ) {
		QString msg = "The number of input values [" + toString( (int)values.size() ) +
		"] and data points [" + toString( (int)data.size() ) + "] do not match.";
		throw IException(IException::Programmer, msg, _FILEINFO_);
		}

		// Check the there is enough data to solve the fit

		// Without continuity conditions, we need 1 data point per coefficient
		int numPointsNeeded = segmentCount * (degree() + 1) * dimensions();

		@@ -254,6 +292,7 @@ namespace Isis {
		throw IException(IException::Programmer, msg, _FILEINFO_);
		}

		// Check that everything is sorted
		for (int i = 0; i < (int)values.size() - 1; i++) {
		if ( values[i + 1] < values[i] ) {
		QString msg = "Input values are not sorted in ascending order.";
		@@ -261,6 +300,7 @@ namespace Isis {
		}
		}

		// Check that each data point is valid
		for (int i = 0; i < (int)data.size(); i++) {
		if ( (int)data[i].size() != dimensions() ) {
		QString msg = "Data point number [" + toString( i + 1 ) +
		@@ -324,6 +364,20 @@ namespace Isis {
		return;
		}

		// If there are only two data points, then use them as the end points and
		// evenly distribute the knots based on value
		if (values.size() == 2) {
		std::vector<double> newKnots(segments + 1);
		newKnots.front() = values.front();
		newKnots.back() = values.back();
		double segmentDistance = ( values.back() - values.front() ) / segments;
		for (int i = 1; i < segments; i++) {
		newKnots[i] = values.front() + segmentDistance * i;
		}
		setKnots(newKnots);
		return;
		}

		// =========================
		// = Setup storage vectors =
		// =========================
		@@ -364,10 +418,23 @@ namespace Isis {
		// ===========================

		double segmentWeight = totalWeight / segments;
		// We always take the first and last times plus padding for the knots.
		// We always take the first and last times for the knots.
		std::vector<double> newKnots(segments + 1);
		newKnots.front() = values.front();
		newKnots.back() = values.back();

		// If there is 0 curvature all along the curve, it's a straight line,
		// then evenly distribute the knots based on the value
		if ( fabs(totalWeight) < 1.0e-15 ) {
		double segmentDistance = ( values.back() - values.front() ) / segments;
		for (int i = 1; i < segments; i++) {
		newKnots[i] = values.front() + segmentDistance * i;
		}
		}

		// Otherwise distribute the knots based on the integral of the curvature
		// over the arc length.
		else {
		for (int i = 1; i < segments; i++) {
		// Find the data point before the knot
		std::vector<double>::iterator timeIt;
		@@ -388,6 +455,7 @@ namespace Isis {
		newKnots[i] = ratio * values[preKnotIndex+1] +
		(1 - ratio) * values[preKnotIndex];
		}
		}

		// Save the new knot locations and reset the polynomials to zero polynomials.
		setKnots(newKnots);
		@@ -736,6 +804,30 @@ namespace Isis {
		}


		/**
		* Convenience method to check if the polynomial is the zero polynomial.
		*
		* @return @b bool If the piecewise polynomial is the zerop polynomial.
		*/
		bool PiecewisePolynomial::isZero() const {
		// Check the coefficients for each segment
		for (int i = 0; i < segments(); i++) {
		std::vector< std::vector<double> > segmentCoeff = coefficients(i);
		// Check the coefficients for each dimension
		for (int j = 0; j < dimensions(); j++) {
		for (int k = 0; k < degree() + 1; k++) {
		if (segmentCoeff[j][k] != 0) {
		return false;
		}
		}
		}
		}

		// Every coefficient has been checked, so return true
		return true;
		}


		/**
		* Sets the degree of the polynomials. All polynomials will be reset to the
		* zero polynomial.

isis/src/base/objs/PiecewisePolynomial/PiecewisePolynomial.h

+5 −0

Original line number	Diff line number	Diff line
		@@ -57,9 +57,12 @@ namespace Isis {

		std::vector<double> evaluate(double value);
		std::vector<double> derivativeVariable(double value);

		void fitPolynomials(const std::vector<double> &values,
		const std::vector< std::vector<double> > &data,
		int segments);
		// TODO Is 100 enough samples to default to? JAM
		void refitPolynomials(int segments, int samples = 100);

		int degree() const;
		std::vector< std::vector<double> > coefficients(int segment) const;
		@@ -67,6 +70,8 @@ namespace Isis {
		const std::vector<double> knots() const;
		int segments() const;

		bool isZero() const;

		void setDegree(int degree);
		void setCoefficients(int segment,
		const std::vector< std::vector<double> > &coefficients);

isis/src/base/objs/PiecewisePolynomial/PiecewisePolynomial.truth

+0 −3

Original line number	Diff line number	Diff line
		@@ -175,9 +175,6 @@ RMS Error: 0.0

		Test error throws

		Attempt to evaluate at value outside of range:
		PROGRAMMER ERROR Failed evaluating piecewise polynomial at value [20.0].
		PROGRAMMER ERROR Value [20.0] is not within valid range [-4.0, 4.0].

		Polynomial fit errors:
		PROGRAMMER ERROR The number of input values [3] and data points [1] do not match.

isis/src/base/objs/PiecewisePolynomial/unitTest.cpp

+0 −8

Original line number	Diff line number	Diff line
		@@ -71,14 +71,6 @@ int main(int argc, char *argv[]) {

		cout << endl << "Test error throws" << endl << endl;

		try {
		cout << "Attempt to evaluate at value outside of range:" << endl;
		testPoly.evaluate(20);
		}
		catch(IException &e) {
		e.print();
		}

		cout << endl << "Polynomial fit errors:" << endl;

		try {

isis/src/base/objs/SpicePosition/SpicePosition.cpp

+3 −0

Original line number	Diff line number	Diff line
		@@ -1699,6 +1699,9 @@ namespace Isis {
		*/
		void SpicePosition::setPolynomialSegments(int segments) {
		m_segments = segments;
		if ( !m_polynomial.isZero() ) {
		m_polynomial.refitPolynomials(segments);
		}
		}