htrdr

Solving radiative transfer in heterogeneous media
git clone git://git.meso-star.fr/htrdr.git
Log | Files | Refs | README | LICENSE

commit 69297c12dc4074c2d6544df1a1ca4076cce1014c
parent ebca6733feb85a75843756e2933cc7d60627c208
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 24 Jul 2018 10:27:59 +0200

Implement the htrdr_draw_radiance_sw functio

Diffstat:
Mcmake/CMakeLists.txt | 1+
Asrc/htrdr_draw_radiance_sw.c | 182+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/htrdr_solve.h | 19++++++++++++++++++-
3 files changed, 201 insertions(+), 1 deletion(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -62,6 +62,7 @@ set(HTRDR_FILES_SRC htrdr_buffer.c htrdr_camera.c htrdr_compute_radiance_sw.c + htrdr_draw_radiance_sw.c htrdr_main.c htrdr_rectangle.c htrdr_sky.c diff --git a/src/htrdr_draw_radiance_sw.c b/src/htrdr_draw_radiance_sw.c @@ -0,0 +1,182 @@ +/* Copyright (C) 2018 Université Paul Sabatier, |Meso|Star> + * + * 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 <http://www.gnu.org/licenses/>. */ + +#include "htrdr.h" +#include "htrdr_buffer.h" +#include "htrdr_camera.h" +#include "htrdr_solve.h" + +#include <star/ssp.h> + +#include <rsys/math.h> + +#define TILE_SIZE 32 /* definition in X & Y of a tile */ +STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2); + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static FINLINE uint16_t +morton2D_decode(const uint32_t u32) +{ + uint32_t x = u32 & 0x55555555; + x = (x | (x >> 1)) & 0x33333333; + x = (x | (x >> 2)) & 0x0F0F0F0F; + x = (x | (x >> 4)) & 0x00FF00FF; + x = (x | (x >> 8)) & 0x0000FFFF; + return (uint16_t)x; +} + +static res_T +draw_tile + (struct htrdr* htrdr, + const size_t tile_org[2], /* Origin of the tile in pixel space */ + const size_t tile_sz[2], /* Definition of the tile */ + const double pix_sz[2], /* Size of a pixel in the normalized image plane */ + const struct htrdr_camera* cam, + const size_t spp, /* #samples per pixel */ + struct ssp_rng* rng, + struct htrdr_buffer* buf) +{ + size_t npixels; + size_t mcode; /* Morton code of tile pixel */ + ASSERT(htrdr && tile_org && tile_sz && pix_sz && cam && spp && buf); + + /* Adjust the #pixels to process them wrt a morton order */ + npixels = round_up_pow2(MMAX(tile_sz[0], tile_sz[1])); + npixels *= npixels; + + FOR_EACH(mcode, 0, npixels) { + struct htrdr_accum* pix_accum; + size_t ipix_tile[2]; /* Pixel coord in the tile */ + size_t ipix[2]; /* Pixel coord in the buffer */ + size_t i; + + ipix_tile[0] = morton2D_decode((uint32_t)(mcode>>0)); + if(ipix_tile[0] >= tile_sz[0]) continue; /* Pixel is out of tile */ + ipix_tile[1] = morton2D_decode((uint32_t)(mcode>>1)); + if(ipix_tile[1] >= tile_sz[1]) continue; /* Pixel is out of tile */ + + /* Compute the pixel coordinate */ + ipix[0] = tile_org[0] + ipix_tile[0]; + ipix[1] = tile_org[1] + ipix_tile[1]; + + /* Fetch and reset the pixel accumulator */ + pix_accum = htrdr_buffer_at(buf, ipix[0], ipix[1]); + *pix_accum = HTRDR_ACCUM_NULL; + + FOR_EACH(i, 0, spp) { + double pix_samp[2]; + double ray_org[3]; + double ray_dir[3]; + double weight; + + /* Sample a position into the pixel, in the normalized image plane */ + pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0]; + pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1]; + + /* Generate a ray starting from the pinhole camera and passing through the + * pixel sample */ + htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir); + + /* Compute the radiance that reach the pixel through the ray */ + weight = htrdr_compute_radiance_sw + (htrdr, rng, ray_org, ray_dir, -1/*FIXME wavelength*/); + ASSERT(weight >= 0); + + /* Update the pixel accumulator */ + pix_accum->sum_weights += weight; + pix_accum->sum_weights_sqr += weight*weight; + pix_accum->nweights += 1; + } + } + + return RES_OK; +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +htrdr_draw_radiance_sw + (struct htrdr* htrdr, + const struct htrdr_camera* cam, + const size_t spp, + struct htrdr_buffer* buf) +{ + struct ssp_rng* rng; + size_t ntiles_x, ntiles_y, ntiles; + int32_t mcode; /* Morton code of the tile */ + struct htrdr_buffer_layout layout; + double pix_sz[2]; /* Pixel size in the normalized image plane */ + res_T res = RES_OK; + ASSERT(htrdr && cam && buf); + + htrdr_buffer_get_layout(buf, &layout); + ASSERT(layout.width || layout.height || layout.elmt_size); + + if(layout.elmt_size != sizeof(struct htrdr_accum) + || layout.alignment != ALIGNOF(struct htrdr_accum)) { + htrdr_log_err(htrdr, + "%s: invalid buffer layout. " + "The pixel size must be the size of an accumulator.\n", + FUNC_NAME); + res = RES_BAD_ARG; + goto error; + } + + res = ssp_rng_create(htrdr->allocator, &ssp_rng_mt19937_64, &rng); + if(res != RES_OK) { + htrdr_log_err(htrdr, "%s: could not create the RNG.\n", FUNC_NAME); + goto error; + } + + ntiles_x = (layout.width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; + ntiles_y = (layout.height+ (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; + ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y)); + ntiles *= ntiles; + + pix_sz[0] = 1.0 / (double)layout.width; + pix_sz[1] = 1.0 / (double)layout.height; + + for(mcode=0; mcode<(int32_t)ntiles; ++mcode) { + size_t tile_org[2]; + size_t tile_sz[2]; + + /* Decode the morton code to retrieve the tile index */ + tile_org[0] = morton2D_decode((uint32_t)(mcode>>0)); + if(tile_org[0] >= ntiles_x) continue; /* Skip border tile */ + tile_org[1] = morton2D_decode((uint32_t)(mcode>>1)); + if(tile_org[1] >= ntiles_y) continue; /* Skip border tile */ + + /* Define the tile origin to pixel space */ + tile_org[0] *= TILE_SIZE; + tile_org[1] *= TILE_SIZE; + + /* Compute the size of the tile clamped by the borders of the buffer */ + tile_sz[0] = MMIN(TILE_SIZE, layout.width - tile_org[0]); + tile_sz[1] = MMIN(TILE_SIZE, layout.height - tile_org[1]); + + res = draw_tile(htrdr, tile_org, tile_sz, pix_sz, cam, spp, rng, buf); + if(res != RES_OK) goto error; + } + +exit: + if(rng) SSP(rng_ref_put(rng)); + return res; +error: + goto exit; +} + diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h @@ -18,8 +18,19 @@ #include <rsys/rsys.h> +/* Monte carlo accumulator */ +struct htrdr_accum { + double sum_weights; /* Sum of Monte-Carlo weights */ + double sum_weights_sqr; /* Sum of Monte-Carlo square weights */ + size_t nweights; /* #accumlated weights */ + size_t nfailures; /* #failures */ +}; +#define HTRDR_ACCUM_NULL__ {0,0,0,0} +static const struct htrdr_accum HTRDR_ACCUM_NULL = HTRDR_ACCUM_NULL__; + /* Forward declarations */ struct htrdr; +struct htrdr_camera; struct ssp_rng; extern LOCAL_SYM res_T @@ -34,5 +45,11 @@ htrdr_compute_radiance_sw const double dir[3], const double wavelength); -#endif /* HTRDR_SOLVE_H */ +extern LOCAL_SYM res_T +htrdr_draw_radiance_sw + (struct htrdr* htrdr, + const struct htrdr_camera* cam, + const size_t spp, /* #samples per pixel, i.e. #realisations */ + struct htrdr_buffer* buf); /* Buffer of struct htrdr_accum */ +#endif /* HTRDR_SOLVE_H */