
// Copyright 2019 Shlomi Steinberg
// 

#include <string>
#include <filesystem>
#include <iostream>
#include <fstream>
#include <cassert>
#include <getopt.h>
#include <vector>
#include <optional>
#include <algorithm>
#include <execution>
#include <chrono>

#include <opencv2/core.hpp>
#include <opencv2/videoio.hpp>

#include <fenv.h>

#include "surface.hpp"
#include "write2d.hpp"

#include "random_generator.hpp"

#include "common.hpp"
#include "XYZ.hpp"
#include "glm.hpp"
#include "erf.hpp"

#include "load_csv.hpp"

#include "thread_pool.hpp"

#include "speckle_field_renderer.hpp"
#include "surface_sampler.hpp"
#include "renderer.hpp"
#include "illuminant.hpp"
#include "refractive_indices.hpp"

#include "model.hpp"
#include "light.hpp"
#include "envmap.hpp"


// Defines basic MC-simulation properties and configurations
class mc_context {
public:
	// Config
	using fp = float;
	using vertex_t = glm::vec<3, fp>;

	// Cheap hack to create different statistics for different materials. 
	static constexpr size_t materials = 1;

	static constexpr size_t M = 1024;
    static constexpr size_t supersampling = 1;
	static constexpr auto incident_eta = std::complex<fp>{ 1,0 };
    static constexpr auto scale = 1e+3f * 1500.f;     // Model scale to um

    static constexpr bool draw_phases_with_toletance = true;
	// /Config
};


template <bool white>
using renderering_context = speckle_field_renderer<
		random_generator,
		mc_context,
		mc_k_correlation_sampler<mc_context>,
		white,
        // Ignore spatial phases as our K-correlation surfaces have very short coherence lengths. Other models might need to take into  
        // account the full term correctly.
        false
	>;


void print_usage() {
	std::cout << "Usage: demo <options>" << std::endl;
	std::cout << "-o <output>, --out=<output>\t\t\tSpecify output path (required)." << std::endl;
	std::cout << "-f <frames>, --frames=<frames>\t\t\tFrame count (required)" << std::endl;
	std::cout << "-a <fps>, --fps=<frames>\t\t\tFrames/second (defaults to 24)" << std::endl;
	std::cout << "-r <resolution>, --resolution=<resolution>\tOutput resolution (defaults to 1024x768)." << std::endl;
	std::cout << "-m <scene>, --scene=<scene>\t\t0 - Bunny, 1 - Dragon, 2 - Buddha, 3 - Teapot (defaults to 0)." << std::endl;
	std::cout << "-v <map,rotation,render_size,render_dy>, --envmap=<map,rotation,render_size,render_dy> Loads an environment map (defaults to no envmap). Rotation can be optionally specified to rotate the map around the y axis." << std::endl;
	std::cout << "-s <speed>, --speed=<speed>\t\t\tFraction of pi the light rotates around the model over the entire video (defaults to 0.25)." << std::endl;
	std::cout << "-b <extend>, --light_extend=<extend>\t\t\tLight source extend in millimetre (defaults to 8mm)." << std::endl;
	std::cout << "-d <distance>, --light_distance=<distance>\t\tLight source distance in millimetre (defaults to 1000mm)." << std::endl;
	std::cout << "-l <wavelength>, --lambda=<wavelength>\t\tLight wavelength, in nanometre. When unspecified assumes a white source." << std::endl;
	std::cout << "-t <threads>, --threads=<threads>\t\tThreads count (defaults to hardware concurrency)." << std::endl;
	std::cout << "-z <name>, --material=<name>\t\tMaterial name (defaults to 'Al', Aluminum)." << std::endl;
	std::cout << "-i <spectrum>, --illuminant=<spectrum>\t\tIlluminant spectrum (defaults to D65)." << std::endl;
	std::cout << "-e <intensity>, --intensity=<intensity>\t\tIlluminant intensity (defaults to 2.0)." << std::endl;
	std::cout << "-w <magnitude,glints>, --artistic=<magnitude,glints>\tArtistic mode speckle control." << std::endl;
	std::cout << "-h <blur>, --blurenv=<blur>\t\t\tBlur the rendered envmap, in radians." << std::endl;
	std::cout << "-q\t\t\t\t\t\tUse lossy x264 codec (defaults to Huffman Lossless Codec otherwise)." << std::endl;
	std::cout << "-x\t\t\t\t\t\tDo not render diffrations, just output mean lobe." << std::endl;

	std::cout << std::endl << std::endl;
	std::cout << "Available materials: " << std::endl;
	std::cout<< "'Al' -- Aluminum" << std::endl;
	std::cout<< "'Ag' -- Silver" << std::endl;
	std::cout<< "'Au' -- Gold" << std::endl;
	std::cout<< "'Cu' -- Copper" << std::endl;
	std::cout<< "'Co' -- Cobalt" << std::endl;
	std::cout<< "'Cr' -- Chromium" << std::endl;
	std::cout<< "'Fe' -- Iron" << std::endl;
	std::cout<< "'Fe3O4' -- Magnetite (Iron oxide)" << std::endl;
	std::cout<< "'Ir' -- Iridium" << std::endl;
	std::cout<< "'Mg' -- Magnesium" << std::endl;
	std::cout<< "'Ni' -- Nickel" << std::endl;
	std::cout<< "'Pb' -- Lead" << std::endl;
	std::cout<< "'Pt' -- Platinum" << std::endl;
	std::cout<< "'Ta' -- Tantalum" << std::endl;
	std::cout<< "'Ti' -- Titanium" << std::endl;
	std::cout<< "'TiN' -- Titanium Nitride" << std::endl;
	std::cout<< "'W'  -- Tungsten" << std::endl;

	std::cout << std::endl << std::endl;
	std::cout << "Available illuminants: " << std::endl;
	std::cout<< "'D50' -- 5000K white" << std::endl;
	std::cout<< "'D65' -- 6500K white" << std::endl;
	std::cout<< "'D75' -- 7500K cold" << std::endl;
	std::cout<< "'F1'  -- F1 white fluorescent" << std::endl;
	std::cout<< "'F2'  -- F2 warm fluorescent" << std::endl;
	std::cout<< "'F3'  -- F3 warm fluorescent" << std::endl;
	std::cout<< "'F9'  -- F9 warm fluorescent" << std::endl;
}

