commit 1f624d506b3978b9681b523d18cec60a260b8d93
parent ac7c9edf63de832c06b6059845f18a2561d98bc4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 8 Apr 2015 11:17:37 +0200
Add the max_<steps|time> parameters to the smc_solve function
Diffstat:
4 files changed, 111 insertions(+), 49 deletions(-)
diff --git a/src/smc.h b/src/smc.h
@@ -33,6 +33,7 @@
#define SMC_H
#include <rsys/rsys.h>
+#include <limits.h>
/* Library symbol management */
#if defined(SMC_SHARED_BUILD) /* Build shared library */
@@ -71,6 +72,10 @@ struct smc_type {
void (*sqrt)(void* result, const void* value);
};
+static const struct smc_type SMC_TYPE_NULL =
+{ NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };
+
+
struct smc_estimator_status {
void* E; /* Expected value */
void* V; /* Variance */
@@ -78,8 +83,15 @@ struct smc_estimator_status {
unsigned long N; /* Samples count */
};
-static const struct smc_type smc_type_null =
-{ NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };
+struct smc_integrator {
+ void (*integrand)(void* value, struct ssp_rng* rng, void* ctx);
+ const struct smc_type* type;
+ unsigned long max_steps; /* Maximum # of experiments */
+ unsigned long max_time; /* Hint on the maximum time in micro seconds */
+};
+
+static const struct smc_integrator SMC_INTEGRATOR_NULL =
+{ NULL, NULL, ULONG_MAX, ULONG_MAX };
/* Forward declaration of opaque types */
struct smc_device; /* Entry point of the library */
@@ -115,8 +127,7 @@ smc_device_ref_put
SMC_API res_T
smc_solve
(struct smc_device* dev,
- void (*integrand)(void* value, struct ssp_rng* rng, void* ctx),
- const struct smc_type* type,
+ struct smc_integrator* integrator,
void* ctx,
struct smc_estimator** estimator);
diff --git a/src/smc_integrator.c b/src/smc_integrator.c
@@ -33,6 +33,7 @@
#include "smc_device_c.h"
#include "smc_type_c.h"
+#include <rsys/clock_time.h>
#include <rsys/mem_allocator.h>
#include <limits.h>
@@ -52,6 +53,17 @@ struct smc_estimator {
/*******************************************************************************
* Helper functions
******************************************************************************/
+static char
+check_integrator(struct smc_integrator* integrator)
+{
+ ASSERT(integrator);
+ return integrator->integrand
+ && integrator->type
+ && integrator->max_steps
+ && integrator->max_time
+ && check_type(integrator->type);
+}
+
static void
estimator_release(ref_T* ref)
{
@@ -81,17 +93,17 @@ estimator_release(ref_T* ref)
res_T
smc_solve
(struct smc_device* dev,
- void (*integrand)(void* value, struct ssp_rng* rng, void* ctx),
- const struct smc_type* type,
+ struct smc_integrator* integrator,
void* ctx,
struct smc_estimator** out_estimator)
{
struct smc_estimator* estimator = NULL;
+ struct time time_begin;
unsigned long i;
void* val = NULL;
res_T res = RES_OK;
- if(!dev || !integrand || !type || !out_estimator || !check_type(type)) {
+ if(!dev || !integrator || !out_estimator || !check_integrator(integrator)) {
res = RES_BAD_ARG;
goto error;
}
@@ -106,12 +118,12 @@ smc_solve
ref_init(&estimator->ref);
#define TYPE_CREATE(Dst) { \
- (Dst) = type->create(dev->allocator); \
+ (Dst) = integrator->type->create(dev->allocator); \
if(!(Dst)) { \
res = RES_MEM_ERR; \
goto error; \
} \
- type->zero((Dst)); \
+ integrator->type->zero((Dst)); \
} (void)0
TYPE_CREATE(estimator->value);
TYPE_CREATE(estimator->square_value);
@@ -122,17 +134,28 @@ smc_solve
#undef TYPE_CREATE
estimator->nsamples = 0;
estimator->status.N = 0;
- estimator->type = *type;
-
- FOR_EACH(i, 0, 16) {
- if(estimator->nsamples == ULONG_MAX) break;
- integrand(val, dev->rng, ctx);
- type->add(estimator->value, estimator->value, val);
- type->mul(val, val, val);
- type->add(estimator->square_value, estimator->square_value, val);
+ estimator->type = *integrator->type;
+
+ if(integrator->max_time != ULONG_MAX)
+ time_current(&time_begin);
+
+ FOR_EACH(i, 0, integrator->max_steps) {
+ struct time elapsed_time;
+
+ integrator->integrand(val, dev->rng, ctx);
+ integrator->type->add(estimator->value, estimator->value, val);
+ integrator->type->mul(val, val, val);
+ integrator->type->add(estimator->square_value, estimator->square_value, val);
++estimator->nsamples;
+
+ if(integrator->max_time != ULONG_MAX) {
+ time_current(&elapsed_time);
+ time_sub(&elapsed_time, &elapsed_time, &time_begin);
+ if((unsigned long)time_val(&elapsed_time, TIME_USEC) > integrator->max_time)
+ break;
+ }
}
- type->destroy(dev->allocator, val);
+ integrator->type->destroy(dev->allocator, val);
exit:
if(out_estimator) *out_estimator = estimator;
diff --git a/src/test_smc_light_path.c b/src/test_smc_light_path.c
@@ -43,10 +43,10 @@
#define IMG_WIDTH 640
#define IMG_HEIGHT 480
-#define LIGHT_PATH_DEPTH 4
-#define SKY_RADIANCE 0.1f
-#define ALBEDO 1.f
-#define EPSILON 1.e-3f
+#define LIGHT_PATH_DEPTH 16
+#define ALBEDO 0.6f
+#define EPSILON 1.e-1f
+#define GAMMA 2.2
struct cbox_desc{
const float* vertices;
@@ -205,13 +205,19 @@ struct integrator_context {
static float
direct_lighting(struct s3d_scene* scn, const float pos[3], const float N[3])
{
- const float light_pos[3] = { 200.f, 0.f, 400.f };
- const float radiance = 20000.f;
+ const float light_pos[3] = { 200.f, 200.f, 400.f };
+ const float flux = 60.0; /* Radiant flux in watt */
+ const float intensity = (float)(flux/(4.0*PI)); /*Radiant intensity in W/sr*/
struct s3d_hit hit;
float len;
float wi[3];
float range[2];
+ NCHECK(scn, NULL);
+ NCHECK(pos, NULL);
+ NCHECK(N, NULL);
+ CHECK(f3_is_normalized(N), 1);
+
f3_sub(wi, light_pos, pos);
len = f3_normalize(wi, wi);
@@ -219,11 +225,25 @@ direct_lighting(struct s3d_scene* scn, const float pos[3], const float N[3])
range[0] = EPSILON;
range[1] = len;
CHECK(s3d_scene_trace_ray(scn, pos, wi, range, &hit), RES_OK);
+ if(!S3D_HIT_NONE(&hit)) return 0.f; /* Light is occluded */
+
+ len *= 0.01f; /* Transform len from centimer to meter */
+ return
+ intensity / (len*len) /* radiance in W/(sr.m^2) */
+ * (float)(ALBEDO / PI) /* BRDF */
+ * f3_dot(wi, N); /* cos(wi,N) */
+}
- if(!S3D_HIT_NONE(&hit))
- return 0.f;
+static float
+skydome_lighting(const float wi[3])
+{
+ const float ground_irradiance = 0.01f;
+ const float sky_irradiance = 12.0f;
+ float u;
- return (radiance * (float)(ALBEDO / PI)) / (len*len) * f3_dot(wi, N);
+ NCHECK(wi, NULL);
+ u = CLAMP(wi[2], 0.f, 0.05f) * 1.f;
+ return u * sky_irradiance + (1.f - u) * ground_irradiance;
}
static void
@@ -267,7 +287,7 @@ light_path_integrator(void* value, struct ssp_rng* rng, void* data)
(ctx->scn, ray_org, ray_dir, ray_range, &hit), RES_OK);
if(S3D_HIT_NONE(&hit)) { /* Skydome lighting */
- L += throughput * SKY_RADIANCE;
+ L += throughput * skydome_lighting(ray_dir);
break;
}
@@ -304,6 +324,7 @@ main(int argc, char** argv)
struct s3d_shape* shape;
struct s3d_vertex_data attribs[2];
struct smc_device* smc;
+ struct smc_integrator integrator = SMC_INTEGRATOR_NULL;
unsigned char* img = NULL;
size_t ix, iy;
float pos[3];
@@ -359,14 +380,17 @@ main(int argc, char** argv)
ctx.pixel_size[0] = 1.f / (float)IMG_WIDTH;
ctx.pixel_size[1] = 1.f / (float)IMG_HEIGHT;
+ integrator.integrand = light_path_integrator;
+ integrator.type = &smc_float;
+ integrator.max_steps = 1;
+
FOR_EACH(iy, 0, IMG_HEIGHT) {
ctx.ipixel[1] = iy;
FOR_EACH(ix, 0, IMG_WIDTH) {
struct smc_estimator* estimator;
struct smc_estimator_status status;
ctx.ipixel[0] = ix;
- CHECK(smc_solve
- (smc, light_path_integrator, &smc_float, &ctx, &estimator), RES_OK);
+ CHECK(smc_solve(smc, &integrator, &ctx, &estimator), RES_OK);
if(img) {
const size_t ipix = (iy*IMG_WIDTH + ix) * 3/*RGB*/;
@@ -374,7 +398,7 @@ main(int argc, char** argv)
unsigned char colu;
CHECK(smc_estimator_get_status(estimator, &status), RES_OK);
- col = (float)pow(SMC_FLOAT(status.E), 1.0/2.2); /* Gamma correction */
+ col = (float)pow(SMC_FLOAT(status.E), 1.0/GAMMA); /* Gamma correction */
colu = (unsigned char)(CLAMP(col, 0.f, 1.f) * 255.f); /* Float to U8 */
img[ipix + 0] = colu;
img[ipix + 1] = colu;
diff --git a/src/test_smc_solve.c b/src/test_smc_solve.c
@@ -66,6 +66,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct smc_device* dev;
+ struct smc_integrator integrator = SMC_INTEGRATOR_NULL;
struct smc_estimator* estimator;
struct smc_estimator_status status;
(void)argc, (void)argv;
@@ -74,23 +75,18 @@ main(int argc, char** argv)
CHECK(smc_device_create(NULL, &allocator, &dev), RES_OK);
- CHECK(smc_solve(NULL, NULL, NULL, NULL, NULL), RES_BAD_ARG);
- CHECK(smc_solve(dev, NULL, NULL, NULL, NULL), RES_BAD_ARG);
- CHECK(smc_solve(NULL, rcp_x, NULL, NULL, NULL), RES_BAD_ARG);
- CHECK(smc_solve(dev, rcp_x, NULL, NULL, NULL), RES_BAD_ARG);
- CHECK(smc_solve(NULL, NULL, &smc_float, NULL, NULL), RES_BAD_ARG);
- CHECK(smc_solve(dev, NULL, &smc_float, NULL, NULL), RES_BAD_ARG);
- CHECK(smc_solve(NULL, rcp_x, &smc_float, NULL, NULL), RES_BAD_ARG);
- CHECK(smc_solve(dev, rcp_x, &smc_float, NULL, NULL), RES_BAD_ARG);
-
- CHECK(smc_solve(NULL, NULL, NULL, NULL, &estimator), RES_BAD_ARG);
- CHECK(smc_solve(dev, NULL, NULL, NULL, &estimator), RES_BAD_ARG);
- CHECK(smc_solve(NULL, rcp_x, NULL, NULL, &estimator), RES_BAD_ARG);
- CHECK(smc_solve(dev, rcp_x, NULL, NULL, &estimator), RES_BAD_ARG);
- CHECK(smc_solve(NULL, NULL, &smc_double, NULL, &estimator), RES_BAD_ARG);
- CHECK(smc_solve(dev, NULL, &smc_double, NULL, &estimator), RES_BAD_ARG);
- CHECK(smc_solve(NULL, rcp_x, &smc_double, NULL, &estimator), RES_BAD_ARG);
- CHECK(smc_solve(dev, rcp_x, &smc_double, NULL, &estimator), RES_OK);
+ integrator.integrand = rcp_x;
+ integrator.type = &smc_double;
+ integrator.max_steps = 4096;
+
+ CHECK(smc_solve(NULL, NULL, NULL, NULL), RES_BAD_ARG);
+ CHECK(smc_solve(dev, NULL, NULL, NULL), RES_BAD_ARG);
+ CHECK(smc_solve(NULL, &integrator, NULL, NULL), RES_BAD_ARG);
+ CHECK(smc_solve(dev, &integrator, NULL, NULL), RES_BAD_ARG);
+ CHECK(smc_solve(NULL, NULL, NULL, &estimator), RES_BAD_ARG);
+ CHECK(smc_solve(dev, NULL, NULL, &estimator), RES_BAD_ARG);
+ CHECK(smc_solve(NULL, &integrator, NULL, &estimator), RES_BAD_ARG);
+ CHECK(smc_solve(dev, &integrator, NULL, &estimator), RES_OK);
CHECK(smc_estimator_get_status(NULL, NULL), RES_BAD_ARG);
CHECK(smc_estimator_get_status(estimator, NULL), RES_BAD_ARG);
@@ -100,7 +96,15 @@ main(int argc, char** argv)
(log(2.0) - log(1.0), SMC_DOUBLE(status.E), SMC_DOUBLE(status.SE)), 1);
CHECK(smc_estimator_ref_put(estimator), RES_OK);
- CHECK(smc_solve(dev, cos_x, &smc_float, (void*)0xC0DE, &estimator), RES_OK);
+ integrator.type = NULL;
+ CHECK(smc_solve(dev, &integrator, NULL, &estimator), RES_BAD_ARG);
+ integrator.type = &smc_double;
+ integrator.integrand = NULL;
+ CHECK(smc_solve(dev, &integrator, NULL, &estimator), RES_BAD_ARG);
+
+ integrator.integrand = cos_x;
+ integrator.type = &smc_float;
+ CHECK(smc_solve(dev, &integrator, (void*)0xC0DE, &estimator), RES_OK);
CHECK(smc_estimator_get_status(estimator, &status), RES_OK);
CHECK(eq_eps
((float)(sin(3.0*PI/4.0) - sin(PI/4.0)),
