Merge branch 'p-liang-barsky' into 'master' (0cd1a6bf) · Commits · XRFLab109 / Allpix Squared

src/core/geometry/PixelDetectorModel.cpp

+7 −4

Original line number	Diff line number	Diff line
		@@ -117,10 +117,13 @@ ROOT::Math::XYZPoint PixelDetectorModel::getSensorIntercept(const ROOT::Math::XY
		// We have to be centered around the sensor box. This means we need to shift by the matrix center
		auto translation_local = ROOT::Math::Translation3D(static_cast<ROOT::Math::XYZVector>(getMatrixCenter()));

		try {
		auto intersection_point =
		LiangBarsky::closestIntersection(direction, translation_local.Inverse()(inside), getSensorSize());

		// Get intersection from Liang-Barsky line clipping and re-transform to local coordinates:
		return translation_local(LiangBarsky(direction, translation_local.Inverse()(inside), getSensorSize()));
		} catch(std::invalid_argument&) {
		if(intersection_point) {
		return translation_local(intersection_point.value());
		} else {
		return inside;
		}
		}

src/tools/liang_barsky.h

+81 −52

Original line number	Diff line number	Diff line
		@@ -2,6 +2,12 @@
		* @file
		* @brief Utility to perform Liang-Barsky clipping checks on volumes
		*
		* @see Liang, Y. D., and Barsky, B., "A New Concept and Method for Line Clipping",
		* ACM Transactions on Graphics, 3(1):1–22 for an in-depth explanation.
		* This method requires the position to be in the coordinate system of the box
		* to be tested for intersections, with the box center at its origin and the box sides aligned with the coordinate
		* axes.
		*
		* @copyright Copyright (c) 2022 CERN and the Allpix Squared authors.
		* This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
		* In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
		@@ -15,25 +21,23 @@
		#include <Math/Point3D.h>
		#include <Math/Vector3D.h>

		#include <optional>
		#include <utility>

		namespace allpix {

		class LiangBarsky {
		public:
		/**
		* @brief Liang–Barsky clipping of a line against faces of a box.
		*
		* See Liang, Y. D., and Barsky, B., "A New Concept and Method for Line Clipping", ACM Transactions on Graphics,
		* 3(1):1–22 for an in-depth explanation. This method requires the position to be in the coordinate system of the box
		* to be tested for intersections, with the box center at its origin and the box sides aligned with the coordinate
		* axes.
		*
		* @param direction Direction vector of the motion
		* @param position Original ("before") position to be considered
		* @brief Check intersection of a line defined by a point and a vector with a box
		* @param direction Defining vector of the line
		* @param position A point on that line
		* @param box Size of the box to calculate the intersections with
		* @return Closest intersection with box in the direction indicated by input vector
		*
		* @throws std::invalid_argument if no intersection of track segment with the box volume can be found in positive
		* direction from the given position.
		* @return Pair of signed distances from `position` to intersection points along the line in units of length of
		* `direction`, with sign of these distances meaning direction w.r.t. line-defining vector or std::nullopt if the
		* line has no intersection with the given box
		*/
		inline ROOT::Math::XYZPoint LiangBarsky(const ROOT::Math::XYZVector& direction,
		static std::optional<std::pair<double, double>> intersectionDistances(const ROOT::Math::XYZVector& direction,
		const ROOT::Math::XYZPoint& position,
		const ROOT::Math::XYZVector& box) {

		@@ -68,8 +72,33 @@ namespace allpix {
		clip(direction.Z(), -position.Z() - box.Z() / 2, t0, t1) &&
		clip(-direction.Z(), position.Z() - box.Z() / 2, t0, t1);

		// The intersection is a point P + t * D. Return closest impact point if positive (i.e. in direction of the motion)
		if(intersect) {
		return std::make_pair(t0, t1);
		}
		return std::nullopt;
		}

		/**
		* @brief Get closest intersection point in positive direction
		* @param direction Direction vector of the motion
		* @param position Original ("before") position to be considered
		* @param box Size of the box to calculate the intersections with
		* @return Closest intersection with box in the direction indicated by input vector
		* or std::nullopt if no intersection of track segment with the box volume can be found in positive
		* direction from the given position.
		*/
		static std::optional<ROOT::Math::XYZPoint> closestIntersection(const ROOT::Math::XYZVector& direction,
		const ROOT::Math::XYZPoint& position,
		const ROOT::Math::XYZVector& box) {

		auto intersect = intersectionDistances(direction, position, box);

		if(!intersect) {
		return std::nullopt;
		}
		// The intersection is a point P + t * D. Return closest impact point if positive (i.e. in direction of the
		// motion)
		auto [t0, t1] = intersect.value();
		if(t0 > 0 && t1 > 0) {
		return (position + std::min(t0, t1) * direction);
		} else if(t0 > 0) {
		@@ -77,11 +106,11 @@ namespace allpix {
		} else if(t1 > 0) {
		return (position + t1 * direction);
		}
		}

		// Otherwise: The line does not intersect the box.
		throw std::invalid_argument("no intersection with volume boundaries found");
		// Otherwise: there is no intersection in positive direction
		return std::nullopt;
		}
		};
		} // namespace allpix

		#endif /* ALLPIX_LIANG_BARSKY_H */