Skip to content

Commit

Permalink
WIP -- MLT: found the heatmap bug
Browse files Browse the repository at this point in the history
  • Loading branch information
w3ntao committed Oct 14, 2024
1 parent 997df6e commit 1eb93b9
Show file tree
Hide file tree
Showing 5 changed files with 63 additions and 50 deletions.
2 changes: 1 addition & 1 deletion src/pbrt/base/sampler.cu
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,7 @@ Point2f Sampler::get_2d() {
}

case (Type::mlt): {
auto converted_ptr = ((MLTSampler *)ptr);
auto converted_ptr = (MLTSampler *)ptr;
return Point2f(converted_ptr->next_sample(), converted_ptr->next_sample());
}

Expand Down
25 changes: 16 additions & 9 deletions src/pbrt/films/grey_scale_film.cu
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,18 @@
#include "pbrt/spectrum_util/color_encoding.h"
#include "pbrt/spectrum_util/rgb.h"
#include "pbrt/util/basic_math.h"
#include <algorithm>

/*
const std::vector<RGB> colors = {
RGB(0, 0, 1), RGB(0, 1, 1), RGB(0, 1, 0), RGB(1, 1, 0), RGB(1, 0, 0),
RGB(0, 0, 0), RGB(0, 0, 1), RGB(0, 1, 1), RGB(0, 1, 0), RGB(1, 1, 0), RGB(1, 0, 0),
};
const std::vector<RGB> colors = {RGB(0, 0, 0), RGB(0, 0, 1), RGB(0, 1, 1), RGB(0, 1, 0),
RGB(1, 1, 0), RGB(1, 0, 0), RGB(1, 1, 1)};
const std::vector<RGB> colors = {
RGB(68, 1, 84) / 255, RGB(59, 82, 139) / 255, RGB(33, 145, 140) / 255,
RGB(94, 201, 98) / 255, RGB(253, 231, 37) / 255,
Expand All @@ -19,22 +23,24 @@ const std::vector<RGB> colors = {

const std::vector<RGB> colors = {RGB(68, 1, 84) / 255, RGB(68, 57, 131) / 255,
RGB(49, 104, 142) / 255, RGB(33, 145, 140) / 255,
RGB(53, 183, 121) / 255, RGB(144, 215, 67) / 255,
RGB(53, 183, 121) / 255, RGB(144, 215, 67) / 251,
RGB(253, 231, 37) / 255};

/*
taken from https://www.andrewnoske.com/wiki/Code_-_heatmaps_and_color_gradients
6 colors rainbow: black, blue, cyan, green, yellow, red
5 colors rainbow: blue, cyan, green, yellow, red
7 colors rainbow: black, blue, cyan, green, yellow, red, white
7 color viridis: https://waldyrious.net/viridis-palette-generator/
*/