template <typename FB>
void write_heatmap(const std::filesystem::path &path, std::size_t res_x, std::size_t res_y, const FB& fb) {
	// Load heatmap table
	std::vector<glm::vec3> colormap;
	std::filesystem::path map_path = "../data/cividis.csv";
	{
		std::ifstream ifs(map_path);
		if (!ifs) {
			std::cerr << "'" << map_path << "' not found!" << std::endl;
			return;
		}
		std::string dummy;
		std::getline(ifs, dummy);
		std::vector<double> data;
		colormap.reserve(300);
		while ((data = load_csv::read_csv_line(ifs)).size())
			colormap.emplace_back(glm::vec3{ static_cast<float>(data[0]), static_cast<float>(data[1]), static_cast<float>(data[2]) });
	}

	std::size_t max{}, min{ std::numeric_limits<std::size_t>::max() };
	for (std::size_t y=0;y<res_y;++y)
	for (std::size_t x=0;x<res_x;++x) {
		if (fb(x,y).c.x<=.0f) continue;
		max = glm::max(max, static_cast<std::size_t>(fb(x,y).c.x));
		min = glm::min(min, static_cast<std::size_t>(fb(x,y).c.x));
	}
	min = glm::min(min,max);

	::surface<> s = { std::make_unique<std::uint8_t[]>(res_x*res_y*3), res_x, res_y, 3 };
	for (std::size_t y=0;y<res_y;++y)
	for (std::size_t x=0;x<res_x;++x) {
		const auto v = fb(x,y).c.x;
		const auto f = (v-min)/static_cast<float>(max-min)*(colormap.size()-1);
		const auto l = glm::clamp<std::size_t>(static_cast<std::size_t>(f), 0, colormap.size()-1);
		const auto h = glm::clamp<std::size_t>(l+1, 0, colormap.size()-1);
		const auto c = v<=0 ? glm::vec3{ .0f } : glm::mix(colormap[l], colormap[h], glm::fract(f));

		const auto packed = glm::packUnorm4x8(glm::vec4{ glm::pow(c, glm::vec3{ 1.f/2.2f }), .0f });
		s.data[(x+y*res_x)*3 + 0] = static_cast<std::uint8_t>((packed >>  0) & 0xFF);
		s.data[(x+y*res_x)*3 + 1] = static_cast<std::uint8_t>((packed >>  8) & 0xFF);
		s.data[(x+y*res_x)*3 + 2] = static_cast<std::uint8_t>((packed >> 16) & 0xFF);
	}

	write_2d(path, s);
}

