simplex.h

Go to the documentation of this file.

 
#pragma once
 
#include "../../algorithm.h"
#include "../matrix.h"
#include "./point.h"
 
#include <algorithm>
#include <array>
#include <vector>
 
namespace emlabcpp
{
 
template < typename Item, std::size_t N >
class simplex
{
public:
        using container = std::array< Item, N + 1 >;
 
private:
        container points_;
 
public:
        using value_type     = Item;
        using iterator       = typename container::iterator;
        using const_iterator = typename container::const_iterator;
 
        constexpr explicit simplex() = default;
 
        constexpr explicit simplex( Item item )
          : points_( { item } )
        {
        }
 
        constexpr simplex( simplex< Item, N - 1 > const& sub, Item item )
          : points_()
        {
                std::copy( sub.begin(), sub.end(), points_.begin() );
                *points_.rbegin() = item;
 
                std::sort( points_.begin(), points_.end() );
        }
 
        constexpr simplex( std::array< Item, N + 1 > data )
          : points_( std::move( data ) )
        {
                std::sort( points_.begin(), points_.end() );
        }
 
        template < typename... Ts >
        constexpr explicit simplex( Ts... ts )
          : points_( { ts... } )
        {
                static_assert(
                    sizeof...( Ts ) == N + 1, "Number of parameters for simplex has to be N+1" );
        }
 
        [[nodiscard]] constexpr const_iterator begin() const
        {
                return points_.begin();
        }
 
        [[nodiscard]] constexpr const_iterator end() const
        {
                return points_.end();
        }
 
        [[nodiscard]] constexpr Item const& operator[]( std::size_t index ) const
        {
                return points_[index];
        }
 
        [[nodiscard]] constexpr std::size_t size() const
        {
                return points_.size();
        }
};
 
template < std::size_t N, std::size_t U >
constexpr point< N > center_of( simplex< point< N >, U > const& s )
{
        vector< N > avg = sum( s, [&]( point< N > const& p ) -> vector< N > {
                return vector_cast( p ) / s.size();
        } );
 
        return point_cast( avg );
}
 
template < std::size_t N >
constexpr float volume_of( simplex< point< N >, N > const& simplex )
{
        matrix< N, N > m;
        for ( std::size_t const i : range( simplex.size() - 1 ) ) {
                auto diff = simplex[i + 1] - simplex[0];
                std::copy( diff.begin(), diff.end(), &m[i][0] );
        }
 
        return std::abs( ( 1 / std::tgammaf( N + 1 ) ) * determinant( m ) );
}
 
template < typename Item, std::size_t N >
constexpr bool operator<( simplex< Item, N > const& lh, simplex< Item, N > const& rh )
{
        // there are actually N+1 items in N simplex;
        std::size_t i = *find_if( range( N ), [&]( std::size_t j ) {
                return lh[j] != rh[j];
        } );
        return lh[i] < rh[i];
}
 
template < typename Item, std::size_t N >
constexpr bool operator>( simplex< Item, N > const& lh, simplex< Item, N > const& rh )
{
        return rh < lh;
}
 
template < typename Item, std::size_t N >
constexpr bool operator<=( simplex< Item, N > const& lh, simplex< Item, N > const& rh )
{
        return !( lh > rh );
}
 
template < typename Item, std::size_t N >
constexpr bool operator>=( simplex< Item, N > const& lh, simplex< Item, N > const& rh )
{
        return !( lh < rh );
}
 
template < typename Item, std::size_t N >
constexpr bool operator==( simplex< Item, N > const& lh, simplex< Item, N > const& rh )
{
        return equal( lh, rh );
}
 
template < typename Item, std::size_t N >
constexpr bool operator!=( simplex< Item, N > const& lh, simplex< Item, N > const& rh )
{
        return !( lh == rh );
}
 
}  // namespace emlabcpp