diff --git a/src/kdtree_splitters.hpp b/src/kdtree_splitters.hpp new file mode 100644 index 0000000..4594298 --- /dev/null +++ b/src/kdtree_splitters.hpp @@ -0,0 +1,71 @@ +#ifndef __KDTREE_SPLITTERS_HPP +#define __KDTREE_SPLITTERS_HPP + +#include + +namespace CosmoTool +{ + + template + struct KD_homogeneous_cell_splitter + { + typedef typename KDDef::KDCoordinates coords; + typedef typename KDDef::CoordType ctype; + + void operator()(KDCell **cells, uint32_t Ncells, uint32_t& split_index, int axis, coords minBound, coords maxBound) + { + ctype midCoord = 0.5*(maxBound[axis]+minBound[axis]); + uint32_t below = 0, above = Ncells-1; + ctype delta_max = std::abs(cells[0]->coord[axis]-midCoord); + uint32_t idx_max = 0; + + while (below < above) + { + ctype delta = cells[below]->coord[axis]-midCoord; + if (delta > 0) + { + if (delta < delta_max) + { + delta_max = delta; + idx_max = above; + } + std::swap(cells[below], cells[above--]); + } + else + { + if (-delta < delta_max) + { + delta_max = -delta; + idx_max = below; + } + below++; + } + } + if (idx_max != above) + { + bool cond1 = cells[idx_max]->coord[axis] > midCoord; + bool cond2 = cells[above]->coord[axis] > midCoord; + if ((cond1 && cond2) || (!cond1 && !cond2)) + { + split_index = above; + std::swap(cells[above], cells[idx_max]); + } + else if (cond2) + { + split_index = above-1; + std::swap(cells[above-1], cells[idx_max]); + } + else + { + split_index = above+1; + std::swap(cells[above+1], cells[idx_max]); + } + } + else split_index = above; + } + }; + + +}; + +#endif