-
Notifications
You must be signed in to change notification settings - Fork 35
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Implement surface translation and transformation (#887)
* Add transposed matrix multiply * Add nodiscard attributes to potentially confusing functions * Implement surface translator * Implement surface transformations * Fix typos in documentation * Add note about transpose policy tag * Move surface transformer/lator to detail
- Loading branch information
Showing
12 changed files
with
913 additions
and
11 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,213 @@ | ||
| //----------------------------------*-C++-*----------------------------------// | ||
| // Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. | ||
| // See the top-level COPYRIGHT file for details. | ||
| // SPDX-License-Identifier: (Apache-2.0 OR MIT) | ||
| //---------------------------------------------------------------------------// | ||
| //! \file orange/surf/detail/SurfaceTransformer.cc | ||
| //---------------------------------------------------------------------------// | ||
| #include "SurfaceTransformer.hh" | ||
|
|
||
| #include "corecel/cont/Range.hh" | ||
| #include "corecel/math/ArrayOperators.hh" | ||
| #include "corecel/math/ArrayUtils.hh" | ||
| #include "orange/MatrixUtils.hh" | ||
|
|
||
| #include "../ConeAligned.hh" | ||
| #include "../CylAligned.hh" | ||
| #include "../CylCentered.hh" | ||
| #include "../GeneralQuadric.hh" | ||
| #include "../Plane.hh" | ||
| #include "../PlaneAligned.hh" | ||
| #include "../SimpleQuadric.hh" | ||
| #include "../Sphere.hh" | ||
| #include "../SphereCentered.hh" | ||
|
|
||
| namespace celeritas | ||
| { | ||
| namespace detail | ||
| { | ||
| namespace | ||
| { | ||
| //---------------------------------------------------------------------------// | ||
| #define ORANGE_INSTANTIATE_OP(OUT, IN) \ | ||
| template OUT SurfaceTransformer::operator()(IN<Axis::x> const&) const; \ | ||
| template OUT SurfaceTransformer::operator()(IN<Axis::y> const&) const; \ | ||
| template OUT SurfaceTransformer::operator()(IN<Axis::z> const&) const | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| } // namespace | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform an axis-aligned plane. | ||
| */ | ||
| template<Axis T> | ||
| Plane SurfaceTransformer::operator()(PlaneAligned<T> const& other) const | ||
| { | ||
| return (*this)(Plane{other}); | ||
| } | ||
|
|
||
| ORANGE_INSTANTIATE_OP(Plane, PlaneAligned); | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform a centered cylinder. | ||
| */ | ||
| template<Axis T> | ||
| GeneralQuadric SurfaceTransformer::operator()(CylCentered<T> const& other) const | ||
| { | ||
| return (*this)(CylAligned<T>{other}); | ||
| } | ||
|
|
||
| ORANGE_INSTANTIATE_OP(GeneralQuadric, CylCentered); | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform a sphere. | ||
| */ | ||
| Sphere SurfaceTransformer::operator()(SphereCentered const& other) const | ||
| { | ||
| return (*this)(Sphere{other}); | ||
| } | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform an axis-aligned cylinder. | ||
| */ | ||
| template<Axis T> | ||
| GeneralQuadric SurfaceTransformer::operator()(CylAligned<T> const& other) const | ||
| { | ||
| return (*this)(SimpleQuadric{other}); | ||
| } | ||
|
|
||
| ORANGE_INSTANTIATE_OP(GeneralQuadric, CylAligned); | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform a plane. | ||
| */ | ||
| Plane SurfaceTransformer::operator()(Plane const& other) const | ||
| { | ||
| // Rotate the normal direction | ||
| Real3 normal = tr_.rotate_up(other.normal()); | ||
|
|
||
| // Transform a point on the original plane | ||
| Real3 point = tr_.transform_up(other.displacement() * other.normal()); | ||
|
|
||
| return Plane{normal, point}; | ||
| } | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform a sphere. | ||
| */ | ||
| Sphere SurfaceTransformer::operator()(Sphere const& other) const | ||
| { | ||
| // Transform origin, keep the same radius | ||
| return Sphere::from_radius_sq(tr_.transform_up(other.origin()), | ||
| other.radius_sq()); | ||
| } | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform a cone. | ||
| */ | ||
| template<Axis T> | ||
| GeneralQuadric SurfaceTransformer::operator()(ConeAligned<T> const& other) const | ||
| { | ||
| return (*this)(SimpleQuadric{other}); | ||
| } | ||
|
|
||
| ORANGE_INSTANTIATE_OP(GeneralQuadric, ConeAligned); | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform a simple quadric. | ||
| */ | ||
| GeneralQuadric SurfaceTransformer::operator()(SimpleQuadric const& other) const | ||
| { | ||
| return (*this)(GeneralQuadric{other}); | ||
| } | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Transform a quadric. | ||
| * | ||
| * See celeritas-doc/nb/geometry/quadric-transform.ipynb . The implementation | ||
| * below is less than optimal because we don't need to explicitly construct the | ||
| * full Q matrix. | ||
| * | ||
| * The inverse transform is: \f[ | ||
| \mathbf{R}^{-1}(\mathbf{x}' - \mathbf{t}) \to \mathbf{x} | ||
| \f] | ||
| * or | ||
| \f[ | ||
| \mathbf{\tilde x} | ||
| \gets | ||
| \begin{bmatrix} | ||
| 1 & 0 \\ | ||
| -\mathbf{R}^{-1} \mathbf{t} & \mathbf{R}^{-1} \\ | ||
| \end{bmatrix} | ||
| \mathbf{\tilde x}' | ||
| = \mathbf{\tilde R}^{-1} \mathbf{\tilde x}' | ||
| \f] | ||
| * | ||
| */ | ||
| GeneralQuadric SurfaceTransformer::operator()(GeneralQuadric const& other) const | ||
| { | ||
| using Vec4 = Array<real_type, 4>; | ||
| using Mat4 = SquareMatrix<real_type, 4>; | ||
|
|
||
| // Build inverse transform matrix | ||
| Mat4 const tr_inv = [this] { | ||
| // Reverse and rotate translation | ||
| Real3 trans = tr_.rotate_down(Real3{0, 0, 0} - tr_.translation()); | ||
|
|
||
| // Combine inverted translation with inverse (transpose) of rotation | ||
| Mat4 tr_inv; | ||
| tr_inv[0][0] = 1; | ||
| for (auto i : range(1, 4)) | ||
| { | ||
| tr_inv[i][0] = trans[i - 1]; | ||
| tr_inv[0][i] = 0; | ||
| for (auto j : range(1, 4)) | ||
| { | ||
| tr_inv[i][j] = tr_.rotation()[j - 1][i - 1]; | ||
| } | ||
| } | ||
| return tr_inv; | ||
| }(); | ||
|
|
||
| auto calc_q = [&other] { | ||
| constexpr auto X = to_int(Axis::x); | ||
| constexpr auto Y = to_int(Axis::y); | ||
| constexpr auto Z = to_int(Axis::z); | ||
|
|
||
| Real3 const second = make_array(other.second()); | ||
| Real3 const cross = make_array(other.cross()) / real_type(2); | ||
| Real3 const first = make_array(other.first()) / real_type(2); | ||
| real_type const zeroth = other.zeroth(); | ||
|
|
||
| return Mat4{Vec4{zeroth, first[X], first[Y], first[Z]}, | ||
| Vec4{first[X], second[X], cross[X], cross[Z]}, | ||
| Vec4{first[Y], cross[X], second[Y], cross[Y]}, | ||
| Vec4{first[Z], cross[Z], cross[Y], second[Z]}}; | ||
| }; | ||
|
|
||
| // Apply original quadric matrix to inverse transforom | ||
| auto qrinv = gemm(calc_q(), tr_inv); | ||
|
|
||
| // Apply transpose of inverse (result should be symmetric) | ||
| auto qprime = gemm(matrix::transpose, tr_inv, qrinv); | ||
|
|
||
| // Extract elements back from the matrix | ||
| return GeneralQuadric( | ||
| {qprime[1][1], qprime[2][2], qprime[3][3]}, | ||
| {2 * qprime[1][2], 2 * qprime[2][3], 2 * qprime[1][3]}, | ||
| {2 * qprime[0][1], 2 * qprime[0][2], 2 * qprime[0][3]}, | ||
| qprime[0][0]); | ||
| } | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| } // namespace detail | ||
| } // namespace celeritas |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,62 @@ | ||
| //----------------------------------*-C++-*----------------------------------// | ||
| // Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. | ||
| // See the top-level COPYRIGHT file for details. | ||
| // SPDX-License-Identifier: (Apache-2.0 OR MIT) | ||
| //---------------------------------------------------------------------------// | ||
| //! \file orange/surf/detail/SurfaceTransformer.hh | ||
| //---------------------------------------------------------------------------// | ||
| #pragma once | ||
|
|
||
| #include "orange/Types.hh" | ||
| #include "orange/transform/Transformation.hh" | ||
|
|
||
| #include "../SurfaceFwd.hh" | ||
|
|
||
| namespace celeritas | ||
| { | ||
| namespace detail | ||
| { | ||
| //---------------------------------------------------------------------------// | ||
| /*! | ||
| * Apply a transformation to a surface to get another surface. | ||
| * | ||
| * The transform gives the new origin and rotation for the surface: rotation is | ||
| * applied first, then translation. | ||
| */ | ||
| class SurfaceTransformer | ||
| { | ||
| public: | ||
| //! Construct with the transformation to apply | ||
| explicit SurfaceTransformer(Transformation const& trans) : tr_{trans} {} | ||
|
|
||
| //// SURFACE FUNCTIONS //// | ||
|
|
||
| template<Axis T> | ||
| Plane operator()(PlaneAligned<T> const&) const; | ||
|
|
||
| template<Axis T> | ||
| GeneralQuadric operator()(CylCentered<T> const&) const; | ||
|
|
||
| Sphere operator()(SphereCentered const&) const; | ||
|
|
||
| template<Axis T> | ||
| GeneralQuadric operator()(CylAligned<T> const&) const; | ||
|
|
||
| Plane operator()(Plane const&) const; | ||
|
|
||
| Sphere operator()(Sphere const&) const; | ||
|
|
||
| template<Axis T> | ||
| GeneralQuadric operator()(ConeAligned<T> const&) const; | ||
|
|
||
| GeneralQuadric operator()(SimpleQuadric const&) const; | ||
|
|
||
| GeneralQuadric operator()(GeneralQuadric const&) const; | ||
|
|
||
| private: | ||
| Transformation tr_; | ||
| }; | ||
|
|
||
| //---------------------------------------------------------------------------// | ||
| } // namespace detail | ||
| } // namespace celeritas |
Oops, something went wrong.