commit a55f64aa902c4c7383cbcc0046bef2d6749e3745
parent 2e5d7bbff2ddd30b25d7e46a096c625dd68c49bf
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 21 Feb 2018 16:08:39 +0100
First draft of the sdis_solve_camera function
Diffstat:
5 files changed, 390 insertions(+), 7 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -154,6 +154,14 @@ struct sdis_interface_shader {
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
+/* Functor use to write estimations performed by sdis_solve_camera */
+typedef res_T
+(*sdis_write_estimations_T)
+ (void* context, /* User data */
+ const size_t origin[2], /* Coordinates of the 1st estimation in image plane */
+ const size_t nestimations[2], /* #estimations in X and Y */
+ const struct sdis_mc* estimations); /* List of row ordered estimations */
+
BEGIN_DECLS
/*******************************************************************************
@@ -369,11 +377,25 @@ sdis_solve_probe
const size_t nrealisations, /* #realisations */
const double position[3], /* Probe position */
const double time, /* Observation time */
- const double fp_to_meter,/* Scale from floating point units to meters */
- const double ambient_radiative_temperature,
- const double reference_temperature,
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double ambient_radiative_temperature, /* In Kelvin */
+ const double reference_temperature, /* In Kelvin */
struct sdis_estimator** estimator);
+SDIS_API res_T
+sdis_solve_camera
+ (struct sdis_scene* scn,
+ const struct sdis_camera* cam, /* Point of view */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double ambient_radiative_temperature, /* In Kelvin */
+ const double reference_temperature, /* In Kelvin */
+ const size_t width, /* Image definition in in X */
+ const size_t height, /* Image definition in Y */
+ const size_t spp, /* #samples per pixel */
+ sdis_write_estimations_T writer,
+ void* writer_data);
+
END_DECLS
#endif /* SDIS_H */
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -52,6 +52,7 @@ device_release(ref_T* ref)
if(dev->s3d) S3D(device_ref_put(dev->s3d));
ASSERT(flist_name_is_empty(&dev->names));
flist_name_release(&dev->names);
+ darray_tile_release(&dev->tiles);
MEM_RM(dev->allocator, dev);
}
@@ -94,11 +95,19 @@ sdis_device_create
dev->nthreads = MMIN(nthreads_hint, (unsigned)omp_get_num_procs());
ref_init(&dev->ref);
flist_name_init(allocator, &dev->names);
+ darray_tile_init(allocator, &dev->tiles);
+
+ res = darray_tile_resize(&dev->tiles, dev->nthreads);
+ if(res != RES_OK) {
+ log_err(dev,
+ "%s: could not allocate the per thread buffer of estimations.\n", FUNC_NAME);
+ goto error;
+ }
res = s2d_device_create(log, allocator, 0, &dev->s2d);
if(res != RES_OK) {
log_err(dev,
- "%s, could not create the Star-2D device on Stardis.\n", FUNC_NAME);
+ "%s: could not create the Star-2D device on Stardis.\n", FUNC_NAME);
}
res = s3d_device_create(log, allocator, 0, &dev->s3d);
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -16,6 +16,7 @@
#ifndef SDIS_DEVICE_C_H
#define SDIS_DEVICE_C_H
+#include <rsys/dynamic_array.h>
#include <rsys/free_list.h>
#include <rsys/ref_count.h>
@@ -23,6 +24,18 @@ struct name { FITEM; };
#define FITEM_TYPE name
#include <rsys/free_list.h>
+#define DARRAY_NAME mc
+#define DARRAY_DATA struct sdis_mc
+#include <rsys/dynamic_array.h>
+
+#define DARRAY_NAME tile
+#define DARRAY_DATA struct darray_mc
+#define DARRAY_FUNCTOR_INIT darray_mc_init
+#define DARRAY_FUNCTOR_RELEASE darray_mc_release
+#define DARRAY_FUNCTOR_COPY darray_mc_copy
+#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_mc_copy_and_release
+#include <rsys/dynamic_array.h>
+
struct sdis_device {
struct logger* logger;
struct mem_allocator* allocator;
@@ -30,6 +43,7 @@ struct sdis_device {
int verbose;
struct flist_name names;
+ struct darray_tile tiles;
struct s2d_device* s2d;
struct s3d_device* s3d;
diff --git a/src/sdis_solve_probe.c b/src/sdis_solve_probe.c
@@ -14,6 +14,7 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "sdis.h"
+#include "sdis_camera.h"
#include "sdis_device_c.h"
#include "sdis_estimator_c.h"
#include "sdis_solve_probe_Xd.h"
@@ -29,6 +30,141 @@
#include <star/ssp.h>
#include <omp.h>
+/*******************************************************************************
+ * 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
+solve_pixel
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* mdm,
+ const struct sdis_camera* cam,
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ const size_t ipix[2], /* Pixel coordinate in the image plane */
+ const size_t nrealisations,
+ const double pix_sz[2], /* Pixel size in the normalized image plane */
+ struct sdis_mc* estimation)
+{
+ double weight = 0;
+ double sqr_weight = 0;
+ size_t N = 0; /* #realisations that do not fail */
+ size_t irealisation;
+ res_T res = RES_OK;
+ ASSERT(scn && mdm && rng && cam && ipix && nrealisations);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+
+ FOR_EACH(irealisation, 0, nrealisations) {
+ double samp[2]; /* Pixel sample */
+ double ray_pos[3];
+ double ray_dir[3];
+ double w = 0;
+
+ /* Generate a sample into the pixel to estimate */
+ samp[0] = (double)ipix[0] + ssp_rng_uniform_double(rng, 0, pix_sz[0]);
+ samp[1] = (double)ipix[1] + ssp_rng_uniform_double(rng, 0, pix_sz[1]);
+
+ /* Generate a ray starting from the camera position and passing through
+ * pixel sample */
+ camera_ray(cam, samp, ray_pos, ray_dir);
+
+ /* Launch the realisation */
+ res = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
+ time, fp_to_meter, Tarad, Tref, &w);
+ if(res != RES_OK) goto error;
+
+ if(res == RES_OK) {
+ weight += w;
+ sqr_weight += w*w;
+ ++N;
+ } else if(res != RES_BAD_OP) {
+ goto error;
+ }
+ }
+
+ if(N == 0) {
+ estimation->E = NaN;
+ estimation->SE = NaN;
+ estimation->V = NaN;
+ } else {
+ estimation->E = weight/(double)N;
+ estimation->V = sqr_weight/(double)N - estimation->E*estimation->E;
+ estimation->SE = sqrt(estimation->V/(double)N);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+solve_tile
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* mdm,
+ const struct sdis_camera* cam,
+ const double time,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ const size_t origin[2], /* Tile origin in image plane */
+ const size_t size[2], /* #pixels in X and Y */
+ const size_t spp, /* #samples per pixel */
+ const double pix_sz[2], /* Pixel size in the normalized image plane */
+ struct sdis_mc* estimations)
+{
+ size_t mcode; /* Morton code of the tile pixel */
+ size_t npixels;
+ res_T res = RES_OK;
+ ASSERT(scn && rng && mdm && cam && spp && origin && estimations);
+ ASSERT(size &&size[0] && size[1]);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+
+ /* Adjust the #pixels to process them wrt a morton order */
+ npixels = round_up_pow2(MMAX(size[0], size[1]));
+ npixels *= npixels;
+
+ FOR_EACH(mcode, 0, npixels) {
+ size_t ipix[2];
+ struct sdis_mc* estimation;
+
+ ipix[0] = morton2D_decode((uint32_t)(mcode>>0));
+ if(ipix[0] >= size[0]) continue;
+ ipix[1] = morton2D_decode((uint32_t)(mcode>>1));
+ if(ipix[1] >= size[1]) continue;
+
+ estimation = estimations + ipix[1]*size[0] + ipix[0];
+ ipix[0] = ipix[0] + origin[0];
+ ipix[1] = ipix[1] + origin[1];
+
+ res = solve_pixel(scn, rng, mdm, cam, time, fp_to_meter, Tarad, Tref, size,
+ spp, pix_sz, estimation);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
res_T
sdis_solve_probe
(struct sdis_scene* scn,
@@ -100,9 +236,7 @@ sdis_solve_probe
(scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w);
}
if(res_local != RES_OK) {
- if(res_local == RES_BAD_OP) {
- ++estimator->nfailures;
- } else {
+ if(res_local != RES_BAD_OP) {
ATOMIC_SET(&res, res_local);
continue;
}
@@ -114,6 +248,7 @@ sdis_solve_probe
}
estimator->nrealisations = N;
+ estimator->nfailures = nrealisations - N;
estimator->temperature.E = weight / (double)N;
estimator->temperature.V =
sqr_weight / (double)N
@@ -138,3 +273,141 @@ error:
goto exit;
}
+res_T
+sdis_solve_camera
+ (struct sdis_scene* scn,
+ const struct sdis_camera* cam,
+ const double time,
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ const size_t width, /* #pixels in X */
+ const size_t height, /* #pixels in Y */
+ const size_t spp, /* #samples per pixel */
+ sdis_write_estimations_T writer,
+ void* writer_data)
+{
+ #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);
+
+ const struct sdis_medium* medium = NULL;
+ struct darray_mc* tiles = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ size_t ntiles_x, ntiles_y, ntiles;
+ double pix_sz[2]; /* Size of a pixel in the normalized image plane */
+ int64_t mcode; /* Morton code of a tile */
+ size_t i;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !cam || time < 0 || fp_to_meter <= 0 || Tref < 0 || !width
+ || !height || !spp || !writer) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(scene_is_2d(scn)) {
+ log_err(scn->dev, "%s: 2D scene are not supported.\n", FUNC_NAME);
+ goto error;
+ }
+
+ /* Retrieve the medium in which the submitted position lies */
+ res = scene_get_medium(scn, cam->position, &medium);
+ if(res != RES_OK) goto error;
+
+ if(medium != SDIS_MEDIUM_FLUID) {
+ log_err(scn->dev, "%s: the camera position `%g %g %g' is not in a fluid.\n",
+ FUNC_NAME, SPLIT3(cam->position));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
+ if(!rngs) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Allocate per thread buffer of estimations */
+ tiles = darray_tile_data_get(&scn->dev->tiles);
+ ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads);
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ const size_t nestimations = TILE_SIZE * TILE_SIZE;
+ res = darray_mc_resize(tiles+i, nestimations);
+ if(res != RES_OK) goto error;
+ }
+
+ ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y));
+ ntiles *= ntiles;
+
+ pix_sz[0] = 1.0 / (double)width;
+ pix_sz[1] = 1.0 / (double)height;
+
+ omp_set_num_threads((int)scn->dev->nthreads);
+ /*#pragma omp parallel for schedule(static)*/
+ for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
+ size_t tile_org[2] = {0, 0};
+ size_t tile_sz[2] = {0, 0};
+ const int ithread = omp_get_thread_num();
+ struct sdis_mc* estimations = NULL;
+ struct ssp_rng* rng = rngs[ithread];
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue;
+
+ tile_org[0] = morton2D_decode((uint32_t)(mcode>>0));
+ if(tile_org[0] >= ntiles_x) continue; /* Discard tile */
+ tile_org[1] = morton2D_decode((uint32_t)(mcode>>1));
+ if(tile_org[1] >= ntiles_y) continue; /* Disaard tile */
+
+ /* Setup the tile coordinates in the image plane */
+ tile_org[0] *= TILE_SIZE;
+ tile_org[1] *= TILE_SIZE;
+ tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
+ tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[1]);
+
+ /* Fetch the estimation buffer */
+ estimations = darray_mc_data_get(tiles+ithread);
+
+ /* Draw the tile */
+ res_local = solve_tile(scn, rng, medium, cam, time, fp_to_meter, Tarad,
+ Tref, tile_org, tile_sz, spp, pix_sz, estimations);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ /* Write the estimations */
+ res_local = writer(writer_data, tile_org, tile_sz, estimations);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ }
+
+exit:
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ return (res_T)res;
+error:
+ goto exit;
+}
+
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -745,6 +745,71 @@ XD(probe_realisation)
return RES_OK;
}
+#if SDIS_SOLVE_PROBE_DIMENSION == 3
+static res_T
+XD(ray_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* medium,
+ const double position[],
+ const double direction[],
+ const double time,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ double* weight)
+{
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ const float range[2] = {0, FLT_MAX};
+ float org[3] = {0, 0, 0};
+ float dir[3] = {0, 0, 0};
+ res_T res = RES_OK;
+ ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight);
+ ASSERT(medium && medium->type == SDIS_MEDIUM_FLUID);
+
+ fX_set_dX(org, position);
+ fX_set_dX(dir, direction);
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &hit, &hit));
+
+ if(SXD_HIT_NONE(&hit)) {
+ if(Tarad >= 0) {
+ *weight = Tarad;
+ } else {
+ log_err(scn->dev,
+"%s: the ray starting from `%g %g %g' and traced along the `%g %g %g' direction\n"
+"reaches an invalid ambient temperature of `%gK'. One has to provide a valid\n"
+"ambient radiative temperature, i.e. it must be greater or equal to 0.\n",
+ FUNC_NAME, SPLIT3(org), SPLIT3(dir), Tarad);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ } else {
+ struct rwalk_context ctx;
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct XD(temperature) T = XD(TEMPERATURE_NULL);
+
+ dX(set)(rwalk.vtx.P, position);
+ XD(move_pos)(rwalk.vtx.P, dir, hit.distance);
+ rwalk.vtx.time = time;
+ rwalk.hit = hit;
+ rwalk.mdm = medium;
+
+ ctx.Tarad = Tarad;
+ ctx.Tref3 = Tref*Tref*Tref;
+
+ T.func = XD(boundary_temperature);
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+
+ *weight = T.value;
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+#endif /* SDIS_SOLVE_PROBE_DIMENSION == 3 */
+
#undef SDIS_SOLVE_PROBE_DIMENSION
#undef DIM
#undef sXd
