Core/Seeding/HoughTransformTest.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/tools/old/interface.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Seeding/HoughTransformUtils.hpp"
0015 #include "Acts/Utilities/Logger.hpp"
0016
0017 #include <array>
0018 #include <format>
0019 #include <vector>
0020
0021 using namespace Acts;
0022
0023 namespace ActsTests {
0024
0025 auto logger = getDefaultLogger("UnitTests", Logging::VERBOSE);
0026
0027 struct DriftCircle {
0028   double y{0.};
0029   double z{0.};
0030   double rDrift{0.};
0031   double rDriftError{0.};
0032
0033   DriftCircle(const double _y, const double _z, const double _r,
0034               const double _rUncert)
0035       : y{_y}, z{_z}, rDrift{_r}, rDriftError{_rUncert} {}
0036 };
0037
0038 BOOST_AUTO_TEST_SUITE(SeedingSuite)
0039
0040 BOOST_AUTO_TEST_CASE(hough_transform_seeder) {
0041   // we are using the slope on yz plane with the y coordinate (hardcoded from
0042   // the csv MuonSimHit data)
0043   std::vector<std::pair<double, double>> simHits = {
0044       {-0.0401472 / 0.994974, -422.612}};
0045
0046   // Define the drift Circles
0047   constexpr double uncert{0.3};
0048   std::array<DriftCircle, 6> driftCircles{
0049       DriftCircle{-427.981, -225.541, 14.5202, uncert},
0050       DriftCircle{-412.964, -199.53, 1.66237, uncert},
0051       DriftCircle{-427.981, -173.519, 12.3176, uncert},
0052       DriftCircle{-427.981, 173.519, 1.5412, uncert},
0053       DriftCircle{-442.999, 199.53, 12.3937, uncert},
0054       DriftCircle{-427.981, 225.541, 3.77967, uncert}
0055
0056   };
0057
0058   // configure the binning of the hough plane
0059   HoughTransformUtils::HoughPlaneConfig planeCfg;
0060   planeCfg.nBinsX = 1000;
0061   planeCfg.nBinsY = 1000;
0062
0063   // instantiate the peak finder
0064   HoughTransformUtils::PeakFinders::IslandsAroundMaxConfig peakFinderCfg;
0065   peakFinderCfg.fractionCutoff = 0.7;
0066   peakFinderCfg.threshold = 3.;
0067   peakFinderCfg.minSpacingBetweenPeaks = {0., 30.};
0068
0069   // and map the hough plane to parameter ranges.
0070   // The first coordinate is tan(theta), the second is z0 in mm
0071   HoughTransformUtils::HoughAxisRanges axisRanges{-3., 3., -2000., 2000.};
0072
0073   // create the functions parametrising the hough space lines for drift circles.
0074   // Note that there are two solutions for each drift circle and angle
0075
0076   // left solution
0077   auto houghParamFromDCleft = [](double tanTheta, const DriftCircle& DC) {
0078     return DC.y - tanTheta * DC.z - DC.rDrift / std::cos(std::atan(tanTheta));
0079   };
0080   // right solution
0081   auto houghParamFromDCright = [](double tanTheta, const DriftCircle& DC) {
0082     return DC.y - tanTheta * DC.z + DC.rDrift / std::cos(std::atan(tanTheta));
0083   };
0084
0085   // create the function parametrising the drift radius uncertainty
0086   auto houghWidthFromDC = [](double, const DriftCircle& DC) {
0087     return std::min(DC.rDriftError * 3.,
0088                     1.0);  // scale reported errors up to at least 1mm or 3
0089                            // times the reported error as drift circle calib not
0090                            // fully reliable at this stage
0091   };
0092
0093   // instantiate the hough plane
0094   HoughTransformUtils::HoughPlane<GeometryIdentifier::Value> houghPlane(
0095       planeCfg);
0096
0097   // also instantiate the peak finder
0098   HoughTransformUtils::PeakFinders::IslandsAroundMax<GeometryIdentifier::Value>
0099       peakFinder(peakFinderCfg);
0100
0101   // loop over the true hits
0102   for (auto& sh : simHits) {
0103     houghPlane.reset();
0104
0105     for (std::size_t k = 0; k < driftCircles.size(); ++k) {
0106       auto dc = driftCircles[k];
0107
0108       houghPlane.fill<DriftCircle>(dc, axisRanges, houghParamFromDCleft,
0109                                    houghWidthFromDC, k);
0110       houghPlane.fill<DriftCircle>(dc, axisRanges, houghParamFromDCright,
0111                                    houghWidthFromDC, k);
0112     }
0113
0114     // now get the peaks
0115     auto maxima = peakFinder.findPeaks(houghPlane, axisRanges);
0116     for (auto& max : maxima) {
0117       // check the Hough Transforms results
0118       BOOST_CHECK_CLOSE(max.x, sh.first, 4.);
0119       BOOST_CHECK_CLOSE(max.y, sh.second, 0.2);
0120     }
0121   }
0122 }
0123
0124 BOOST_AUTO_TEST_CASE(hough_transform_sliding_window) {
0125   // Create a simple 10x10 Hough space and fill it with two triangular
0126   // (pyramid-shaped) distributions. Each distribution has its maximum
0127   // value at the specified bin coordinates and decreases linearly with
0128   // Manhattan distance. When distributions overlap we take the maximum
0129   // of the two contributions so that the global maximum remains the
0130   // requested peak value.
0131   const std::size_t nX = 10;
0132   const std::size_t nY = 10;
0133   HoughTransformUtils::HoughPlaneConfig config{nX, nY};
0134   HoughTransformUtils::HoughPlane<int> plane(config);
0135
0136   auto addTriangular = [&](const std::vector<std::array<std::size_t, 2>>& peaks,
0137                            HoughTransformUtils::YieldType peak) {
0138     for (std::size_t x = 0; x < nX; ++x) {
0139       for (std::size_t y = 0; y < nY; ++y) {
0140         HoughTransformUtils::YieldType val = 0.0f;
0141         for (auto [cx, cy] : peaks) {
0142           int dist = (x >= cx ? x - cx : cx - x) +
0143                      (y >= cy ? y - cy : cy - y);  // Manhattan distance
0144           val = std::max(
0145               val, peak - static_cast<HoughTransformUtils::YieldType>(dist));
0146         }
0147         plane.fillBin(x, y, x, y, val);
0148       }
0149     }
0150   };
0151
0152   // Add 3 triangular peaks with maxima 4 at specified locations
0153   std::vector<std::array<std::size_t, 2>> testPeaks({{4, 4}, {4, 5}, {2, 6}});
0154   addTriangular(testPeaks, 4.0f);
0155
0156   // Verify that the maxima at the requested locations are exactly 4
0157   BOOST_CHECK_EQUAL(plane.nHits(4, 4), 4.0f);
0158   BOOST_CHECK_EQUAL(plane.nHits(4, 5), 4.0f);
0159
0160   // config for unisolated max finding
0161   HoughTransformUtils::PeakFinders::SlidingWindowConfig cfg1{4,     0, 0,
0162                                                              false, 0, 0};
0163   auto peaks1 = slidingWindowPeaks(plane, cfg1);
0164   BOOST_CHECK_EQUAL(peaks1[0][0], 2);
0165   BOOST_CHECK_EQUAL(peaks1[0][1], 6);
0166   BOOST_CHECK_EQUAL(peaks1[1][0], 4);
0167   BOOST_CHECK_EQUAL(peaks1[1][1], 4);
0168   BOOST_CHECK_EQUAL(peaks1[2][0], 4);
0169   BOOST_CHECK_EQUAL(peaks1[2][1], 5);
0170
0171   // config for removing duplicates, no recentering
0172   HoughTransformUtils::PeakFinders::SlidingWindowConfig cfg2{4,     2, 2,
0173                                                              false, 0, 0};
0174   auto peaks2 = slidingWindowPeaks(plane, cfg2);
0175   BOOST_CHECK_EQUAL(peaks2.size(), 1);
0176   BOOST_CHECK_EQUAL(peaks2[0][0], 4);
0177   BOOST_CHECK_EQUAL(peaks2[0][1], 4);
0178
0179   // config for removing duplicates, with recentering
0180   HoughTransformUtils::PeakFinders::SlidingWindowConfig cfg3{4,    2, 2,
0181                                                              true, 2, 2};
0182   auto peaks3 = slidingWindowPeaks(plane, cfg3);
0183   BOOST_CHECK_EQUAL(peaks3.size(), 1);
0184   BOOST_CHECK_EQUAL(peaks3[0][0], 3);
0185   BOOST_CHECK_EQUAL(peaks3[0][1], 5);
0186
0187   // test behaviour at the edge of the plane (safety check)
0188   plane.reset();
0189   addTriangular({{1, 1}, {5, 5}, {8, 9}}, 4.0f);
0190   peaks3 = slidingWindowPeaks(plane, cfg3);
0191   BOOST_CHECK_EQUAL(peaks3.size(), 1);
0192
0193   {
0194     auto img1 =
0195         HoughTransformUtils::PeakFinders::hitsCountImage(plane, {1, 1}, 4, 4);
0196     // the image should reflect this content (peak at 1, 1)
0197     //      0 1 2 3 4 5 - indices
0198     //   +--------+
0199     //   |        |
0200     // 0 |  2 3 2 |
0201     // 1 |  3 4 3 |
0202     // 2 |  2 3 2 |
0203     //   +________+
0204     // content outside of valid plane indices is padded with zeros
0205     std::vector<unsigned char> expected(
0206         {0, 0, 0, 0, 0, 2, 3, 2, 0, 3, 4, 3, 0, 2, 3, 2});
0207     BOOST_CHECK_EQUAL(expected.size(), img1.size());
0208     BOOST_CHECK_EQUAL_COLLECTIONS(img1.begin(), img1.end(), expected.begin(),
0209                                   expected.end());
0210   }
0211   {
0212     // customized summary function: gets layers bitmask
0213     auto maskGetter = [](const HoughTransformUtils::HoughPlane<int>& p, int x,
0214                          int y) -> unsigned {
0215       unsigned mask = 0;
0216       for (unsigned l : p.layers(x, y)) {
0217         mask |= 1 << std::min(l, 16u);
0218         std::cout << x << " " << y << " " << l << "\n";
0219       }
0220       return mask;
0221     };
0222
0223     plane.reset();
0224     plane.fillBin(1, 1, 0, 1, 1);
0225     plane.fillBin(1, 1, 1, 2, 1);
0226     plane.fillBin(1, 1, 2, 3, 1);
0227     plane.fillBin(2, 2, 3, 3, 1);
0228
0229     auto img = HoughTransformUtils::PeakFinders::hitsCountImage<int, unsigned>(
0230         plane, {1, 1}, 4, 4, maskGetter);
0231
0232     std::vector<unsigned char> expected({0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0233                                          1 << 0x1 | 1 << 0x2 | 1 << 0x3, 0, 0,
0234                                          0, 0, 1 << 0x3});
0235
0236     BOOST_CHECK_EQUAL(expected.size(), img.size());
0237     BOOST_CHECK_EQUAL_COLLECTIONS(img.begin(), img.end(), expected.begin(),
0238                                   expected.end());
0239   }
0240 }
0241
0242 BOOST_AUTO_TEST_CASE(hough_plane_layers_hits) {
0243   // Create 1x1 Hough plane and fill it with hits. For each layer add one more
0244   // hit than to the previous one.
0245
0246   const std::size_t nX = 1;
0247   const std::size_t nY = 1;
0248   HoughTransformUtils::HoughPlaneConfig config{nX, nY};
0249   HoughTransformUtils::HoughPlane<std::uint8_t> plane(config);
0250
0251   std::uint8_t nHits = 0;
0252   static constexpr std::uint8_t nLayers = 10;
0253   for (std::uint8_t layer = 1; layer <= nLayers; ++layer) {
0254     // Add hits equal to the layer number
0255     for (std::uint8_t hit = 0; hit < layer; ++hit) {
0256       plane.fillBin(0, 0, nHits++, layer);
0257     }
0258   }
0259
0260   BOOST_CHECK_EQUAL(nLayers, plane.nLayers(0, 0));
0261   BOOST_CHECK_EQUAL(nHits, plane.nHits(0, 0));
0262 }
0263
0264 BOOST_AUTO_TEST_SUITE_END()
0265
0266 }  // namespace ActsTests