diff --git a/include/genericls/UsgsAstroLsSensorModel.h b/include/genericls/UsgsAstroLsSensorModel.h index 4c860981fc00f402f30d22029eb9ce9539762038..ec262a7cacba5369ccc2ffc5a99b7093df2f5851 100644 --- a/include/genericls/UsgsAstroLsSensorModel.h +++ b/include/genericls/UsgsAstroLsSensorModel.h @@ -930,6 +930,14 @@ private: double& vy, double& vz) const; + // Computes the imaging locus that would view a ground point at a specific + // time. Computationally, this is the opposite of losToEcf. + csm::ImageCoord computeViewingPixel( + const double& time, // The time to use the EO at + const csm::EcefCoord& groundPoint, // The ground coordinate + const std::vector& adj // Parameter Adjustments for partials + ) const; + // The linear approximation for the sensor model is used as the starting point // for iterative rigorous calculations. void computeLinearApproximation( diff --git a/src/UsgsAstroLsSensorModel.cpp b/src/UsgsAstroLsSensorModel.cpp index 7a2d392746d2513b644fceb47ef5ea2c7150c4cc..1d4cb5f1aab199db08a066fb405267fc8fdf85c7 100644 --- a/src/UsgsAstroLsSensorModel.cpp +++ b/src/UsgsAstroLsSensorModel.cpp @@ -200,170 +200,120 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage( double* achieved_precision, csm::WarningList* warnings) const { - //# func_description - // Computes line and sample given the ground coordinates in ECF cs. - // The solution is iterative and repeatedly calls the routine - // imageToPlane. If convergence is not achieved, a warning is issued. - // This method uses Newton-Raphson method on planes to iterate. - - // Initialize variables - const int MKTR = 20; - const double DELTA_IMAGE = 0.1; - double preSquare = desired_precision * desired_precision; - int mode = -1; - - // Initialize line, sample to center, CSM convention - csm::ImageCoord ip; - computeLinearApproximation(ground_pt, ip); - double lineTemp = ip.line; - double sampleTemp = ip.samp; - - // Compute ground elevation - - double height, aPrec; - double x = ground_pt.x; - double y = ground_pt.y; - double z = ground_pt.z; - computeElevation(x, y, z, height, aPrec, desired_precision); - - // Compute ground delta for the seed line and sample. - - double xSeed = x; - double ySeed = y; - double zSeed = z; - imageToPlane(lineTemp, sampleTemp, height, adj, - xSeed, ySeed, zSeed, mode); - double dx = x - xSeed; - double dy = y - ySeed; - double dz = z - zSeed; - double len = dx * dx + dy * dy + dz * dz; - - // Iterate until ground delta is acceptable. - - double xLine, yLine, zLine; - double xPLine, yPLine, zPLine; - double xSample, ySample, zSample; - double xPSamp, yPSamp, zPSamp; - double dLine, dSamp, det; - int ktr = 0; - while (preSquare < len && ktr < MKTR) - { - // Check for convergence - mode = -1; - ktr++; + // Search for the line, sample coordinate that viewed a given ground point. + // This method uses an iterative bisection method to search for the image + // line. + // + // For a given search window, this routine involves projecting the + // ground point onto the focal plane based on the instrument orientation + // at the start and end of the search window. Then, it computes the focal + // plane intersection at a mid-point of the search window. Then, it reduces + // the search window based on the signs of the intersected line offsets from + // the center of the ccd. For example, if the line offset is -145 at the + // start of the window, 10 at the mid point, and 35 at the end of the search + // window, the window will be reduced to the start of the old window to the + // middle of the old window. + // + // In order to achieve faster convergence, the mid point is calculated + // using the False Position Method instead of simple bisection. This method + // uses the zero of the line between the two ends of the search window for + // the mid point instead of a simple bisection. In most cases, this will + // converge significantly faster, but it can be slower than simple bisection + // in some situations. + + // Start bisection search on the image lines + double sampCtr = _data.m_TotalSamples / 2.0; + double firstTime = getImageTime(csm::ImageCoord(0.0, sampCtr)); + double lastTime = getImageTime(csm::ImageCoord(_data.m_TotalLines, sampCtr)); + double firstOffset = computeViewingPixel(firstTime, ground_pt, adj).line - 0.5; + double lastOffset = computeViewingPixel(lastTime, ground_pt, adj).line - 0.5; + + // Check if both offsets have the same sign. + // This means there is not guaranteed to be a zero. + if ((firstOffset > 0) != (lastOffset < 0)) { + throw csm::Error( + csm::Error::ALGORITHM, + "Ground point is not viewed by the image.", + "UsgsAstroLsSensorModel::groundToImage"); + } - // Compute the approximate Jacobian using finite differences on planes. - - // Compute partial of ground coordinates w.r.t. line - - xLine = xSeed; - yLine = ySeed; - zLine = zSeed; - mode = -1; - imageToPlane(lineTemp + DELTA_IMAGE, sampleTemp, height, adj, - xLine, yLine, zLine, mode); - xPLine = xLine - xSeed; - yPLine = yLine - ySeed; - zPLine = zLine - zSeed; - - // Compute partial of ground coordinates w.r.t. sample - - xSample = xSeed; - ySample = ySeed; - zSample = zSeed; - mode = -1; - imageToPlane(lineTemp, sampleTemp + DELTA_IMAGE, height, adj, - xSample, ySample, zSample, mode); - xPSamp = xSample - xSeed; - yPSamp = ySample - ySeed; - zPSamp = zSample - zSeed; - - // Compute delta line and sample by fixing one coordinate - // axis. This axis has the largest image ray component and - // not updated in the imageToPlane() routine. Since - // imageToPlane() updates only two of the three coordinates. - // Therefore, one of dx, dy, dz must be = 0.0 exactly. - // The following if else if should work just fine. - - // Compute the adjustment step by multiplying the ground delta by the - // inverse of the (approximate) Jacobian. - - if (0.0 == dx) - { - det = yPLine * zPSamp - yPSamp * zPLine; - dLine = (zPSamp * dy - yPSamp * dz); - dSamp = (yPLine * dz - zPLine * dy); - } - else if (0.0 == dy) - { - det = xPLine * zPSamp - xPSamp * zPLine; - dLine = (zPSamp * dx - xPSamp * dz); - dSamp = (xPLine * dz - zPLine * dx); + // 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 < _data.m_TotalLines) { + ++approxNextPoint.line; + } + else { + --approxNextPoint.line; + } + csm::EcefCoord approxIntersect = imageToGround(approxPoint, _data.m_RefElevation); + csm::EcefCoord approxNextIntersect = imageToGround(approxNextPoint, _data.m_RefElevation); + 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 bisection search for zero + for (int it = 0; it < 30; it++) { + double nextTime = ((firstTime * lastOffset) - (lastTime * firstOffset)) + / (lastOffset - firstOffset); + double nextOffset = computeViewingPixel(nextTime, ground_pt, adj).line - 0.5; + // We're looking for a zero, so check that either firstLine and middleLine have + // opposite signs, or middleLine and lastLine have opposite signs. + if ((firstOffset > 0) == (nextOffset < 0)) { + lastTime = nextTime; + lastOffset = nextOffset; } - else if (0.0 == dz) - { - det = xPLine * yPSamp - xPSamp * yPLine; - dLine = (yPSamp * dx - xPSamp * dy); - dSamp = (xPLine * dy - yPLine * dx); + else { + firstTime = nextTime; + firstOffset = nextOffset; } - else - { - throw csm::Error( - csm::Error::ALGORITHM, - "Undefined case.", - "UsgsAstroLsSensorModel::groundToImage"); - } - - // Update line and sample estimates - - if (0.0 == det) - { - throw csm::Error( - csm::Error::ALGORITHM, - "Divide by zero.", - "UsgsAstroLsSensorModel::groundToImage"); + if (fabs(lastOffset - firstOffset) < pixelPrec) { + break; } + } - // We have to divide by the determinant of the Jacobian here as part of - // the inverse calculation. - // The multiplication by DELTA_IMAGE is because the calculation of the - // ground partials with respect to sample and line does not divide by - // DELTA_IMAGE. This also means the determinant of the Jacobian should - // be divided by DELTA_IMAGE^2. So, dLine and dSamp should be multiplied - // by DELTA_IMAGE^2/DELTA_IMAGE, which is just DELTA_IMAGE. - lineTemp += (dLine / det * DELTA_IMAGE); - sampleTemp += (dSamp / det * DELTA_IMAGE); - - // Update ground delta - - xSeed = x; - ySeed = y; - zSeed = z; - mode = -1; - imageToPlane(lineTemp, sampleTemp, height, adj, - xSeed, ySeed, zSeed, mode); - dx = x - xSeed; - dy = y - ySeed; - dz = z - zSeed; - len = dx * dx + dy * dy + dz * dz; - - } // while (preSquare < len) + // Check that the desired precision was met + + double computedTime = ((firstTime * lastOffset) - (lastTime * firstOffset)) + / (lastOffset - firstOffset); + csm::ImageCoord calculatedPixel = computeViewingPixel(computedTime, ground_pt, adj); + // The computed viewing line is the detector line, so we need to convert that to image lines + auto referenceTimeIt = std::upper_bound(_data.m_IntTimeStartTimes.begin(), + _data.m_IntTimeStartTimes.end(), + computedTime); + if (referenceTimeIt != _data.m_IntTimeStartTimes.begin()) { + --referenceTimeIt; + } + size_t referenceIndex = std::distance(_data.m_IntTimeStartTimes.begin(), referenceTimeIt); + calculatedPixel.line += _data.m_IntTimeLines[referenceIndex] - 1 + + (computedTime - _data.m_IntTimeStartTimes[referenceIndex]) + / _data.m_IntTimes[referenceIndex]; + csm::EcefCoord calculatedPoint = imageToGround(calculatedPixel, _data.m_RefElevation); + 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; // If the final correction is greater than 10 meters, // the solution is not valid enough to report even with a warning - if (len > 100.0) - { + if (len > 100.0) { throw csm::Error( csm::Error::ALGORITHM, "Did not converge.", "UsgsAstroLsSensorModel::groundToImage"); } - if (achieved_precision) + if (achieved_precision) { *achieved_precision = sqrt(len); + } - if (warnings && (desired_precision > 0.0) && (preSquare < len)) - { + double preSquare = desired_precision * desired_precision; + if (warnings && (desired_precision > 0.0) && (preSquare < len)) { std::stringstream msg; msg << "Desired precision not achieved. "; msg << len << " " << preSquare << "\n"; @@ -373,7 +323,7 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage( "UsgsAstroLsSensorModel::groundToImage()")); } - return csm::ImageCoord(lineTemp, sampleTemp); + return calculatedPixel; } //*************************************************************************** @@ -2026,6 +1976,135 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel( } +//*************************************************************************** +// UsgsAstroLineScannerSensorModel::computeViewingPixel +//*************************************************************************** +csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel( + const double& time, + const csm::EcefCoord& groundPoint, + const std::vector& adj) const +{ + // Get the exterior orientation + double xc, yc, zc, vx, vy, vz; + getAdjSensorPosVel(time, adj, xc, yc, zc, vx, vy, vz); + + // Compute the look vector + double bodyLookX = groundPoint.x - xc; + double bodyLookY = groundPoint.y - yc; + double bodyLookZ = groundPoint.z - zc; + + // Rotate the look vector into the camera reference frame + int nOrder = 8; + if (_data.m_PlatformFlag == 0) + nOrder = 4; + int nOrderQuat = nOrder; + if (_data.m_NumQuaternions < 6 && nOrder == 8) + nOrderQuat = 4; + double q[4]; + lagrangeInterp( + _data.m_NumQuaternions, &_data.m_Quaternions[0], _data.m_T0Quat, _data.m_DtQuat, + time, 4, nOrderQuat, q); + double norm = sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]); + // Divide by the negative norm for 0 through 2 to invert the quaternion + q[0] /= -norm; + q[1] /= -norm; + q[2] /= -norm; + q[3] /= norm; + double bodyToCamera[9]; + bodyToCamera[0] = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3]; + bodyToCamera[1] = 2 * (q[0] * q[1] - q[2] * q[3]); + bodyToCamera[2] = 2 * (q[0] * q[2] + q[1] * q[3]); + bodyToCamera[3] = 2 * (q[0] * q[1] + q[2] * q[3]); + bodyToCamera[4] = -q[0] * q[0] + q[1] * q[1] - q[2] * q[2] + q[3] * q[3]; + bodyToCamera[5] = 2 * (q[1] * q[2] - q[0] * q[3]); + bodyToCamera[6] = 2 * (q[0] * q[2] - q[1] * q[3]); + bodyToCamera[7] = 2 * (q[1] * q[2] + q[0] * q[3]); + bodyToCamera[8] = -q[0] * q[0] - q[1] * q[1] + q[2] * q[2] + q[3] * q[3]; + double cameraLookX = bodyToCamera[0] * bodyLookX + + bodyToCamera[1] * bodyLookY + + bodyToCamera[2] * bodyLookZ; + double cameraLookY = bodyToCamera[3] * bodyLookX + + bodyToCamera[4] * bodyLookY + + bodyToCamera[5] * bodyLookZ; + double cameraLookZ = bodyToCamera[6] * bodyLookX + + bodyToCamera[7] * bodyLookY + + bodyToCamera[8] * bodyLookZ; + + // Invert the attitude correction + double aTime = time - _data.m_T0Quat; + double euler[3]; + double nTime = aTime / _data.m_HalfTime; + double nTime2 = nTime * nTime; + euler[0] = + (getValue(6, adj) + getValue(9, adj)* nTime + getValue(12, adj)* nTime2) / _data.m_FlyingHeight; + euler[1] = + (getValue(7, adj) + getValue(10, adj)* nTime + getValue(13, adj)* nTime2) / _data.m_FlyingHeight; + euler[2] = + (getValue(8, adj) + getValue(11, adj)* nTime + getValue(14, adj)* nTime2) / _data.m_HalfSwath; + double cos_a = cos(euler[0]); + double sin_a = sin(euler[0]); + double cos_b = cos(euler[1]); + double sin_b = sin(euler[1]); + double cos_c = cos(euler[2]); + double sin_c = sin(euler[2]); + double attCorr[9]; + attCorr[0] = cos_b * cos_c; + attCorr[1] = -cos_a * sin_c + sin_a * sin_b * cos_c; + attCorr[2] = sin_a * sin_c + cos_a * sin_b * cos_c; + attCorr[3] = cos_b * sin_c; + attCorr[4] = cos_a * cos_c + sin_a * sin_b * sin_c; + attCorr[5] = -sin_a * cos_c + cos_a * sin_b * sin_c; + attCorr[6] = -sin_b; + attCorr[7] = sin_a * cos_b; + attCorr[8] = cos_a * cos_b; + double adjustedLookX = attCorr[0] * cameraLookX + + attCorr[3] * cameraLookY + + attCorr[6] * cameraLookZ; + double adjustedLookY = attCorr[1] * cameraLookX + + attCorr[4] * cameraLookY + + attCorr[7] * cameraLookZ; + double adjustedLookZ = attCorr[2] * cameraLookX + + attCorr[5] * cameraLookY + + attCorr[8] * cameraLookZ; + + // Invert the boresight correction + double correctedLookX = _data.m_MountingMatrix[0] * adjustedLookX + + _data.m_MountingMatrix[3] * adjustedLookY + + _data.m_MountingMatrix[6] * adjustedLookZ; + double correctedLookY = _data.m_MountingMatrix[1] * adjustedLookX + + _data.m_MountingMatrix[4] * adjustedLookY + + _data.m_MountingMatrix[7] * adjustedLookZ; + double correctedLookZ = _data.m_MountingMatrix[2] * adjustedLookX + + _data.m_MountingMatrix[5] * adjustedLookY + + _data.m_MountingMatrix[8] * adjustedLookZ; + + // Convert to focal plane coordinate + double lookScale = _data.m_Focal / correctedLookZ; + double focalX = correctedLookX * lookScale; + double focalY = correctedLookY * lookScale; + + // TODO invert distortion here + // We probably only want to try and invert the distortion if we are + // reasonably close to the actual CCD because the distortion equations are + // sometimes only well behaved close to the CCD. + + // Convert to detector line and sample + double detectorLine = _data.m_ITransL[0] + + _data.m_ITransL[1] * focalX + + _data.m_ITransL[2] * focalY; + double detectorSample = _data.m_ITransS[0] + + _data.m_ITransS[1] * focalX + + _data.m_ITransS[2] * focalY; + + // Convert to image sample line + double line = detectorLine + _data.m_DetectorLineOrigin - _data.m_DetectorLineOffset + - _data.m_OffsetLines + 0.5; + double sample = (detectorSample + _data.m_DetectorSampleOrigin - _data.m_StartingSample) + / _data.m_DetectorSampleSumming - _data.m_OffsetSamples + 0.5; + return csm::ImageCoord(line, sample); +} + + //*************************************************************************** // UsgsAstroLineScannerSensorModel::computeLinearApproximation //*************************************************************************** diff --git a/tests/test_genericls.py b/tests/test_genericls.py index ba4d28045e6c5b6aaba65fb5f951d3cda8207987..4f44b752db84bbf369a7b3e65d03491971bb892b 100644 --- a/tests/test_genericls.py +++ b/tests/test_genericls.py @@ -58,89 +58,119 @@ class TestCTX: class TestHRSCNadir: @pytest.mark.parametrize('image, ground',[ - ((0.5, 0.5, 0), (985377.89225802, 3243484.7261571, 206009.0045894)), - ((0.5, 5175.5, 0), (849592.95667544, 3281524.2988248, 208281.94748872)), - ((25503.5, 0.5, 0), (995812.30843397, 3161653.5178064, 734845.5368816)), - ((25503.5, 5175.5, 0), (817255.33536118, 3211512.2358805, 738854.37876771)), - ((12751.5, 2587.5, 0), (914645.41695902, 3237864.2204448, 459626.50243137)) + # ((0.5, 0.5, 0), (985377.89225802, 3243484.7261571, 206009.0045894)), + # ((0.5, 5175.5, 0), (849592.95667544, 3281524.2988248, 208281.94748872)), + # ((25503.5, 0.5, 0), (995812.30843397, 3161653.5178064, 734845.5368816)), + # ((25503.5, 5175.5, 0), (817255.33536118, 3211512.2358805, 738854.37876771)), + # ((12751.5, 2587.5, 0), (914645.41695902, 3237864.2204448, 459626.50243137)) + ((12751.5, 2587.5, 0), (914645.4170054727, 3237864.2203280977, 459626.5031521894)), + ((0.5, 0.5, 0), (985377.8921162762, 3243484.7262333957, 206009.00407163633)), + ((0.5, 5175.5, 0), (849592.9565413876, 3281524.298893377, 208281.94696064907)), + ((25503.5, 0.5, 0), (995811.943316327, 3161654.168981308, 734843.2570738876)), + ((25503.5, 5175.5, 0), (817254.9608542051, 3211512.854057241, 738852.1327129134)) ]) def test_image_to_ground(self, hrsc_nadir_model, image, ground): gx, gy, gz = ground x, y, z = hrsc_nadir_model.imageToGround(*image) - assert x == pytest.approx(gx, abs=20) - assert y == pytest.approx(gy, abs=20) - assert z == pytest.approx(gz, abs=20) + assert x == pytest.approx(gx) + assert y == pytest.approx(gy) + assert z == pytest.approx(gz) @pytest.mark.parametrize('image, ground',[ - ((0.5, 0.5, 0), (985377.89225802, 3243484.7261571, 206009.0045894)), - ((0.5, 5175.5, 0), (849592.95667544, 3281524.2988248, 208281.94748872)), - ((25503.5, 0.5, 0), (995812.30843397, 3161653.5178064, 734845.5368816)), - ((25503.5, 5175.5, 0), (817255.33536118, 3211512.2358805, 738854.37876771)), - ((12751.5, 2587.5, 0), (914645.41695902, 3237864.2204448, 459626.50243137)) + # ((0.5, 0.5, 0), (985377.89225802, 3243484.7261571, 206009.0045894)), + # ((0.5, 5175.5, 0), (849592.95667544, 3281524.2988248, 208281.94748872)), + # ((25503.5, 0.5, 0), (995812.30843397, 3161653.5178064, 734845.5368816)), + # ((25503.5, 5175.5, 0), (817255.33536118, 3211512.2358805, 738854.37876771)), + # ((12751.5, 2587.5, 0), (914645.41695902, 3237864.2204448, 459626.50243137)) + ((12751.5, 2587.5, 0), (914645.4170054727, 3237864.2203280977, 459626.5031521894)), + ((0.5, 0.5, 0), (985377.8921162762, 3243484.7262333957, 206009.00407163633)), + ((0.5, 5175.5, 0), (849592.9565413876, 3281524.298893377, 208281.94696064907)), + ((25503.5, 0.5, 0), (995811.943316327, 3161654.168981308, 734843.2570738876)), + ((25503.5, 5175.5, 0), (817254.9608542051, 3211512.854057241, 738852.1327129134)) ]) def test_ground_to_image(self, hrsc_nadir_model, image, ground): y, x = hrsc_nadir_model.groundToImage(*ground) iy, ix, _ = image - assert x == pytest.approx(ix, abs=1) - assert y == pytest.approx(iy, abs=1) + assert x == pytest.approx(ix, abs=0.001) + assert y == pytest.approx(iy, abs=0.001) class TestHRSCStereo1: @pytest.mark.parametrize('image, ground',[ - ((0.5, 0.5, 0), (979590.30174957, 3242179.3919958, 249088.83886381)), - ((0.5, 2583.5, 0), (855179.96698611, 3277243.42159, 248427.66616907)), - ((12503.5, 0.5, 0), (984596.8817312, 3152386.5706555, 787256.94481508)), - ((12503.5, 2583.5, 0), (824462.92414568, 3197824.047306, 787981.15939371)), - ((6251.5, 1291.5, 0), (913762.62459356, 3231516.1914323, 503425.73182682)) + # ((0.5, 0.5, 0), (979590.30174957, 3242179.3919958, 249088.83886381)), + # ((0.5, 2583.5, 0), (855179.96698611, 3277243.42159, 248427.66616907)), + # ((12503.5, 0.5, 0), (984596.8817312, 3152386.5706555, 787256.94481508)), + # ((12503.5, 2583.5, 0), (824462.92414568, 3197824.047306, 787981.15939371)), + # ((6251.5, 1291.5, 0), (913762.62459356, 3231516.1914323, 503425.73182682)) + ((6251.5, 1291.5, 0), (913776.1812336871, 3231512.357457191, 503425.7355612734)), + ((0.5, 0.5, 0), (979602.4581277388, 3242175.7082415046, 249088.97789060202)), + ((0.5, 2583.5, 0), (855191.9410074502, 3277240.29787523, 248427.65484663288)), + ((12503.5, 0.5, 0), (984612.2803597612, 3152382.0160443066, 787255.9359414595)), + ((12503.5, 2583.5, 0), (824478.2804830034, 3197820.3642705963, 787980.0517643926)) ]) def test_image_to_ground(self, hrsc_stereo_1_model, image, ground): gx, gy, gz = ground x, y, z = hrsc_stereo_1_model.imageToGround(*image) - assert x == pytest.approx(gx, abs=100) - assert y == pytest.approx(gy, abs=100) - assert z == pytest.approx(gz, abs=100) + assert x == pytest.approx(gx) + assert y == pytest.approx(gy) + assert z == pytest.approx(gz) @pytest.mark.parametrize('image, ground',[ - ((0.5, 0.5, 0), (979590.30174957, 3242179.3919958, 249088.83886381)), - ((0.5, 2583.5, 0), (855179.96698611, 3277243.42159, 248427.66616907)), - ((12503.5, 0.5, 0), (984596.8817312, 3152386.5706555, 787256.94481508)), - ((12503.5, 2583.5, 0), (824462.92414568, 3197824.047306, 787981.15939371)), - ((6251.5, 1291.5, 0), (913762.62459356, 3231516.1914323, 503425.73182682)) + # ((0.5, 0.5, 0), (979590.30174957, 3242179.3919958, 249088.83886381)), + # ((0.5, 2583.5, 0), (855179.96698611, 3277243.42159, 248427.66616907)), + # ((12503.5, 0.5, 0), (984596.8817312, 3152386.5706555, 787256.94481508)), + # ((12503.5, 2583.5, 0), (824462.92414568, 3197824.047306, 787981.15939371)), + # ((6251.5, 1291.5, 0), (913762.62459356, 3231516.1914323, 503425.73182682)) + ((6251.5, 1291.5, 0), (913776.1812336871, 3231512.357457191, 503425.7355612734)), + ((0.5, 0.5, 0), (979602.4581277388, 3242175.7082415046, 249088.97789060202)), + ((0.5, 2583.5, 0), (855191.9410074502, 3277240.29787523, 248427.65484663288)), + ((12503.5, 0.5, 0), (984612.2803597612, 3152382.0160443066, 787255.9359414595)), + ((12503.5, 2583.5, 0), (824478.2804830034, 3197820.3642705963, 787980.0517643926)) ]) def test_ground_to_image(self, hrsc_stereo_1_model, image, ground): y, x = hrsc_stereo_1_model.groundToImage(*ground) iy, ix, _ = image - assert x == pytest.approx(ix, abs=1) - assert y == pytest.approx(iy, abs=1) + assert x == pytest.approx(ix, abs=0.001) + assert y == pytest.approx(iy, abs=0.001) class TestHRSCStereo2: @pytest.mark.parametrize('image, ground',[ - ((0.5, 0.5, 0), (994968.78141471, 3243624.373581, 150910.83677465)), - ((0.5, 2583.5, 0), (840154.182443, 3286852.5036334, 156704.92657188)), - ((14263.5, 0.5, 0), (1015683.4735728, 3170802.9252193, 665761.34487006)), - ((14263.5, 2583.5, 0), (802816.11887483, 3229529.8393814, 674042.87178595)), - ((7131.5, 1291.5, 0), (915544.01372502, 3243112.89855, 419540.13516932)) + # ((0.5, 0.5, 0), (994968.78141471, 3243624.373581, 150910.83677465)), + # ((0.5, 2583.5, 0), (840154.182443, 3286852.5036334, 156704.92657188)), + # ((14263.5, 0.5, 0), (1015683.4735728, 3170802.9252193, 665761.34487006)), + # ((14263.5, 2583.5, 0), (802816.11887483, 3229529.8393814, 674042.87178595)), + # ((7131.5, 1291.5, 0), (915544.01372502, 3243112.89855, 419540.13516932)) + ((7131.5, 1291.5, 0), (915561.6002628063, 3243108.0221202467, 419539.4600493548)), + ((0.5, 0.5, 0), (994983.8770025282, 3243619.7787126494, 150910.07885617757)), + ((0.5, 2583.5, 0), (840169.1316321647, 3286848.7002639165, 156704.55683200277)), + ((14263.5, 0.5, 0), (1015704.1721659054, 3170796.569969528, 665760.0501533676)), + ((14263.5, 2583.5, 0), (802836.7715344077, 3229524.77329203, 674042.5500676908)) ]) def test_image_to_ground(self, hrsc_stereo_2_model, image, ground): gx, gy, gz = ground x, y, z = hrsc_stereo_2_model.imageToGround(*image) - assert x == pytest.approx(gx, abs=100) - assert y == pytest.approx(gy, abs=100) - assert z == pytest.approx(gz, abs=100) + assert x == pytest.approx(gx) + assert y == pytest.approx(gy) + assert z == pytest.approx(gz) @pytest.mark.parametrize('image, ground',[ - ((0.5, 0.5, 0), (994968.78141471, 3243624.373581, 150910.83677465)), - ((0.5, 2583.5, 0), (840154.182443, 3286852.5036334, 156704.92657188)), - ((14263.5, 0.5, 0), (1015683.4735728, 3170802.9252193, 665761.34487006)), - ((14263.5, 2583.5, 0), (802816.11887483, 3229529.8393814, 674042.87178595)), - ((7131.5, 1291.5, 0), (915544.01372502, 3243112.89855, 419540.13516932)) + # ((0.5, 0.5, 0), (994968.78141471, 3243624.373581, 150910.83677465)), + # ((0.5, 2583.5, 0), (840154.182443, 3286852.5036334, 156704.92657188)), + # ((14263.5, 0.5, 0), (1015683.4735728, 3170802.9252193, 665761.34487006)), + # ((14263.5, 2583.5, 0), (802816.11887483, 3229529.8393814, 674042.87178595)), + # ((7131.5, 1291.5, 0), (915544.01372502, 3243112.89855, 419540.13516932)) + ((7131.5, 1291.5, 0), (915561.6002628063, 3243108.0221202467, 419539.4600493548)), + ((0.5, 0.5, 0), (994983.8770025282, 3243619.7787126494, 150910.07885617757)), + ((0.5, 2583.5, 0), (840169.1316321647, 3286848.7002639165, 156704.55683200277)), + ((14263.5, 0.5, 0), (1015704.1721659054, 3170796.569969528, 665760.0501533676)), + ((14263.5, 2583.5, 0), (802836.7715344077, 3229524.77329203, 674042.5500676908)) ]) def test_ground_to_image(self, hrsc_stereo_2_model, image, ground): y, x = hrsc_stereo_2_model.groundToImage(*ground) iy, ix, _ = image - assert x == pytest.approx(ix, abs=1) - assert y == pytest.approx(iy, abs=1) + assert x == pytest.approx(ix, abs=0.001) + assert y == pytest.approx(iy, abs=0.001)