template <typename R>
void render_video(const std::filesystem::path &outpath, 
                  const std::filesystem::path &outpath_speckle1,
                  const std::filesystem::path &outpath_speckle2,
                  const std::filesystem::path &outpath_mean, 
				  const bool render_diffractions,
                  const renderer::model_t &model,
                  const glm::mat4 &modelmat, const glm::mat4 &viewmat,
                  const std::vector<light> &lights, 
				  const envmap *env,
				  bool lossy,
				  std::size_t res_x, std::size_t res_y, 
				  std::optional<std::size_t> threads, 
				  std::size_t frame_count, 
				  std::size_t fps,
				  R renderer, 
				  thread_pool &tp) {
	print_progress(.0f);

	// Enqueue all frames
	using future_t = std::future<std::unique_ptr<renderer::render_output_t>>;
	std::vector<future_t> futures;
	futures.reserve(frame_count);
	for (std::size_t idx=0; idx<frame_count; ++idx)
		futures.emplace_back(tp.enqueue([idx, res_x, res_y, &renderer, &model, &modelmat, &viewmat, &lights, env, frame_count]() {
			const auto t = (idx+.5f)/static_cast<float>(frame_count);
			std::vector<float> coherence_statistics_per_light; coherence_statistics_per_light.resize(lights.size(), 0.f);
			return ::renderer::template do_render(res_x, res_y, renderer, 
												  model, modelmat, viewmat, lights, env, mc_context::materials,
												  t, mc_context::scale);
		}));

	// Write video out frame-by-frame
    const auto fourcc = lossy ? cv::VideoWriter::fourcc('X','2','6','4') : cv::VideoWriter::fourcc('H','F','Y','U');
    cv::VideoWriter outputVideo0, outputVideo1, outputVideo2, outputVideo3;
	if (render_diffractions) {
		outputVideo0.open(outpath.string().c_str(), fourcc, fps, cv::Size(res_x, res_y), true);
		outputVideo1.open(outpath_speckle1.string().c_str(), fourcc, fps, cv::Size(res_x, res_y), true);
		outputVideo2.open(outpath_speckle2.string().c_str(), fourcc, fps, cv::Size(res_x, res_y), true);
	}
    outputVideo3.open(outpath_mean.string().c_str(), fourcc, fps, cv::Size(res_x, res_y), true);
	std::size_t prog = 0;
	for (auto&& frame_future : futures) {
		while (frame_future.wait_for(std::chrono::seconds(1)) != std::future_status::ready)
			print_progress((prog+1e-10f)/static_cast<float>(frame_count));
		const auto f = frame_future.get();
		
		cv::Mat frame0(res_y, res_x, CV_8UC3, cv::Scalar(0, 0, 0));
		cv::Mat frame1(res_y, res_x, CV_8UC3, cv::Scalar(0, 0, 0));
		cv::Mat frame2(res_y, res_x, CV_8UC3, cv::Scalar(0, 0, 0));
		cv::Mat frame3(res_y, res_x, CV_8UC3, cv::Scalar(0, 0, 0));
		for (std::size_t y=0; y<res_y; ++y) {
			auto pixels0 = reinterpret_cast<std::uint8_t*>(frame0.data) + frame0.step1()*y;
			auto pixels1 = reinterpret_cast<std::uint8_t*>(frame1.data) + frame1.step1()*y;
			auto pixels2 = reinterpret_cast<std::uint8_t*>(frame2.data) + frame2.step1()*y;
			auto pixels3 = reinterpret_cast<std::uint8_t*>(frame3.data) + frame3.step1()*y;
			for (std::size_t x=0; x<res_x; ++x) {
				// Gamma correct and pack
				const auto packed0 = glm::packUnorm4x8(glm::sqrt(glm::max(glm::vec4{ .0f }, (*f)(x,y).c0 * 0.9f)));
				const auto packed1 = glm::packUnorm4x8(glm::sqrt(glm::max(glm::vec4{ .0f }, (*f)(x,y).c1 * 0.9f)));
				const auto packed2 = glm::packUnorm4x8(glm::sqrt(glm::max(glm::vec4{ .0f },-(*f)(x,y).c1 * 0.9f)));
				const auto packed3 = glm::packUnorm4x8(glm::sqrt(glm::max(glm::vec4{ .0f }, (*f)(x,y).c2 * 0.9f)));
				pixels0[x*3 + 2] = static_cast<std::uint8_t>((packed0 >>  0) & 0xFF);
				pixels0[x*3 + 1] = static_cast<std::uint8_t>((packed0 >>  8) & 0xFF);
				pixels0[x*3 + 0] = static_cast<std::uint8_t>((packed0 >> 16) & 0xFF);
				pixels1[x*3 + 2] = static_cast<std::uint8_t>((packed1 >>  0) & 0xFF);
				pixels1[x*3 + 1] = static_cast<std::uint8_t>((packed1 >>  8) & 0xFF);
				pixels1[x*3 + 0] = static_cast<std::uint8_t>((packed1 >> 16) & 0xFF);
				pixels2[x*3 + 2] = static_cast<std::uint8_t>((packed2 >>  0) & 0xFF);
				pixels2[x*3 + 1] = static_cast<std::uint8_t>((packed2 >>  8) & 0xFF);
				pixels2[x*3 + 0] = static_cast<std::uint8_t>((packed2 >> 16) & 0xFF);
				pixels3[x*3 + 2] = static_cast<std::uint8_t>((packed3 >>  0) & 0xFF);
				pixels3[x*3 + 1] = static_cast<std::uint8_t>((packed3 >>  8) & 0xFF);
				pixels3[x*3 + 0] = static_cast<std::uint8_t>((packed3 >> 16) & 0xFF);
			}
		}

		if (render_diffractions) {
			outputVideo0 << frame0;
			outputVideo1 << frame1;
			outputVideo2 << frame2;
		}
		outputVideo3 << frame3;

		print_progress((++prog+1e-10f)/static_cast<float>(frame_count));
	}

	print_progress(1.f);
}

