SAH kd tree construction CUDA

#pragma once

/* SAH kd-tree construction algorithm implementation
 *
 * Copyright (c) 2019-2020, Anatoliy V. Tomilov
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following condition is met:
 * Redistributions of source code must retain the above copyright notice, this condition and the following disclaimer.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include "sah_kd_tree.h"

#include <thrust/tuple.h>
#include <thrust/functional.h>
#include <thrust/advance.h>
#include <thrust/distance.h>
#include <thrust/sequence.h>
#include <thrust/fill.h>
#include <thrust/copy.h>
#include <thrust/partition.h>
#include <thrust/remove.h>
#include <thrust/gather.h>
#include <thrust/scatter.h>
#include <thrust/scan.h>
#include <thrust/transform_scan.h>
#include <thrust/transform.h>
#include <thrust/transform_reduce.h>
#include <thrust/reduce.h>
#include <thrust/extrema.h>
#include <thrust/count.h>
#include <thrust/sort.h>
#include <thrust/unique.h>

#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/discard_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/reverse_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/transform_output_iterator.h>
#include <thrust/iterator/zip_iterator.h>

#include <type_traits>
#include <utility>
#include <limits>
#include <iterator>

#include <cassert>

//#ifdef DEBUG
#include <thrust/equal.h>
#include <thrust/logical.h>

#include <utility>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <chrono>
#include <string>

struct timepoint
{
    std::string prefix = "time delta";
    std::chrono::time_point< std::chrono::system_clock > start = std::chrono::system_clock::now();

    void reset()
    {
        start = std::chrono::system_clock::now();
    }

    void print(std::string what)
    {
        auto stop = std::chrono::system_clock::now();
        std::cout << prefix << ' ' << what << ' ' << std::chrono::duration_cast< std::chrono::milliseconds >(stop - std::exchange(start, stop)).count() << '\n';
    }
};
//#endif

#ifndef __CUDACC_EXTENDED_LAMBDA__
#error "nvcc --expt-extended-lambda"
#endif
#ifndef __CUDACC_RELAXED_CONSTEXPR__
#error "nvcc --expt-relaxed-constexpr"
#endif

namespace sah_kd_tree
{

template< typename T >
struct add_leaf_const_reference
{
    using type = const T &; // I, U, F
};

template< typename T >
using add_leaf_const_reference_t = typename add_leaf_const_reference< T >::type;

template<>
struct add_leaf_const_reference< thrust::null_type >
{
    using type = thrust::null_type;
};

template< typename ...Ts >
struct add_leaf_const_reference< thrust::tuple< Ts... > >
{
    using type = thrust::tuple< add_leaf_const_reference_t< Ts >... >;
};

template< typename iterator >
using input_iterator_value_type = add_leaf_const_reference_t< typename iterator::value_type >;

template< typename iterator >
using output_iterator_value_type = typename iterator::value_type;

template< typename policy, template< typename ... > class container >
struct builder
{

    policy p;

    template< I dimension >
    struct projection
    {

        policy p;

        struct
        {
            container< F > a, b, c;
        } triangle;

        struct
        {
            container< F > split_cost;
            container< U > split_event;
            container< F > split_pos;
            struct
            {
                container< U > l, r;
            } polygon_count;
        } layer;

        struct
        {
            container< F > min, max;
        } node;

        struct
        {
            container< F > min, max;
        } polygon;

        struct
        {
            container< U > node;
            container< F > pos;
            container< I > kind;
            container< U > polygon;

            container< U > l, r;
        } event;

        template< typename a_iterator, typename b_iterator, typename c_iterator >
        void set_triangle(a_iterator a_begin, a_iterator a_end,
                          b_iterator b_begin, b_iterator b_end,
                          c_iterator c_begin, c_iterator c_end)
        {
            auto triangle_count = U(thrust::distance(a_begin, a_end));
            assert(triangle_count == U(thrust::distance(b_begin, b_end)));
            assert(triangle_count == U(thrust::distance(c_begin, c_end)));

            triangle.a.assign(a_begin, a_end);
            triangle.b.assign(b_begin, b_end);
            triangle.c.assign(c_begin, c_end);
        }

        void calculate_triangle_bbox()
        {
            auto triangle_count = U(triangle.a.size());
            assert(triangle_count == U(triangle.b.size()));
            assert(triangle_count == U(triangle.c.size()));

            std::cout << "triangle_count = " << triangle_count << "\n";

            polygon.min.resize(triangle_count);
            polygon.max.resize(triangle_count);

            auto triangle_begin = thrust::make_zip_iterator(thrust::make_tuple(triangle.a.cbegin(), triangle.b.cbegin(), triangle.c.cbegin()));
            using triangle_type = input_iterator_value_type< decltype(triangle_begin) >;
            auto polygon_bbox_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon.min.begin(), polygon.max.begin()));
            using polygon_bbox_type = output_iterator_value_type< decltype(polygon_bbox_begin) >;
            auto to_triangle_bbox = [] __host__ __device__ (triangle_type triangle) -> polygon_bbox_type
            {
                F a = thrust::get< 0 >(triangle);
                F b = thrust::get< 1 >(triangle);
                F c = thrust::get< 2 >(triangle);
                return {thrust::min(a, thrust::min(b, c)), thrust::max(a, thrust::max(b, c))};
            };
            thrust::transform(p, triangle_begin, thrust::next(triangle_begin, triangle_count), polygon_bbox_begin, to_triangle_bbox);
        }

        void caluculate_root_node_bbox()
        {
            auto root_bbox_min_begin = thrust::min_element(p, polygon.min.cbegin(), polygon.min.cend());
            node.min.assign(root_bbox_min_begin, thrust::next(root_bbox_min_begin));

            auto root_bbox_max_begin = thrust::max_element(p, polygon.max.cbegin(), polygon.max.cend());
            node.max.assign(root_bbox_max_begin, thrust::next(root_bbox_max_begin));
        }

        template< typename type >
        struct doubler
        {
            __host__ __device__ auto operator () (type value) -> thrust::tuple< type, type >
            {
                return {value, value};
            }
        };

        void generate_initial_event()
        {
            auto triangle_count = U(triangle.a.size());

            auto triangle_bbox_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon.min.cbegin(), polygon.max.cbegin()));
            using bbox_type = input_iterator_value_type< decltype(triangle_bbox_begin) >;
            auto is_planar_event = [] __host__ __device__ (bbox_type bbox) -> bool { return !(thrust::get< 0 >(bbox) < thrust::get< 1 >(bbox)); };
            auto planar_event_count = U(thrust::count_if(p, triangle_bbox_begin, thrust::next(triangle_bbox_begin, triangle_count), is_planar_event));

            std::cout << "planar_event_count = " << planar_event_count << "\n";

            auto event_count = triangle_count - planar_event_count + triangle_count;

            event.pos.resize(event_count);
            event.kind.resize(event_count, I(0));
            event.polygon.resize(event_count);

            auto event_kind_lr_begin = thrust::make_zip_iterator(thrust::make_tuple(event.kind.begin(), event.kind.rbegin()));
            auto planar_event_kind = thrust::fill_n(p, event_kind_lr_begin, triangle_count - planar_event_count, thrust::make_tuple< I, I >(+1, -1)); // right event sequenced before left event if positions are equal
            //thrust::fill_n(p, thrust::get< 0 >(planar_event_kind.get_iterator_tuple()), planar_event_count, I(0));

            auto triangle_begin = thrust::make_counting_iterator< U >(0);
            auto planar_event_begin = thrust::next(event.polygon.begin(), triangle_count - planar_event_count);
            auto event_pair_begin = thrust::make_zip_iterator(thrust::make_tuple(event.polygon.begin(), event.polygon.rbegin()));
            auto solid_event_begin = thrust::make_transform_output_iterator(event_pair_begin, doubler< U >{});
            thrust::partition_copy(p, triangle_begin, thrust::next(triangle_begin, triangle_count), triangle_bbox_begin, planar_event_begin, solid_event_begin, is_planar_event);

            auto event_polygon_bbox_begin = thrust::make_permutation_iterator(triangle_bbox_begin, event.polygon.cbegin());
            auto to_event_pos = [] __host__ __device__ (I event_kind, bbox_type bbox) -> F
            {
                return (event_kind < 0) ? thrust::get< 1 >(bbox) : thrust::get< 0 >(bbox);
            };
            thrust::transform(p, event.kind.cbegin(), event.kind.cend(), event_polygon_bbox_begin, event.pos.begin(), to_event_pos);

            auto event_begin = thrust::make_zip_iterator(thrust::make_tuple(event.pos.begin(), event.kind.begin(), event.polygon.begin()));
            thrust::sort(p, event_begin, thrust::next(event_begin, event_count));

            event.node.resize(event_count, U(0));
        }

        using Y = projection< (dimension + 1) % 3 >;
        using Z = projection< (dimension + 2) % 3 >;

        void find_perfect_split(const params & sah, U node_count, const container< U > & layer_node_offset, const Y & y, const Z & z)
        {
            auto event_count = U(event.kind.size());

            { // calculate l and r polygon count
                event.l.resize(event_count);
                event.r.resize(event_count);

                auto l_triangle_count_begin = thrust::make_transform_iterator(event.kind.cbegin(), [] __host__ __device__ (I event_kind) -> U { return (event_kind < 0) ? 0 : 1; });
                thrust::exclusive_scan_by_key(p, event.node.cbegin(), event.node.cend(), l_triangle_count_begin, event.l.begin());

                auto r_triangle_count_begin = thrust::make_transform_iterator(event.kind.crbegin(), [] __host__ __device__ (I event_kind) -> U { return (0 < event_kind) ? 0 : 1; });
                thrust::exclusive_scan_by_key(p, event.node.crbegin(), event.node.crend(), r_triangle_count_begin, event.r.rbegin());
            }

            layer.split_cost.resize(node_count);
            layer.split_event.resize(node_count);
            layer.split_pos.resize(node_count);

            layer.polygon_count.l.resize(node_count);
            layer.polygon_count.r.resize(node_count);

            auto node_limits_begin = thrust::make_zip_iterator(thrust::make_tuple(node.min.cbegin(), node.max.cbegin(), y.node.min.cbegin(), y.node.max.cbegin(), z.node.min.cbegin(), z.node.max.cbegin()));
            auto node_bbox_begin = thrust::make_permutation_iterator(node_limits_begin, event.node.cbegin());
            auto split_event_begin = thrust::make_counting_iterator< U >(0);
            auto perfect_split_input_begin = thrust::make_zip_iterator(thrust::make_tuple(node_bbox_begin, event.pos.cbegin(), event.kind.cbegin(), split_event_begin, event.l.cbegin(), event.r.cbegin()));
            using perfect_split_input_type = input_iterator_value_type< decltype(perfect_split_input_begin) >;
            auto perfect_split_begin = thrust::make_zip_iterator(thrust::make_tuple(layer.split_cost.begin(), layer.split_event.begin(), layer.split_pos.begin(), layer.polygon_count.l.begin(), layer.polygon_count.r.begin()));
            using perfect_split_type = output_iterator_value_type< decltype(perfect_split_begin) >;
            auto to_split = [sah] __host__ __device__ (perfect_split_input_type perfect_split_value) -> perfect_split_type
            {
                const auto & node_bbox = thrust::get< 0 >(perfect_split_value);
                F split_pos = thrust::get< 1 >(perfect_split_value);
                I event_kind = thrust::get< 2 >(perfect_split_value);
                U split_event = thrust::get< 3 >(perfect_split_value);
                U l = thrust::get< 4 >(perfect_split_value), r = thrust::get< 5 >(perfect_split_value);
                F min = thrust::get< 0 >(node_bbox), max = thrust::get< 1 >(node_bbox);
                F L = split_pos - min;
                F R = max - split_pos;
                if (event_kind < 0) {
                    ++split_event;
                } else if (event_kind == 0) {
                    if (L < R) {
                        ++l;
                        ++split_event;
                    } else {
                        ++r;
                    }
                }
                F x = max - min;
                F y = thrust::get< 3 >(node_bbox) - thrust::get< 2 >(node_bbox);
                F z = thrust::get< 5 >(node_bbox) - thrust::get< 4 >(node_bbox);
                F perimeter = y + z;
                F square = y * z;
                F split_cost = (l * (square + perimeter * L) + r * (square + perimeter * R)) / (square + perimeter * x);
                split_cost *= sah.intersection_cost;
                split_cost += sah.traversal_cost;
                if ((l == 0) || (r == 0)) {
                    split_cost *= sah.emptiness_factor;
                }
                return {split_cost, split_event, split_pos, l, r};
            };
            auto perfect_split_value_begin = thrust::make_transform_iterator(perfect_split_input_begin, to_split);
            auto ends = thrust::reduce_by_key(p, event.node.cbegin(), event.node.cend(), perfect_split_value_begin, thrust::make_discard_iterator(), thrust::make_permutation_iterator(perfect_split_begin, layer_node_offset.cbegin()), thrust::equal_to< U >{}, thrust::minimum< perfect_split_type >{});
            assert(ends.first == thrust::make_discard_iterator(layer_node_offset.size()));
        }

        void calculate_polygon_side(const container< I > & node_split_dimension,
                                    U base_node,
                                    container< U > & polygon_event_l, container< U > & polygon_event_r,
                                    container< I > & polygon_side)
        {
            auto event_count = U(event.kind.size());

            auto split_dimension_begin = thrust::make_permutation_iterator(node_split_dimension.cbegin(), event.node.cbegin());

            { // find event counterpart
                auto event_begin = thrust::make_counting_iterator< U >(0);
                auto event_end = thrust::next(event_begin, event_count);

                auto event_side_begin = thrust::make_zip_iterator(thrust::make_tuple(split_dimension_begin, event.kind.cbegin()));
                using event_side_type = input_iterator_value_type< decltype(event_side_begin) >;

                auto is_not_r_event = [] __host__ __device__ (event_side_type event_side) -> bool
                {
                    U node_split_dimension = thrust::get< 0 >(event_side);
                    if (node_split_dimension != dimension) {
                        return false;
                    }
                    I event_kind = thrust::get< 1 >(event_side);
                    return !(event_kind < 0);
                };
                thrust::scatter_if(p, event_begin, event_end, event.polygon.cbegin(), event_side_begin, polygon_event_l.begin(), is_not_r_event);

                // TODO: optimize out one of (polygon_event_l, polygon_event_r)
                auto is_not_l_event = [] __host__ __device__ (event_side_type event_side) -> bool
                {
                    U node_split_dimension = thrust::get< 0 >(event_side);
                    if (node_split_dimension != dimension) {
                        return false;
                    }
                    I event_kind = thrust::get< 1 >(event_side);
                    return !(0 < event_kind);
                };
                thrust::scatter_if(p, event_begin, event_end, event.polygon.cbegin(), event_side_begin, polygon_event_r.begin(), is_not_l_event);
            }

            //std::cout << "dimension " << dimension << ": event count " << event_count << std::endl;

            auto split_event_begin = thrust::make_permutation_iterator(thrust::prev(layer.split_event.cbegin(), base_node), event.node.cbegin());
            auto polygon_event_lr_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon_event_l.cbegin(), polygon_event_r.cbegin()));
            auto event_lr_begin = thrust::make_permutation_iterator(polygon_event_lr_begin, event.polygon.cbegin());
            using event_lr_type = input_iterator_value_type< decltype(event_lr_begin) >;
            auto polygon_side_begin = thrust::make_permutation_iterator(polygon_side.begin(), event.polygon.cbegin());
            auto to_polygon_side = [] __host__ __device__ (event_lr_type event_lr, U split_event) -> I
            {
                U l = thrust::get< 0 >(event_lr);
                U r = thrust::get< 1 >(event_lr);
                assert(l != U(~0));
                assert(r != U(~0));
                if (r < l) {
                    printf("event1: %u %u %u\n", l, r, split_event);
                }
                assert(!(r < l));
                if (r < split_event) {
                    return -1; // goes to left node
                } else if (l < split_event) {
                    return  0; // goes to both left child node and right child node (splitted), assert(event_kind != 0)
                } else {
                    return +1; // goes to right node
                }
            };
            auto is_x = [] __host__ __device__ (I node_split_dimension) -> bool { return node_split_dimension == dimension; };
            thrust::transform_if(p, event_lr_begin, thrust::next(event_lr_begin, event_count), split_event_begin, split_dimension_begin, polygon_side_begin, to_polygon_side, is_x);
        }

        void decouple_lr_event(const container< I > & node_split_dimension, const container< I > & polygon_side)
        {
            auto event_count = U(event.kind.size());

            auto event_begin = thrust::make_counting_iterator< U >(0);

            auto split_dimension_begin = thrust::make_permutation_iterator(node_split_dimension.cbegin(), event.node.cbegin());
            auto polygon_side_begin = thrust::make_permutation_iterator(polygon_side.cbegin(), event.polygon.cbegin());
            auto polygon_stencil_begin = thrust::make_zip_iterator(thrust::make_tuple(split_dimension_begin, polygon_side_begin));
            using polygon_stencil_type = input_iterator_value_type< decltype(polygon_stencil_begin) >;

            auto is_l_polygon = [] __host__ __device__ (polygon_stencil_type polygon_stencil) -> bool
            {
                I split_dimension = thrust::get< 0 >(polygon_stencil);
                if (split_dimension < 0) {
                    return false;
                }
                I polygon_side = thrust::get< 1 >(polygon_stencil);
                assert(polygon_side != 101325);
                return polygon_side < 0;
            };
            auto event_l_end = thrust::copy_if(p, event_begin, thrust::next(event_begin, event_count), polygon_stencil_begin, event.l.begin(), is_l_polygon);
            event.l.erase(event_l_end, event.l.end());

            auto is_r_polygon = [] __host__ __device__ (polygon_stencil_type polygon_stencil) -> bool
            {
                I split_dimension = thrust::get< 0 >(polygon_stencil);
                if (split_dimension < 0) {
                    return false;
                }
                I polygon_side = thrust::get< 1 >(polygon_stencil);
                assert(polygon_side != 101325);
                return 0 < polygon_side;
            };
            auto event_r_end = thrust::copy_if(p, event_begin, thrust::next(event_begin, event_count), polygon_stencil_begin, event.r.begin(), is_r_polygon);
            event.r.erase(event_r_end, event.r.end());
        }

        void split_polygon(const container< I > & node_split_dimension, const container< F > & node_split_pos,
                           const container< U > & polygon_triangle, const container< U > & polygon_node,
                           U polygon_count, U splitted_polygon_count,
                           const container< U > & splitted_polygon,
                           const Y & y, const Z & z)
        {
            // node of right part of splitted polygon (starting from polygon_count) is still node from previous layer

            polygon.min.resize(polygon_count + splitted_polygon_count);
            polygon.max.resize(polygon_count + splitted_polygon_count);

            auto polygon_bbox_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon.min.begin(), polygon.max.begin()));
            auto polygon_l_bbox_begin = thrust::make_permutation_iterator(polygon_bbox_begin, splitted_polygon.cbegin());
            using polygon_bbox_input_type = input_iterator_value_type< decltype(polygon_l_bbox_begin) >;
            auto node_begin = thrust::make_zip_iterator(thrust::make_tuple(node_split_dimension.cbegin(), node_split_pos.cbegin()));
            auto polygon_node_begin = thrust::make_permutation_iterator(node_begin, polygon_node.cbegin());
            auto triangle_begin = thrust::make_zip_iterator(thrust::make_tuple(triangle.a.cbegin(), triangle.b.cbegin(), triangle.c.cbegin(), y.triangle.a.cbegin(), y.triangle.b.cbegin(), y.triangle.c.cbegin(), z.triangle.a.cbegin(), z.triangle.b.cbegin(), z.triangle.c.cbegin()));
            auto polygon_triangle_begin = thrust::make_permutation_iterator(triangle_begin, polygon_triangle.cbegin());
            auto polygon_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon_node_begin, polygon_triangle_begin));
            using polygon_type = input_iterator_value_type< decltype(polygon_begin) >;
            auto polygon_r_bbox_begin = thrust::next(polygon_bbox_begin, polygon_count);
            auto splitted_polygon_bbox_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon_l_bbox_begin, polygon_r_bbox_begin));
            using splitted_polygon_bbox_type = output_iterator_value_type< decltype(splitted_polygon_bbox_begin) >;
            auto to_splitted_polygon = [] __host__ __device__ (polygon_bbox_input_type bbox, polygon_type polygon) -> splitted_polygon_bbox_type
            {
                F min = thrust::get< 0 >(bbox), max = thrust::get< 1 >(bbox);
                assert(!(max < min));
                const auto & polygon_node = thrust::get< 0 >(polygon);
                I node_split_dimension = thrust::get< 0 >(polygon_node);
                F node_split_pos = thrust::get< 1 >(polygon_node);
                if (node_split_dimension == dimension) {
                    assert(!(node_split_pos < min) && !(max < node_split_pos));
                    return {{min, node_split_pos}, {node_split_pos, max}};
                } else if (!(min < max)) {
                    return {bbox, bbox};
                }

                const auto & triangle = thrust::get< 1 >(polygon);
                F a, b, c;
                if (node_split_dimension == ((dimension + 1) % 3)) {
                    a = thrust::get< 3 >(triangle);
                    b = thrust::get< 4 >(triangle);
                    c = thrust::get< 5 >(triangle);
                } else {
                    assert(node_split_dimension == ((dimension + 2) % 3));
                    a = thrust::get< 6 >(triangle);
                    b = thrust::get< 7 >(triangle);
                    c = thrust::get< 8 >(triangle);
                }

                bool a_side = (a < node_split_pos);
                bool b_side = (b < node_split_pos);
                bool c_side = (c < node_split_pos);

                F ax = thrust::get< 0 >(triangle);
                F bx = thrust::get< 1 >(triangle);
                F cx = thrust::get< 2 >(triangle);

                if (a_side == b_side) {
                    assert(a_side != c_side);
                    thrust::swap(a_side, c_side);
                    thrust::swap(ax, cx);
                    thrust::swap(a, c);
                } else if (a_side == c_side) {
                    assert(a_side != b_side);
                    thrust::swap(a_side, b_side);
                    thrust::swap(ax, bx);
                    thrust::swap(a, b);
                }
                assert(b_side == c_side);

                if (a_side != ((bx - ax) * (c - a) < (b - a) * (cx - ax))) {
                    thrust::swap(bx, cx);
                    thrust::swap(b, c);
                }

                assert((b < a) || (a < b));
                F lpos = ax + (bx - ax) * (node_split_pos - a) / (b - a);
                assert((c < a) || (a < c));
                F rpos = ax + (cx - ax) * (node_split_pos - a) / (c - a);

                if (std::is_floating_point_v< F >) {
                    if (max < lpos) {
                        lpos = max;
                    }
                    if (rpos < min) {
                        rpos = min;
                    }

                    if (rpos < lpos) {
                        rpos = lpos = (lpos + rpos) / F(2);
                    }
                }

                assert(!(rpos < lpos));

                F lmin, rmin;
                if (min < lpos) {
                    if ((ax < bx) == (a < b)) {
                        lmin = min;
                        rmin = lpos;
                    } else {
                        lmin = lpos;
                        rmin = min;
                    }
                } else {
                    lmin = min;
                    rmin = min;
                }
                F lmax, rmax;
                if (rpos < max) {
                    if ((ax < cx) == (a < c)) {
                        lmax = rpos;
                        rmax = max;
                    } else {
                        lmax = max;
                        rmax = rpos;
                    }
                } else {
                    lmax = max;
                    rmax = max;
                }
                assert(!(lmax < lmin));
                assert(!(rmax < rmin));
                return {{lmin, lmax}, {rmin, rmax}};
            };
            thrust::transform(p, polygon_l_bbox_begin, thrust::next(polygon_l_bbox_begin, splitted_polygon_count), thrust::next(polygon_begin, polygon_count), splitted_polygon_bbox_begin, to_splitted_polygon);
        }

        void merge_event(U polygon_count, const container< U > & polygon_node, U splitted_polygon_count, const container< U > & splitted_polygon)
        {
            auto event_count = U(event.kind.size());

            auto polygon_bbox_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon.min.cbegin(), polygon.max.cbegin()));
            using bbox_type = input_iterator_value_type< decltype(polygon_bbox_begin) >;
            auto is_planar_polygon = [] __host__ __device__ (bbox_type bbox) -> bool { return !(thrust::get< 0 >(bbox) < thrust::get< 1 >(bbox)); };

            auto polygon_l_bbox_begin = thrust::make_permutation_iterator(polygon_bbox_begin, splitted_polygon.cbegin());
            auto polygon_r_bbox_begin = thrust::next(polygon_bbox_begin, polygon_count);

            auto planar_polygon_l_count = U(thrust::count_if(p, polygon_l_bbox_begin, thrust::next(polygon_l_bbox_begin, splitted_polygon_count), is_planar_polygon));
            auto planar_polygon_r_count = U(thrust::count_if(p, polygon_r_bbox_begin, thrust::next(polygon_r_bbox_begin, splitted_polygon_count), is_planar_polygon));

            U splitted_event_l_count = splitted_polygon_count + splitted_polygon_count - planar_polygon_l_count;
            U splitted_event_r_count = splitted_polygon_count + splitted_polygon_count - planar_polygon_r_count;

            //std::cout << "splitted_event_l_count = " << splitted_event_l_count << "\n";
            //std::cout << "splitted_event_r_count = " << splitted_event_r_count << "\n";

            U splitted_event_count = splitted_event_l_count + splitted_event_r_count;

            auto l = U(event.l.size());
            auto r = U(event.r.size());

            U event_storage_size = std::exchange(event_count, l + r + splitted_event_count);
            if (event_storage_size < event_count) {
                event_storage_size = event_count;
            }

            // grant additional storage for all merge operations
            event.node.resize(event_storage_size + event_count, 101325u);
            event.pos.resize(event_storage_size + event_count);
            event.kind.resize(event_storage_size + event_count, I(0));
            event.polygon.resize(event_storage_size + event_count);

            // merge l and r event
            auto event_node_begin = thrust::make_permutation_iterator(polygon_node.cbegin(), event.polygon.cbegin());
            auto event_value_begin = thrust::make_zip_iterator(thrust::make_tuple(event.pos.begin(), event.kind.begin(), event.polygon.begin()));

            auto event_key_l_begin = thrust::make_permutation_iterator(event_node_begin, event.l.cbegin());
            auto event_value_l_begin = thrust::make_permutation_iterator(event_value_begin, event.l.cbegin());

            auto event_key_r_begin = thrust::make_permutation_iterator(event_node_begin, event.r.cbegin());
            auto event_value_r_begin = thrust::make_permutation_iterator(event_value_begin, event.r.cbegin());

            if ((false)) {
                auto check_node = [] __host__ __device__ (U node)
                {
                    if (node == 0 || node == ~0u) {
                        printf("node node %u\n", node);
                    }
                };
                thrust::for_each_n(p, event_key_l_begin, l, check_node);
                thrust::for_each_n(p, event_key_r_begin, r, check_node);
            }

            auto event_lr_key_begin = thrust::next(event.node.begin(), event_storage_size);
            auto event_lr_value_begin = thrust::next(event_value_begin, event_storage_size);
            thrust::merge_by_key(p, event_key_l_begin, thrust::next(event_key_l_begin, l), event_key_r_begin, thrust::next(event_key_r_begin, r), event_value_l_begin, event_value_r_begin, event_lr_key_begin, event_lr_value_begin);

            auto splitted_event_offset = event_storage_size + l + r;

            // calculate event kind for event of splitted polygon
            auto event_l_kind_l_begin = thrust::next(event.kind.begin(), splitted_event_offset);
            auto event_r_kind_l_begin = thrust::next(event_l_kind_l_begin, splitted_event_l_count);
            auto event_l_kind_r_begin = thrust::make_reverse_iterator(event_r_kind_l_begin);
            auto event_r_kind_r_begin = thrust::prev(event_l_kind_r_begin, splitted_event_r_count);

            auto event_l_kind_lr_begin = thrust::make_zip_iterator(thrust::make_tuple(event_l_kind_l_begin, event_l_kind_r_begin));
            auto event_r_kind_lr_begin = thrust::make_zip_iterator(thrust::make_tuple(event_r_kind_l_begin, event_r_kind_r_begin));

            auto kind_lr = thrust::make_tuple(+1, -1);
            thrust::fill_n(p, event_l_kind_lr_begin, splitted_polygon_count - planar_polygon_l_count, kind_lr);
            thrust::fill_n(p, event_r_kind_lr_begin, splitted_polygon_count - planar_polygon_r_count, kind_lr);

            // calculate l and r polygon for event of splitted polygon
            auto event_l_polygon_l_begin = thrust::next(event.polygon.begin(), splitted_event_offset);
            auto event_r_polygon_l_begin = thrust::next(event_l_polygon_l_begin, splitted_event_l_count);
            auto event_l_polygon_r_begin = thrust::make_reverse_iterator(event_r_polygon_l_begin);
            auto event_r_polygon_r_begin = thrust::prev(event_l_polygon_r_begin, splitted_event_r_count);

            auto event_l_polygon_lr_begin = thrust::make_zip_iterator(thrust::make_tuple(event_l_polygon_l_begin, event_l_polygon_r_begin));
            auto event_r_polygon_lr_begin = thrust::make_zip_iterator(thrust::make_tuple(event_r_polygon_l_begin, event_r_polygon_r_begin));

            auto event_l_polygon_lr_output_begin = thrust::make_transform_output_iterator(event_l_polygon_lr_begin, doubler< U >{});
            auto event_r_polygon_lr_output_begin = thrust::make_transform_output_iterator(event_r_polygon_lr_begin, doubler< U >{});

            auto event_l_polygon_planar_begin = thrust::next(event_l_polygon_l_begin, splitted_polygon_count - planar_polygon_l_count);
            auto event_r_polygon_planar_begin = thrust::next(event_r_polygon_l_begin, splitted_polygon_count - planar_polygon_r_count);

            auto polygon_l_begin = splitted_polygon.cbegin();
            auto polygon_r_begin = thrust::make_counting_iterator< U >(polygon_count);

            auto ends_l = thrust::partition_copy(p, polygon_l_begin, thrust::next(polygon_l_begin, splitted_polygon_count), polygon_l_bbox_begin, event_l_polygon_planar_begin, event_l_polygon_lr_output_begin, is_planar_polygon);
            assert(thrust::next(event_l_polygon_planar_begin, planar_polygon_l_count) == ends_l.first);
            auto ends_r = thrust::partition_copy(p, polygon_r_begin, thrust::next(polygon_r_begin, splitted_polygon_count), polygon_r_bbox_begin, event_r_polygon_planar_begin, event_r_polygon_lr_output_begin, is_planar_polygon);
            assert(thrust::next(event_r_polygon_planar_begin, planar_polygon_r_count) == ends_r.first);

            // calculate event pos
#if 1
            auto event_polygon_begin = event_l_polygon_l_begin;
            auto event_polygon_end = thrust::next(event_l_polygon_l_begin, splitted_polygon_count);
            auto event_pos_begin = thrust::next(event.pos.begin(), splitted_event_offset);
            event_pos_begin = thrust::gather(p, std::exchange(event_polygon_begin, event_polygon_end), event_polygon_end, polygon.min.cbegin(), event_pos_begin);
            event_polygon_end = thrust::next(event_polygon_begin, splitted_polygon_count - planar_polygon_l_count);
            event_pos_begin = thrust::gather(p, std::exchange(event_polygon_begin, event_polygon_end), event_polygon_end, polygon.max.cbegin(), event_pos_begin);
            event_polygon_end = thrust::next(event_polygon_begin, splitted_polygon_count);
            event_pos_begin = thrust::gather(p, std::exchange(event_polygon_begin, event_polygon_end), event_polygon_end, polygon.min.cbegin(), event_pos_begin);
            assert(event.polygon.end() == thrust::next(event_polygon_begin, splitted_polygon_count - planar_polygon_r_count));
            event_pos_begin = thrust::gather(p, event_polygon_begin, event.polygon.end(), polygon.max.cbegin(), event_pos_begin);
            assert(event_pos_begin == event.pos.end());
#else
            auto event_polygon_bbox_begin = thrust::make_permutation_iterator(polygon_bbox_begin, event_l_polygon_l_begin);
            auto to_event_pos = [] __host__ __device__ (I event_kind, bbox_type bbox) -> F
            {
                return (event_kind < 0) ? thrust::get< 1 >(bbox) : thrust::get< 0 >(bbox);
            };
            thrust::transform(p, event_l_kind_l_begin, event.kind.end(), event_polygon_bbox_begin, thrust::next(event.pos.begin(), splitted_event_offset), to_event_pos);
#endif

            // calculate event node
            thrust::gather(p, event_l_polygon_l_begin, event.polygon.end(), polygon_node.cbegin(), thrust::next(event.node.begin(), splitted_event_offset));

            auto event_begin = thrust::make_zip_iterator(thrust::make_tuple(event.node.begin(), event.pos.begin(), event.kind.begin(), event.polygon.begin()));

            // sort splitted event
            auto splitted_event_begin = thrust::next(event_begin, splitted_event_offset);
            auto splitted_event_end = thrust::next(splitted_event_begin, splitted_event_count);
            thrust::sort(p, splitted_event_begin, splitted_event_end);

            // cleanup repeating planar events
            if (std::is_floating_point_v< F >) {
                auto clean_splitted_event_end = thrust::unique(p, splitted_event_begin, splitted_event_end);
                //std::cout << "CLEAN: " << thrust::distance(clean_splitted_event_end, splitted_event_end) << std::endl;
                event_count -= U(thrust::distance(clean_splitted_event_end, splitted_event_end));
                splitted_event_end = clean_splitted_event_end;
            }

            if ((false)) {
                decltype(event.node) node{thrust::get< 0 >(splitted_event_begin.get_iterator_tuple()), thrust::get< 0 >(splitted_event_end.get_iterator_tuple())};
                decltype(event.polygon) polygon{thrust::get< 3 >(splitted_event_begin.get_iterator_tuple()), thrust::get< 3 >(splitted_event_end.get_iterator_tuple())};
                auto it = thrust::make_zip_iterator(thrust::make_tuple(polygon.begin(), node.begin()));
                thrust::sort_by_key(p, polygon.begin(), polygon.end(), node.begin());
                auto d = thrust::make_discard_iterator();
                using P = thrust::tuple< U, U, U >;
                container< P > uniq_nodes{node.size()};
                auto transform = [] __host__ __device__ (U node) -> P
                {
                    return {node, 0, 1};
                };
                auto node_begin = thrust::make_transform_iterator(node.begin(), transform);
                auto sum = [] __host__ __device__ (const P & l, const P & r) -> P
                {
                    U node_l = thrust::get< 0 >(l);
                    U node_r = thrust::get< 0 >(r);
                    bool b = node_l == node_r;
                    if (!b) {
                        printf("nodes1 %u %u\n", node_l, node_r);
                    }
                    U eq = thrust::get< 1 >(l) + thrust::get< 1 >(r) + (b ? 0 : 1);
                    return {node_r, eq, thrust::get< 2 >(l) + thrust::get< 2 >(r)};
                };
                auto ends = thrust::reduce_by_key(p, polygon.begin(), polygon.end(), node_begin, d, uniq_nodes.begin(), thrust::equal_to< U >{}, sum);
                uniq_nodes.erase(ends.second, uniq_nodes.end());
                auto check = [] __host__ __device__ (const P & x)
                {
                    if (!((thrust::get< 2 >(x) == 2) || (thrust::get< 2 >(x) == 1))) {
                        printf("not two: %u %u\n", thrust::get< 2 >(x), thrust::get< 0 >(x));
                    }
                    assert(thrust::get< 1 >(x) == 0);
                    assert((thrust::get< 2 >(x) == 2) || (thrust::get< 2 >(x) == 1));
                };
                thrust::for_each(p, uniq_nodes.begin(), uniq_nodes.end(), check);
            }

            // merge splitted event w/ lr event
            auto event_lr_begin = thrust::next(event_begin, event_storage_size);
            auto event_end = thrust::merge(p, event_lr_begin, splitted_event_begin, splitted_event_begin, splitted_event_end, event_begin);
            assert(thrust::next(event_begin, event_count) == event_end);

            assert(thrust::is_sorted(p, event_begin, event_end));

            if ((false)) {
                auto event_polygon_bbox_begin1 = thrust::make_zip_iterator(thrust::make_tuple(event.kind.cbegin(), thrust::make_permutation_iterator(polygon_bbox_begin, event.polygon.cbegin())));
                using event_polygon_bbox_type1 = input_iterator_value_type< decltype(event_polygon_bbox_begin1) >;
                auto to_event_pos1 = [] __host__ __device__ (event_polygon_bbox_type1 x) -> F
                {
                    I event_kind = thrust::get< 0 >(x);
                    const bbox_type & bbox = thrust::get< 1 >(x);
                    return (event_kind < 0) ? thrust::get< 1 >(bbox) : thrust::get< 0 >(bbox);
                };
                assert(thrust::equal(p, event.pos.cbegin(), thrust::next(event.pos.cbegin(), event_count), thrust::make_transform_iterator(event_polygon_bbox_begin1, to_event_pos1)));
            }

            // crop
            event.node.resize(event_count);
            event.pos.resize(event_count);
            event.kind.resize(event_count);
            event.polygon.resize(event_count);

            assert(thrust::equal(p, event.node.cbegin(), event.node.cend(), thrust::make_permutation_iterator(polygon_node.cbegin(), event.polygon.cbegin())));

            if ((false)) {
                auto node = event.node;
                auto polygon = event.polygon;
                auto it = thrust::make_zip_iterator(thrust::make_tuple(polygon.begin(), node.begin()));
                thrust::sort_by_key(p, polygon.begin(), polygon.end(), node.begin());
                auto d = thrust::make_discard_iterator();
                using P = thrust::tuple< U, U, U >;
                container< P > uniq_nodes{node.size()};
                auto transform = [] __host__ __device__ (U node) -> P
                {
                    return {node, 0, 1};
                };
                auto node_begin = thrust::make_transform_iterator(node.begin(), transform);
                auto sum = [] __host__ __device__ (const P & l, const P & r) -> P
                {
                    U node_l = thrust::get< 0 >(l);
                    U node_r = thrust::get< 0 >(r);
                    bool b = node_l == node_r;
                    if (!b) {
                        printf("nodes %u %u\n", node_l, node_r);
                    }
                    U eq = thrust::get< 1 >(l) + thrust::get< 1 >(r) + (b ? 0 : 1);
                    return {node_r, eq, thrust::get< 2 >(l) + thrust::get< 2 >(r)};
                };
                auto ends = thrust::reduce_by_key(p, polygon.begin(), polygon.end(), node_begin, d, uniq_nodes.begin(), thrust::equal_to< U >{}, sum);
                uniq_nodes.erase(ends.second, uniq_nodes.end());
                auto check = [] __host__ __device__ (const P & x)
                {
                    if (!((thrust::get< 2 >(x) == 2) || (thrust::get< 2 >(x) == 1))) {
                        printf("not two: %u %u\n", thrust::get< 2 >(x), thrust::get< 0 >(x));
                    }
                    assert(thrust::get< 1 >(x) == 0);
                    assert((thrust::get< 2 >(x) == 2) || (thrust::get< 2 >(x) == 1));
                };
                thrust::for_each(p, uniq_nodes.begin(), uniq_nodes.end(), check);
            }
        }

        void set_node_count(U node_count)
        {
            node.min.resize(node_count);
            node.max.resize(node_count);
        }

        void split_node(U prev_base_node, U base_node,
                        const container< I > & node_split_dimension, const container< U > & node_split_pos,
                        const container< U > & node_l, const container< U > & node_r)
        {
            auto node_split_pos_begin = thrust::next(node_split_pos.cbegin(), prev_base_node);
            auto node_split_pos_end = thrust::next(node_split_pos.cbegin(), base_node);
            auto node_split_dimension_begin = thrust::next(node_split_dimension.cbegin(), prev_base_node);
            auto is_x = [] __host__ __device__ (I node_split_dimension) -> bool { return node_split_dimension == dimension; };
            thrust::scatter_if(p, node_split_pos_begin, node_split_pos_end, thrust::next(node_l.cbegin(), prev_base_node), node_split_dimension_begin, node.max.begin(), is_x);
            thrust::scatter_if(p, node_split_pos_begin, node_split_pos_end, thrust::next(node_r.cbegin(), prev_base_node), node_split_dimension_begin, node.min.begin(), is_x);
        }

    };

    projection< 0 > x{p};
    projection< 1 > y{p};
    projection< 2 > z{p};

    struct
    {
        container< U > triangle;
        container< U > node;
        container< I > side;
        struct
        {
            container< U > l, r; // corresponding left and right event in different best dimensions
        } event;
    } polygon;

    struct
    {
        container< I > split_dimension;
        container< F > split_pos;
        container< U > l, r; // left child node and right child node if not leaf, polygon range otherwise
        struct
        {
            container< U > node, l, r; // unique polygon count in current node, its left child node and right child node
        } polygon_count;
    } node; // TODO: optimize out node.l

    container< U > layer_node_offset;

    container< U > splitted_polygon;

    void calculate_node_best_split(const params & sah, U base_node, U node_count)
    {
        auto node_split_cost_begin = thrust::make_zip_iterator(thrust::make_tuple(x.layer.split_cost.cbegin(), y.layer.split_cost.cbegin(), z.layer.split_cost.cbegin()));
        auto node_split_pos_begin = thrust::make_zip_iterator(thrust::make_tuple(x.layer.split_pos.cbegin(), y.layer.split_pos.cbegin(), z.layer.split_pos.cbegin()));
        auto node_l_polygon_count_begin = thrust::make_zip_iterator(thrust::make_tuple(x.layer.polygon_count.l.cbegin(), y.layer.polygon_count.l.cbegin(), z.layer.polygon_count.l.cbegin()));
        auto node_r_polygon_count_begin = thrust::make_zip_iterator(thrust::make_tuple(x.layer.polygon_count.r.cbegin(), y.layer.polygon_count.r.cbegin(), z.layer.polygon_count.r.cbegin()));
        auto node_split_begin = thrust::make_zip_iterator(thrust::make_tuple(node_split_cost_begin, node_split_pos_begin, node_l_polygon_count_begin, node_r_polygon_count_begin));
        using node_split_type = input_iterator_value_type< decltype(node_split_begin) >;
        auto node_polygon_count_begin = thrust::next(node.polygon_count.node.cbegin(), base_node);
        auto node_best_split_begin = thrust::next(thrust::make_zip_iterator(thrust::make_tuple(node.split_dimension.begin(), node.split_pos.begin(), node.polygon_count.l.begin(), node.polygon_count.r.begin())), base_node);
        using node_best_split_type = output_iterator_value_type< decltype(node_best_split_begin) >;
        auto to_node_best_split = [sah] __host__ __device__ (node_split_type node_split, U node_polygon_count) -> node_best_split_type
        {
            assert(node_polygon_count != 0);
            const auto & node_split_cost = thrust::get< 0 >(node_split);
            const auto & node_split_pos = thrust::get< 1 >(node_split);
            const auto & node_l_polygon_count = thrust::get< 2 >(node_split);
            const auto & node_r_polygon_count = thrust::get< 3 >(node_split);
            F x = thrust::get< 0 >(node_split_cost);
            F y = thrust::get< 1 >(node_split_cost);
            F z = thrust::get< 2 >(node_split_cost);
            F best_node_split_cost = thrust::min(sah.intersection_cost * node_polygon_count, thrust::min(x, thrust::min(y, z)));
            if (!(best_node_split_cost < x)) {
                return {0, thrust::get< 0 >(node_split_pos), thrust::get< 0 >(node_l_polygon_count), thrust::get< 0 >(node_r_polygon_count)};
            } else if (!(best_node_split_cost < y)) {
                return {1, thrust::get< 1 >(node_split_pos), thrust::get< 1 >(node_l_polygon_count), thrust::get< 1 >(node_r_polygon_count)};
            } else if (!(best_node_split_cost < z)) {
                return {2, thrust::get< 2 >(node_split_pos), thrust::get< 2 >(node_l_polygon_count), thrust::get< 2 >(node_r_polygon_count)};
            } else {
                return {-1};
            }
        };
        auto is_node_not_empty = [] __host__ __device__ (U node_polygon_count) -> bool { return node_polygon_count != 0; };
        thrust::transform_if(p, node_split_begin, thrust::next(node_split_begin, node_count), node_polygon_count_begin, node_polygon_count_begin, node_best_split_begin, to_node_best_split, is_node_not_empty);
    }

    template< typename iterator >
    void print(iterator first, U count, const char * const what) const
    {
        std::copy(first, thrust::next(first, count), std::ostream_iterator< typename thrust::iterator_traits< iterator >::value_type >(std::cout << what, ", "));
        std::cout << "\n";
    }

    template< typename iterator >
    void print(iterator first, iterator last, const char * const what) const
    {
        std::copy(first, last, std::ostream_iterator< typename thrust::iterator_traits< iterator >::value_type >(std::cout << what, ", "));
        std::cout << "\n";
    }

    U get_splitted_polygon_count(U base_node, U node_count) const
    {
#if 0
        std::cout << "node_size " << node.split_pos.size() << "\n";
        std::cout << "base_node " << base_node << "\n";
        std::cout << "node_count " << node_count << "\n";
        print(x.layer.polygon_count.l.cbegin(), x.layer.polygon_count.l.cend(), "x.layer.polygon_count.l: ");
        print(x.layer.polygon_count.r.cbegin(), x.layer.polygon_count.r.cend(), "x.layer.polygon_count.r: ");
        print(y.layer.polygon_count.l.cbegin(), y.layer.polygon_count.l.cend(), "y.layer.polygon_count.l: ");
        print(y.layer.polygon_count.r.cbegin(), y.layer.polygon_count.r.cend(), "y.layer.polygon_count.r: ");
        print(z.layer.polygon_count.l.cbegin(), z.layer.polygon_count.l.cend(), "z.layer.polygon_count.l: ");
        print(z.layer.polygon_count.r.cbegin(), z.layer.polygon_count.r.cend(), "z.layer.polygon_count.r: ");
        print(thrust::next(node.split_pos.cbegin(), base_node), node_count, "node.split_pos: ");
        print(x.layer.split_pos.cbegin(), node_count, "x.layer.split_pos: ");
        print(y.layer.split_pos.cbegin(), node_count, "y.layer.split_pos: ");
        print(z.layer.split_pos.cbegin(), node_count, "z.layer.split_pos: ");
        print(thrust::next(node.split_dimension.cbegin(), base_node), node_count, "node.split_dimension: ");
        print(thrust::next(node.polygon_count.l.cbegin(), base_node), node_count, "node.polygon_count.l: ");
        print(thrust::next(node.polygon_count.r.cbegin(), base_node), node_count, "node.polygon_count.r: ");
        print(thrust::next(node.polygon_count.node.cbegin(), base_node), node_count, "node.polygon_count.node: ");
        std::cout << std::endl;
#endif

        auto node_polygon_count_begin = thrust::make_zip_iterator(thrust::make_tuple(node.split_dimension.cbegin(), node.polygon_count.l.cbegin(), node.polygon_count.r.cbegin(), node.polygon_count.node.cbegin()));
        using node_polygon_count_type = input_iterator_value_type< decltype(node_polygon_count_begin) >;
        auto to_splitted_polygon_count = [] __host__ __device__ (node_polygon_count_type node_polygon_count) -> U
        {
            I split_dimension = thrust::get< 0 >(node_polygon_count);
            if (split_dimension < 0) {
                return 0;
            }
            U polygon_count_l = thrust::get< 1 >(node_polygon_count);
            U polygon_count_r = thrust::get< 2 >(node_polygon_count);
            U polygon_count = thrust::get< 3 >(node_polygon_count);
            return polygon_count_l + polygon_count_r - polygon_count;
        };
        auto splitted_polygon_count_begin = thrust::next(node_polygon_count_begin, base_node);
        return thrust::transform_reduce(p, splitted_polygon_count_begin, thrust::next(splitted_polygon_count_begin, node_count), to_splitted_polygon_count, U(0), thrust::plus< U >{});
    }

    void separate_splitted_polygon(U base_node, U polygon_count, U splitted_polygon_count)
    {
        polygon.triangle.resize(polygon_count + splitted_polygon_count);
        polygon.node.resize(polygon_count + splitted_polygon_count);
        splitted_polygon.resize(splitted_polygon_count);

        auto polygon_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon.triangle.begin(), polygon.node.begin()));
        auto indexed_polygon_begin = thrust::make_zip_iterator(thrust::make_tuple(thrust::make_counting_iterator< U >(0), polygon_begin));
        auto split_dimension_begin = thrust::make_permutation_iterator(node.split_dimension.cbegin(), polygon.node.cbegin());
        auto splitted_polygon_stencil_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon.node.cbegin(), split_dimension_begin, polygon.side.cbegin()));
        using splitted_polygon_stencil_type = input_iterator_value_type< decltype(splitted_polygon_stencil_begin) >;
        auto splitted_polygon_begin = thrust::make_zip_iterator(thrust::make_tuple(splitted_polygon.begin(), thrust::next(polygon_begin, polygon_count)));
        auto is_splitted_polygon = [base_node] __host__ __device__ (splitted_polygon_stencil_type polygon_stencil) -> bool
        {
            U polygon_node = thrust::get< 0 >(polygon_stencil);
            if (polygon_node < base_node) {
                return false;
            }
            I split_dimension = thrust::get< 1 >(polygon_stencil);
            if (split_dimension < 0) {
                return false;
            }
            I polygon_side = thrust::get< 2 >(polygon_stencil);
            assert(polygon_side != 101325);
            if (polygon_side != 0) {
                return false;
            }
            return true;
        };
        auto splitted_polygon_end = thrust::copy_if(p, indexed_polygon_begin, thrust::next(indexed_polygon_begin, polygon_count), splitted_polygon_stencil_begin, splitted_polygon_begin, is_splitted_polygon);
        //std::cout << "DIST: " << thrust::distance(splitted_polygon_begin, splitted_polygon_end) << std::endl;
        assert(thrust::next(splitted_polygon_begin, splitted_polygon_count) == splitted_polygon_end);
    }

    void update_polygon_node(U base_node, U polygon_count)
    {
        auto node_lr_begin = thrust::make_zip_iterator(thrust::make_tuple(node.l.cbegin(), node.r.cbegin()));
        auto polygon_node_lr_begin = thrust::make_permutation_iterator(node_lr_begin, polygon.node.cbegin());
        using polygon_node_lr_type = input_iterator_value_type< decltype(polygon_node_lr_begin) >;
        auto to_polygon_node = [] __host__ __device__ (I polygon_side, polygon_node_lr_type polygon_node_lr) -> U
        {
            return (0 < polygon_side) ? thrust::get< 1 >(polygon_node_lr) : thrust::get< 0 >(polygon_node_lr); // splitted event assigned to l node
        };
        auto split_dimension_begin = thrust::make_permutation_iterator(node.split_dimension.cbegin(), polygon.node.cbegin());
        auto node_stencil_begin = thrust::make_zip_iterator(thrust::make_tuple(polygon.node.cbegin(), split_dimension_begin));
        using node_stencil_type = input_iterator_value_type< decltype(node_stencil_begin) >;
        auto is_current_layer = [base_node] __host__ __device__ (node_stencil_type splitted_node_stencil) -> bool
        {
            U polygon_node = thrust::get< 0 >(splitted_node_stencil);
            if (polygon_node < base_node) {
                return false;
            }
            I split_dimension = thrust::get< 1 >(splitted_node_stencil);
            if (split_dimension < 0) {
                return false;
            }
            return true;
        };
        thrust::transform_if(p, polygon.side.cbegin(), thrust::next(polygon.side.cbegin(), polygon_count), polygon_node_lr_begin, node_stencil_begin, polygon.node.begin(), to_polygon_node, is_current_layer);
    }

    void update_splitted_polygon_node(U polygon_count, U splitted_polygon_count)
    {
        auto splitted_polygon_node_begin = thrust::next(polygon.node.begin(), polygon_count);
        thrust::gather(p, splitted_polygon_node_begin, thrust::next(splitted_polygon_node_begin, splitted_polygon_count), node.r.cbegin(), splitted_polygon_node_begin);
    }

    U d = 0;

    tree build(const params & sah)
    {
        auto triangle_count = U(x.triangle.a.size());
        assert(triangle_count > 0);
        assert(triangle_count == U(y.triangle.a.size()));
        assert(triangle_count == U(z.triangle.a.size()));

        polygon.triangle.resize(triangle_count);
        thrust::sequence(p, polygon.triangle.begin(), polygon.triangle.end());

        x.calculate_triangle_bbox();
        y.calculate_triangle_bbox();
        z.calculate_triangle_bbox();

        x.caluculate_root_node_bbox();
        y.caluculate_root_node_bbox();
        z.caluculate_root_node_bbox();

        polygon.node.assign(triangle_count, U(0));

        x.generate_initial_event();
        y.generate_initial_event();
        z.generate_initial_event();

        // layer
        U base_node = 0;
        U node_count = 1;

        node.split_dimension.resize(1);
        node.split_pos.resize(1);
        node.l.resize(1);
        node.r.resize(1);
        node.polygon_count.node.assign(1, triangle_count);
        node.polygon_count.l.resize(1);
        node.polygon_count.r.resize(1);

        for (U depth = 0; (sah.max_depth == 0) || (depth < sah.max_depth); ++depth) {
            d = depth;
            std::cout << "\n\n\ndepth = " << depth << "\n";

            assert(thrust::is_sorted(p, x.event.node.cbegin(), x.event.node.cend()));
            assert(thrust::is_sorted(p, y.event.node.cbegin(), y.event.node.cend()));
            assert(thrust::is_sorted(p, z.event.node.cbegin(), z.event.node.cend()));

            assert(thrust::none_of(p, x.event.node.cbegin(), x.event.node.cend(), thrust::placeholders::_1 < base_node));
            assert(thrust::none_of(p, y.event.node.cbegin(), y.event.node.cend(), thrust::placeholders::_1 < base_node));
            assert(thrust::none_of(p, z.event.node.cbegin(), z.event.node.cend(), thrust::placeholders::_1 < base_node));

            assert(thrust::equal(p, x.event.node.cbegin(), x.event.node.cend(), thrust::make_permutation_iterator(polygon.node.cbegin(), x.event.polygon.cbegin())));
            assert(thrust::equal(p, y.event.node.cbegin(), y.event.node.cend(), thrust::make_permutation_iterator(polygon.node.cbegin(), y.event.polygon.cbegin())));
            assert(thrust::equal(p, z.event.node.cbegin(), z.event.node.cend(), thrust::make_permutation_iterator(polygon.node.cbegin(), z.event.polygon.cbegin())));

            {
                layer_node_offset.resize(node_count);

                auto is_node_not_empty = [] __host__ __device__ (U node_polygon_count) -> bool { return node_polygon_count != 0; };
                auto layer_node_begin = thrust::make_counting_iterator< U >(0);
                auto layer_node_offset_end = thrust::copy_if(p, layer_node_begin, thrust::next(layer_node_begin, node_count), thrust::next(node.polygon_count.node.cbegin(), base_node), layer_node_offset.begin(), is_node_not_empty);
                layer_node_offset.erase(layer_node_offset_end, layer_node_offset.end());

                std::cout << "node_count = " << node_count << std::endl;
                std::cout << "layer_node_offset_count = " << layer_node_offset.size() << std::endl;
            }
#if 0
            std::copy(x.event.node.cbegin(), x.event.node.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;
            std::copy(y.event.node.cbegin(), y.event.node.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;
            std::copy(z.event.node.cbegin(), z.event.node.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;

            std::copy(x.event.polygon.cbegin(), x.event.polygon.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;
            std::copy(y.event.polygon.cbegin(), y.event.polygon.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;
            std::copy(z.event.polygon.cbegin(), z.event.polygon.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;

            std::copy(polygon.node.cbegin(), polygon.node.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;
            std::copy(polygon.node.cbegin(), polygon.node.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;
            std::copy(polygon.node.cbegin(), polygon.node.cend(), std::ostream_iterator< U >(std::cout, ", "));
            std::cout << std::endl;
#endif

            x.find_perfect_split(sah, node_count, layer_node_offset, y, z);
            y.find_perfect_split(sah, node_count, layer_node_offset, z, x);
            z.find_perfect_split(sah, node_count, layer_node_offset, x, y);

            calculate_node_best_split(sah, base_node, node_count);

            auto node_split_dimension_begin = thrust::next(node.split_dimension.cbegin(), base_node);
            auto completed_node_count = U(thrust::count(p, node_split_dimension_begin, thrust::next(node_split_dimension_begin, node_count), I(-1)));
            std::cout << "completed node count " << completed_node_count << ' ' << node_count << std::endl;
            if (completed_node_count == node_count) {
                break;
            }

            auto polygon_count = U(polygon.triangle.size());

#ifdef NDEBUG
            polygon.side.resize(polygon_count);
            polygon.event.l.resize(polygon_count);
            polygon.event.r.resize(polygon_count);
#else
            polygon.side.assign(polygon_count, I(101325));
            polygon.event.l.assign(polygon_count, U(~0u));
            polygon.event.r.resize(polygon_count, U(~0u));
#endif

            x.calculate_polygon_side(node.split_dimension, base_node, polygon.event.l, polygon.event.r, polygon.side);
            y.calculate_polygon_side(node.split_dimension, base_node, polygon.event.l, polygon.event.r, polygon.side);
            z.calculate_polygon_side(node.split_dimension, base_node, polygon.event.l, polygon.event.r, polygon.side);

            U splitted_polygon_count = get_splitted_polygon_count(base_node, node_count);

            std::cout << "splitted_polygon_count = " << splitted_polygon_count << "\n";

            { // generate index for child node
                auto node_l_begin = thrust::next(node.l.begin(), base_node);
                auto to_node_count = [] __host__ __device__ (I node_split_dimension) -> U { return (node_split_dimension < 0) ? 0 : 2; };
                auto node_l_end = thrust::transform_exclusive_scan(p, node_split_dimension_begin, thrust::next(node_split_dimension_begin, node_count), node_l_begin, to_node_count, base_node + node_count, thrust::plus< U >{});

                auto node_r_begin = thrust::next(node.r.begin(), base_node);
                auto to_node_r = [] __host__ __device__ (U node_l) { return node_l + 1; };
                thrust::transform(p, node_l_begin, node_l_end, node_r_begin, to_node_r);
            }

            separate_splitted_polygon(base_node, polygon_count, splitted_polygon_count);
            // filled polygon.triangle and polygon.node (copied) for splitted polygons

            x.decouple_lr_event(node.split_dimension, polygon.side);
            y.decouple_lr_event(node.split_dimension, polygon.side);
            z.decouple_lr_event(node.split_dimension, polygon.side);
            // l and r event indices are copied

            update_polygon_node(base_node, polygon_count);
            // l and r polygon nodes updated, splitted polygon node turned into corresponding l polygon node

            x.split_polygon(node.split_dimension, node.split_pos, polygon.triangle, polygon.node, polygon_count, splitted_polygon_count, splitted_polygon, y, z);
            y.split_polygon(node.split_dimension, node.split_pos, polygon.triangle, polygon.node, polygon_count, splitted_polygon_count, splitted_polygon, z, x);
            z.split_polygon(node.split_dimension, node.split_pos, polygon.triangle, polygon.node, polygon_count, splitted_polygon_count, splitted_polygon, x, y);
            // bboxes of polygones resulting after split are updated (unsing polygon.node from splitted polygones part)

            update_splitted_polygon_node(polygon_count, splitted_polygon_count);

            x.merge_event(polygon_count, polygon.node, splitted_polygon_count, splitted_polygon);
            y.merge_event(polygon_count, polygon.node, splitted_polygon_count, splitted_polygon);
            z.merge_event(polygon_count, polygon.node, splitted_polygon_count, splitted_polygon);

            U prev_base_node = base_node;
            base_node += node_count;
            node_count -= completed_node_count;
            node_count += node_count;

            node.polygon_count.node.resize(base_node + node_count);

            auto is_not_leaf = [] __host__ __device__ (I node_split_dimension) -> bool { return !(node_split_dimension < 0); };

            thrust::scatter_if(p, thrust::next(node.polygon_count.l.cbegin(), prev_base_node), thrust::next(node.polygon_count.l.cbegin(), base_node), thrust::next(node.l.cbegin(), prev_base_node), node_split_dimension_begin, node.polygon_count.node.begin(), is_not_leaf);
            thrust::scatter_if(p, thrust::next(node.polygon_count.r.cbegin(), prev_base_node), thrust::next(node.polygon_count.r.cbegin(), base_node), thrust::next(node.r.cbegin(), prev_base_node), node_split_dimension_begin, node.polygon_count.node.begin(), is_not_leaf);

            x.set_node_count(base_node + node_count);
            y.set_node_count(base_node + node_count);
            z.set_node_count(base_node + node_count);

            // at first copy parent bbox to children one
            auto node_bbox_begin = thrust::make_zip_iterator(thrust::make_tuple(x.node.min.begin(), x.node.max.begin(), y.node.min.begin(), y.node.max.begin(), z.node.min.begin(), z.node.max.begin()));
            auto layer_bbox_begin = thrust::next(node_bbox_begin, prev_base_node);
            auto layer_bbox_end = thrust::next(node_bbox_begin, base_node);
            thrust::scatter_if(p, layer_bbox_begin, layer_bbox_end, thrust::next(node.l.cbegin(), prev_base_node), node_split_dimension_begin, node_bbox_begin, is_not_leaf);
            thrust::scatter_if(p, layer_bbox_begin, layer_bbox_end, thrust::next(node.r.cbegin(), prev_base_node), node_split_dimension_begin, node_bbox_begin, is_not_leaf);

            // then move wall
            x.split_node(prev_base_node, base_node, node.split_dimension, node.split_pos, node.l, node.r);
            y.split_node(prev_base_node, base_node, node.split_dimension, node.split_pos, node.l, node.r);
            z.split_node(prev_base_node, base_node, node.split_dimension, node.split_pos, node.l, node.r);

            node.split_dimension.resize(base_node + node_count, I(-1));
            node.split_pos.resize(base_node + node_count);
            node.l.resize(base_node + node_count);
            node.r.resize(base_node + node_count);
            node.polygon_count.l.resize(base_node + node_count);
            node.polygon_count.r.resize(base_node + node_count);
        }
        // calculate node parent
        // sort value (polygon) by key (polygon.node)
        // reduce value (counter, 1) by operation (project1st, plus) and key (node) to (key (node), value (offset, count))
        // scatter value (offset, count) to (node.l, node.r) at key (node)
        return {};
    }

};

template< template< typename ... > class container, typename policy >
builder< policy, container > make_builder(const policy & p)
{
    return {p};
}

}