diff --git a/six/modules/c++/scene/source/ProjectionModel.cpp b/six/modules/c++/scene/source/ProjectionModel.cpp index fec40d54bd..988b7420b9 100644 --- a/six/modules/c++/scene/source/ProjectionModel.cpp +++ b/six/modules/c++/scene/source/ProjectionModel.cpp @@ -459,22 +459,24 @@ math::linear::MatrixMxN<2, 2> ProjectionModel::slantToImagePartials( mTimeCOAPoly(imageGridPoint.row, imageGridPoint.col); const Vector3 rARP = mARPPoly(timeCOA); const Vector3 vARP = mARPVelPoly(timeCOA); - const Vector3 imageGridPointECEF = imageGridToECEF(imageGridPoint); - Vector3 slantRange = imageGridPointECEF - rARP; + Vector3 slantRange = mSCP - rARP; slantRange.normalize(); - Vector3 slantNormal = math::linear::cross(slantRange, vARP); + Vector3 slantNormal = math::linear::cross(vARP, slantRange) * mLookDir; slantNormal.normalize(); Vector3 slantAzimuth = math::linear::cross(slantNormal, slantRange); slantAzimuth.normalize(); + + // Second, map image grid point to the slant plane and compute finite differences + const Vector3 refPoint = imageToScene(imageGridPoint, mSCP, slantNormal); const types::RowCol rangePerturb = - sceneToImage(imageGridPointECEF + delta * slantRange); + sceneToImage(refPoint + delta * slantRange); const types::RowCol azimuthPerturb = - sceneToImage(imageGridPointECEF + delta * slantAzimuth); + sceneToImage(refPoint + delta * slantAzimuth); math::linear::MatrixMxN<2, 2> partials(0.0); - partials[0][0] = (imageGridPoint.row - rangePerturb.row) / delta; - partials[0][1] = (imageGridPoint.row - azimuthPerturb.row) / delta; - partials[1][0] = (imageGridPoint.col - rangePerturb.col) / delta; - partials[1][1] = (imageGridPoint.col - azimuthPerturb.col) / delta; + partials[0][0] = (rangePerturb.row - imageGridPoint.row) / delta; + partials[0][1] = (azimuthPerturb.row - imageGridPoint.row) / delta; + partials[1][0] = (rangePerturb.col - imageGridPoint.col) / delta; + partials[1][1] = (azimuthPerturb.col - imageGridPoint.col) / delta; return partials; }