<commit_before>#include "kdtree.h"
#include "clipping.h"

#include <stack>


namespace {

// Cf. [WH06], 5.2, Table 1
static constexpr int COST_TRAVERSAL = 15;
static constexpr int COST_INTERSECTION = 20;


// Cost function bias
inline float lambda(size_t num_ltris, size_t num_rtris) {
    if (num_ltris == 0 || num_rtris == 0) {
        return 0.8f;
    }
    return 1;
}

//
// Cost function of splitting box b at a given plane p
//
// Args:
//   {l,r}area_ratio - ratio of the surface area of the left resp. right
//     box over the box
//   num_{l,r}tris - number of triangles in the left resp. right box
//
// Return:
//   cost to split the box
//
inline float cost(
    float larea_ratio, float rarea_ratio, size_t num_ltris, size_t num_rtris)
{
    return lambda(num_ltris, num_rtris) *
        (COST_TRAVERSAL + COST_INTERSECTION *
            (larea_ratio * num_ltris + rarea_ratio * num_rtris));
}

enum class Dir { LEFT, RIGHT };

//
// SAH function
//
// Args:
//   ax, pos: splitting plane
//   box: AABB to split
//   num_{l,r}tris: number of triangles in the left resp. right box produces
//     by splitting `box` at `p`
//   num_ptris: number of triangles lying in the plane `p`
//
// Return:
//   cost to split the box + wether the planar triangles should be appended to
//   the lhs or rhs of the box
//
inline std::pair<float /*cost*/, Dir> surface_area_heuristics(
    Axis ax, float pos, const Box& box,
    size_t num_ltris, size_t num_rtris, size_t num_ptris)
{
    Box lbox, rbox;
    std::tie(lbox, rbox) = box.split(ax, pos);
    float area = box.surface_area();
    float larea_ratio = lbox.surface_area()/area;
    float rarea_ratio = rbox.surface_area()/area;

    auto left_planar_cost = cost(
        larea_ratio, rarea_ratio,
        num_ltris + num_ptris, num_rtris);
    auto right_planar_cost = cost(
        larea_ratio, rarea_ratio,
        num_ltris, num_ptris + num_rtris);

    if (left_planar_cost < right_planar_cost) {
        return {left_planar_cost, Dir::LEFT};
    } else {
        return {right_planar_cost, Dir::RIGHT};
    }
}

}


KDTree::KDTree(Triangles tris)
    : tris_(std::move(tris))
{
    assert(tris_.size() > 0);
    assert(tris_.size() <= Node::MAX_ID);

    std::vector<TriangleId> ids(tris_.size(), 0);

    box_ = tris_.front().bbox();
    for (size_t i = 1; i < tris_.size(); ++i) {
        box_ = box_ + tris_[i].bbox();
        ids[i] = i;
    }

    root_ = std::unique_ptr<Node>(build(std::move(ids), box_));
}

KDTree::Node* KDTree::build(KDTree::TriangleIds tris, const Box& box) {
    if (tris.size() == 0) {
        return nullptr;
    }

    // to few triangles -> terminate
    if (tris.size() <= 3) {
        return new Node(tris);
    }

    // box too small -> no need to split further -> terminate
    if (box.surface_area() == 0) {
        return new Node(tris);
    }

    float min_cost;
    Axis plane_ax; float plane_pos;  // plane
    TriangleIds ltris, rtris;
    std::tie(min_cost, plane_ax, plane_pos, ltris, rtris) =
        find_plane_and_classify(tris, box);

    // automatic termination
    // remove lambda factor from cost again, otherwise we may stuck in an
    // empty space split forever
    if (COST_INTERSECTION * tris.size()
        * lambda(ltris.size(), rtris.size()) < min_cost)
    {
        return new Node(tris);
    }

    Box lbox, rbox;
    std::tie(lbox, rbox) = box.split(plane_ax, plane_pos);

    return new Node(plane_ax, plane_pos,
        build(std::move(ltris), lbox),
        build(std::move(rtris), rbox));
}

std::tuple<
    float /*cost*/, Axis /* plane axis */, float /* plane pos */,
    KDTree::TriangleIds /*left*/, KDTree::TriangleIds /*right*/>
