commit eab1bdb5d874df902a054389dfb25ac5ee4d286c
parent a70d7d2a0b5c62df4951a037d96da83292dcec51
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 4 Jun 2018 09:43:20 +0200
Draft of a null collision algorithm to estimate the transmission
Diffstat:
5 files changed, 115 insertions(+), 46 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -151,11 +151,10 @@ dump_buffer(struct htrdr* htrdr)
(unsigned long)buf_layout.height);
FOR_EACH(y, 0, buf_layout.height) {
FOR_EACH(x, 0, buf_layout.width) {
- if(*((double*)htrdr_buffer_at(htrdr->buf, x, y)) < 1.0) {
- fprintf(htrdr->output, "0 0 0\n");
- } else {
- fprintf(htrdr->output, "255 255 255\n");
- }
+ double val = *((double*)htrdr_buffer_at(htrdr->buf, x, y));
+ int i;
+ i = (int)(val * 255);
+ fprintf(htrdr->output, "%d %d %d\n", i, i, i);
}
}
}
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -19,6 +19,12 @@
#include <rsys/logger.h>
#include <rsys/ref_count.h>
+enum htrdr_voxel_data {
+ HTRDR_K_EXT_MIN,
+ HTRDR_K_EXT_MAX,
+ HTRDR_VOXEL_DATA_COUNT__
+};
+
/* Forward declaration */
struct htrdr_args;
struct htrdr_buffer;
diff --git a/src/htrdr_clouds.c b/src/htrdr_clouds.c
@@ -19,6 +19,7 @@
#include <star/svx.h>
#include <high_tune/htcp.h>
+#include <rsys/double3.h>
#include <rsys/dynamic_array_double.h>
#include <rsys/dynamic_array_size_t.h>
#include <rsys/hash_table.h>
@@ -46,6 +47,12 @@ vertex_eq(const struct vertex* v0, const struct vertex* v1)
#define HTABLE_KEY_FUNCTOR_EQ vertex_eq
#include <rsys/hash_table.h>
+struct build_octree_context {
+ struct htcp_desc* htcp_desc;
+ double dst_max; /* Max traversal distance */
+ double tau_threshold; /* Threshold criteria for the merge process */
+};
+
struct octree_data {
struct htable_vertex vertex2id; /* Map a coordinate to its vertex id */
struct darray_double vertices; /* Array of unique vertices */
@@ -84,10 +91,13 @@ register_leaf
{
struct octree_data* ctx = context;
struct vertex v[8];
+ const double* data;
int i;
ASSERT(leaf && ctx);
(void)ileaf;
+ data = leaf->data;
+
/* Compute the leaf vertices */
v[0].x = leaf->lower[0]; v[0].y = leaf->lower[1]; v[0].z = leaf->lower[2];
v[1].x = leaf->upper[0]; v[1].y = leaf->lower[1]; v[1].z = leaf->lower[2];
@@ -113,39 +123,71 @@ register_leaf
/* Add the vertex id to the leaf cell */
CHK(RES_OK == darray_size_t_push_back(&ctx->cells, &id));
}
- CHK(RES_OK == darray_double_push_back(&ctx->data, leaf->data));
+ FOR_EACH(i, 0, HTRDR_VOXEL_DATA_COUNT__) {
+ CHK(RES_OK == darray_double_push_back(&ctx->data, data+i));
+ }
}
static void
-vox_get(const size_t xyz[3], void* dst, void* ctx)
+vox_get(const size_t xyz[3], void* dst, void* context)
{
- struct htcp_desc* desc = ctx;
+ struct build_octree_context* ctx = context;
+ const double c_ext = 6.e-10; /* Extinction cross section in m^2.particle^-1 */
+ const double rho_air = 1.293; /* Volumic mass of terrestrial air in kg.m^-3 */
+ const double rho_h2o = 1000; /* Volumic mass of water in kg.m^-3 */
+ const double sigma = 0.1; /* Std deviation of the `r' parameter */
+ const double r = 9.76617; /* Modal radius of water particles in um */
+ double ql; /* Mass of the liquid water in the voxel in kg.m^-3 */
+ double v; /* typical volume of a particle */
+ double nv; /* Number density */
+ double k_ext; /* Extinction coefficient in m^-1 */
double* pdbl = dst;
ASSERT(xyz && dst && ctx);
- *pdbl = htcp_desc_RCT_at(desc, xyz[0], xyz[1], xyz[2], 0/*time*/);
+
+ ql = htcp_desc_RCT_at(ctx->htcp_desc, xyz[0], xyz[1], xyz[2], 0/*time*/);
+ v = 4*PI*(r*r*r)*exp(4.5*sigma*sigma) / 3.0;
+ nv = 1.e18 * ql * rho_air / (v * rho_h2o);
+ k_ext = c_ext*nv;
+
+ pdbl[HTRDR_K_EXT_MIN] = k_ext;
+ pdbl[HTRDR_K_EXT_MAX] = k_ext;
}
static void
vox_merge(void* dst, const void* voxels[], const size_t nvoxs, void* ctx)
{
double* pdbl = dst;
- double dbl_max = -DBL_MAX;
- size_t i;
+ double k_ext_min = DBL_MAX;
+ double k_ext_max =-DBL_MAX;
+ size_t ivox;
ASSERT(dst && voxels && nvoxs);
(void)ctx;
- FOR_EACH(i, 0, nvoxs) dbl_max = MMAX(dbl_max, *(double*)(voxels[i]));
- *pdbl = dbl_max;
+
+ FOR_EACH(ivox, 0, nvoxs) {
+ const double* voxel_data = voxels[ivox];
+ ASSERT(voxel_data[HTRDR_K_EXT_MIN] <= voxel_data[HTRDR_K_EXT_MAX]);
+ k_ext_min = MMIN(k_ext_min, voxel_data[HTRDR_K_EXT_MIN]);
+ k_ext_max = MMAX(k_ext_max, voxel_data[HTRDR_K_EXT_MAX]);
+ }
+ pdbl[HTRDR_K_EXT_MIN] = k_ext_min;
+ pdbl[HTRDR_K_EXT_MAX] = k_ext_max;
}
static int
-vox_challenge_merge(const void* voxels[], const size_t nvoxs, void* ctx)
+vox_challenge_merge(const void* voxels[], const size_t nvoxs, void* context)
{
- size_t i;
- (void)nvoxs, (void)ctx;
- FOR_EACH(i, 0, nvoxs) {
- if(*(double*)(voxels[i]) != 0) break;
+ struct build_octree_context* ctx = context;
+ double k_ext_min = DBL_MAX;
+ double k_ext_max =-DBL_MAX;
+ size_t ivox;
+ ASSERT(voxels && nvoxs);
+
+ FOR_EACH(ivox, 0, nvoxs) {
+ const double* voxel_data = voxels[ivox];
+ k_ext_min = MMIN(k_ext_min, voxel_data[HTRDR_K_EXT_MIN]);
+ k_ext_max = MMAX(k_ext_max, voxel_data[HTRDR_K_EXT_MAX]);
}
- return i >= nvoxs;
+ return (k_ext_max - k_ext_min)*ctx->dst_max <= ctx->tau_threshold;
}
/*******************************************************************************
@@ -158,8 +200,9 @@ clouds_load(struct htrdr* htrdr, const int verbose, const char* filename)
struct svx_tree* clouds = NULL;
struct htcp_desc htcp_desc = HTCP_DESC_NULL;
struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL;
+ struct build_octree_context ctx;
size_t nvoxs[3];
- double low[3], upp[3];
+ double low[3], upp[3], sz[3];
res_T res = RES_OK;
ASSERT(htrdr && filename);
@@ -196,12 +239,22 @@ clouds_load(struct htrdr* htrdr, const int verbose, const char* filename)
upp[1] = low[1] + (double)nvoxs[1] * htcp_desc.vxsz_y;
upp[2] = low[2] + (double)nvoxs[2] * htcp_desc.vxsz_z[0];
+ /* Compute the size of of the AABB */
+ sz[0] = upp[0] - low[0];
+ sz[1] = upp[1] - low[1];
+ sz[2] = upp[2] - low[2];
+
+ /* Setup the build context */
+ ctx.htcp_desc = &htcp_desc;
+ ctx.dst_max = sz[2];
+ ctx.tau_threshold = 0.0;
+
/* Setup the voxel descriptor */
vox_desc.get = vox_get;
vox_desc.merge = vox_merge;
vox_desc.challenge_merge = vox_challenge_merge;
- vox_desc.context = &htcp_desc;
- vox_desc.size = sizeof(double);
+ vox_desc.context = &ctx;
+ vox_desc.size = sizeof(double)*HTRDR_VOXEL_DATA_COUNT__;
/* Create the octree */
res = svx_octree_create(htrdr->svx, low, upp, nvoxs, &vox_desc, &clouds);
@@ -228,6 +281,7 @@ clouds_dump_vtk(struct htrdr* htrdr, FILE* stream)
{
struct svx_tree_desc desc;
struct octree_data data;
+ const double* leaf_data;
size_t nvertices;
size_t ncells;
size_t i;
@@ -277,13 +331,17 @@ clouds_dump_vtk(struct htrdr* htrdr, FILE* stream)
FOR_EACH(i, 0, ncells) fprintf(stream, "11\n");
/* Write the cell data */
+ leaf_data = darray_double_cdata_get(&data.data);
fprintf(stream, "CELL_DATA %lu\n", (unsigned long)ncells);
- fprintf(stream, "SCALARS Val double 1\n");
+ fprintf(stream, "SCALARS Val double %d\n", HTRDR_VOXEL_DATA_COUNT__);
fprintf(stream, "LOOKUP_TABLE default\n");
FOR_EACH(i, 0, ncells) {
- fprintf(stream, "%g\n", darray_double_cdata_get(&data.data)[i]);
+ size_t idata;
+ FOR_EACH(idata, 0, HTRDR_VOXEL_DATA_COUNT__) {
+ fprintf(stream, "%g ", leaf_data[i*HTRDR_VOXEL_DATA_COUNT__ + idata]);
+ }
+ fprintf(stream, "\n");
}
-
octree_data_release(&data);
return RES_OK;
}
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -22,13 +22,6 @@
struct htrdr;
extern LOCAL_SYM res_T
-htrdr_solve_transmission
- (struct htrdr* htrdr,
- const double pos[3],
- const double dir[3],
- double *val);
-
-extern LOCAL_SYM res_T
htrdr_solve_transmission_buffer
(struct htrdr* htrdr);
diff --git a/src/htrdr_transmission.c b/src/htrdr_transmission.c
@@ -25,9 +25,11 @@
#include <omp.h>
struct transmit_context {
- double tau;
+ double tau_sampled;
+ double tau_max_min;
+ double tau_min;
};
-static const struct transmit_context TRANSMIT_CONTEXT_NULL = {0};
+static const struct transmit_context TRANSMIT_CONTEXT_NULL = {0,0,0};
/*******************************************************************************
* Helper functions
@@ -51,25 +53,29 @@ discard_hit
const double range[2],
void* context)
{
+ const double* vox_data = NULL;
struct transmit_context* ctx = context;
double dst;
- double k;
+ double k_ext_min;
+ double k_ext_max;
ASSERT(hit && ctx && !SVX_HIT_NONE(hit));
(void)org, (void)dir, (void)range;
- k = *((double*)hit->voxel.data);
+ vox_data = hit->voxel.data;
+ k_ext_min = vox_data[HTRDR_K_EXT_MIN];
+ k_ext_max = vox_data[HTRDR_K_EXT_MAX];
+
dst = hit->distance[1] - hit->distance[0];
ASSERT(dst >= 0);
- ctx->tau += k*dst;
- return 1;
+ ctx->tau_min += k_ext_min*dst;
+ ctx->tau_max_min += (k_ext_max - k_ext_min)*dst;
+ return ctx->tau_max_min < ctx->tau_sampled;
}
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_solve_transmission
+static res_T
+transmission_realisation
(struct htrdr* htrdr,
+ struct ssp_rng* rng,
const double pos[3],
const double dir[3],
double *val)
@@ -80,6 +86,7 @@ htrdr_solve_transmission
res_T res = RES_OK;
ASSERT(htrdr && pos && dir && val);
+ ctx.tau_sampled = ssp_ran_exp(rng, 1.0);
res = svx_octree_trace_ray(htrdr->clouds, pos, dir, range, NULL,
discard_hit, &ctx, &hit);
if(res != RES_OK) {
@@ -89,8 +96,11 @@ htrdr_solve_transmission
SPLIT3(pos), SPLIT3(dir));
goto error;
}
- ASSERT(SVX_HIT_NONE(&hit));
- *val = ctx.tau ? exp(-ctx.tau) : 1.0;
+ if(!SVX_HIT_NONE(&hit)) {
+ *val = 0;
+ } else {
+ *val = ctx.tau_min ? exp(-ctx.tau_min) : 1.0;
+ }
exit:
return res;
@@ -98,6 +108,9 @@ error:
goto exit;
}
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
res_T
htrdr_solve_transmission_buffer(struct htrdr* htrdr)
{
@@ -165,7 +178,7 @@ htrdr_solve_transmission_buffer(struct htrdr* htrdr)
htrdr_rectangle_sample_pos(htrdr->rect, pixsamp, pos);
/* Solve the transmission for the sample position wrt the main dir */
- res_local = htrdr_solve_transmission(htrdr, pos, htrdr->main_dir, &T);
+ res_local = transmission_realisation(htrdr, rng, pos, htrdr->main_dir, &T);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
continue;
