GCC Code Coverage Report

Directory: gridformat/
File: gridformat/grid/image_grid.hpp
Date: 2024-11-10 16:24:00
Exec Total Coverage
Lines: 99 100 99.0%
Functions: 67 67 100.0%
Branches: 49 94 52.1%
Line Branch Exec Source
1 // SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2 // SPDX-License-Identifier: MIT
3 /*!
4 * \file
5 * \ingroup Grid
6 * \brief Predefined image grid implementation.
7 */
8 #ifndef GRIDFORMAT_GRID_IMAGE_GRID_HPP_
9 #define GRIDFORMAT_GRID_IMAGE_GRID_HPP_
10
11 #include <cmath>
12 #include <array>
13 #include <ranges>
14 #include <cassert>
15 #include <utility>
16 #include <numeric>
17 #include <algorithm>
18 #include <functional>
19 #include <type_traits>
20
21 #include <gridformat/common/ranges.hpp>
22 #include <gridformat/common/concepts.hpp>
23 #include <gridformat/common/exceptions.hpp>
24 #include <gridformat/common/type_traits.hpp>
25 #include <gridformat/common/iterator_facades.hpp>
26 #include <gridformat/common/flat_index_mapper.hpp>
27 #include <gridformat/common/md_index.hpp>
28
29 #include <gridformat/grid/cell_type.hpp>
30 #include <gridformat/grid/traits.hpp>
31
32 namespace GridFormat {
33
34 /*!
35 * \ingroup Grid
36 * \brief Predefined grid implementation that represents a structured, equispaced grid.
37 * \tparam dim The dimension of the grid (1 <= dim <= 3)
38 * \tparam CoordinateType The type used to represent coordinates (e.g. `double`)
39 */
40 template<std::size_t dim, Concepts::Scalar CoordinateType>
41 class ImageGrid {
42 static_assert(dim > 0 && dim <= 3);
43
44 template<int codim>
45 struct Entity { std::array<std::size_t, dim> location; };
46
47 public:
48 using Cell = Entity<0>; //!< The type used for grid cells
49 using Point = Entity<dim>; //!< The type used for grid points
50
51 //! Constructor overload for general range types
52 template<Concepts::StaticallySizedMDRange<1> Size,
53 Concepts::StaticallySizedMDRange<1> Cells>
54 4 ImageGrid(Size&& size, Cells&& cells)
55 : ImageGrid(
56
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 Ranges::filled_array<dim>(CoordinateType{0}),
57 std::forward<Size>(size),
58 std::forward<Cells>(cells)
59
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 8 )
60 4 {}
61
62 //! Constructor overload for general range types
63 template<Concepts::StaticallySizedMDRange<1> Origin,
64 Concepts::StaticallySizedMDRange<1> Size,
65 Concepts::StaticallySizedMDRange<1> Cells>
66 requires(static_size<Origin> == dim and
67 static_size<Size> == dim and
68 static_size<Cells> == dim and
69 std::integral<std::ranges::range_value_t<Cells>>)
70 4 ImageGrid(Origin&& origin, Size&& size, Cells&& cells)
71 : ImageGrid(
72 Ranges::to_array<dim, CoordinateType>(origin),
73 Ranges::to_array<dim, CoordinateType>(size),
74 Ranges::to_array<dim, std::size_t>(cells)
75 4 )
76 4 {}
77
78 /*!
79 * \brief Construct a grid with the given size and discretization.
80 * \param size The size of the grid
81 * \param cells The number of cells in all directions
82 */
83 3 ImageGrid(std::array<CoordinateType, dim> size,
84 std::array<std::size_t, dim> cells)
85
2/4
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 3 : ImageGrid(Ranges::filled_array<dim>(CoordinateType{0}), std::move(size), std::move(cells))
86 3 {}
87
88 /*!
89 * \brief Construct a grid with the given size and discretization.
90 * \param origin The origin of the grid (e.g. the lower-left corner in 2d)
91 * \param size The size of the grid
92 * \param cells The number of cells in all directions
93 */
94 14 ImageGrid(std::array<CoordinateType, dim> origin,
95 std::array<CoordinateType, dim> size,
96 std::array<std::size_t, dim> cells)
97 14 : _lower_right{std::move(origin)}
98
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 14 , _upper_right{Ranges::apply_pairwise<CoordinateType>(std::plus{}, _lower_right, size)}
99 14 , _spacing{_compute_spacing(cells)}
100
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 14 , _cell_index_tuples{MDLayout{cells}}
101
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 14 , _point_index_tuples{MDLayout{Ranges::incremented(cells, 1)}}
102
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 14 , _cell_point_offsets{MDLayout{Ranges::filled_array<dim, std::size_t>(2)}}
103
7/14
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 12 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 12 times.
✗ Branch 16 not taken.
✓ Branch 19 taken 12 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 12 times.
✗ Branch 23 not taken.
 56 , _point_mapper{Ranges::incremented(cells, 1)} {
104
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 12 times.
 39 if (std::ranges::any_of(cells, [] (const std::integral auto i) { return i == 0; }))
105 ✗ throw ValueError("Number of cells in each direction must be > 0");
106 14 }
107
108 1 std::size_t number_of_cells() const { return _cell_index_tuples.size(); }
109 1 std::size_t number_of_points() const { return _point_index_tuples.size(); }
110 2 std::size_t number_of_cells(unsigned int direction) const { return _cell_index_tuples.size(direction); }
111 64 std::size_t number_of_points(unsigned int direction) const { return _point_index_tuples.size(direction); }
112 69 const auto& origin() const { return _lower_right; }
113 115 const auto& spacing() const { return _spacing; }
114
115 //! Return an array containing the number of cells in all directions
116 409 auto extents() const {
117
1/2
✓ Branch 1 taken 381 times.
✗ Branch 2 not taken.
 409 return Ranges::to_array<dim, std::size_t>(
118
3/6
✓ Branch 1 taken 381 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 381 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 381 times.
✗ Branch 8 not taken.
 887 std::views::iota(std::size_t{0}, dim) | std::views::transform([&] (std::size_t i) {
119 775 return _cell_index_tuples.size(i);
120 })
121 818 );
122 }
123
124 //! Return the ordinates of the grid along the given direction
125 60 auto ordinates(int direction) const {
126
2/4
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
 60 std::vector<CoordinateType> result(number_of_points(direction));
127
2/2
✓ Branch 1 taken 678 times.
✓ Branch 2 taken 55 times.
 789 for (unsigned int i = 0; i < result.size(); ++i)
128 729 result[i] = _ordinate_at(i, direction);
129 60 return result;
130 }
131
132 //! Return the physical position of the given grid point
133 181006 auto position(const Point& p) const {
134
1/2
✓ Branch 1 taken 148878 times.
✗ Branch 2 not taken.
 181006 return Ranges::to_array<dim, CoordinateType>(
135
1/2
✓ Branch 1 taken 148878 times.
✗ Branch 2 not taken.
 181006 std::views::iota(std::size_t{0}, dim)
136
2/4
✓ Branch 1 taken 148878 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 148878 times.
✗ Branch 5 not taken.
 635858 | std::views::transform([&] (std::size_t direction) {
137 326340 return _ordinate_at(p.location[direction], direction);
138 })
139 362012 );
140 }
141
142 //! Return the physical center position of the given grid cell
143 22620 auto center(const Cell& c) const {
144 22620 unsigned corner_count = 0;
145 22620 auto result = Ranges::filled_array<dim>(CoordinateType{0});
146
8/14
✓ Branch 1 taken 19620 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19620 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19620 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 88080 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 88080 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 107700 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 88080 times.
✓ Branch 19 taken 19620 times.
 241980 for (const auto& p : points(c)) {
147 109680 corner_count++;
148
1/2
✓ Branch 1 taken 88080 times.
✗ Branch 2 not taken.
 109680 const auto ppos = position(p);
149
2/2
✓ Branch 0 taken 195360 times.
✓ Branch 1 taken 88080 times.
 367440 for (unsigned i = 0; i < dim; ++i)
150 257760 result[i] += ppos[i];
151 }
152
1/2
✓ Branch 1 taken 19620 times.
✗ Branch 2 not taken.
 31020 std::ranges::for_each(result, [&] (auto& v) { v /= static_cast<CoordinateType>(corner_count); });
153 27420 return result;
154 }
155
156 //! Return a unique id for the given point
157 79736 std::size_t id(const Point& p) const {
158 79736 return _point_mapper.map(p.location);
159 }
160
161 //! Return a range over all cells of the grid
162 2764 std::ranges::range auto cells() const {
163
2/4
✓ Branch 1 taken 2603 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2603 times.
✗ Branch 5 not taken.
 33129 return _cell_index_tuples | std::views::transform([] (const auto& indices) {
164 116620 return Cell{Ranges::to_array<dim>(indices)};
165 5528 });
166 }
167
168 //! Return a range over all points of the grid
169 3288 std::ranges::range auto points() const {
170
2/4
✓ Branch 1 taken 3130 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3130 times.
✗ Branch 5 not taken.
 14936 return _point_index_tuples | std::views::transform([] (const auto& indices) {
171 61434 return Point{Ranges::to_array<dim>(indices)};
172 6576 });
173 }
174
175 //! Return a range over all points in the given grid cell
176 82080 std::ranges::range auto points(const Cell& cell) const {
177 82080 return _cell_point_offsets
178
2/4
✓ Branch 1 taken 66840 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 66840 times.
✗ Branch 5 not taken.
 205920 | std::views::transform([loc = cell.location] (const auto& offset) {
179 158880 auto point_loc = loc;
180 158880 std::transform(
181 point_loc.begin(), point_loc.end(), offset.begin(), point_loc.begin(),
182 356160 [] (std::size_t p, std::size_t o) { return p + o; }
183 );
184 158880 return Point{std::move(point_loc)};
185 164160 });
186 }
187
188 private:
189 419909 CoordinateType _ordinate_at(std::size_t i, int direction) const {
190 419909 return _lower_right[direction] + _spacing[direction]*static_cast<CoordinateType>(i);
191 }
192
193 template<Concepts::StaticallySizedMDRange<1> Cells>
194 14 auto _compute_spacing(Cells&& cells) const {
195
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 14 return Ranges::apply_pairwise(
196 std::divides{},
197
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 28 Ranges::apply_pairwise(std::minus{}, _upper_right, _lower_right),
198 cells
199 28 );
200 }
201
202 std::array<unsigned int, dim> _directions_in_descending_size() const {
203 std::array<unsigned int, dim> result;
204 for (int i = 0; i < dim; ++i)
205 result[i] = i;
206 std::ranges::sort(result, [&] (unsigned int dir1, unsigned int dir2) {
207 return _cell_index_tuples.size(dir1) < _cell_index_tuples.size(dir2);
208 });
209 return result;
210 }
211
212 std::array<CoordinateType, dim> _lower_right;
213 std::array<CoordinateType, dim> _upper_right;
214 std::array<CoordinateType, dim> _spacing;
215 MDIndexRange _cell_index_tuples;
216 MDIndexRange _point_index_tuples;
217 MDIndexRange _cell_point_offsets;
218 FlatIndexMapper<dim> _point_mapper;
219 };
220
221 template<Concepts::StaticallySizedMDRange<1> S, Concepts::StaticallySizedMDRange<1> C>
222 requires(static_size<S> == static_size<C>)
223 ImageGrid(S&&, C&&) -> ImageGrid<static_size<S>, std::ranges::range_value_t<S>>;
224
225
226 #ifndef DOXYGEN
227 namespace Traits {
228
229 template<std::size_t dim, typename CT>
230 struct Points<ImageGrid<dim, CT>> {
231 3288 static std::ranges::range auto get(const ImageGrid<dim, CT>& grid) {
232 3288 return grid.points();
233 }
234 };
235
236 template<std::size_t dim, typename CT>
237 struct Cells<ImageGrid<dim, CT>> {
238 2764 static std::ranges::range auto get(const ImageGrid<dim, CT>& grid) {
239 2764 return grid.cells();
240 }
241 };
242
243 template<std::size_t dim, typename CT>
244 struct Origin<ImageGrid<dim, CT>> {
245 69 static const auto& get(const ImageGrid<dim, CT>& grid) {
246 69 return grid.origin();
247 }
248 };
249
250 template<std::size_t dim, typename CT>
251 struct Spacing<ImageGrid<dim, CT>> {
252 111 static const auto& get(const ImageGrid<dim, CT>& grid) {
253 111 return grid.spacing();
254 }
255 };
256
257 template<std::size_t dim, typename CT>
258 struct Extents<ImageGrid<dim, CT>> {
259 405 static auto get(const ImageGrid<dim, CT>& grid) {
260 405 return grid.extents();
261 }
262 };
263
264 template<std::size_t dim, typename CT, typename Entity>
265 struct Location<ImageGrid<dim, CT>, Entity> {
266 static_assert(
267 std::is_same_v<Entity, typename ImageGrid<dim, CT>::Cell> ||
268 std::is_same_v<Entity, typename ImageGrid<dim, CT>::Point>
269 );
270 62156 static auto get(const ImageGrid<dim, CT>&, const Entity& e) {
271 62156 return e.location;
272 }
273 };
274
275 template<std::size_t dim, typename CT>
276 struct Ordinates<ImageGrid<dim, CT>> {
277 60 static auto get(const ImageGrid<dim, CT>& grid, unsigned int dir) {
278 60 return grid.ordinates(dir);
279 }
280 };
281
282 // required for structured grid concept
283 template<std::size_t dim, typename CT>
284 struct PointCoordinates<ImageGrid<dim, CT>, typename ImageGrid<dim, CT>::Point> {
285 43230 static decltype(auto) get(const ImageGrid<dim, CT>& grid, const typename ImageGrid<dim, CT>::Point& p) {
286 43230 return grid.position(p);
287 }
288 };
289
290 // register as unstructured grid as well
291 template<std::size_t dim, typename CT>
292 struct CellPoints<ImageGrid<dim, CT>, typename ImageGrid<dim, CT>::Cell> {
293 59460 static auto get(const ImageGrid<dim, CT>& grid, const typename ImageGrid<dim, CT>::Cell& c) {
294 59460 return grid.points(c);
295 }
296 };
297
298 template<std::size_t dim, typename CT>
299 struct PointId<ImageGrid<dim, CT>, typename ImageGrid<dim, CT>::Point> {
300 79736 static auto get(const ImageGrid<dim, CT>& grid, const typename ImageGrid<dim, CT>::Point& p) {
301 79736 return grid.id(p);
302 }
303 };
304
305 template<std::size_t dim, typename CT>
306 struct CellType<ImageGrid<dim, CT>, typename ImageGrid<dim, CT>::Cell> {
307 78440 static auto get(const ImageGrid<dim, CT>&, const typename ImageGrid<dim, CT>::Cell&) {
308 if constexpr (dim == 1)
309 return GridFormat::CellType::segment;
310 if constexpr (dim == 2)
311 57320 return GridFormat::CellType::pixel;
312 if constexpr (dim == 3)
313 21120 return GridFormat::CellType::voxel;
314 }
315 };
316
317 } // namespace Traits
318 #endif // DOXYGEN
319
320 } // namespace GridFormat
321
322 #endif // GRIDFORMAT_GRID_IMAGE_GRID_HPP_
323