KDTree::find_plane_and_classify(const TriangleIds& tris, const Box& box) const
{
    // The box should have some surface, otherwise the surface area heuristics
    // does not make any sense.
    assert(box.surface_area() != 0);

    // auxiliary event structure
    static constexpr int STARTING = 2;
    static constexpr int ENDING = 0;
    static constexpr int PLANAR = 1;

    struct Event {
        TriangleId id;
        float point;
        // auxiliary point
        // for type == ENDING, point_aux is the starting point
        // for type == STARTNG, point_aux is the ending point
        // for type == PLANAR, points_aux is the same point as `point`
        float point_aux;
        int type;
    };

    std::vector<Event> event_lists[AXES.size()];

    // generate events
    size_t num_tris = 0;
    for (const auto& id : tris) {
        auto clipped_box = clip_triangle_at_aabb(tris_[id], box);
        if (clipped_box.is_trivial()) {
            continue;
        }
        num_tris += 1;

        for (auto ax : AXES) {
            auto& events = event_lists[static_cast<int>(ax)];
            if (clipped_box.is_planar(ax)) {
                events.emplace_back(
                    Event{id, clipped_box.min[ax], clipped_box.min[ax],
                    PLANAR});
            } else {
                events.emplace_back(
                    Event{id, clipped_box.min[ax], clipped_box.max[ax],
                    STARTING});
                events.emplace_back(
                    Event{id, clipped_box.max[ax], clipped_box.min[ax],
                    ENDING});
            }
        }
    }

    // all clipped?
    if (num_tris == 0) {
        return std::make_tuple(0, Axis::X, 0, TriangleIds(), TriangleIds());
    }

    // sweep for min_cost
    float min_cost = std::numeric_limits<float>::max();
    Axis min_plane_ax; float min_plane_pos;
    Dir min_side;
    // we also store those for asserts below
    float min_ltris, min_rtris, min_ptris;

    for (auto ax : AXES) {
        auto& events = event_lists[static_cast<int>(ax)];
        std::sort(events.begin(), events.end(),
            [](const Event& e1, const Event& e2) {
                return e1.point < e2.point ||
                    (e1.point == e2.point && e1.type < e2.type);
            });

        int num_ltris = 0, num_ptris = 0, num_rtris = num_tris;

        for (size_t i = 0; i < events.size(); ) {
            auto& event = events[i];

            auto p = event.point;
            int point_starting = 0;
            int point_ending = 0;
            int point_planar = 0;
<commit_msg>Fix kdtree: Return max cost in case all triangles were clipped.<commit_after>#include "kdtree.h"
#include "clipping.h"

#include <stack>


namespace {

// Cf. [WH06], 5.2, Table 1
static constexpr int COST_TRAVERSAL = 15;
static constexpr int COST_INTERSECTION = 20;


// Cost function bias
inline float lambda(size_t num_ltris, size_t num_rtris) {
    if (num_ltris == 0 || num_rtris == 0) {
        return 0.8f;
    }
    return 1;
}

//
// Cost function of splitting box b at a given plane p
//
// Args:
//   {l,r}area_ratio - ratio of the surface area of the left resp. right
//     box over the box
//   num_{l,r}tris - number of triangles in the left resp. right box
//
// Return:
//   cost to split the box
//
inline float cost(
    float larea_ratio, float rarea_ratio, size_t num_ltris, size_t num_rtris)
{
    return lambda(num_ltris, num_rtris) *
        (COST_TRAVERSAL + COST_INTERSECTION *
            (larea_ratio * num_ltris + rarea_ratio * num_rtris));
}

enum class Dir { LEFT, RIGHT };

//
// SAH function
//
// Args:
//   ax, pos: splitting plane
//   box: AABB to split
//   num_{l,r}tris: number of triangles in the left resp. right box produces
//     by splitting `box` at `p`
//   num_ptris: number of triangles lying in the plane `p`
//
// Return:
//   cost to split the box + wether the planar triangles should be appended to
//   the lhs or rhs of the box
//
inline std::pair<float /*cost*/, Dir> surface_area_heuristics(
    Axis ax, float pos, const Box& box,
    size_t num_ltris, size_t num_rtris, size_t num_ptris)
{
    Box lbox, rbox;
    std::tie(lbox, rbox) = box.split(ax, pos);
    float area = box.surface_area();
    float larea_ratio = lbox.surface_area()/area;
    float rarea_ratio = rbox.surface_area()/area;

    auto left_planar_cost = cost(
        larea_ratio, rarea_ratio,
        num_ltris + num_ptris, num_rtris);
    auto right_planar_cost = cost(
        larea_ratio, rarea_ratio,
        num_ltris, num_ptris + num_rtris);

    if (left_planar_cost < right_planar_cost) {
        return {left_planar_cost, Dir::LEFT};
    } else {
        return {right_planar_cost, Dir::RIGHT};
    }
}

}


KDTree::KDTree(Triangles tris)
    : tris_(std::move(tris))
{
    assert(tris_.size() > 0);
    assert(tris_.size() <= Node::MAX_ID);

    std::vector<TriangleId> ids(tris_.size(), 0);

    box_ = tris_.front().bbox();
    for (size_t i = 1; i < tris_.size(); ++i) {
        box_ = box_ + tris_[i].bbox();
        ids[i] = i;
    }

    root_ = std::unique_ptr<Node>(build(std::move(ids), box_));
}

KDTree::Node* KDTree::build(KDTree::TriangleIds tris, const Box& box) {
    if (tris.size() == 0) {
        return nullptr;
    }

    // to few triangles -> terminate
    if (tris.size() <= 3) {
        return new Node(tris);
    }

    // box too small -> no need to split further -> terminate
    if (box.surface_area() == 0) {
        return new Node(tris);
    }

    float min_cost;
    Axis plane_ax; float plane_pos;  // plane
    TriangleIds ltris, rtris;
    std::tie(min_cost, plane_ax, plane_pos, ltris, rtris) =
        find_plane_and_classify(tris, box);

    // automatic termination
    // remove lambda factor from cost again, otherwise we may stuck in an
    // empty space split forever
    if (COST_INTERSECTION * tris.size()
        * lambda(ltris.size(), rtris.size()) < min_cost)
    {
        return new Node(tris);
    }

    Box lbox, rbox;
    std::tie(lbox, rbox) = box.split(plane_ax, plane_pos);

    return new Node(plane_ax, plane_pos,
        build(std::move(ltris), lbox),
        build(std::move(rtris), rbox));
}

std::tuple<
    float /*cost*/, Axis /* plane axis */, float /* plane pos */,
    KDTree::TriangleIds /*left*/, KDTree::TriangleIds /*right*/>
KDTree::find_plane_and_classify(const TriangleIds& tris, const Box& box) const
{
    // The box should have some surface, otherwise the surface area heuristics
    // does not make any sense.
    assert(box.surface_area() != 0);

    // auxiliary event structure
    static constexpr int STARTING = 2;
    static constexpr int ENDING = 0;
    static constexpr int PLANAR = 1;

    struct Event {
        TriangleId id;
        float point;
        // auxiliary point
        // for type == ENDING, point_aux is the starting point
        // for type == STARTNG, point_aux is the ending point
        // for type == PLANAR, points_aux is the same point as `point`
        float point_aux;
        int type;
    };

    std::vector<Event> event_lists[AXES.size()];

    // generate events
    size_t num_tris = 0;
    for (const auto& id : tris) {
        auto clipped_box = clip_triangle_at_aabb(tris_[id], box);
        if (clipped_box.is_trivial()) {
            continue;
        }
        num_tris += 1;

        for (auto ax : AXES) {
            auto& events = event_lists[static_cast<int>(ax)];
            if (clipped_box.is_planar(ax)) {
                events.emplace_back(
                    Event{id, clipped_box.min[ax], clipped_box.min[ax],
                    PLANAR});
            } else {
                events.emplace_back(
                    Event{id, clipped_box.min[ax], clipped_box.max[ax],
                    STARTING});
                events.emplace_back(
                    Event{id, clipped_box.max[ax], clipped_box.min[ax],
                    ENDING});
            }
        }
    }

    // all clipped?
    if (num_tris == 0) {
        return std::make_tuple(std::numeric_limits<float>::max(), Axis::X, 0,
                               TriangleIds(), TriangleIds());
    }

    // sweep for min_cost
    float min_cost = std::numeric_limits<float>::max();
    Axis min_plane_ax; float min_plane_pos;
    Dir min_side;
    // we also store those for asserts below
    float min_ltris, min_rtris, min_ptris;

    for (auto ax : AXES) {
        auto& events = event_lists[static_cast<int>(ax)];
        std::sort(events.begin(), events.end(),
            [](const Event& e1, const Event& e2) {
                return e1.point < e2.point ||
                    (e1.point == e2.point && e1.type < e2.type);
            });

        int num_ltris = 0, num_ptris = 0, num_rtris = num_tris;

        for (size_t i = 0; i < events.size(); ) {
            auto& event = events[i];

            auto p = event.point;
            int point_starting = 0;
            int point_ending = 0;
            int point_planar = 0;

{'fraction_non_alphanumeric': 0.12483210871454531, 'fraction_numerical': 0.006715651418187564, 'mean_word_length': 2.2784252784252783, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4773468', 'n_tokens_mistral': 4501, 'n_tokens_neox': 4279, 'n_words': 1152}