static RGB convert_to_heatmap_rgb(double linear) {
// linear = linear * 2;

linear = clamp<double>(linear, 0, 1);

/*
const auto gamma = 2.0;
linear = pow(linear, 1.0 / gamma);
*/

const auto gap = 1.0 / (colors.size() - 1);

Expand All @@ -48,11 +54,12 @@ static RGB convert_to_heatmap_rgb(double linear) {
return colors[colors.size() - 1];
}

void GreyScaleFilm::write_to_png(const std::string &filename) {
FloatType max_intensity = 0.0;
for (float pixel : pixels) {
max_intensity = std::max(max_intensity, pixel);
}
void GreyScaleFilm::write_to_png(const std::string &filename) const {
auto sorted_value = pixels;
std::sort(sorted_value.begin(), sorted_value.end(),std::greater{});

const auto ratio = 0.01;
const double max_intensity = sorted_value[pixels.size() * ratio];

SRGBColorEncoding srgb_encoding;

Expand Down
4 changes: 2 additions & 2 deletions src/pbrt/films/grey_scale_film.h
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,10 @@ class GreyScaleFilm {
pixels(std::vector<FloatType>(_resolution.x * _resolution.y, 0)) {}

void add_sample(const Point2i &coord, FloatType val) {
pixels[coord.x + resolution.x * coord.y] += val;
pixels[coord.y * resolution.x + coord.x] += val;
}

void write_to_png(const std::string &filename);
void write_to_png(const std::string &filename) const;

private:
Point2i resolution;
Expand Down
74 changes: 38 additions & 36 deletions src/pbrt/integrators/mlt_path.cu
Original file line number Diff line number Diff line change
Expand Up @@ -10,17 +10,18 @@
#include "pbrt/scene/parameter_dictionary.h"
#include "pbrt/spectrum_util/global_spectra.h"

const size_t NUM_MLT_SAMPLERS = 1024 * 4;
// for unknown reason bugs triggered when NUM_MLT_SAMPLERS set too large (> 1024 * 4)
const size_t NUM_MLT_SAMPLERS = 1024 * 256;

struct BatchSample {
struct MLTSample {
PathSample path_sample;
FloatType weight;
FloatType sampling_density;

PBRT_GPU
BatchSample(FloatType _x, FloatType _y, const SampledSpectrum &_radiance,
const SampledWavelengths &_lambda, FloatType _weight)
: path_sample(PathSample(_x, _y, _radiance, _lambda)), weight(_weight) {}
MLTSample(FloatType _x, FloatType _y, const SampledSpectrum &_radiance,
const SampledWavelengths &_lambda, FloatType _weight, FloatType _sampling_density)
: path_sample(PathSample(_x, _y, _radiance, _lambda)), weight(_weight),
sampling_density(_sampling_density) {}
};

__global__ void build_seed_path(PathSample *initial_paths, PathSample *seed_path_candidates,
Expand Down Expand Up @@ -79,7 +80,7 @@ __global__ void build_seed_path(PathSample *initial_paths, PathSample *seed_path
seed_path_candidates[worker_idx * num_seed_paths_per_worker + selected_path_id];
}

__global__ void wavefront_render(BatchSample *batch_samples, PathSample *path_samples,
__global__ void wavefront_render(MLTSample *mlt_samples, PathSample *path_samples,
MLTPathIntegrator *mlt_integrator, RNG *rngs, const Filter *filter,
double brightness) {
const uint worker_idx = blockIdx.x * blockDim.x + threadIdx.x;
Expand Down Expand Up @@ -107,17 +108,18 @@ __global__ void wavefront_render(BatchSample *batch_samples, PathSample *path_sa
mlt_integrator->compute_luminance(new_path.radiance, new_path.lambda);
const auto path_luminance = mlt_integrator->compute_luminance(path.radiance, path.lambda);

double accept_probability = std::min<FloatType>(1.0, new_path_luminance / path_luminance);
FloatType accept_probability = std::min<FloatType>(1.0, new_path_luminance / path_luminance);

const double new_path_weight = (accept_probability + mlt_sampler->large_step) /
(new_path_luminance / brightness + p_large);
const double path_weight = (1.0 - accept_probability) / (path_luminance / brightness + p_large);
const FloatType new_path_weight = (accept_probability + mlt_sampler->large_step) /
(new_path_luminance / brightness + p_large);
const FloatType path_weight =
(1.0 - accept_probability) / (path_luminance / brightness + p_large);

batch_samples[offset + 0] =
BatchSample(path.x, path.y, path.radiance, path.lambda, path_weight);
mlt_samples[offset + 0] =
MLTSample(path.x, path.y, path.radiance, path.lambda, path_weight, accept_probability);

batch_samples[offset + 1] =
BatchSample(new_path.x, new_path.y, new_path.radiance, new_path.lambda, new_path_weight);
mlt_samples[offset + 1] = MLTSample(new_path.x, new_path.y, new_path.radiance, new_path.lambda,
new_path_weight, 1 - accept_probability);

if (rng->uniform<FloatType>() < accept_probability) {
path = new_path;
Expand Down Expand Up @@ -172,17 +174,12 @@ MLTPathIntegrator *MLTPathIntegrator::create(std::optional<int> samples_per_pixe
PBRT_GPU
PathSample MLTPathIntegrator::generate_new_path(Sampler *sampler, const Filter *filter) const {
// TODO: move filter into IntegratorBase
auto pixel_pos = sampler->get_2d();
auto u = sampler->get_2d();

auto pixel_x = pixel_pos.x * film_dimension.x;
auto pixel_y = pixel_pos.y * film_dimension.y;
auto pixel_x = u.x * film_dimension.x;
auto pixel_y = u.y * film_dimension.y;

auto pixel_coord = Point2i(std::floor(pixel_x), std::floor(pixel_y));
auto pixel_offset = Point2f(pixel_x - pixel_coord.x, pixel_y - pixel_coord.y);

auto fs = filter->sample(pixel_offset);

auto camera_sample = CameraSample(Point2f(pixel_x, pixel_y) + Vector2f(0.5, 0.5), fs.weight);
auto camera_sample = CameraSample(Point2f(pixel_x, pixel_y), 1);

auto lu = sampler->get_1d();
auto lambda = SampledWavelengths::sample_visible(lu);
Expand All @@ -191,7 +188,7 @@ PathSample MLTPathIntegrator::generate_new_path(Sampler *sampler, const Filter *

auto radiance = ray.weight * PathIntegrator::eval_li(ray.ray, lambda, base, sampler, max_depth);

return PathSample(pixel_pos.x, pixel_pos.y, radiance, lambda);
return PathSample(pixel_x, pixel_y, radiance, lambda);
}

void MLTPathIntegrator::render(Film *film, GreyScaleFilm &heat_map, const Filter *filter) {
Expand Down Expand Up @@ -234,38 +231,43 @@ void MLTPathIntegrator::render(Film *film, GreyScaleFilm &heat_map, const Filter
CHECK_CUDA_ERROR(cudaFree(ptr));
}

BatchSample *batch_samples;
CHECK_CUDA_ERROR(cudaMallocManaged(&batch_samples, sizeof(BatchSample) * 2 * NUM_MLT_SAMPLERS));
MLTSample *mlt_samples;
CHECK_CUDA_ERROR(cudaMallocManaged(&mlt_samples, sizeof(MLTSample) * 2 * NUM_MLT_SAMPLERS));

auto total_pass =
(long long)(mutation_per_pixel)*film_dimension.x * film_dimension.y / NUM_MLT_SAMPLERS;

auto gap = total_pass / 10;

for (uint pass = 0; pass < total_pass; ++pass) {
if (pass % 50 == 0 || pass == total_pass - 1) {
if (pass % gap == 0 || pass == total_pass - 1) {
printf("mutation: %u, pass %d/%llu\n", mutation_per_pixel, pass + 1, total_pass);
}

// TODO: render preview with OpenGL
wavefront_render<<<blocks, threads>>>(batch_samples, path_samples, this, rngs, filter,
wavefront_render<<<blocks, threads>>>(mlt_samples, path_samples, this, rngs, filter,
brightness);
CHECK_CUDA_ERROR(cudaGetLastError());
CHECK_CUDA_ERROR(cudaDeviceSynchronize());

for (uint idx = 0; idx < 2 * NUM_MLT_SAMPLERS; ++idx) {
auto path_sample = &batch_samples[idx].path_sample;
auto weight = batch_samples[idx].weight;
auto path_sample = &mlt_samples[idx].path_sample;
auto weight = mlt_samples[idx].weight;
auto sampling_density = mlt_samples[idx].sampling_density;

auto pixel_x = clamp<int>(std::floor(path_sample->x), 0, film_dimension.x - 1);
auto pixel_y = clamp<int>(std::floor(path_sample->y), 0, film_dimension.y - 1);

auto pixel_coord =
Point2i(path_sample->x * film_dimension.x, path_sample->y * film_dimension.y);
auto pixel_coord = Point2i(pixel_x, pixel_y);

film->add_sample(pixel_coord, path_sample->radiance, path_sample->lambda, weight);
// TODO: change film->add_sample() to Film::AddSplat() from PBRT-v4

heat_map.add_sample(pixel_coord, weight);
heat_map.add_sample(pixel_coord, sampling_density);
}
}

for (auto ptr : std::vector<void *>({path_samples, batch_samples, rngs})) {
for (auto ptr : std::vector<void *>({path_samples, mlt_samples, rngs})) {
CHECK_CUDA_ERROR(cudaFree(ptr));
}
}
8 changes: 6 additions & 2 deletions src/pbrt/samplers/mlt.h
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,17 @@ class MLTSampler {
constexpr FloatType s2 = 1.0 / 16.0;
const FloatType r = rng.uniform<FloatType>();
const FloatType dx = s1 / (s1 / s2 + std::abs(2.0 * r - 1.0)) - s1 / (s1 / s2 + 1.0);

FloatType mutated_val = NAN;
if (r < 0.5) {
const FloatType x1 = x + dx;
return (x1 < 1.0) ? x1 : x1 - 1.0;
mutated_val = (x1 >= 1.0) ? x1 - 1.0 : x1;
} else {
const FloatType x1 = x - dx;
return (x1 < 0.0) ? x1 + 1.0 : x1;
mutated_val = (x1 < 0.0) ? x1 + 1.0 : x1;
}

return clamp<FloatType>(mutated_val, 0, OneMinusEpsilon);
}

public:
Expand Down

0 comments on commit 1eb93b9

Please sign in to comment.