Merge pull request #718 from lsst/tickets/DM-41345

DM-41345: handle non-finite positions when calculating covariance
lsst · Dec 6, 2023 · 2ea0927 · 2ea0927
2 parents 2c44ac4 + 2e4cf33
commit 2ea0927
Showing 1 changed file with 13 additions and 13 deletions.
diff --git a/src/table/wcsUtils.cc b/src/table/wcsUtils.cc
@@ -22,6 +22,7 @@
  * see <http://www.lsstcorp.org/LegalNotices/>.
  */
 
+#include <cmath>
 #include <memory>
 #include <vector>
 
@@ -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())) {
+        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);
+    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()) {
@@ -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