viltrum

VILTRUM: Varied Integration Layouts for arbiTRary integrals in a Unified Manner - A C++17 header-only library that provides a set of numerical integration routines. This library was generated during research for our paper:

Primary-Space Adaptive Control Variates using Piecewise-Polynomial Approximations

ACM Transactions on Graphics - 2021
Miguel Crespo · Adrian Jarabo · Adolfo Muñoz

Installation

viltrum is a header-only library and has no external dependencies, so installation is rather simple (it does not require any compilation process). You just need to download the library in a specific folder:

folder> git clone https://github.com/adolfomunoz/viltrum.git

Then make sure folder is within the included directories (-Ifolder parameter in g++, include_directories("folder") in CMake) and include it from C++.

#include <viltrum/viltrum.h>

That would be enough to use all the features of the library. There are other alternatives that might be more comfortable for you, see them here. There is a CMake-based building system for several example and test executable files that automatically downloads external dependencies and compiles all executables but it is not needed for the libray's usage.

Usage

Integrating a function in a specific n-dimensional range is rather simple. You need the following information:

• An integrator, a numerical integration method, for which there are several to choose from.
• An integrand, a function to be integrated. It's only parameter has to be a std::array<F,N>, where F is a floating point number and N is the number of dimensions. There are several ways in which such integrand can be defined. The integrand can return any numeric value, or, in general, any data type that supports addition and multiplication by a scalar (tested with Eigen arrays ).
• A range, the integration domain, that is composed on two std::array<F,N> marking the limits of the potentially multidimensional integration domain, which can be defined in different ways.

Example:

float sphere(const std::array<float,3>& x) {
    return (x[0]*x[0]+x[1]*x[1]+x[2]*x[2])<=1.0f?1.0f:0.0f;
}
int main() {
    unsigned long samples = 1024;
    auto method = viltrum::integrator_monte_carlo_uniform(samples);
    auto range = viltrum::range(std::array<float,3>{-1.0f,-1.0f,-1.0f},std::array<float,3>{1.0f,1.0f,1.0f});
    std::cout<<"Sphere volume = "<<method.integrate(sphere,range)<<"\n";
}

The return type of the integrate method will be the same return type of the integrand (float in the example above).

Integrating into bins

It is also possible to obtain not only a single integral value, but to obtain integrals in a regular n-dimensional grid of cells or bins. For this purpose, this library provides bin integrators, that is, integrators that are capable of integrating into bins. For using them, in addition to the integrand and range you need:

• A bin accesor, that is, a function that given an n-dimensional array (where n is the dimensionality of the bin regular structure) of positions (of type std::size_t) gives read and write access to the binning structure in that specific position.
• A resolution, a n-dimensional array of std::size_t that gives the resolution of the bin structure for each dimension n.

Given these common structure:Bin integrators can be used as follows:

float slope(const std::array<float,2>& x) { return (x[1]<x[0])?1.0f:0.0f; }

int main(int argc, char **argv) {
    auto integrator_bins = viltrum::integrator_bins_monte_carlo_uniform(1024);
    auto range = viltrum::range(std::array<float,2>{0,0},std::array<float,2>{1,1});
    float output_array[16];
    //...
}

bin integrators can be used as follows (1D example):

    auto output_array_access = [&output_array] (const std::array<std::size_t,1>& i) -> float& { return output_array[i[0]]; }; 
    integrator_bins.integrate(output_array_access,std::array<std::size_t,1>{16},slope,range);
    for (float f : output_array) std::cout<<f<<" "; std::cout<<std::endl;

and can be expanded to two dimensions as follows:

    auto output_array_access_2d = [&output_array] (const std::array<std::size_t,2>& i) -> float& { return output_array[4*i[0]+i[1]]; }; 
    integrator_bins.integrate(output_array_access_2d,std::array<std::size_t,2>{4,4},slope,range);

Additionally, the viltrum library offers a functional version of the same method:

    integrate_bins(integrator_bins,output_array_access,std::array<std::size_t,1>{16},slope,range);
    integrate_bins(integrator_bins,output_array_access_2d,std::array<std::size_t,2>{4,4},slope,range);

This functional version is able to deduce the bin accesor as well as the resolution for std::vector data types, in which there is only one parameter defining the binning structure:

    std::vector<float> output_vector(16);
    integrate_bins(integrator_bins, output_vector, slope, range);
    for (float f : output_vector) std::cout<<f<<" ";
    std::cout<<std::endl;
    
    std::vector<std::vector<float>> output_matrix(4,std::vector<float>(4));
    integrate_bins(integrator_bins, output_matrix, slope, range);
    std::cout<<std::endl;
    for (const std::vector<float>& row : output_matrix) {
        for (float f : row) std::cout<<f<<" ";
        std::cout<<std::endl;
    }
    std::cout<<std::endl;

If you want to define your own automatic bin accesors you can take inspiration from our source code for the case of vectors but we do not provide documentation for that (you can always use the general version by explicitly indicating bin accesor and resolution independently as illustrated above).

Last, there are several bin integrators to choose from. However, if you want to use exactly the same integrator with the default parameters than we use in our paper define the bin integrator as follows:

unsigned long spp_cv = std::max(1UL,(unsigned long)(spp*(1.0/16.0)));
auto integrator = integrator_optimized_perpixel_adaptive_stratified_control_variates(
    viltrum::nested(viltrum::simpson,viltrum::trapezoidal), // nested rule, order 2 polynomials
    viltrum::error_single_dimension_size(1.e-5), // error heuristic
    spp_cv*bins/(2*std::pow(3, dim-1)), // number of adaptive iterations calculated from the spps
    std::max(1UL,spp-spp_cv) // number of spps for the residual
);

where:

• spp is the number of samples per pixel.
• dim is the number of dimensions to be explored.
• bins is the total number of bins (pixels in the case of images) in all dimensions.

License

This code is released under the GPL v3. Additionally, if you are using this code in academic research, we would be grateful if you cited our paper, for which we generated with this source code:

@article{crespo21primary,
    title = {Primary-Space Adaptive Control Variates using Piecewise-Polynomial Approximations.},
    year = {2021},
    journal = {ACM Transactions on Graphics},
    author = {Crespo, Miguel and Jarabo, Adrian and Mu\~{n}oz, Adolfo},
    volume = {40},
    number = {3},
    issn = {0730-0301},
    url = {https://doi.org/10.1145/3450627},
    doi = {10.1145/3450627},
    issue_date = {July 2021},
    month = jul,
    articleno = {25},
    numpages = {15},
    keywords = {adaptive quadrature, numerical integration, control variates, rendering, piecewise polynomial}
}