auto random_vec() {
	using fp = typename mc_context::fp;
    using Rand = random_generator;
    return glm::vec3{ std::normal_distribution<fp>{ 0,10.f }(Rand{}()), std::normal_distribution<fp>{ 0,10.f }(Rand{}()), std::normal_distribution<fp>{ 0,10.f }(Rand{}()) };
}
int main(int argc, char* argv[]) {
	feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);

	using fp = typename mc_context::fp;

	// Output opts
	std::filesystem::path outpath{}, envmap_path;
	std::size_t res_x=1024, res_y=768, frames_count=0, fps = 24, scene = 0;
	std::optional<size_t> lambda, threads;
	std::string material_name{}, illuminant_name{};
	float light_distance_mm = 1000.f, light_extend_mm = 8.f, speed = 0.25f, light_intensity = 2.f;
	float envmap_rot = .0f, blurenvmap = .0f, envmap_dy = .0f, envmap_size = 1.f;
	bool lossy_codec = false;
    bool artistic = false;
	bool render_diffractions = true;
    fp artistic_glints, artistic_magnitude;

	// Parse intput
	{
		option long_options[] = {
			{"out", required_argument, nullptr, 'o'},
			{"resolution", required_argument, nullptr, 'r'},
			{"lambda", required_argument, nullptr, 'l'},
			{"threads", required_argument, nullptr, 't'},
			{"frames", required_argument, nullptr, 'f'},
			{"fps", required_argument, nullptr, 'a'},
			{"light_extend", required_argument, nullptr, 'b'},
			{"light_distance", required_argument, nullptr, 'd'},
			{"speed", required_argument, nullptr, 's'},
			{"scene", required_argument, nullptr, 'm'},
			{"material", required_argument, nullptr, 'z'},
			{"envmap", required_argument, nullptr, 'v'},
			{"illuminant", required_argument, nullptr, 'i'},
			{"intensity", required_argument, nullptr, 'e'},
			{"artistic", required_argument, nullptr, 'w'},
			{"blurenv", required_argument, nullptr, 'h'},
			{"", no_argument, nullptr, 'q'},
			{"", no_argument, nullptr, 'x'},
			{nullptr, 0, nullptr, 0}
		};
		char ch;
		while ((ch = getopt_long(argc, argv, "v:w:z:i:e:o:r:a:m:l:t:f:s:b:d:h:qx", long_options, nullptr)) != -1) {
			switch (ch) {
			case 'o':
				outpath = optarg;
				break;
			case 'r': {
				res_x = res_y = atoi(optarg);
				const auto *x = strstr(optarg,"x");
				if (x)
					res_y = atoi(++x);
				break;
			}
			case 'l':
				lambda = atoi(optarg);
				break;
			case 't':
				threads = atoi(optarg);
				break;
			case 'f':
				frames_count = atoi(optarg);
				break;
			case 'a':
				fps = atoi(optarg);
				break;
			case 'm':
				scene = atoi(optarg);
				break;
			case 's':
				speed = atof(optarg);
				break;
			case 'h':
				blurenvmap = atof(optarg);
				break;
			case 'b':
				light_extend_mm = atof(optarg);
				break;
			case 'd':
				light_distance_mm = atof(optarg);
				break;
			case 'q':
				lossy_codec = true;
				break;
			case 'x':
				render_diffractions = false;
				break;
			case 'z':
				material_name = optarg;
				break;
			case 'i':
				illuminant_name = optarg;
				break;
			case 'e':
				light_intensity = atof(optarg);
				break;
            case 'v': {
                auto nextarg = std::strstr(optarg,",");
				envmap_path = nextarg ? std::string(optarg, nextarg) : std::string(optarg);
				if (!nextarg) break;
				envmap_rot = atof(nextarg+1);
                nextarg = std::strstr(nextarg+1,",");
				if (!nextarg) break;
				envmap_size = atof(nextarg+1);
                nextarg = std::strstr(nextarg+1,",");
				if (!nextarg) break;
				envmap_dy = atof(nextarg+1);
                break;
            }
            case 'w': {
                artistic = true;
                artistic_magnitude = atof(optarg);
                const auto nextarg = std::strstr(optarg,",");
                if (nextarg==nullptr) {
                    print_usage();
                    return 1;
                }
                artistic_glints = atof(nextarg+1);
                break;
            }
			default:
				std::cerr << "Unknown option!" << std::endl << std::endl;
				print_usage();
				return 1;
			}
		}

		if (!outpath.string().length() || res_x<=0 || res_y<=0 || threads.value_or(1)<1 || frames_count<1 || scene>3) {
			print_usage();
			return 1;
		}
	}

    const auto extension = outpath.extension().string();
    auto speckle_path1 = outpath, speckle_path2 = outpath, mean_path = outpath;
    const auto filename = speckle_path1.replace_extension({}).filename().string();
    speckle_path1.replace_filename({});
    speckle_path2.replace_filename({});
    mean_path.replace_filename({});
    speckle_path1 = speckle_path1 / (filename + "_speckle_p" + extension);
    speckle_path2 = speckle_path2 / (filename + "_speckle_m" + extension);
    mean_path = mean_path / (filename + "_mean" + extension);

	using renderer_white = renderering_context<true>;
	using renderer_single_wavelength = renderering_context<false>;

    // Load model
    const auto model_path = scene == 0 ? "../data/bunny/bun_zipper.ply" : 
							scene == 1 ? "../data/dragon_recon/dragon_vrip_res2.ply" :
							scene == 2 ? "../data/happy_recon/happy_vrip_res2.ply" :
							"../data/teapot.ply";
	const auto model_rescale = scene == 3 ? 19e-2f : 1.f;
    const auto model = renderer::load_model(model_path, model_rescale);
    {
        decltype(model.verts[0].x) minx{},maxx{},miny{},maxy{};
        for (const auto &v : model.verts) {
            minx = glm::min(minx,v.x);
            maxx = glm::max(maxx,v.x);
            miny = glm::min(miny,v.y);
            maxy = glm::max(maxy,v.y);
        }
        std::cout << "Model: " << model_path << std::endl;
        std::cout << "Model size (W x H):  " << std::to_string((maxx-minx)*mc_context::scale*1e-3f) << "mm x " 
                                             << std::to_string((maxy-miny)*mc_context::scale*1e-3f) << "mm" << std::endl; 
    }

    // Scene parameters
    auto modelmat=glm::mat4(1);
    auto viewmat=glm::mat4(1);
    auto light_vec = glm::vec3{ 15,-20,-80 };
    auto light_rotate_vec = glm::normalize(glm::vec3{ 2,-5,-3 });
    auto light_rotate_start = -0.025*glm::pi<float>();
    {
        modelmat[1][1]=-1.0;
        if (scene == 0) {
            modelmat[0][0]*=12.5 * res_y/float(res_x);
            modelmat[1][1]*=12.5;
            modelmat[2][2]*=12.5;
            modelmat[0][3]=0.25;
            modelmat[1][3]=1.4;
            modelmat[2][3]=0.1;
            const auto s=glm::sin(-0.1f);
            const auto c=glm::cos(-0.1f);
            viewmat[0][0]=viewmat[2][2]=c;
            viewmat[0][2]=s;
            viewmat[2][0]=-s;
            viewmat[3][3]=viewmat[1][1]=1.0;
            light_vec = glm::vec3{ -30,-50,-80 };
            light_rotate_vec = glm::normalize(glm::vec3{ 20,-5,10 });
            light_rotate_start = 0.07*glm::pi<float>();
        }
        if (scene == 1) {
            modelmat[0][0]*=9;
            modelmat[1][1]*=9;
            modelmat[2][2]*=9;
            modelmat[0][3]=-0.15;
            modelmat[1][3]=1;
            modelmat[2][3]=0.4;
            viewmat = glm::rotate(glm::mat4(1), -glm::pi<float>()/8, glm::normalize(glm::vec3{ -4,5,0 }));
            light_vec = glm::vec3{ 20,-175,-80 };
        }
        if (scene == 2) {
            modelmat[0][0]*=10 * res_y/float(res_x);
            modelmat[1][1]*=10;
            modelmat[2][2]*=10;
            modelmat[0][3]=0;
            modelmat[1][3]=1.5;
            modelmat[2][3]=0.4;
            viewmat = glm::rotate(glm::mat4(1), -glm::pi<float>()/8, glm::normalize(glm::vec3{ 1,5,0 }));
        }
		if (scene == 3) {
			modelmat[2][2]*=1.7;
			modelmat[2][3]=-0.5;
            viewmat = glm::rotate(glm::mat4(1), glm::pi<float>()/2, glm::normalize(glm::vec3{ -1,0.01,0.1 }));
            light_vec = glm::vec3{ 20,40,20 };
            light_rotate_vec = glm::normalize(glm::vec3{ -20,-5,10 });
            light_rotate_start = 0.07*glm::pi<float>();
		}
    }
    auto light_rotate_end = light_rotate_start + speed * glm::pi<float>();

	// Material
	std::array<std::array<std::complex<fp>, 8>, mc_context::materials> materials;
	materials[0] = ior::template aluminum<fp>;
	if (material_name.length()) {
		if (material_name == "Al")			materials[0] = ior::template aluminum<fp>;
		else if (material_name == "Ag")		materials[0] = ior::template silver<fp>;
		else if (material_name == "Au")		materials[0] = ior::template gold<fp>;
		else if (material_name == "Cu")		materials[0] = ior::template copper<fp>;
		else if (material_name == "Co")		materials[0] = ior::template cobalt<fp>;
		else if (material_name == "Cr")		materials[0] = ior::template chromium<fp>;
		else if (material_name == "Fe")		materials[0] = ior::template iron<fp>;
		else if (material_name == "Ir")		materials[0] = ior::template iridium<fp>;
		else if (material_name == "Mg")		materials[0] = ior::template magnesium<fp>;
		else if (material_name == "Ni")		materials[0] = ior::template nickel<fp>;
		else if (material_name == "Fe3O4")	materials[0] = ior::template magnetite<fp>;
		else if (material_name == "Pb")		materials[0] = ior::template lead<fp>;
		else if (material_name == "Pt")		materials[0] = ior::template platinum<fp>;
		else if (material_name == "Ta")		materials[0] = ior::template tantalum<fp>;
		else if (material_name == "Ti")		materials[0] = ior::template titanium<fp>;
		else if (material_name == "TiN")	materials[0] = ior::template titanium_nitride<fp>;
		else if (material_name == "W")		materials[0] = ior::template tungsten<fp>;
		else {
			std::cerr << "Unknown material!" << std::endl;
			print_usage();
			return 1;
		}
	}
	if constexpr (mc_context::materials>1) {
		// for Buddha
		materials[1] = ior::template gold<fp>;
	}

	// Illuminants
	std::vector<light> lights;
	{
		std::unique_ptr<illuminant_base> illuminant = std::make_unique<illuminant_d65>();
		if (illuminant_name.length()) {
			if (illuminant_name == "D65")		illuminant = std::make_unique<illuminant_d65>();
			else if (illuminant_name == "D50")	illuminant = std::make_unique<illuminant_d50>();
			else if (illuminant_name == "D75")	illuminant = std::make_unique<illuminant_d75>();
			else if (illuminant_name == "F1")	illuminant = std::make_unique<illuminant_f1>();
			else if (illuminant_name == "F2")	illuminant = std::make_unique<illuminant_f2>();
			else if (illuminant_name == "F3")	illuminant = std::make_unique<illuminant_f3>();
			else if (illuminant_name == "F9")	illuminant = std::make_unique<illuminant_f9>();
			else {
				std::cerr << "Unknown illuminant!" << std::endl;
				print_usage();
				return 1;
			}
		}
		lights.emplace_back(std::move(illuminant),
						    glm::normalize(light_vec), 
							light_distance_mm*1e+3f, light_distance_mm*1e+3f,
							light_rotate_vec, light_rotate_start, light_rotate_end, light_intensity, 
							light_extend_mm*1e+3f);
	}
	if constexpr (mc_context::materials>1) {
		// for Buddha
		light_vec = glm::vec3{ -50,2,-30 };
		light_rotate_vec = glm::normalize(glm::vec3{ 0,1,0 });
		light_extend_mm = 1000;
		light_distance_mm = 1000;
		lights.emplace_back(std::make_unique<illuminant_d65>(),
						    glm::normalize(light_vec), 
							light_distance_mm*1e+3f, 75*light_distance_mm*1e+3f,
							light_rotate_vec, fp(0), fp(0), fp(1.35), 
							light_extend_mm*1e+3f,true);
	}

	// Env map
	std::vector<envmap> envmaps;
	if (envmap_path.string().length()) {
		std::cout << "Generating env maps..." << std::endl;
		for (std::size_t m=0; m<mc_context::materials; ++m) {
			const auto blur_angle = glm::pi<float>() / glm::max(1.f, 22.f*glm::pow(mc_k_correlation_sampler<mc_context>::l_cor(m),1.30f));
			envmaps.emplace_back(envmap_path, blur_angle, blurenvmap, envmap_dy, envmap_size, res_x/(float)res_y, envmap_rot);
		}
	}
	const auto *env = envmaps.size() ? envmaps.data() : nullptr;

	// Start thread pool
	thread_pool tp(threads.value_or(std::thread::hardware_concurrency()));

	// Do very fancy renderings
	{
        std::cout << "Output: " << outpath << std::endl;
        std::cout << "Output (speckle): " << speckle_path1 << " & " << speckle_path2 << std::endl;
        std::cout << "Output (mean): " << mean_path << std::endl;
		std::cout << "Rendering " << std::to_string(res_x) << "x" << std::to_string(res_y) << " (x" << std::to_string(sqr(mc_context::supersampling)) << " SS), " << std::to_string(frames_count) 
				  << " frames, " << std::to_string(mc_context::M) << " MC samples."
				  << " Light rotation of " << std::to_string(speed) << "*pi (" 
				  << std::to_string(speed/frames_count*fps*180.f/glm::pi<float>()) << "degrees/sec)..." << std::endl;
        if (artistic) {
            std::cout << "Artistic speckle control (" << std::to_string(artistic_magnitude) << ", " << std::to_string(artistic_glints) << ")" << std::endl;
        }
        else {
            const auto cr = lights[0].coherence_radius(lights[0].solid_angle(0,{0,0,0}));
            std::cout << "Average coherence radius at normal incidence: ~" << std::to_string(cr) << "um" << std::endl;
        }

		lambda ? 
			render_video(outpath, speckle_path1, speckle_path2, mean_path, render_diffractions, 
                         model, modelmat, viewmat, lights, env,
                         lossy_codec, res_x, res_y, 
						 threads, frames_count, fps,
						 renderer_single_wavelength{ render_diffractions, res_x*res_y, std::move(materials), artistic, artistic_glints, artistic_magnitude, lambda.value() }, tp)
			: 
			render_video(outpath, speckle_path1, speckle_path2, mean_path, render_diffractions,
                         model, modelmat, viewmat, lights, env,
                         lossy_codec, res_x, res_y,
						 threads, frames_count, fps,
						 renderer_white{ render_diffractions, res_x*res_y, std::move(materials), artistic, artistic_glints, artistic_magnitude }, tp);
	}

	std::cout << std::endl;
	return 0;
}
