Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-41345: handle non-finite positions when calculating covariance #718

Merged
merged 3 commits into from
Nov 2, 2023
Merged

DM-41345: handle non-finite positions when calculating covariance #718

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6c621bc
clang-format pass
parejkoj Oct 27, 2023
59d70a0
Make constants const and static
parejkoj Oct 27, 2023
2e4cf33
Properly handle non-finite position
parejkoj Oct 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 13 additions & 13 deletions src/table/wcsUtils.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
* see <http://www.lsstcorp.org/LegalNotices/>.
*/

#include <cmath>
#include <memory>
#include <vector>

Expand Down Expand Up @@ -91,35 +92,35 @@ void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList) {
}
}

Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const& wcs, lsst::geom::Point2D center, Eigen::Matrix2f err) {
Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const &wcs, lsst::geom::Point2D center,
Eigen::Matrix2f err) {
if (!std::isfinite(center.getX()) || !std::isfinite(center.getY())) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this would be the end result anyway, but I agree it's better to just skip going forward in this situation.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, if the center is non-finite, the computed skyCenter will be NaN and the call to geom::makeSkyWcs will raise because it's not a valid WCS.

return Eigen::Matrix2f::Constant(NAN);
}
// Get the derivative of the pixel-to-sky transformation, then use it to
// propagate the centroid uncertainty to coordinate uncertainty. Note that
// the calculation is done in arcseconds, then converted to radians in
// order to achieve higher precision.
double scale = 1.0 / 3600.0;
Eigen::Matrix2d cdMatrix { {scale, 0}, {0, scale}};
const static double scale = 1.0 / 3600.0;
const static Eigen::Matrix2d cdMatrix{{scale, 0}, {0, scale}};

lsst::geom::SpherePoint skyCenter = wcs.pixelToSky(center);

auto localGnomonicWcs = geom::makeSkyWcs(center, skyCenter, cdMatrix);
auto measurementToLocalGnomonic = wcs.getTransform()->then(
*localGnomonicWcs->getTransform()->inverted()
);

auto measurementToLocalGnomonic = wcs.getTransform()->then(*localGnomonicWcs->getTransform()->inverted());

Eigen::Matrix2d localMatrix = measurementToLocalGnomonic->getJacobian(center);
Eigen::Matrix2f d = localMatrix.cast < float > () * scale * (lsst::geom::PI / 180.0);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know how those spaces got merged! Oops.

Eigen::Matrix2f d = localMatrix.cast<float>() * scale * (lsst::geom::PI / 180.0);

Eigen::Matrix2f skyCov = d * err * d.transpose();

// Multiply by declination correction matrix in order to get sigma(RA) * cos(Dec) for the uncertainty
// in RA, and cov(RA, Dec) * cos(Dec) for the RA/Dec covariance:
float cosDec = std::cos(skyCenter.getDec().asRadians());
Eigen::Matrix2f decCorr{ {cosDec, 0}, {0, 1.0}};
Eigen::Matrix2f decCorr{{cosDec, 0}, {0, 1.0}};
Eigen::Matrix2f skyCovCorr = decCorr * skyCov * decCorr;
return skyCovCorr;
}


template <typename SourceCollection>
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList, bool include_covariance) {
if (sourceList.empty()) {
Expand Down Expand Up @@ -164,8 +165,7 @@ template void updateRefCentroids(geom::SkyWcs const &, SimpleCatalog &);

template void updateSourceCoords(geom::SkyWcs const &, std::vector<std::shared_ptr<SourceRecord>> &,
bool include_covariance);
template void updateSourceCoords(geom::SkyWcs const &, SourceCatalog &,
bool include_covariance);
template void updateSourceCoords(geom::SkyWcs const &, SourceCatalog &, bool include_covariance);
/// @endcond

} // namespace table
Expand Down