geocel/rasterize/Raytracer.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file geocel/rasterize/Raytracer.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/math/Algorithms.hh"
0010
0011 #include "ImageLineView.hh"
0012
0013 namespace celeritas
0014 {
0015 //---------------------------------------------------------------------------//
0016 /*!
0017  * Trace each pixel along a line.
0018  *
0019  * \tparam GTV Geometry Track View
0020  * \tparam F   Calculate the pixel value (result_type) given GTV
0021  */
0022 template<class GTV, class F>
0023 class Raytracer
0024 {
0025   public:
0026     //!@{
0027     //! \name Type aliases
0028     using result_type = int;
0029     //!@}
0030
0031   public:
0032     // Construct from geo track view, function to calculate ID, and image
0033     inline CELER_FUNCTION
0034     Raytracer(GTV&& geo, F&& calc_id, ImageLineView const& image);
0035
0036     // Calculate the value for the next pixel
0037     inline CELER_FUNCTION result_type operator()(size_type pix);
0038
0039   private:
0040     //// DATA ////
0041
0042     GTV geo_;
0043     F calc_id_;
0044     ImageLineView const& image_;
0045
0046     size_type pixel_{invalid_pixel()};  //!< Current pixel
0047     real_type distance_{0};  //!< Distance to next boundary
0048     int cur_id_{-1};  //!< Current ID
0049
0050     //// HELPER FUNCTIONS ////
0051
0052     // Initialize the geometry at this pixel index
0053     inline CELER_FUNCTION void initialize_at_pixel(size_type pix);
0054
0055     // Find the next step
0056     inline CELER_FUNCTION void find_next_step();
0057
0058     //! Sentinel value for needing the pixel to be reset
0059     static CELER_CONSTEXPR_FUNCTION size_type invalid_pixel()
0060     {
0061         return static_cast<size_type>(-1);
0062     }
0063
0064     //! Stop tracking the pixel if this many crossings occur
0065     static CELER_CONSTEXPR_FUNCTION short int max_crossings_per_pixel()
0066     {
0067         return 32;
0068     }
0069 };
0070
0071 //---------------------------------------------------------------------------//
0072 // DEDUCTION GUIDES
0073 //---------------------------------------------------------------------------//
0074
0075 template<class GTV, class F>
0076 CELER_FUNCTION Raytracer(GTV&&, F&&, ImageLineView const&) -> Raytracer<GTV, F>;
0077
0078 //---------------------------------------------------------------------------//
0079 // INLINE DEFINITIONS
0080 //---------------------------------------------------------------------------//
0081 /*!
0082  * Construct from geo track view, function to calculate ID, and image.
0083  */
0084 template<class GTV, class F>
0085 CELER_FUNCTION
0086 Raytracer<GTV, F>::Raytracer(GTV&& geo, F&& calc_id, ImageLineView const& image)
0087     : geo_{celeritas::forward<GTV>(geo)}
0088     , calc_id_{celeritas::forward<F>(calc_id)}
0089     , image_{image}
0090 {
0091 }
0092
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * Calculate the value for a single pixel.
0096  */
0097 template<class GTV, class F>
0098 CELER_FUNCTION auto Raytracer<GTV, F>::operator()(size_type pix) -> result_type
0099 {
0100     if (pix != pixel_)
0101     {
0102         this->initialize_at_pixel(pix);
0103         if (pixel_ == invalid_pixel())
0104         {
0105             // Pixel starts outside the geometry
0106             return cur_id_;
0107         }
0108     }
0109
0110     // Remaining distance to the end of the pixel
0111     real_type pix_dist = image_.pixel_width();
0112     // Distance associated with the largest segment so far
0113     real_type max_dist = 0;
0114     int max_id = cur_id_;
0115     // Countdown to limit stuck tracks or overdetailed pixels
0116     auto abort_counter = this->max_crossings_per_pixel();
0117
0118     while (cur_id_ >= 0 && distance_ < pix_dist)
0119     {
0120         // Move to geometry boundary
0121         if (max_id == cur_id_)
0122         {
0123             max_dist += distance_;
0124         }
0125         else if (distance_ > max_dist)
0126         {
0127             max_dist = distance_;
0128             max_id = cur_id_;
0129         }
0130
0131         // Update pixel and boundary distance
0132         pix_dist -= distance_;
0133         distance_ = 0;
0134
0135         // Cross surface and update post-crossing ID
0136         geo_.move_to_boundary();
0137         geo_.cross_boundary();
0138
0139         if (--abort_counter == 0)
0140         {
0141             // Give up and move to end of pixel, reinitialize on the next call
0142             pix_dist = 0;
0143         }
0144         if (pix_dist > 0)
0145         {
0146             // Update next distance and current ID for the remaining pixel
0147             this->find_next_step();
0148         }
0149     }
0150
0151     if (pix_dist > 0)
0152     {
0153         // Move to pixel boundary
0154         distance_ -= pix_dist;
0155         if (pix_dist > max_dist)
0156         {
0157             max_dist = pix_dist;
0158             max_id = cur_id_;
0159         }
0160     }
0161
0162     if (abort_counter != 0)
0163     {
0164         ++pixel_;
0165     }
0166     else
0167     {
0168         pixel_ = invalid_pixel();
0169     }
0170
0171     return max_id;
0172 }
0173
0174 //---------------------------------------------------------------------------//
0175 /*!
0176  * Initialize the geometry at the given pixel position.
0177  */
0178 template<class GTV, class F>
0179 CELER_FUNCTION void Raytracer<GTV, F>::initialize_at_pixel(size_type pix)
0180 {
0181     CELER_EXPECT(pix < image_.max_index());
0182
0183     GeoTrackInitializer init{image_.start_pos(), image_.start_dir()};
0184     axpy(pix * image_.pixel_width(), init.dir, &init.pos);
0185     geo_ = init;
0186     pixel_ = pix;
0187     distance_ = 0;
0188     this->find_next_step();
0189 }
0190
0191 //---------------------------------------------------------------------------//
0192 /*!
0193  * Find the distance to the next step.
0194  */
0195 template<class GTV, class F>
0196 CELER_FUNCTION void Raytracer<GTV, F>::find_next_step()
0197 {
0198     CELER_EXPECT(pixel_ != invalid_pixel());
0199     CELER_EXPECT(distance_ <= 0);
0200     if (geo_.is_outside())
0201     {
0202         // Skip this pixel since not all navigation engines can trace from
0203         // outside the geometry
0204         pixel_ = invalid_pixel();
0205         distance_ = numeric_limits<real_type>::infinity();
0206     }
0207     else
0208     {
0209         // Search to distances just past the end of the saved pixels (so the
0210         // last point isn't coincident with a boundary)
0211         distance_ = geo_.find_next_step((image_.max_index() + 1 - pixel_)
0212                                         * image_.pixel_width())
0213                         .distance;
0214     }
0215     cur_id_ = this->calc_id_(geo_);
0216
0217     CELER_ENSURE(distance_ > 0);
0218 }
0219
0220 //---------------------------------------------------------------------------//
0221 }  // namespace celeritas