diff --git a/include/usgscsm/UsgsAstroLsSensorModel.h b/include/usgscsm/UsgsAstroLsSensorModel.h index 8a6c3e054cc36c270b1323fe3d083018495c6fde..b2d0a9b4b62b84f27cc1fdbe443cceb4764f7a20 100644 --- a/include/usgscsm/UsgsAstroLsSensorModel.h +++ b/include/usgscsm/UsgsAstroLsSensorModel.h @@ -108,11 +108,6 @@ public: std::vector<double> m_positions; std::vector<double> m_velocities; std::vector<double> m_quaternions; - std::vector<int> m_detectorNodes; - std::vector<double> m_detectorXCoords; - std::vector<double> m_detectorYCoords; - std::vector<double> m_detectorLineCoeffs; - double m_averageDetectorSize; std::vector<double> m_currentParameterValue; std::vector<csm::param::Type> m_parameterType; csm::EcefCoord m_referencePointXyz; @@ -1043,11 +1038,10 @@ private: // Computes the imaging locus that would view a ground point at a specific // time. Computationally, this is the opposite of losToEcf. - csm::ImageCoord computeViewingPixel( + std::vector<double> computeDetectorView( const double& time, // The time to use the EO at const csm::EcefCoord& groundPoint, // The ground coordinate - const std::vector<double>& adj, // Parameter Adjustments for partials - const double& desiredPrecision // Desired precision for distortion inversion + const std::vector<double>& adj // Parameter Adjustments for partials ) const; // The linear approximation for the sensor model is used as the starting point diff --git a/include/usgscsm/Utilities.h b/include/usgscsm/Utilities.h index c5250aa935852d42ce823c1af919efee118a50b0..8110e3e28797bbb9d0a878400c09140a0dfa646e 100644 --- a/include/usgscsm/Utilities.h +++ b/include/usgscsm/Utilities.h @@ -11,25 +11,6 @@ #include <json/json.hpp> #include <Warning.h> - -double distanceToLine(double x, double y, - double a, double b, double c); - -double distanceToPlane(double x, double y, double z, - double a, double b, double c, double d); - -void line(double x1, double y1, double x2, double y2, - double &a, double &b, double &c); - -void plane(double x0, double y0, double z0, - double v1x, double v1y, double v1z, - double v2x, double v2y, double v2z, - double &a, double &b, double &c, double &d); - -std::vector<int> fitLinearApproximation(const std::vector<double> &x, - const std::vector<double> &y, - double tolerance); - // methods pulled out of los2ecf and computeViewingPixel // for now, put everything in here. diff --git a/src/UsgsAstroLsSensorModel.cpp b/src/UsgsAstroLsSensorModel.cpp index 2bf46b9c60bdc6d5bef44356e55b5b89312d361e..5fc142b7a3257a3111448c0873204e4e35718723 100644 --- a/src/UsgsAstroLsSensorModel.cpp +++ b/src/UsgsAstroLsSensorModel.cpp @@ -516,11 +516,6 @@ void UsgsAstroLsSensorModel::reset() m_positions.clear(); // 42 m_velocities.clear(); // 43 m_quaternions.clear(); // 44 - m_detectorNodes.clear(); - m_detectorXCoords.clear(); - m_detectorYCoords.clear(); - m_detectorLineCoeffs.clear(); - m_averageDetectorSize = 0.0; m_currentParameterValue.assign(NUM_PARAMETERS,0.0); m_parameterType.assign(NUM_PARAMETERS,csm::param::REAL); @@ -614,37 +609,6 @@ void UsgsAstroLsSensorModel::updateState() MESSAGE_LOG(m_logger, "updateState: half time duration set to {}", m_halfTime) - // compute linearized detector coordinates - m_detectorXCoords.resize(m_nSamples); - m_detectorYCoords.resize(m_nSamples); - - for (int i = 0; i < m_nSamples; i++) { - computeDistortedFocalPlaneCoordinates( - 0.5, i, - m_detectorSampleOrigin, m_detectorLineOrigin, - m_detectorSampleSumming, 1.0, - m_startingSample, 0.0, - m_iTransS, m_iTransL, - m_detectorXCoords[i], m_detectorYCoords[i] - ); - } - m_averageDetectorSize = 0; - for (int i = 1; i < m_nSamples; i++) { - double xDist = m_detectorXCoords[i] - m_detectorXCoords[i-1]; - double yDist = m_detectorYCoords[i] - m_detectorYCoords[i-1]; - m_averageDetectorSize += sqrt(xDist*xDist + yDist*yDist); - } - m_averageDetectorSize /= (m_nSamples-1); - m_detectorNodes = fitLinearApproximation(m_detectorXCoords, m_detectorYCoords, - m_averageDetectorSize); - - m_detectorLineCoeffs.resize(3 * (m_detectorNodes.size() - 1)); - for (size_t i = 0; i + 1 < m_detectorNodes.size(); i++) { - line(m_detectorXCoords[m_detectorNodes[i]], m_detectorYCoords[m_detectorNodes[i]], - m_detectorXCoords[m_detectorNodes[i+1]], m_detectorYCoords[m_detectorNodes[i+1]], - m_detectorLineCoeffs[3*i], m_detectorLineCoeffs[3*i+1], m_detectorLineCoeffs[3*i+2]); - } - // Parameter covariance, hardcoded accuracy values int num_params = NUM_PARAMETERS; int num_paramsSquare = num_params * num_params; @@ -687,129 +651,76 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage( // UsgsAstroLsSensorModel::groundToImage (internal version) //*************************************************************************** csm::ImageCoord UsgsAstroLsSensorModel::groundToImage( - const csm::EcefCoord& ground_pt, + const csm::EcefCoord& groundPt, const std::vector<double>& adj, - double desired_precision, - double* achieved_precision, + double desiredPrecision, + double* achievedPrecision, csm::WarningList* warnings) const { - MESSAGE_LOG(m_logger, "Computing groundToImage (with adjustments) for {}, {}, {}, with desired precision {}", - ground_pt.x, ground_pt.y, ground_pt.z, desired_precision); // Search for the line, sample coordinate that viewed a given ground point. - // This method uses an iterative secant method to search for the image - // line. Since this is looking for a zero we subtract 0.5 each time the offsets - // are calculated. This allows it to converge on the center of the pixel where - // there is only one correct answer instead of the top or bottom of the pixel - // where there are two correct answers. - - // Convert the ground precision to pixel precision so we can - // check for convergence without re-intersecting - csm::ImageCoord approxPoint; - computeLinearApproximation(ground_pt, approxPoint); - csm::ImageCoord approxNextPoint = approxPoint; - if (approxNextPoint.line + 1 < m_nLines) { - ++approxNextPoint.line; - } - else { - --approxNextPoint.line; - } - double height, aPrec; - computeElevation(ground_pt.x, ground_pt.y, ground_pt.z, height, aPrec, desired_precision); - csm::EcefCoord approxIntersect = imageToGround(approxPoint, height); - csm::EcefCoord approxNextIntersect = imageToGround(approxNextPoint, height); - double lineDX = approxNextIntersect.x - approxIntersect.x; - double lineDY = approxNextIntersect.y - approxIntersect.y; - double lineDZ = approxNextIntersect.z - approxIntersect.z; - double approxLineRes = sqrt(lineDX * lineDX + lineDY * lineDY + lineDZ * lineDZ); - // Increase the precision by a small amount to ensure the desired precision is met - double pixelPrec = desired_precision / approxLineRes * 0.9; - - // Start secant method search on the image lines - double sampCtr = m_nSamples / 2.0; - double firstTime = getImageTime(csm::ImageCoord(0.0, sampCtr)); - double lastTime = getImageTime(csm::ImageCoord(m_nLines, sampCtr)); - // See comment above for why (- 0.5) - double firstOffset = computeViewingPixel(firstTime, ground_pt, adj, pixelPrec/2).line - 0.5; - double lastOffset = computeViewingPixel(lastTime, ground_pt, adj, pixelPrec/2).line - 0.5; - MESSAGE_LOG(m_logger, "groundToImage: Initial firstOffset {}, lastOffset {}", firstOffset, lastOffset) - - double closestLine; - csm::ImageCoord detectorCoord; - - // Start secant method search - for (int it = 0; it < 30; it++) { - double nextTime = ((firstTime * lastOffset) - (lastTime * firstOffset)) - / (lastOffset - firstOffset); - // Because time across the image is not continuous, find the exposure closest - // to the computed nextTime and use that. - - // I.E. if the computed nextTime is 0.3, and the middle exposure times for - // lines are 0.07, 0.17, 0.27, 0.37, and 0.47; then use 0.27 because it is - // the closest middle exposure time. - auto referenceTimeIt = std::upper_bound(m_intTimeStartTimes.begin(), - m_intTimeStartTimes.end(), - nextTime); - if (referenceTimeIt != m_intTimeStartTimes.begin()) { - --referenceTimeIt; - } + // This method first uses a linear approximation to get an initial point. + // Then the detector offset for the line is continuously computed and + // applied to the line until we achieve the desired precision. - size_t referenceIndex = std::distance(m_intTimeStartTimes.begin(), referenceTimeIt); - MESSAGE_LOG(m_logger, "groundToImage: Find reference time, line number {}, start time {}, duration {} ", - m_intTimeLines[referenceIndex], m_intTimeStartTimes[referenceIndex], m_intTimes[referenceIndex]) - double computedLine = (nextTime - m_intTimeStartTimes[referenceIndex]) / m_intTimes[referenceIndex] - + m_intTimeLines[referenceIndex]; - // Remove -0.5 once ISIS linescanrate is fixed - closestLine = floor(computedLine - 0.5); - nextTime = getImageTime(csm::ImageCoord(closestLine, sampCtr)); - - detectorCoord = computeViewingPixel(nextTime, ground_pt, adj, pixelPrec/2); - double nextOffset = detectorCoord.line - 0.5; - - // remove the farthest away node - if (fabs(firstTime - nextTime) > fabs(lastTime - nextTime)) { - firstTime = nextTime; - firstOffset = nextOffset; - } - else { - lastTime = nextTime; - lastOffset = nextOffset; - } - MESSAGE_LOG(m_logger, "groundToImage: Loop firstOffset {}, firstTime {}, lastOffset {}, lastTime {}, nextOffset {}, nextTime {}", - firstOffset, firstTime, lastOffset, lastTime, nextOffset, nextTime) + csm::ImageCoord approxPt; + computeLinearApproximation(groundPt, approxPt); - if (fabs(lastOffset - firstOffset) < pixelPrec) { - MESSAGE_LOG(m_logger, "groundToImage: Final firstOffset {}, lastOffset {}, pixelPre {}", firstOffset, lastOffset, pixelPrec) - break; - } + std::vector<double> detectorView; + double detectorLine = m_nLines; + double count = 0; + double timei; + + while(abs(detectorLine) > desiredPrecision && ++count < 15) { + timei = getImageTime(approxPt); + detectorView = computeDetectorView(timei, groundPt, adj); + + // Convert to detector line + detectorLine = m_iTransL[0] + + m_iTransL[1] * detectorView[0] + + m_iTransL[2] * detectorView[1]; + + // Convert to image line + approxPt.line += detectorLine; } - // The computed pixel is the detector pixel, so we need to convert that to image lines - csm::ImageCoord calculatedPixel(detectorCoord.line + closestLine, detectorCoord.samp); + timei = getImageTime(approxPt); + detectorView = computeDetectorView(timei, groundPt, adj); + + // Invert distortion + double distortedFocalX, distortedFocalY; + applyDistortion(detectorView[0], detectorView[1], distortedFocalX, distortedFocalY, + m_opticalDistCoeffs, m_distortionType, desiredPrecision); + + // Convert to detector line and sample + detectorLine = m_iTransL[0] + + m_iTransL[1] * distortedFocalX + + m_iTransL[2] * distortedFocalY; + double detectorSample = m_iTransS[0] + + m_iTransS[1] * distortedFocalX + + m_iTransS[2] * distortedFocalY; - // Reintersect to ensure the image point actually views the ground point. - csm::EcefCoord calculatedPoint = imageToGround(calculatedPixel, height); - double dx = ground_pt.x - calculatedPoint.x; - double dy = ground_pt.y - calculatedPoint.y; - double dz = ground_pt.z - calculatedPoint.z; - double len = dx * dx + dy * dy + dz * dz; + // Convert to image sample line + approxPt.line += detectorLine; + approxPt.samp = (detectorSample + m_detectorSampleOrigin - m_startingSample) + / m_detectorSampleSumming; - if (achieved_precision) { - *achieved_precision = sqrt(len); + if (achievedPrecision) { + *achievedPrecision = detectorLine; } - double preSquare = desired_precision * desired_precision; - if (warnings && (desired_precision > 0.0) && (preSquare < len)) { - MESSAGE_LOG(m_logger, "groundToImage: Desired precision not achieved {}", preSquare) - std::stringstream msg; - msg << "Desired precision not achieved. "; - msg << len << " " << preSquare << "\n"; - warnings->push_back(csm::Warning(csm::Warning::PRECISION_NOT_MET, - msg.str().c_str(), - "UsgsAstroLsSensorModel::groundToImage()")); + MESSAGE_LOG(m_logger, "groundToImage: image line sample {} {}", + approxPt.line, approxPt.samp) + + if (warnings && (desiredPrecision > 0.0) && (abs(detectorLine) > desiredPrecision)) + { + warnings->push_back( + csm::Warning( + csm::Warning::PRECISION_NOT_MET, + "Desired precision not achieved.", + "UsgsAstroLsSensorModel::groundToImage()")); } - MESSAGE_LOG(m_logger, "groundToImage: Final pixel line {}, sample {}", calculatedPixel.line, calculatedPixel.samp) - return calculatedPixel; + return approxPt; } //*************************************************************************** @@ -1356,7 +1267,7 @@ double UsgsAstroLsSensorModel::getImageTime( const csm::ImageCoord& image_pt) const { // Remove 0.5 after ISIS dependency in the linescanrate is gone - double lineFull = floor(image_pt.line) + 0.5; + double lineFull = image_pt.line; auto referenceLineIt = std::upper_bound(m_intTimeLines.begin(), m_intTimeLines.end(), @@ -1969,12 +1880,11 @@ void UsgsAstroLsSensorModel::losToEcf( // USGS image convention: UL pixel center == (1.0, 1.0) double sampleCSMFull = sample; double sampleUSGSFull = sampleCSMFull; - double fractionalLine = line - floor(line); // Compute distorted image coordinates in mm (sample, line on image (pixels) -> focal plane double distortedFocalPlaneX, distortedFocalPlaneY; computeDistortedFocalPlaneCoordinates( - fractionalLine, sampleUSGSFull, + m_detectorLineOrigin, sampleUSGSFull, m_detectorSampleOrigin, m_detectorLineOrigin, m_detectorSampleSumming, 1.0, m_startingSample, 0.0, @@ -2369,6 +2279,7 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel( double sensVelNom[3]; lagrangeInterp(m_numPositions/3, &m_velocities[0], m_t0Ephem, m_dtEphem, time, 3, nOrder, sensVelNom); + MESSAGE_LOG(m_logger, "getAdjSensorPosVel: using {} order Lagrange", nOrder) @@ -2437,18 +2348,18 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel( //*************************************************************************** -// UsgsAstroLineScannerSensorModel::computeViewingPixel +// UsgsAstroLineScannerSensorModel::computeDetectorView //*************************************************************************** -csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel( +std::vector<double> UsgsAstroLsSensorModel::computeDetectorView( const double& time, const csm::EcefCoord& groundPoint, - const std::vector<double>& adj, - const double& desiredPrecision) const + const std::vector<double>& adj) const { - MESSAGE_LOG(m_logger, "Computing computeViewingPixel (with adjusments)" + MESSAGE_LOG(m_logger, "Computing computeDetectorView (with adjusments)" "for ground point {} {} {} at time {} ", groundPoint.x, groundPoint.y, groundPoint.z, time) + // Helper function to compute the CCD pixel that views a ground point based // on the exterior orientation at a given time. @@ -2460,7 +2371,7 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel( double bodyLookX = groundPoint.x - xc; double bodyLookY = groundPoint.y - yc; double bodyLookZ = groundPoint.z - zc; - MESSAGE_LOG(m_logger, "computeViewingPixel: look vector {} {} {}", + MESSAGE_LOG(m_logger, "computeDetectorView: look vector {} {} {}", bodyLookX, bodyLookY, bodyLookZ) // Rotate the look vector into the camera reference frame @@ -2479,7 +2390,7 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel( double cameraLookZ = bodyToCamera[2] * bodyLookX + bodyToCamera[5] * bodyLookY + bodyToCamera[8] * bodyLookZ; - MESSAGE_LOG(m_logger, "computeViewingPixel: look vector (camrea ref frame)" + MESSAGE_LOG(m_logger, "computeDetectorView: look vector (camrea ref frame)" "{} {} {}", cameraLookX, cameraLookY, cameraLookZ) @@ -2497,7 +2408,7 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel( double adjustedLookZ = attCorr[2] * cameraLookX + attCorr[5] * cameraLookY + attCorr[8] * cameraLookZ; - MESSAGE_LOG(m_logger, "computeViewingPixel: adjusted look vector" + MESSAGE_LOG(m_logger, "computeDetectorView: adjusted look vector" "{} {} {}", adjustedLookX, adjustedLookY, adjustedLookZ) @@ -2505,31 +2416,12 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel( double lookScale = (m_focalLength + getValue(15, adj)) / adjustedLookZ; double focalX = adjustedLookX * lookScale; double focalY = adjustedLookY * lookScale; - double distortedFocalX, distortedFocalY; - MESSAGE_LOG(m_logger, "computeViewingPixel: focal plane coordinates" - "x:{} y:{} scale:{}", - focalX, focalY, lookScale) - // Invert distortion - applyDistortion(focalX, focalY, distortedFocalX, distortedFocalY, - m_opticalDistCoeffs, m_distortionType, desiredPrecision); - - // Convert to detector line and sample - double detectorLine = m_iTransL[0] - + m_iTransL[1] * distortedFocalX - + m_iTransL[2] * distortedFocalY; - double detectorSample = m_iTransS[0] - + m_iTransS[1] * distortedFocalX - + m_iTransS[2] * distortedFocalY; - - // Convert to image sample line - double line = detectorLine + m_detectorLineOrigin; - double sample = (detectorSample + m_detectorSampleOrigin - m_startingSample) - / m_detectorSampleSumming; - MESSAGE_LOG(m_logger, "computeViewingPixel: image line sample {} {}", - line, sample) + MESSAGE_LOG(m_logger, "computeDetectorView: focal plane coordinates" + "x:{} y:{} scale:{}", + focalX, focalY, lookScale) - return csm::ImageCoord(line, sample); + return std::vector<double> {focalX, focalY}; } diff --git a/src/Utilities.cpp b/src/Utilities.cpp index 1f0b5907e9f9d8bc17dd7be1ffd4310041f2bfb4..39ce6ed2ac2c2b56f994e37b23b06201a9951c8f 100644 --- a/src/Utilities.cpp +++ b/src/Utilities.cpp @@ -7,102 +7,6 @@ using json = nlohmann::json; -// Calculate the signed distance from a 2D point to a line given in standard form -// -// --NOTE-- This function assumes that the coefficients of the line have been reduced -// so that sqrt(a^2 + b^2) = 1 -double distanceToLine(double x, double y, - double a, double b, double c) { - return a*x + b*y + c; -} - -// Calculate the distance from a 3D point to a plane given in standard form -// -// --NOTE-- This function assumes that the coefficients of the plane have been reduced -// so that sqrt(a^2 + b^2 + c^2) = 1 -double distanceToPlane(double x, double y, double z, - double a, double b, double c, double d) { - return std::abs(a*x + b*y + c*z + d); -} - -// Calculate the standard form coefficients of a line that passes through two points -// -// --NOTE-- The coefficients have been reduced such that sqrt(a^2 + b^2) = 1 -void line(double x1, double y1, double x2, double y2, - double &a, double &b, double &c) { - a = y1 - y2; - b = x2 - x1; - c = x1*y2 - x2*y1; - double scale = sqrt(a*a + b*b); - a /= scale; - b /= scale; - c /= scale; -} - -// Calculate the standard form coefficients of a plane that passes through a point -// and contains two vectors. -// -// --NOTE-- The coefficients have been reduced such that sqrt(a^2 + b^2 + c^2) = 1 -void plane(double x0, double y0, double z0, - double v1x, double v1y, double v1z, - double v2x, double v2y, double v2z, - double &a, double &b, double &c, double &d) { - // Compute normal vector and scale - a = v1y*v2z - v1z*v2y; - b = v1z*v2x - v1x*v2z; - c = v1x*v2y - v1y*v2x; - double scale = sqrt(a*a + b*b + c*c); - a /= scale; - b /= scale; - c /= scale; - // Shift to point - d = - (a*x0 + b*y0 + c*z0); -} - - - -// Fit a piecewise-linear approximations to 2D data within a tolerance -// -// Returns a vector of node indices -std::vector<int> fitLinearApproximation(const std::vector<double> &x, - const std::vector<double> &y, - double tolerance) { - std::vector<int> nodes; - nodes.push_back(0); - - std::stack<std::pair<int, int>> workStack; - workStack.push(std::make_pair(0, x.size() - 1)); - while (!workStack.empty()) { - std::pair<int, int> range = workStack.top(); - workStack.pop(); - double a, b, c; - line(x[range.first], y[range.first], x[range.second], y[range.second], a, b, c); - double maxError = 0; - int maxIndex = (range.second + range.first) / 2; - for (int i = range.first + 1; i < range.second; i++) { - double error = std::abs(distanceToLine(x[i], y[i], a, b, c)); - if (error > maxError) { - maxError = error; - maxIndex = i; - } - } - // If the max error is greater than the tolerance, split at the largest error - // Do not split if the range only contains two nodes - if (maxError > tolerance && range.second - range.first > 1) { - // Append the second range and then the first range - // so that nodes are added in the same order they came in - workStack.push(std::make_pair(maxIndex, range.second)); - workStack.push(std::make_pair(range.first, maxIndex)); - } - else { - // segment is good so append last point to nodes - nodes.push_back(range.second); - } - } - - return nodes; -} - // Calculates a rotation matrix from Euler angles // in - euler[3] // out - rotationMatrix[9] diff --git a/tests/Fixtures.h b/tests/Fixtures.h index 48164b335a16f7b66e4254d59eb2565d11a738fe..5b89d9d7d33aebc3b4eb0b161bcf97b5d5838263 100644 --- a/tests/Fixtures.h +++ b/tests/Fixtures.h @@ -240,33 +240,6 @@ class ConstVelocityLineScanSensorModel : public ::testing::Test { } }; -class ConstAngularVelocityLineScanSensorModel : public ::testing::Test { - protected: - csm::Isd isd; - UsgsAstroLsSensorModel *sensorModel; - - void SetUp() override { - sensorModel = NULL; - - isd.setFilename("data/constAngularVelocityLineScan.img"); - UsgsAstroPlugin cameraPlugin; - - csm::Model *model = cameraPlugin.constructModelFromISD( - isd, - "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL"); - sensorModel = dynamic_cast<UsgsAstroLsSensorModel *>(model); - - ASSERT_NE(sensorModel, nullptr); - } - - void TearDown() override { - if (sensorModel) { - delete sensorModel; - sensorModel = NULL; - } - } -}; - class OrbitalLineScanSensorModel : public ::testing::Test { protected: csm::Isd isd; diff --git a/tests/LineScanCameraTests.cpp b/tests/LineScanCameraTests.cpp index f4a07e68b9322319683477ef7d8b1da915a37a67..5daec7f73159d2bbd69853178cf2dc1ab46e0387 100644 --- a/tests/LineScanCameraTests.cpp +++ b/tests/LineScanCameraTests.cpp @@ -34,12 +34,14 @@ TEST_F(ConstVelocityLineScanSensorModel, Center) { } TEST_F(ConstVelocityLineScanSensorModel, Inversion) { + double achievedPrecision; csm::ImageCoord imagePt(8.5, 8); csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0); - csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt); + csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt, 0.001, &achievedPrecision); - EXPECT_DOUBLE_EQ(imagePt.line, imageReprojPt.line); - EXPECT_DOUBLE_EQ(imagePt.samp, imageReprojPt.samp); + EXPECT_LT(achievedPrecision, 0.001); + EXPECT_NEAR(imagePt.line, imageReprojPt.line, 1e-3); + EXPECT_NEAR(imagePt.samp, imageReprojPt.samp, 1e-3); } TEST_F(ConstVelocityLineScanSensorModel, OffBody) { @@ -73,39 +75,12 @@ TEST_F(ConstVelocityLineScanSensorModel, RemoteImageLocus) { EXPECT_NEAR(locus.point.z, 0.0, 1e-12); } -// Pan tests - -TEST_F(ConstAngularVelocityLineScanSensorModel, Center) { - csm::ImageCoord imagePt(8.5, 8.0); - csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0); - EXPECT_DOUBLE_EQ(groundPt.x, 10.0); - EXPECT_DOUBLE_EQ(groundPt.y, 0.0); - EXPECT_DOUBLE_EQ(groundPt.z, 0.0); -} - -TEST_F(ConstAngularVelocityLineScanSensorModel, Inversion) { - csm::ImageCoord imagePt(8.5, 8); - csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0); - csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt); - - EXPECT_DOUBLE_EQ(imagePt.line, imageReprojPt.line); - EXPECT_DOUBLE_EQ(imagePt.samp, imageReprojPt.samp); -} - -TEST_F(ConstAngularVelocityLineScanSensorModel, OffBody) { - csm::ImageCoord imagePt(4.5, 4.0); - csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0); - EXPECT_NEAR(groundPt.x, 4.15414478, 1e-8); - EXPECT_NEAR(groundPt.y, -7.98311067, 1e-8); - EXPECT_NEAR(groundPt.z, 63.82129588, 1e-8); -} - TEST_F(ConstVelocityLineScanSensorModel, calculateAttitudeCorrection) { std::vector<double> adj; double attCorr[9]; adj.resize(15, 0); // Pi/2 with simply compensating for member variable m_flyingHeight in UsgsAstroLsSensorModel - adj[7] = (M_PI / 2) * 990.0496255790623081338708; + adj[7] = (M_PI / 2) * sensorModel->m_flyingHeight; sensorModel->calculateAttitudeCorrection(999.5, adj, attCorr); // EXPECT_NEARs are used here instead of EXPECT_DOUBLE_EQs because index 0 and 8 of the matrix @@ -138,9 +113,9 @@ TEST_F(OrbitalLineScanSensorModel, getIlluminationDirectionStationary) { // Calculate expected sun direction // These are the ground point coordinates minus constant sun positions. - double expected_x = 999899.680000017; - double expected_y = -100; - double expected_z = -899.99991466668735; + double expected_x = groundPt.x - sensorModel->m_sunPosition[0]; + double expected_y = groundPt.y - sensorModel->m_sunPosition[1]; + double expected_z = groundPt.z - sensorModel->m_sunPosition[2]; //normalize double scale = sqrt((expected_x * expected_x) + (expected_y * expected_y) + (expected_z * expected_z)); @@ -194,19 +169,21 @@ TEST_F(OrbitalLineScanSensorModel, getSunPositionStationary){ TEST_F(OrbitalLineScanSensorModel, Center) { csm::ImageCoord imagePt(8.5, 8.0); csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0); - EXPECT_DOUBLE_EQ(groundPt.x, 999999.680000017); + EXPECT_DOUBLE_EQ(groundPt.x, 999999.70975015126); EXPECT_DOUBLE_EQ(groundPt.y, 0.0); - EXPECT_DOUBLE_EQ(groundPt.z, -799.99991466668735); + EXPECT_DOUBLE_EQ(groundPt.z, -761.90525203960692); } TEST_F(OrbitalLineScanSensorModel, Inversion) { + double achievedPrecision; for (double line = 0.5; line < 16; line++) { csm::ImageCoord imagePt(line, 8); csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0); - csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt); + csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt, 0.001, &achievedPrecision); // groundToImage has a default precision of 0.001m and each pixel is 100m // so we should be within 0.1 pixels + EXPECT_LT(achievedPrecision, 0.001); EXPECT_NEAR(imagePt.line, imageReprojPt.line, 0.1); EXPECT_NEAR(imagePt.samp, imageReprojPt.samp, 0.1); } @@ -238,7 +215,7 @@ TEST_F(OrbitalLineScanSensorModel, InversionHeight) { TEST_F(OrbitalLineScanSensorModel, InversionReallyHigh) { for (double line = 0.5; line < 16; line++) { csm::ImageCoord imagePt(line, 8); - csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 49000.0); + csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 4900.0); csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt); // groundToImage has a default precision of 0.001m and each pixel is 2m diff --git a/tests/UtilitiesTests.cpp b/tests/UtilitiesTests.cpp index d6fdca569ba8617f5b011db8ff6c0c5f2a28d455..f837808e0dcf28c4091f5d212b2a63031b6629c0 100644 --- a/tests/UtilitiesTests.cpp +++ b/tests/UtilitiesTests.cpp @@ -11,82 +11,6 @@ using json = nlohmann::json; -TEST(UtilitiesTests, distanceToLine) { - // Line passing through (0, 1) and (2,2) - double a = -1.0 / sqrt(5); - double b = 2.0 / sqrt(5); - double c = -2.0 / sqrt(5); - // Point on line - EXPECT_NEAR(distanceToLine(4.0, 3.0, a, b, c), 0, 1e-12); - // Point above line - EXPECT_DOUBLE_EQ(distanceToLine(1.0, 4.0, a, b, c), sqrt(5)); - // Point below line - EXPECT_DOUBLE_EQ(distanceToLine(4.5, 2.0, a, b, c), -sqrt(5) / 2.0); -} - -TEST(UtilitiesTests, line) { - double x1 = 0.0; - double y1 = 1.0; - double x2 = 2.0; - double y2 = 2.0; - double a, b, c; - line(x1, y1, x2, y2, a, b, c); - EXPECT_DOUBLE_EQ(a, -1.0 / sqrt(5)); - EXPECT_DOUBLE_EQ(b, 2.0 / sqrt(5)); - EXPECT_DOUBLE_EQ(c, -2.0 / sqrt(5)); -} - -TEST(UtilitiesTests, distanceToPlane) { - // Plane passing through (0, 0, 1), (0, 2, 0) and (3, 0, 0) - double a = 2.0 / 7.0; - double b = 3.0 / 7.0; - double c = 6.0 / 7.0; - double d = -12.0 / 7.0; - // Points on plane - EXPECT_NEAR(distanceToPlane(0.0, 0.0, 2.0, a, b, c, d), 0.0, 1e-12); - EXPECT_NEAR(distanceToPlane(0.0, 4.0, 0.0, a, b, c, d), 0.0, 1e-12); - EXPECT_NEAR(distanceToPlane(6.0, 0.0, 0.0, a, b, c, d), 0.0, 1e-12); - // Points off plane - EXPECT_DOUBLE_EQ(distanceToPlane(0.0, 0.0, 0.0, a, b, c, d), 12.0 / 7.0); - EXPECT_DOUBLE_EQ(distanceToPlane(0.0, 0.0, 1.0, a, b, c, d), 6.0 / 7.0); - EXPECT_DOUBLE_EQ(distanceToPlane(0.0, 2.0, 0.0, a, b, c, d), 6.0 / 7.0); - EXPECT_DOUBLE_EQ(distanceToPlane(3.0, 0.0, 0.0, a, b, c, d), 6.0 / 7.0); -} - -TEST(UtilitiesTests, plane) { - double x0 = 0.0; - double y0 = 0.0; - double z0 = 2.0; - double v1x = 0.0; - double v1y = 4.0; - double v1z = -2.0; - double v2x = 6.0; - double v2y = 0.0; - double v2z = -2.0; - double a, b, c, d; - plane(x0, y0, z0, v1x, v1y, v1z, v2x, v2y, v2z, a, b, c, d); - EXPECT_DOUBLE_EQ(a, -2.0 / 7.0); - EXPECT_DOUBLE_EQ(b, -3.0 / 7.0); - EXPECT_DOUBLE_EQ(c, -6.0 / 7.0); - EXPECT_DOUBLE_EQ(d, 12.0 / 7.0); -} - -TEST(UtilitiesTests, fitLinearApproximation) { - std::vector<double> x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; - std::vector<double> y = {1, 0, 5, 7, -2, 1, 2, 2, 1, 0, 1}; - - std::vector<int> nodes = fitLinearApproximation(x, y, 1.0); - ASSERT_EQ(nodes.size(), 7); - EXPECT_EQ(nodes[0], 0); - EXPECT_EQ(nodes[1], 1); - EXPECT_EQ(nodes[2], 3); - EXPECT_EQ(nodes[3], 4); - EXPECT_EQ(nodes[4], 6); - EXPECT_EQ(nodes[5], 9); - EXPECT_EQ(nodes[6], 10); -} - - TEST(UtilitiesTests, calculateRotationMatrixFromEuler) { double euler[3], rotationMatrix[9]; euler[0] = 0; diff --git a/tests/data/constAngularVelocityLineScan.json b/tests/data/constAngularVelocityLineScan.json deleted file mode 100644 index 69cb6b7716407f9933f2b1192545c132a1fd5aa1..0000000000000000000000000000000000000000 --- a/tests/data/constAngularVelocityLineScan.json +++ /dev/null @@ -1,92 +0,0 @@ -{ - "name_model" : "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL", - "name_platform" : "TEST_PLATFORM", - "name_sensor" : "TEST_SENSOR", - "center_ephemeris_time": 1000.0, - "starting_ephemeris_time": 999.2, - "line_scan_rate": [ - [0.5, -0.85, 0.1] - ], - "detector_sample_summing": 1, - "detector_line_summing": 1, - "t0_ephemeris": -0.8, - "dt_ephemeris": 0.2, - "t0_quaternion": -0.8, - "dt_quaternion": 0.2, - "focal2pixel_lines": [0.0, 10.0, 0.0], - "focal2pixel_samples": [0.0, 0.0, 10.0], - "focal_length_model": { - "focal_length": 50 - }, - "image_lines": 16, - "image_samples": 16, - "detector_center" : { - "line" : 0.5, - "sample" : 8.0 - }, - "interpolation_method": "lagrange", - "optical_distortion": { - "radial": { - "coefficients": [0.0, 0.0, 0.0] - } - }, - "radii": { - "semimajor": 10, - "semiminor": 10, - "unit":"m" - }, - "reference_height": { - "maxheight": 1, - "minheight": -1, - "unit": "m" - }, - "sensor_position": { - "unit": "m", - "positions": [ - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0], - [1000.0, 0.0, 0.0] - ], - "velocities": [ - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0], - [0.0, 10.0, 0.0] - ] - }, - "sensor_orientation": { - "quaternions": [ - [0.0, -0.660435174378928, 0, 0.750883067090392], - [0.0, -0.672364256957188, 0, 0.740220444169444], - [0.0, -0.68412121795764, 0, 0.729368328857344], - [0.0, -0.695703047662478, 0, 0.718329499236346], - [0.0, -0.707106781186547, 0.0, 0.707106781186547], - [0.0, -0.718329499236346, 0, 0.695703047662478], - [0.0, -0.729368328857344, 0, 0.68412121795764], - [0.0, -0.740220444169444, 0, 0.672364256957188], - [0.0, -0.750883067090392, 0, 0.660435174378928] - ] - }, - "starting_detector_line": 0, - "starting_detector_sample": 0, - "sun_position": { - "positions": [ - [100.0, 100.0, 0.0] - ], - "velocities": [ - [0.0, 0.0, 0.0] - ], - "unit": "m" - } -} diff --git a/tests/data/orbitalLineScan.json b/tests/data/orbitalLineScan.json index f1e00928eae6c7dba80d78610d29e7a4a908e6de..f101451cedfaee24b47248435c9518f387f9ddb7 100644 --- a/tests/data/orbitalLineScan.json +++ b/tests/data/orbitalLineScan.json @@ -43,60 +43,60 @@ "sensor_position": { "unit": "km", "positions": [ - [1050.0, 0.0, 0.0], - [1049.99999475, 0.0, -0.104999999825], - [1049.9999790000002, 0.0, -0.2099999986], - [1049.9999527500004, 0.0, -0.3149999952750001], - [1049.9999160000013, 0.0, -0.4199999888000002], - [1049.9998687500029, 0.0, -0.5249999781250003], - [1049.9998110000056, 0.0, -0.6299999622000008], - [1049.9997427500105, 0.0, -0.7349999399750016], - [1049.999664000018, 0.0, -0.839999910400003], - [1049.9995747500286, 0.0, -0.9449998724250054], - [1049.999475000044, 0.0, -1.049999825000009], - [1049.999364750064, 0.0, -1.1549997670750143], - [1049.9992440000908, 0.0, -1.259999697600022], - [1049.9991127501248, 0.0, -1.3649996155250328], - [1049.9989710001682, 0.0, -1.4699995198000473], - [1049.9988187502213, 0.0, -1.5749994093750666] + [1050.0, 0.0, -0.0], + [1049.99999524, 0.0, -0.1], + [1049.99998095, 0.0, -0.2], + [1049.99995714, 0.0, -0.3], + [1049.99992381, 0.0, -0.39999999], + [1049.99988095, 0.0, -0.49999998], + [1049.99982857, 0.0, -0.59999997], + [1049.99976667, 0.0, -0.69999995], + [1049.99969524, 0.0, -0.79999992], + [1049.99961429, 0.0, -0.89999989], + [1049.99952381, 0.0, -0.99999985], + [1049.99942381, 0.0, -1.0999998], + [1049.99931429, 0.0, -1.19999974], + [1049.99919524, 0.0, -1.29999967], + [1049.99906667, 0.0, -1.39999959], + [1049.99892857, 0.0, -1.49999949] ], "velocities": [ - [ 0.0, 0.0, -1.0], - [-0.0000999999998333, 0.0, -0.999999995], - [-0.0001999999986667, 0.0, -0.9999999800000001], - [-0.0002999999955, 0.0, -0.9999999550000004], - [-0.0003999999893333, 0.0, -0.9999999200000013], - [-0.0004999999791667, 0.0, -0.9999998750000026], - [-0.000599999964, 0.0, -0.9999998200000054], - [-0.0006999999428333, 0.0, -0.9999997550000099], - [-0.0007999999146667, 0.0, -0.999999680000017], - [-0.0008999998785, 0.0, -0.9999995950000273], - [-0.0009999998333333, 0.0, -0.9999995000000418], - [-0.0010999997781667, 0.0, -0.999999395000061], - [-0.001199999712, 0.0, -0.9999992800000864], - [-0.0012999996338334, 0.0, -0.9999991550001189], - [-0.0013999995426667, 0.0, -0.9999990200001602], - [-0.0014999994375001, 0.0, -0.9999988750002108] + [-0.0, 0.0, -1050.0], + [-0.1, 0.0, -1049.99999524], + [-0.2, 0.0, -1049.99998095], + [-0.3, 0.0, -1049.99995714], + [-0.39999999, 0.0 ,-1049.99992381], + [-0.49999998, 0.0 ,-1049.99988095], + [-0.59999997, 0.0 ,-1049.99982857], + [-0.69999995, 0.0 ,-1049.99976667], + [-0.79999992, 0.0 ,-1049.99969524], + [-0.89999989, 0.0 ,-1049.99961429], + [-0.99999985, 0.0 ,-1049.99952381], + [-1.0999998, 0.0 ,-1049.99942381], + [-1.19999974, 0.0 ,-1049.99931429], + [-1.29999967, 0.0 ,-1049.99919524], + [-1.39999959, 0.0 ,-1049.99906667], + [-1.49999949, 0.0, -1049.99892857] ] }, "sensor_orientation": { "quaternions": [ - [0.0, 0.707106781186547, 0, 0.707106781186547], - [0.0, 0.70707142496362, 0, -0.707142135641709], - [0.0, 0.707036066973013, 0, -0.707177488329014], - [0.0, 0.707000707214816, 0, -0.707212839248377], - [0.0, 0.706965345689117, 0, -0.707248188399706], - [0.0, 0.706929982396005, 0, -0.707283535782916], - [0.0, 0.706894617335569, 0, -0.707318881397917], - [0.0, 0.706859250507895, 0, -0.70735422524462], - [0.0, 0.706823881913074, 0, -0.707389567322938], - [0.0, 0.706788511551192, 0, -0.707424907632782], - [0.0, 0.70675313942234, 0, -0.707460246174064], - [0.0, 0.706717765526604, 0, -0.707495582946695], - [0.0, 0.706682389864075, 0, -0.707530917950587], - [0.0, 0.706647012434839, 0, -0.707566251185652], - [0.0, 0.706611633238985, 0, -0.707601582651801], - [0.0, 0.706576252276603, 0, -0.707636912348946] + [0.0, -0.70710678, 0.0, 0.70710678], + [0.0, -0.70707311, 0.0, 0.70714045], + [0.0, -0.70703943, 0.0, 0.70717412], + [0.0, -0.70700576, 0.0, 0.70720779], + [0.0, -0.70697208, 0.0, 0.70724146], + [0.0, -0.7069384, 0.0, 0.70727512], + [0.0, -0.70690472, 0.0, 0.70730878], + [0.0, -0.70687104, 0.0, 0.70734244], + [0.0, -0.70683736, 0.0, 0.7073761], + [0.0, -0.70680367, 0.0, 0.70740976], + [0.0, -0.70676998, 0.0, 0.70744342], + [0.0, -0.70673629, 0.0, 0.70747707], + [0.0, -0.7067026, 0.0, 0.70751073], + [0.0, -0.70666891, 0.0, 0.70754438], + [0.0, -0.70663522, 0.0, 0.70757803], + [0.0, -0.70660152, 0.0, 0.70761168] ] }, "starting_detector_line": 0,