BioDynaMo  v1.03.58-27764645
morton_order.cc
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 
16 #include <morton/morton.h> // NOLINT
17 #include <cmath>
18 
19 namespace bdm {
20 
21 // -----------------------------------------------------------------------------
22 template <typename T>
23 class Stack {
24  public:
25  Stack(uint64_t max_elements) { data_.resize(max_elements); }
26 
27  void Push(const T& el) {
28  data_[++top_] = el;
29  assert(top_ < static_cast<int64_t>(data_.size()));
30  }
31 
32  const T& Pop() {
33  assert(top_ < static_cast<int64_t>(data_.size()));
34  return data_[top_--];
35  }
36 
37  T& Top() { return data_[top_]; }
38  uint64_t Size() const { return static_cast<uint64_t>(top_ + 1); }
39 
40  private:
41  std::vector<T> data_;
42  int64_t top_ = -1;
43 };
44 
45 // -----------------------------------------------------------------------------
46 struct Aux {
47  uint64_t index = 0;
48  uint64_t length = 0;
49  uint64_t xpos = 0;
50  uint64_t ypos = 0;
51  uint64_t zpos = 0;
52 };
53 
54 // -----------------------------------------------------------------------------
55 void MortonOrder::Update(const std::array<uint64_t, 3>& num_boxes_axis) {
56  num_boxes_ = num_boxes_axis[0] * num_boxes_axis[1] * num_boxes_axis[2];
57  std::array<uint64_t, 3> max_dimensions{
58  num_boxes_axis[0] - 1, num_boxes_axis[1] - 1, num_boxes_axis[2] - 1};
59  auto max_dim = std::max(num_boxes_axis[0],
60  std::max(num_boxes_axis[1], num_boxes_axis[2]));
61  auto max_depth =
62  static_cast<uint64_t>(std::ceil(std::log(max_dim) / std::log(2)));
63 
64  offset_index_.clear();
65  offset_index_.reserve(2048);
66  offset_index_.push_back({0, 0});
67 
68  if (max_depth == 0) {
69  return;
70  }
71 
72  auto length = static_cast<uint64_t>(std::pow(2, max_depth - 1));
73 
74  Stack<Aux> s(max_depth);
75  s.Push({0, length, length, length, length});
76 
77  uint64_t box_cnt = 0;
78  uint64_t offset = 0;
79  bool record = false;
80 
81  while (s.Size() != 0) {
82  auto& c = s.Top();
83  // determine borders
84  auto xmax = c.xpos;
85  if (c.index == 1 || c.index == 3 || c.index == 5 || c.index == 7) {
86  xmax += c.length;
87  }
88  auto ymax = c.ypos;
89  if ((c.index > 1 && c.index < 4) || c.index > 5) {
90  ymax += c.length;
91  }
92  auto zmax = c.zpos;
93  if (c.index > 3) {
94  zmax += c.length;
95  }
96  auto xmin = xmax - c.length;
97  auto ymin = ymax - c.length;
98  auto zmin = zmax - c.length;
99 
100  if (max_dimensions[0] >= xmax - 1 && max_dimensions[1] >= ymax - 1 &&
101  max_dimensions[2] >= zmax - 1) {
102  // full
103  if (record) {
104  offset_index_.push_back({box_cnt, offset});
105  record = false;
106  }
107  box_cnt += c.length * c.length * c.length;
108  } else if (max_dimensions[0] < xmin || max_dimensions[1] < ymin ||
109  max_dimensions[2] < zmin) {
110  // empty
111  record = true;
112  auto elements = c.length * c.length * c.length;
113  offset += elements;
114  } else {
115  auto next_length = c.length >> 1;
116  s.Push({0, next_length, xmin + next_length, ymin + next_length,
117  zmin + next_length});
118  continue;
119  }
120 
121  // evaluation was full or empty
122  while (s.Size() > 0 && ++(s.Top().index) == 8) {
123  s.Pop();
124  }
125  }
126 }
127 
128 // -----------------------------------------------------------------------------
129 template <typename T>
130 std::pair<uint64_t, uint64_t> BinarySearch(uint64_t search_val,
131  const T& container, uint64_t from,
132  uint64_t to) {
133  if (container.size() == 0) {
134  return {0, 0};
135  }
136  if (to <= from) {
137  if (container[from].first > search_val && from > 0) {
138  from--;
139  }
140  return {container[from].second, from};
141  }
142 
143  auto m = (from + to) / 2;
144  if (container[m].first == search_val) {
145  return {container[m].second, m};
146  } else if (container[m].first > search_val) {
147  return BinarySearch(search_val, container, from, m);
148  } else {
149  return BinarySearch(search_val, container, m + 1, to);
150  }
151 }
152 
153 // -----------------------------------------------------------------------------
154 uint64_t MortonOrder::GetMortonCode(uint64_t box_index) const {
155  return BinarySearch(box_index, offset_index_, 0, offset_index_.size() - 1)
156  .first +
157  box_index;
158 }
159 
160 // -----------------------------------------------------------------------------
161 class MortonIterator : public Iterator<uint64_t> {
162  public:
163  MortonIterator(uint64_t start_index, uint64_t last_index,
164  uint64_t start_morton_code, uint64_t offset_pos,
165  const std::vector<std::pair<uint64_t, uint64_t>>& offset_index)
166  : next_index_(start_index),
167  last_index_(last_index),
168  next_morton_code_(start_morton_code),
169  offset_pos_(offset_pos),
170  offset_index_(offset_index) {}
171  virtual ~MortonIterator() = default;
172 
173  bool HasNext() const override { return next_index_ <= last_index_; }
174 
175  uint64_t Next() override {
176  if (!HasNext()) {
177  return 0;
178  }
179  auto ret = next_morton_code_;
180  next_index_++;
181  if (HasNext()) {
182  if (offset_pos_ + 1 < offset_index_.size() &&
183  offset_index_[offset_pos_ + 1].first == next_index_) {
184  offset_pos_++;
185  }
187  }
188  return ret;
189  }
190 
191  private:
192  uint64_t next_index_;
193  uint64_t last_index_;
194  uint64_t next_morton_code_ = 0;
195  uint64_t offset_pos_ = 0;
196  const std::vector<std::pair<uint64_t, uint64_t>>& offset_index_;
197 };
198 
199 // -----------------------------------------------------------------------------
201  uint64_t start_index, uint64_t end_index,
202  Functor<void, Iterator<uint64_t>*>& f) const {
203  auto result =
204  BinarySearch(start_index, offset_index_, 0, offset_index_.size() - 1);
205  MortonIterator it(start_index, end_index, start_index + result.first,
206  result.second, offset_index_);
207  f(&it);
208 }
209 
210 } // namespace bdm
bdm::Stack::top_
int64_t top_
Definition: morton_order.cc:42
bdm::Stack
Definition: morton_order.cc:23
bdm::MortonOrder::num_boxes_
uint64_t num_boxes_
Definition: morton_order.h:36
bdm::MortonIterator::offset_pos_
uint64_t offset_pos_
Definition: morton_order.cc:195
bdm::Aux::zpos
uint64_t zpos
Definition: morton_order.cc:51
bdm::Aux
Definition: morton_order.cc:46
bdm
Definition: agent.cc:39
bdm::Iterator
Definition: iterator.h:21
bdm::MortonOrder::offset_index_
std::vector< std::pair< uint64_t, uint64_t > > offset_index_
Definition: morton_order.h:37
bdm::Stack::Top
T & Top()
Definition: morton_order.cc:37
bdm::MortonIterator::offset_index_
const std::vector< std::pair< uint64_t, uint64_t > > & offset_index_
Definition: morton_order.cc:196
bdm::Stack::Size
uint64_t Size() const
Definition: morton_order.cc:38
bdm::MortonOrder::GetMortonCode
uint64_t GetMortonCode(uint64_t box_index) const
Runtime O(log(num_boxes))
Definition: morton_order.cc:154
bdm::MortonIterator
Definition: morton_order.cc:161
bdm::MortonIterator::next_morton_code_
uint64_t next_morton_code_
Definition: morton_order.cc:194
bdm::Stack::Push
void Push(const T &el)
Definition: morton_order.cc:27
bdm::MortonIterator::~MortonIterator
virtual ~MortonIterator()=default
bdm::MortonIterator::last_index_
uint64_t last_index_
Definition: morton_order.cc:193
bdm::Functor
Definition: functor.h:24
bdm::Aux::ypos
uint64_t ypos
Definition: morton_order.cc:50
bdm::MortonOrder::Update
void Update(const std::array< uint64_t, 3 > &num_boxes_axis)
Definition: morton_order.cc:55
bdm::MortonIterator::Next
uint64_t Next() override
Definition: morton_order.cc:175
bdm::BinarySearch
uint64_t BinarySearch(const TSearch &search_val, const TContainer &container, uint64_t from, uint64_t to)
Definition: algorithm.h:76
bdm::Stack::data_
std::vector< T > data_
Definition: morton_order.cc:41
bdm::Aux::xpos
uint64_t xpos
Definition: morton_order.cc:49
bdm::Stack::Stack
Stack(uint64_t max_elements)
Definition: morton_order.cc:25
bdm::Aux::length
uint64_t length
Definition: morton_order.cc:48
bdm::MortonIterator::HasNext
bool HasNext() const override
Definition: morton_order.cc:173
bdm::MortonIterator::next_index_
uint64_t next_index_
Definition: morton_order.cc:192
bdm::MortonIterator::MortonIterator
MortonIterator(uint64_t start_index, uint64_t last_index, uint64_t start_morton_code, uint64_t offset_pos, const std::vector< std::pair< uint64_t, uint64_t >> &offset_index)
Definition: morton_order.cc:163
bdm::Stack::Pop
const T & Pop()
Definition: morton_order.cc:32
morton_order.h
bdm::MortonOrder::CallMortonIteratorConsumer
void CallMortonIteratorConsumer(uint64_t start_index, uint64_t end_index, Functor< void, Iterator< uint64_t > * > &f) const
Definition: morton_order.cc:200
bdm::Aux::index
uint64_t index
Definition: morton_order.cc:47