// COMPILE_LINE: g++ -o dali22_code dali22_code.cpp -lm -lstdc++ -lgomp -fopenmp
// Supplementary material code for the paper Soler, Molazem, Subr'2022,
// A Theoretical Analysis of Compactness of the Light Transport Operator
//
// Copyright (C) 2022 Cyril Soler (cyril.soler@inria.fr)
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
#include
#include
#include
#include
#include
#include
void test_compactness(double p_distance);
int main(int,char *[])
{
double p = 2.0;
test_compactness(p);
return 0;
}
// Computes the form factor of patch to a disk. The disk is tangent to the support plane of
// the patch and orthogonal to it. Parameters are:
// r : radius of the disk
// x,y: coordinates of the patch in its plane, assuming the tangent point of the disk is (0,0),
// x is the distance orthogonal to the axis of the disk
// y is the distance along the axis of the disk
//
// We use the formula of Eq.2 in
// NARAGHI, M. H. N. (1988).
// Radiation view factors from differential plane sources to disks - A general formulation.
// Journal of Thermophysics and Heat Transfer, 2(3), 271–274. doi:10.2514/3.96
//
// using:
// theta=pi/2, h=y, b=x, a=r, phi=pi/2
double point_to_tangent_disk_form_factor(double x,double y,double r)
{
double alpha = 2*r*r+x*x+y*y;
return 0.5*r*y/(r*r+x*x) * (alpha/sqrt(alpha*alpha - 4*r*r*(r*r+x*x))-1);
}
// Monte-Calo integration version of the previous function, optionaly used for testing.
//
double point_to_tangent_disk_form_factor_MC(double x,double y,double r)
{
// uniform sample the disk
uint32_t NR = 100;
uint32_t NT = 300;
double z = 0.0;
double res = 0.0;
for(uint32_t i=0;i static std::string number(T x,int decimals=0,bool leading_zeros=false)
{
std::ostringstream o ;
if(leading_zeros)
{
o << std::setw(decimals);
o << std::setiosflags(std::ios::fixed);
o << std::setfill('0') ;
}
if(decimals > 0) o << std::setprecision(decimals) ;
o << x ;
return o.str() ;
}
bool test_point_to_tangent_disk_form_factor()
{
double r=1;
bool succeed = true;
std::cerr << "Testing analytical point-to-disk form factor formula..." << std::endl;
for(int i=0;i<10;++i)
{
double x = 10*(drand48()-0.5);
double y = 10*drand48();
double ff_analytical = point_to_tangent_disk_form_factor(x,y,r);
double ff_test_MC = point_to_tangent_disk_form_factor_MC(x,y,r);
double error = 100.0*fabs(ff_analytical-ff_test_MC)/std::max(ff_test_MC,ff_analytical);
succeed = succeed && (error<0.01);
std::cerr << " x=" << number(x,8,true) << ",\t y=" << number(y,8,true) << "\t FF_analytical=" << number(ff_analytical,8,true) << ", FF_MC=" << number(ff_test_MC,8,true) << " error=" << number(error,8,true) << "% : " << ((error<0.01)?"OK":"FAILED") << std::endl;
}
return succeed;
}
void compute_pairwise_distances(double r1,double r2,int N,double& Lp_dist,double p,double& Lp_v)
{
double res1 = 0.0;
double res3 = 0.0;
double S_n = M_PI*r1*r1;
double S_np1 = M_PI*r2*r2;
// make || L_n ||_p = 1
double L_n = 1.0/pow(M_PI*S_n ,1.0/p);
double L_np1 = 1.0/pow(M_PI*S_np1,1.0/p);
for(int i=0;i& a,int N,int P,std::ostream& o,int max_size=0)
{
if(max_size == 0)
max_size = INT_MAX;
for(size_t i=0;i < max_size && i<(size_t)N;++i)
{
o << "[ " ;
for(size_t j=0;j < max_size && j<(size_t)P;++j)
{
o.width(12) ;
o << a[i+N*j] ;
o.width(1) ;
o << "," ;
}
if(max_size < P)
o << " ..." ;
o << " ]" << std::endl;
}
if(max_size < N)
o << "..." << std::endl;
return o ;
}
void test_compactness(double p_distance)
{
// This test considers a scene made of two squares connected by an edge, making an angle of pi/2
// We form a bounded sequence of functions {u_n}_{n>0} defined by
//
// u_n (x) = (\pi r^2)^{-1/p} delta(C_n)(x)
//
// ...where C_n is the circle of radius r=1/1.2^n included in the support plan of one of the two squares,
// and tangent to the common edge at the middle point. For all n, we have ||u_n||_p = 1, so the series is bounded.
//
// We compute the illumination function on the other square (i.e. T u_n) using explicit form-factor formulas
// and verify whether or not the series is Cauchy, meaning that the distance between e.g. consecutive values converges to
// 0. If not, the series is not convergent since the space of square-integrable functions over the scene is Hilbert and therefore
// complete (meaning all Cauchy series are convergent).
test_point_to_tangent_disk_form_factor();
// First, draw the shape of the transported radiant exitance on the receiver (Fig.3, middle)
std::cerr << std::endl;
std::cerr << "Generating data for Figure 3.b to file fig3b_data.gnu" << std::endl;
std::cerr << "Plot with gnuplot, using: splot 'fig3b_data.gnu' matrix w li" << std::endl;
std::ofstream of("fig3b_data.gnu");
const int N = 80; // number of points on the receiving surface
double r_n=1.0/pow(1.2,7);
double S_n=M_PI*r_n*r_n;
for(int i=0;i Lp_distances(n_layers*n_layers,0.0) ;
std::vector Lp_v(n_layers,0.0);
for(uint32_t n=2;n