commit b42b3747d155954a32888a81bc7b2b14209363ad
parent d6f2e64f78e01381df6241e70594cabe69bee1d4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 20 May 2019 09:49:32 +0200
Add an estimator of the overall image to the camera solver
Add support of path registration to the camera solver
Diffstat:
7 files changed, 161 insertions(+), 30 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -926,8 +926,11 @@ sdis_solve_camera
const size_t width, /* Image definition in in X */
const size_t height, /* Image definition in Y */
const size_t spp, /* #samples per pixel */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
sdis_write_accums_T writer,
- void* writer_data);
+ void* writer_data,
+ /* Estimator of the whole image. May be NULL */
+ struct sdis_estimator** estimator);
SDIS_API res_T
sdis_solve_medium
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -32,6 +32,7 @@ ray_realisation_3d
const double fp_to_meter,
const double Tarad,
const double Tref,
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -50,9 +51,14 @@ ray_realisation_3d
ctx.Tarad = Tarad;
ctx.Tref3 = Tref*Tref*Tref;
+ ctx.heat_path = heat_path;
f3_set_d3(dir, direction);
+ /* Register the starting position against the heat path */
+ res = register_heat_vertex(heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE);
+ if(res != RES_OK) goto error;
+
res = trace_radiative_path_3d(scn, dir, fp_to_meter, &ctx, &rwalk, rng, &T);
if(res != RES_OK) goto error;
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -46,8 +46,8 @@ probe_realisation_2d
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -61,8 +61,8 @@ probe_realisation_3d
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
/*******************************************************************************
@@ -79,8 +79,8 @@ boundary_realisation_2d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -94,8 +94,8 @@ boundary_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -140,6 +140,7 @@ ray_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
#endif /* SDIS_REALISATION_H */
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -121,8 +121,8 @@ XD(probe_realisation)
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -213,8 +213,8 @@ XD(boundary_realisation)
const double fp_to_meter,
const double Tarad,
const double Tref,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -72,8 +72,12 @@ solve_pixel
const double Tref, /* In Kelvin */
const size_t ipix[2], /* Pixel coordinate in the image plane */
const size_t nrealisations,
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_accum* accum)
+ struct accum* acc_temp,
+ struct accum* acc_time,
+ struct sdis_accum* accum,
+ struct sdis_estimator* estimator)
{
double sum_weights = 0;
double sum_weights_sqr = 0;
@@ -82,12 +86,25 @@ solve_pixel
res_T res = RES_OK;
ASSERT(scn && mdm && rng && cam && ipix && nrealisations && Tref >= 0);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+ ASSERT(acc_time && acc_temp && accum && estimator);
FOR_EACH(irealisation, 0, nrealisations) {
+ struct time t0, t1;
double samp[2]; /* Pixel sample */
double ray_pos[3];
double ray_dir[3];
double w = 0;
+ struct sdis_heat_path* pheat_path = NULL;
+ struct sdis_heat_path heat_path;
+ res_T res_simul = RES_OK;
+
+ /* Begin time registration */
+ time_current(&t0);
+
+ if(register_paths) {
+ heat_path_init(scn->dev->allocator, &heat_path);
+ pheat_path = &heat_path;
+ }
/* Generate a sample into the pixel to estimate */
samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
@@ -98,19 +115,46 @@ solve_pixel
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) {
+ res_simul = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
+ time, fp_to_meter, Tarad, Tref, pheat_path, &w);
+
+ /* Handle fatal error */
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ res = res_simul;
+ goto error;
+ }
+
+ if(pheat_path) {
+ pheat_path->status = res_simul == RES_OK
+ ? SDIS_HEAT_PATH_SUCCEED
+ : SDIS_HEAT_PATH_FAILED;
+
+ /* Check if the path must be saved regarding the register_paths mask */
+ if(!(register_paths & (int)pheat_path->status)) {
+ heat_path_release(pheat_path);
+ } else { /* Register the sampled path */
+ res = estimator_add_and_release_heat_path(estimator, pheat_path);
+ if(res != RES_OK) goto error;
+ }
+ }
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ if(res_simul == RES_OK) {
+ /* Update global accumulators */
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+ acc_temp->sum += w; acc_temp->sum2 += w*w; ++acc_temp->count;
+ acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count;
+
+ /* Update per pixel accumulator */
sum_weights += w;
sum_weights_sqr += w*w;
++N;
- } else if(res == RES_BAD_OP) {
- res = RES_OK;
- } else {
- goto error;
}
}
+ /* Setup per pixel accumulator */
accum->sum_weights = sum_weights;
accum->sum_weights_sqr = sum_weights_sqr;
accum->nweights = N;
@@ -135,14 +179,18 @@ solve_tile
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 int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_accum* accums)
+ struct accum* acc_temp,
+ struct accum* acc_time,
+ struct sdis_accum* accums,
+ struct sdis_estimator* estimator)
{
size_t mcode; /* Morton code of the tile pixel */
size_t npixels;
res_T res = RES_OK;
ASSERT(scn && rng && mdm && cam && spp && origin && accums && Tref >= 0);
- ASSERT(size &&size[0] && size[1]);
+ ASSERT(size &&size[0] && size[1] && acc_temp && acc_time && estimator);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
/* Adjust the #pixels to process them wrt a morton order */
@@ -163,7 +211,7 @@ solve_tile
ipix[1] = ipix[1] + origin[1];
res = solve_pixel(scn, rng, mdm, cam, time, fp_to_meter, Tarad, Tref, ipix,
- spp, pix_sz, accum);
+ spp, register_paths, pix_sz, acc_temp, acc_time, accum, estimator);
if(res != RES_OK) goto error;
}
@@ -367,14 +415,19 @@ sdis_solve_camera
const size_t width, /* #pixels in X */
const size_t height, /* #pixels in Y */
const size_t spp, /* #samples per pixel */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
sdis_write_accums_T writer,
- void* writer_data)
+ void* writer_data,
+ struct sdis_estimator** out_estimator)
{
#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);
+ struct sdis_estimator* estimator = NULL;
struct sdis_medium* medium = NULL;
struct darray_accum* tiles = NULL;
+ struct accum* acc_temps = NULL;
+ struct accum* acc_times = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
size_t ntiles_x, ntiles_y, ntiles;
@@ -384,11 +437,10 @@ sdis_solve_camera
ATOMIC res = RES_OK;
if(!scn || !cam || time < 0 || fp_to_meter <= 0 || Tref < 0 || !width
- || !height || !spp || !writer) {
+ || !height || !spp || !writer || !out_estimator) {
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;
@@ -422,6 +474,14 @@ sdis_solve_camera
if(res != RES_OK) goto error;
}
+ /* Create the per thread accumulators */
+ acc_temps = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps));
+ if(!acc_temps) { res = RES_MEM_ERR; goto error; }
+ acc_times = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_times));
+ if(!acc_times) { res = RES_MEM_ERR; goto error; }
+
/* Allocate per thread buffer of accumulations */
tiles = darray_tile_data_get(&scn->dev->tiles);
ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads);
@@ -439,6 +499,10 @@ sdis_solve_camera
pix_sz[0] = 1.0 / (double)width;
pix_sz[1] = 1.0 / (double)height;
+ /* Create the estimator */
+ res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator);
+ if(res != RES_OK) goto error;
+
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static, 1/*chunk size*/)
for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
@@ -447,6 +511,8 @@ sdis_solve_camera
const int ithread = omp_get_thread_num();
struct sdis_accum* accums = NULL;
struct ssp_rng* rng = rngs[ithread];
+ struct accum* acc_temp = &acc_temps[ithread];
+ struct accum* acc_time = &acc_times[ithread];
res_T res_local = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue;
@@ -467,7 +533,8 @@ sdis_solve_camera
/* Draw the tile */
res_local = solve_tile(scn, rng, medium, cam, time, fp_to_meter, Tarad,
- Tref, tile_org, tile_sz, spp, pix_sz, accums);
+ Tref, tile_org, tile_sz, spp, register_paths, pix_sz, acc_temp, acc_time,
+ accums, estimator);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
continue;
@@ -481,6 +548,22 @@ sdis_solve_camera
}
}
+ /* Setup the estimated temperature and per realisation time for the whole
+ * image */
+ if(out_estimator) {
+ const size_t nrealisations = width*height*spp;
+ struct accum acc_temp;
+ struct accum acc_time;
+
+ sum_accums(acc_temps, scn->dev->nthreads, &acc_temp);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ ASSERT(acc_temp.count == acc_time.count);
+
+ estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count);
+ estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
+ }
+
exit:
if(rngs) {
FOR_EACH(i, 0, scn->dev->nthreads) {
@@ -488,7 +571,10 @@ exit:
}
MEM_RM(scn->dev->allocator, rngs);
}
+ if(acc_temps) MEM_RM(scn->dev->allocator, acc_temps);
+ if(acc_times) MEM_RM(scn->dev->allocator, acc_times);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_estimator) *out_estimator = estimator;
return (res_T)res;
error:
goto exit;
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -90,7 +90,7 @@ XD(solve_probe)
if(res != RES_OK) goto error;
}
- /* Create the per thread accumulator */
+ /* Create the per thread accumulators */
acc_temps = MEM_CALLOC
(scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps));
if(!acc_temps) { res = RES_MEM_ERR; goto error; }
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -509,9 +509,12 @@ main(int argc, char** argv)
struct geometry geom = GEOMETRY_NULL;
struct s3dut_mesh* msh = NULL;
struct s3dut_mesh_data msh_data;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_accum_buffer* buf = NULL;
struct sdis_camera* cam = NULL;
struct sdis_device* dev = NULL;
+ struct sdis_estimator* estimator = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid0 = NULL;
struct sdis_medium* fluid1 = NULL;
@@ -522,6 +525,7 @@ main(int argc, char** argv)
struct solid solid_param = SOLID_NULL;
struct interf interface_param = INTERF_NULL;
size_t ntris, npos;
+ size_t nreals, nfails;
double pos[3];
double tgt[3];
double up[3];
@@ -606,15 +610,46 @@ main(int argc, char** argv)
dump_mesh(stdout, geom.positions, npos, geom.indices, ntris);
exit(0);
#endif
+ BA(sdis_solve_camera(NULL, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+ BA(sdis_solve_camera(scn, NULL, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+ BA(sdis_solve_camera(scn, cam, INF, 0, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+ BA(sdis_solve_camera(scn, cam, INF, 1, 300, -1, IMG_WIDTH, IMG_HEIGHT, SPP,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+ BA(sdis_solve_camera(scn, cam, INF, 1, 300, 300, 0, IMG_HEIGHT, SPP,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+ BA(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, 0, SPP,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+ BA(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, 0,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+ BA(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
+ SDIS_HEAT_PATH_NONE, NULL, buf, &estimator));
+ BA(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, NULL));
/* Launch the simulation */
OK(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
- sdis_accum_buffer_write, buf));
+ SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ CHK(nreals + nfails == IMG_WIDTH*IMG_HEIGHT*SPP);
+
+ fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE);
+ fprintf(stderr, "Time per realisation (in usec) ~ %g +/- %g\n", time.E, time.SE);
+ fprintf(stderr, "#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP));
/* Write the image */
dump_image(buf);
/* Release memory */
+ OK(sdis_estimator_ref_put(estimator));
OK(sdis_scene_ref_put(scn));
OK(sdis_camera_ref_put(cam));
OK(sdis_interface_ref_put(interf0));
