reverb_time.h

#pragma once

#include "core/geo/triangle_vec.h"
#include "core/scene_data.h"

#include "utilities/aligned/set.h"

namespace wayverb {
namespace core {


template <typename Vertex, typename Surface>
double area(const generic_scene_data<Vertex, Surface>& scene,
            size_t surface_index) {
    //  TODO this is dumb because of the linear search.
    //  If this becomes a bottleneck, maybe scene_data could store surfaces
    //  pre-sorted by material.
    return std::accumulate(
            begin(scene.get_triangles()),
            end(scene.get_triangles()),
            0.0,
            [&](auto running_total, auto tri) {
                return running_total +
                       (tri.surface == surface_index
                                ? geo::area(geo::get_triangle_vec3(
                                          tri, scene.get_vertices().data()))
                                : 0.0);
            });
}

template <typename Vertex, typename Surface>
double area(const generic_scene_data<Vertex, Surface>& scene) {
    //  This is OK - we have to look at every triangle anyway.
    return std::accumulate(
            begin(scene.get_triangles()),
            end(scene.get_triangles()),
            0.0,
            [&](auto running_total, auto tri) {
                return running_total +
                       geo::area(geo::get_triangle_vec3(
                               tri, scene.get_vertices().data()));
            });
}

template <typename Vertex, typename Surface>
auto absorption_area(const generic_scene_data<Vertex, Surface>& scene,
                     size_t surface_index) {
    return scene.get_surfaces()[surface_index].absorption *
           area(scene, surface_index);
}

template <typename Vertex, typename Surface>
auto equivalent_absorption_area(
        const generic_scene_data<Vertex, Surface>& scene) {
    const auto num_surfaces = scene.get_surfaces().size();
    decltype(absorption_area(scene, 0)) running_total{};
    for (auto i = 0u; i != num_surfaces; ++i) {
        running_total += absorption_area(scene, i);
    }
    return running_total;
}


std::array<std::pair<cl_uint, cl_uint>, 3> get_index_pairs(const triangle& t);

template <typename It>  
bool triangles_are_oriented(It begin, It end) {
    //  The scene is consistently oriented if no triangles in the scene share
    //  an edge with the same winding direction.
    util::aligned::set<std::pair<cl_uint, cl_uint>> table;

    //  For each triangle.
    for (; begin != end; ++begin) {
        //  For each pair of vertices.
        for (const auto& pair : get_index_pairs(*begin)) {
            //  See whether this pair of vertices has already been found.
            if (!table.insert(pair).second) {
                //  The pair already existed.
                return false;
            }
        }
    }
    //  No triangles share a pair of vertices.
    return true;
}

float six_times_tetrahedron_volume(const geo::triangle_vec3& t);

template <typename Vertex, typename Surface>
float estimate_room_volume(const generic_scene_data<Vertex, Surface>& scene) {
    return std::abs(std::accumulate(
                   begin(scene.get_triangles()),
                   end(scene.get_triangles()),
                   0.0f,
                   [&](auto running_total, auto tri) {
                       return running_total +
                              six_times_tetrahedron_volume(
                                      geo::get_triangle_vec3(
                                              tri,
                                              scene.get_vertices().data()));
                   })) /
           6;
}


float estimate_air_intensity_absorption(float frequency, float humidity);

template <size_t... Ix>
auto estimate_air_intensity_absorption(
        const std::array<float, sizeof...(Ix)>& band_centres,
        float humidity,
        std::index_sequence<Ix...>) {
    return std::array<float, sizeof...(Ix)>{{estimate_air_intensity_absorption(
            std::get<Ix>(band_centres), humidity)...}};
}

template <size_t Bands>
auto estimate_air_intensity_absorption(
        const std::array<float, Bands>& band_centres, float humidity) {
    return estimate_air_intensity_absorption(
            band_centres, humidity, std::make_index_sequence<Bands>{});
}


//  USE THESE ONES

template <typename Absorption, typename Coeff>
auto sabine_reverb_time(double room_volume,
                        Absorption absorption_area,
                        Coeff air_coefficient) {
    if (room_volume <= 0) {
        throw std::runtime_error{"Room volume is non-positive"};
    }
    if (any(absorption_area <= 0)) {
        throw std::runtime_error{"Absorption area is non-positive"};
    }

    const auto numerator = 0.161f * room_volume;
    const auto denominator =
            absorption_area + (4 * room_volume * air_coefficient);
    return numerator / denominator;
}

template <typename Vertex, typename Surface, typename Coeff>
auto sabine_reverb_time(const generic_scene_data<Vertex, Surface>& scene,
                        Coeff air_coefficient) {
    return sabine_reverb_time(estimate_room_volume(scene),
                              equivalent_absorption_area(scene),
                              air_coefficient);
}

template <typename Absorption, typename Coeff>
auto eyring_reverb_time(double room_volume,
                        Absorption absorption_area,
                        double full_area,
                        Coeff air_coefficient) {
    if (room_volume <= 0) {
        throw std::runtime_error{"Room volume is non-positive"};
    }
    if (any(absorption_area <= 0)) {
        throw std::runtime_error{"Absorption area is non-positive"};
    }

    const auto numerator = 0.161f * room_volume;
    const auto denominator =
            -full_area * log(1 - (absorption_area / full_area)) +
            (4 * room_volume * air_coefficient);
    return numerator / denominator;
}

template <typename Vertex, typename Surface, typename Coeff>
auto eyring_reverb_time(const generic_scene_data<Vertex, Surface>& scene,
                        Coeff air_coefficient) {
    return eyring_reverb_time(estimate_room_volume(scene),
                              equivalent_absorption_area(scene),
                              area(scene),
                              air_coefficient);
}

}  // namespace core
}  // namespace wayverb