solstice-solver

Solver library of the solstice app
git clone git://git.meso-star.com/solstice-solver.git
Log | Files | Refs | README | LICENSE

commit 48fed335d585e7a9add505e3b7dac14b099a4086
parent 64ab6b264087b556c6e454a73d685756d9ff7975
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  2 Jun 2017 15:27:31 +0200

Merge branch 'release-0.2'

Diffstat:
MREADME.md | 10++++++++++
Mcmake/CMakeLists.txt | 18+++++++++++-------
Msrc/ssol.h | 185++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------
Msrc/ssol_atmosphere.c | 76++++++++++++++++++++++++++++++++++------------------------------------------
Msrc/ssol_atmosphere_c.h | 18+-----------------
Asrc/ssol_data.c | 99+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/ssol_device.c | 2+-
Msrc/ssol_draw_draft.c | 41+++++++++++++++++++++++++++++++++--------
Msrc/ssol_draw_pt.c | 99+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------
Msrc/ssol_estimator.c | 47+++++++++++++++++++++++++++++++++++++++++------
Msrc/ssol_estimator_c.h | 38++++++++++++++++----------------------
Msrc/ssol_image.c | 79++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
Msrc/ssol_instance.c | 2+-
Msrc/ssol_material.c | 334+++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------
Msrc/ssol_material_c.h | 55+++++++++++++++++++++++++++----------------------------
Msrc/ssol_mc_receiver.c | 40+++++++++++++++++++++++++++++++++++++---
Msrc/ssol_object.c | 7++++---
Msrc/ssol_param_buffer.c | 43+++++++++++++++++++++++++++++++++----------
Msrc/ssol_ranst_sun_wl.c | 18++++++++++--------
Msrc/ssol_scene.c | 48++++++++++++++++++++++++++++++++++++------------
Msrc/ssol_scene_c.h | 5+++++
Msrc/ssol_shape.c | 539+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------
Msrc/ssol_shape_c.h | 53+++++++++++++++++++++++++++++++++++++++++++++--------
Msrc/ssol_solver.c | 350++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
Msrc/ssol_spectrum.c | 85++++++++++++++++++++++++++++++++++++++++++++-----------------------------------
Msrc/ssol_spectrum_c.h | 13+++++++------
Msrc/ssol_sun.c | 85++++++++++++++++++++++++++++++++++++++++++-------------------------------------
Msrc/ssol_sun_c.h | 8++++++--
Msrc/test_ssol_atmosphere.c | 48++++++++++++++++++++++++++++++++++--------------
Asrc/test_ssol_data.c | 85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_ssol_draw.c | 57+++++++++++++++++++++++++++++++++++++++++++++++++++------
Msrc/test_ssol_image.c | 98++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---
Msrc/test_ssol_instance.c | 6+++++-
Msrc/test_ssol_material.c | 32++++++++++++++++----------------
Msrc/test_ssol_materials.h | 64+++++++++++++++-------------------------------------------------
Msrc/test_ssol_object.c | 5+++++
Msrc/test_ssol_param_buffer.c | 94++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------------
Msrc/test_ssol_scene.c | 26+++++++-------------------
Msrc/test_ssol_shape.c | 32+++++++++++++++++++++++++-------
Msrc/test_ssol_solver1.c | 143++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
Msrc/test_ssol_solver2.c | 11+++++++----
Msrc/test_ssol_solver2b.c | 9++++++---
Msrc/test_ssol_solver3.c | 12+++++++++---
Msrc/test_ssol_solver4.c | 19+++++++++++--------
Msrc/test_ssol_solver5.c | 8+++++---
Msrc/test_ssol_solver6.c | 4++++
Msrc/test_ssol_solver7.c | 29+++++++++++++++--------------
Msrc/test_ssol_solver8.c | 8++++----
Msrc/test_ssol_solver9.c | 8+++-----
Msrc/test_ssol_spectrum.c | 8++++++++
50 files changed, 2250 insertions(+), 953 deletions(-)

diff --git a/README.md b/README.md @@ -24,6 +24,16 @@ well as all the aforementioned prerequisites. Finally generate the project from the `cmake/CMakeLists.txt` file by appending to the `CMAKE_PREFIX_PATH` variable the install directories of its dependencies. +## Release notes + +### Version 0.2 + +- Add normal maps to describe spatially varying normals in the tangent space of + the surface. +- Add support of spectral data to the atmosphere and the materials. +- Fix the per primitive irradiance estimate by dividing the result by the area + of the primitive in order to have watts per square meter. + ## Licenses Solstice-Solver is developed by [|Meso|Star>](http://www.meso-star.com) for the diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -23,10 +23,11 @@ option(NO_TEST "Do not build tests" OFF) ################################################################################ # Check dependencies ################################################################################ -find_package(RCMake 0.2.3 REQUIRED) +find_package(RCMake 0.3 REQUIRED) find_package(RSys 0.4 REQUIRED) -find_package(Star3D 0.4 REQUIRED) -find_package(StarCPR REQUIRED) +find_package(Star3D 0.4.1 REQUIRED) +find_package(Star3DUT 0.2 REQUIRED) +find_package(StarCPR 0.1 REQUIRED) find_package(StarSF 0.1 REQUIRED) find_package(StarSP 0.4 REQUIRED) find_package(OpenMP 1.2 REQUIRED) @@ -38,23 +39,25 @@ include(rcmake_runtime) include_directories( ${RSys_INCLUDE_DIR} ${Star3D_INCLUDE_DIR} + ${Star3DUT_INCLUDE_DIR} ${StarCPR_INCLUDE_DIR} ${StarSF_INCLUDE_DIR} ${StarSP_INCLUDE_DIR}) -rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarCPR StarSF StarSP) +rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D Star3DUT StarCPR StarSF StarSP) ################################################################################ # Configure and define targets ################################################################################ set(VERSION_MAJOR 0) -set(VERSION_MINOR 1) -set(VERSION_PATCH 1) +set(VERSION_MINOR 2) +set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(SSOL_FILES_SRC ssol_atmosphere.c ssol_camera.c + ssol_data.c ssol_device.c ssol_draw.c ssol_draw_pt.c @@ -104,7 +107,7 @@ rcmake_prepend_path(SSOL_FILES_INC_API ${SSOL_SOURCE_DIR}) rcmake_prepend_path(SSOL_FILES_DOC ${PROJECT_SOURCE_DIR}/../) add_library(solstice-solver SHARED ${SSOL_FILES_SRC} ${SSOL_FILES_INC} ${SSOL_FILES_INC_API}) -target_link_libraries(solstice-solver RSys Star3D StarCPR StarSF StarSP) +target_link_libraries(solstice-solver RSys Star3D Star3DUT StarCPR StarSF StarSP) if(CMAKE_COMPILER_IS_GNUCC) target_link_libraries(solstice-solver m) @@ -152,6 +155,7 @@ if(NOT NO_TEST) new_test(test_ssol_atmosphere) new_test(test_ssol_by_receiver_integration) new_test(test_ssol_camera) + new_test(test_ssol_data) new_test(test_ssol_device) new_test(test_ssol_image) new_test(test_ssol_material) diff --git a/src/ssol.h b/src/ssol.h @@ -92,9 +92,14 @@ enum ssol_pixel_format { SSOL_PIXEL_FORMATS_COUNT__ }; -enum ssol_parametrization_type { - SSOL_PARAMETRIZATION_TEXCOORD, /* Map from 3D to 2D with texcoord */ - SSOL_PARAMETRIZATION_PRIMITIVE_ID /* Map from 3D to 1D with primitive id */ +enum ssol_filter_mode { + SSOL_FILTER_LINEAR, + SSOL_FILTER_NEAREST +}; + +enum ssol_address_mode { + SSOL_ADDRESS_CLAMP, + SSOL_ADDRESS_REPEAT }; enum ssol_quadric_type { @@ -102,6 +107,7 @@ enum ssol_quadric_type { SSOL_QUADRIC_PARABOL, SSOL_QUADRIC_HYPERBOL, SSOL_QUADRIC_PARABOLIC_CYLINDER, + SSOL_QUADRIC_HEMISPHERE, SSOL_QUADRIC_TYPE_COUNT__ }; @@ -113,6 +119,11 @@ enum ssol_attrib_usage { SSOL_ATTRIBS_COUNT__ }; +enum ssol_data_type { + SSOL_DATA_REAL, + SSOL_DATA_SPECTRUM +}; + /* Describe a vertex data */ struct ssol_vertex_data { enum ssol_attrib_usage usage; /* Semantic of the data */ @@ -166,7 +177,7 @@ static const struct ssol_quadric_parabol SSOL_QUADRIC_PARABOL_NULL = SSOL_QUADRIC_PARABOL_NULL__; struct ssol_quadric_hyperbol { - /* Define (x^2 + y^2) / a^2 - (z - 1/2)^2 / b^2 + 1 = 0 + /* Define (x^2 + y^2) / a^2 - (z - 1/2)^2 / b^2 + 1 = 0; z > 0 * with a^2 = f - f^2; b = f -1/2; f = real_focal / (img_focal + real_focal) */ double img_focal, real_focal; }; @@ -181,6 +192,14 @@ struct ssol_quadric_parabolic_cylinder { static const struct ssol_quadric_parabolic_cylinder SSOL_QUADRIC_PARABOLIC_CYLINDER_NULL = SSOL_QUADRIC_PARABOLIC_CYLINDER_NULL__; +struct ssol_quadric_hemisphere { + /* Define x^2 + y^2 + (z-radius)^2 - radius^2 = 0 with z <= r */ + double radius; +}; +#define SSOL_QUADRIC_HEMISPHERE_NULL__ { -1.0 } +static const struct ssol_quadric_hemisphere SSOL_QUADRIC_HEMISPHERE_NULL = +SSOL_QUADRIC_HEMISPHERE_NULL__; + struct ssol_quadric { enum ssol_quadric_type type; union { @@ -188,12 +207,13 @@ struct ssol_quadric { struct ssol_quadric_parabol parabol; struct ssol_quadric_hyperbol hyperbol; struct ssol_quadric_parabolic_cylinder parabolic_cylinder; + struct ssol_quadric_hemisphere hemisphere; } data; /* 3x4 column major transformation of the quadric in object space */ double transform[12]; - /* Hint on the how to discretised */ + /* Hint on how to discretise */ size_t slices_count_hint; }; @@ -227,23 +247,44 @@ struct ssol_punched_surface { static const struct ssol_punched_surface SSOL_PUNCHED_SURFACE_NULL = SSOL_PUNCHED_SURFACE_NULL__; +struct ssol_data { + enum ssol_data_type type; + union { + double real; + struct ssol_spectrum* spectrum; + } value; +}; +#define SSOL_DATA_NULL__ {SSOL_DATA_REAL, {0.0}} +static const struct ssol_data SSOL_DATA_NULL = SSOL_DATA_NULL__; + struct ssol_medium { - double absorptivity; - double refractive_index; + struct ssol_data absorption; + struct ssol_data refractive_index; }; -#define SSOL_MEDIUM_VACUUM__ { 0, 1 } +#define SSOL_MEDIUM_VACUUM__ {{SSOL_DATA_REAL, {0}}, {SSOL_DATA_REAL, {1}}} static const struct ssol_medium SSOL_MEDIUM_VACUUM = SSOL_MEDIUM_VACUUM__; +struct ssol_surface_fragment { + double dir[3]; /* World space incoming direction. Point forward the surface */ + double P[3]; /* World space position */ + double Ng[3]; /* Normalized world space geometry normal */ + double Ns[3]; /* Normalized world space shading normal */ + double uv[2]; /* Texture coordinates */ + double dPdu[3]; /* Partial derivative of the position in u */ + double dPdv[3]; /* Partial derivative of the position in v */ +}; + +#define SSOL_SURFACE_FRAGMENT_NULL__ \ + {{0,0,0}, {0,0,0}, {0,0,0}, {0,0,0}, {0,0}, {0,0,0}, {0,0,0}} +static const struct ssol_surface_fragment SSOL_SURFACE_FRAGMENT_NULL = + SSOL_SURFACE_FRAGMENT_NULL__; + typedef void (*ssol_shader_getter_T) (struct ssol_device* dev, struct ssol_param_buffer* buf, - const double wavelength, /* In nanometer */ - const double P[3], /* World space position */ - const double Ng[3], /* World space geometry normal */ - const double Ns[3], /* World space shading normal */ - const double uv[2], /* Texture coordinates */ - const double w[3], /* Incoming direction. Point toward the surface */ + const double wavelength, + const struct ssol_surface_fragment* fragment, double* val); /* Returned value */ /* Dielectric material shader */ @@ -285,7 +326,6 @@ SSOL_THIN_DIELECTRIC_SHADER_NULL = SSOL_THIN_DIELECTRIC_SHADER_NULL__; * FILE* argument */ struct ssol_receiver_data { uint64_t realization_id; - int64_t date; uint32_t segment_id; /* Its absolute value is the identifier of an SSOL instance. A negative @@ -347,22 +387,28 @@ struct ssol_mc_result { static const struct ssol_mc_result SSOL_MC_RESULT_NULL = SSOL_MC_RESULT_NULL__; struct ssol_mc_global { - struct ssol_mc_result cos_loss; /* In W */ + struct ssol_mc_result cos_factor; /* [0 1] */ + struct ssol_mc_result absorbed; /* In W */ struct ssol_mc_result shadowed; /* In W */ struct ssol_mc_result missing; /* In W */ + struct ssol_mc_result atmosphere; /* In W */ + struct ssol_mc_result reflectivity; /* In W */ }; #define SSOL_MC_GLOBAL_NULL__ { \ SSOL_MC_RESULT_NULL__, \ SSOL_MC_RESULT_NULL__, \ + SSOL_MC_RESULT_NULL__, \ + SSOL_MC_RESULT_NULL__, \ + SSOL_MC_RESULT_NULL__, \ SSOL_MC_RESULT_NULL__ \ } static const struct ssol_mc_global SSOL_MC_GLOBAL_NULL = SSOL_MC_GLOBAL_NULL__; struct ssol_mc_receiver { struct ssol_mc_result integrated_irradiance; /* In W */ + struct ssol_mc_result integrated_absorbed_irradiance; /* In W */ struct ssol_mc_result absorptivity_loss; /* In W */ struct ssol_mc_result reflectivity_loss; /* In W */ - struct ssol_mc_result cos_loss; /* In W TODO remove this */ /* Internal data */ size_t N__; @@ -379,6 +425,14 @@ struct ssol_mc_receiver { static const struct ssol_mc_receiver SSOL_MC_RECEIVER_NULL = SSOL_MC_RECEIVER_NULL__; +#define MC_RCV_NONE__ { \ + { -1, -1, -1 }, \ + { -1, -1, -1 }, \ + { -1, -1, -1 }, \ + { -1, -1, -1 }, \ + 0, NULL, NULL \ +} + struct ssol_mc_shape { /* Internal data */ size_t N__; @@ -388,11 +442,17 @@ struct ssol_mc_shape { #define SSOL_MC_SHAPE_NULL__ { 0, NULL, NULL } static const struct ssol_mc_shape SSOL_MC_SHAPE_NULL = SSOL_MC_SHAPE_NULL__; +struct ssol_mc_sampled { + struct ssol_mc_result cos_factor; /* [0 1] */ + struct ssol_mc_result shadowed; + size_t nb_samples; +}; + struct ssol_mc_primitive { struct ssol_mc_result integrated_irradiance; /* In W */ + struct ssol_mc_result integrated_absorbed_irradiance; /* In W */ struct ssol_mc_result absorptivity_loss; /* In W */ struct ssol_mc_result reflectivity_loss; /* In W */ - struct ssol_mc_result cos_loss; /* In W TODO remove this */ }; #define SSOL_MC_PRIMITIVE_NULL__ { \ SSOL_MC_RESULT_NULL__, \ @@ -516,12 +576,21 @@ ssol_image_get_layout SSOL_API res_T ssol_image_map (const struct ssol_image* image, - void** memory); + char** memory); SSOL_API res_T ssol_image_unmap (const struct ssol_image* image); +SSOL_API res_T +ssol_image_sample + (const struct ssol_image* image, + const enum ssol_filter_mode filter, + const enum ssol_address_mode address_u, + const enum ssol_address_mode address_v, + const double uv[2], + void* val); + /* Helper function that matches the `ssol_write_pixels_T' functor type */ SSOL_API res_T ssol_image_write @@ -567,7 +636,7 @@ ssol_scene_compute_aabb /* Detach all the instances from the scene and release the reference that the * scene takes onto them. - * Also detach the attached sun if any. */ + * Also detach the attached sun and atmosphere if any. */ SSOL_API res_T ssol_scene_clear (struct ssol_scene* scn); @@ -863,7 +932,10 @@ SSOL_API void* ssol_param_buffer_allocate (struct ssol_param_buffer* buf, const size_t size, - const size_t alignment); /* Power of 2 in [1, 64] */ + const size_t alignment, /* Power of 2 in [1, 64] */ + /* Functor to invoke on the allocated memory priorly to its destruction. + * May be NULL */ + void (*release)(void*)); /* Retrieve the address of the first allocated parameter */ SSOL_API void* @@ -971,7 +1043,7 @@ ssol_sun_set_buie_param ******************************************************************************/ /* The atmosphere describes absorption along the light paths */ SSOL_API res_T -ssol_atmosphere_create_uniform +ssol_atmosphere_create (struct ssol_device* dev, struct ssol_atmosphere** atmosphere); @@ -983,11 +1055,10 @@ SSOL_API res_T ssol_atmosphere_ref_put (struct ssol_atmosphere* atmosphere); -/* List of per wavelength power of the sun */ SSOL_API res_T -ssol_atmosphere_set_uniform_absorption +ssol_atmosphere_set_absorption (struct ssol_atmosphere* atmosphere, - struct ssol_spectrum* spectrum); + struct ssol_data* absorption); /******************************************************************************* * Estimator API - Describe the state of a simulation. @@ -1014,7 +1085,7 @@ ssol_estimator_get_mc_sampled_x_receiver struct ssol_mc_receiver* rcv); SSOL_API res_T -ssol_estimator_get_count +ssol_estimator_get_realisation_count (const struct ssol_estimator* estimator, size_t* count); @@ -1029,6 +1100,17 @@ ssol_estimator_get_sampled_area (const struct ssol_estimator* estimator, double* area); +SSOL_API res_T +ssol_estimator_get_sampled_count + (const struct ssol_estimator* estimator, + size_t* count); + +SSOL_API res_T +ssol_estimator_get_mc_sampled + (struct ssol_estimator* estimator, + const struct ssol_instance* samp_instance, + struct ssol_mc_sampled* sampled); + /******************************************************************************* * Tracked paths ******************************************************************************/ @@ -1110,9 +1192,60 @@ ssol_draw_pt const size_t width, /* #pixels in X */ const size_t height, /* #pixels in Y */ const size_t spp, + const double up[3], /* Direction toward the top of the skydome */ ssol_write_pixels_T writer, void* writer_data); +/******************************************************************************* + * Data API + ******************************************************************************/ +SSOL_API struct ssol_data* +ssol_data_set_real + (struct ssol_data* data, + const double real); + +/* Get a reference onto the submitted spectrum */ +SSOL_API struct ssol_data* +ssol_data_set_spectrum + (struct ssol_data* data, + struct ssol_spectrum* spectrum); + +/* Release the reference on its associated spectrum, if defined */ +SSOL_API struct ssol_data* +ssol_data_clear + (struct ssol_data* data); + +SSOL_API struct ssol_data* +ssol_data_copy + (struct ssol_data* dst, + const struct ssol_data* src); + +SSOL_API double +ssol_data_get_value + (const struct ssol_data* data, + const double wavelength); + +/******************************************************************************* + * Medium API + ******************************************************************************/ +static FINLINE struct ssol_medium* +ssol_medium_clear(struct ssol_medium* medium) +{ + ASSERT(medium); + ssol_data_clear(&medium->absorption); + ssol_data_clear(&medium->refractive_index); + return medium; +} + +static FINLINE struct ssol_medium* +ssol_medium_copy(struct ssol_medium* dst, const struct ssol_medium* src) +{ + ASSERT(dst && src); + ssol_data_copy(&dst->absorption, &src->absorption); + ssol_data_copy(&dst->refractive_index, &src->refractive_index); + return dst; +} + END_DECLS #endif /* SSOL_H */ diff --git a/src/ssol_atmosphere.c b/src/ssol_atmosphere.c @@ -33,48 +33,63 @@ atmosphere_release(ref_T* ref) ASSERT(ref); dev = atmosphere->dev; ASSERT(dev && dev->allocator); - switch (atmosphere->type) { - case ATMOS_UNIFORM: - if (atmosphere->data.uniform.spectrum) - SSOL(spectrum_ref_put(atmosphere->data.uniform.spectrum)); - break; - default: FATAL("Unreachable code\n"); break; - } + ssol_data_clear(&atmosphere->absorption); MEM_RM(dev->allocator, atmosphere); SSOL(device_ref_put(dev)); } +static INLINE int +check_absorption(const struct ssol_data* absorption) +{ + if(!absorption) return 0; + + /* Check absorptivity in [0, INF) */ + switch(absorption->type) { + case SSOL_DATA_REAL: + if(absorption->value.real < 0 || absorption->value.real > 1) + return 0; + break; + case SSOL_DATA_SPECTRUM: + if(!absorption->value.spectrum + || !spectrum_check_data(absorption->value.spectrum, 0, 1)) + return 0; + break; + default: FATAL("Unreachable code\n"); break; + } + + return 1; +} + /******************************************************************************* * Exported ssol_atmosphere functions ******************************************************************************/ res_T -ssol_atmosphere_create_uniform +ssol_atmosphere_create (struct ssol_device* dev, struct ssol_atmosphere** out_atmosphere) { struct ssol_atmosphere* atmosphere = NULL; res_T res = RES_OK; - if (!dev || !out_atmosphere) { + if(!dev || !out_atmosphere) { return RES_BAD_ARG; } atmosphere = (struct ssol_atmosphere*)MEM_CALLOC (dev->allocator, 1, sizeof(struct ssol_atmosphere)); - if (!atmosphere) { + if(!atmosphere) { res = RES_MEM_ERR; goto error; } SSOL(device_ref_get(dev)); atmosphere->dev = dev; - atmosphere->type = ATMOS_UNIFORM; ref_init(&atmosphere->ref); exit: - if (out_atmosphere) *out_atmosphere = atmosphere; + if(out_atmosphere) *out_atmosphere = atmosphere; return res; error: - if (atmosphere) { + if(atmosphere) { SSOL(atmosphere_ref_put(atmosphere)); atmosphere = NULL; } @@ -85,8 +100,7 @@ res_T ssol_atmosphere_ref_get (struct ssol_atmosphere* atmosphere) { - if (!atmosphere) - return RES_BAD_ARG; + if(!atmosphere) return RES_BAD_ARG; ref_get(&atmosphere->ref); return RES_OK; } @@ -95,41 +109,19 @@ res_T ssol_atmosphere_ref_put (struct ssol_atmosphere* atmosphere) { - if (!atmosphere) - return RES_BAD_ARG; + if(!atmosphere) return RES_BAD_ARG; ref_put(&atmosphere->ref, atmosphere_release); return RES_OK; } res_T -ssol_atmosphere_set_uniform_absorption +ssol_atmosphere_set_absorption (struct ssol_atmosphere* atmosphere, - struct ssol_spectrum* spectrum) + struct ssol_data* absorption) { - struct atm_uniform* uni; - if (!atmosphere || !spectrum || atmosphere->type != ATMOS_UNIFORM) + if(!atmosphere || !absorption || !check_absorption(absorption)) return RES_BAD_ARG; - uni = &atmosphere->data.uniform; - if (spectrum == uni->spectrum) /* no change */ - return RES_OK; - if (uni->spectrum) - SSOL(spectrum_ref_put(uni->spectrum)); - SSOL(spectrum_ref_get(spectrum)); - uni->spectrum = spectrum; + ssol_data_copy(&atmosphere->absorption, absorption); return RES_OK; } -/******************************************************************************* - * Local functions - ******************************************************************************/ -double -atmosphere_uniform_get_absorption - (const struct ssol_atmosphere* atmosphere, - const double wavelength) -{ - const struct ssol_spectrum* spectrum; - ASSERT(atmosphere && atmosphere->type == ATMOS_UNIFORM && wavelength >= 0); - spectrum = atmosphere->data.uniform.spectrum; - return spectrum_interpolate(spectrum, wavelength); -} - diff --git a/src/ssol_atmosphere_c.h b/src/ssol_atmosphere_c.h @@ -21,29 +21,13 @@ struct ssol_scene; -enum atmosphere_type { - ATMOS_UNIFORM, - ATMOS_TYPES_COUNT__ -}; - -struct atm_uniform { - struct ssol_spectrum* spectrum; -}; - struct ssol_atmosphere { - enum atmosphere_type type; struct ssol_scene* scene_attachment; - union { - struct atm_uniform uniform; - } data; + struct ssol_data absorption; struct ssol_device* dev; ref_T ref; }; -extern LOCAL_SYM double -atmosphere_uniform_get_absorption - (const struct ssol_atmosphere* atmosphere, - const double wavelength); #endif /* SSOL_ATMOSPHERE_C_H */ diff --git a/src/ssol_data.c b/src/ssol_data.c @@ -0,0 +1,99 @@ +/* Copyright (C) CNRS 2016-2017 + * + * 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/>. */ + +#ifndef SSOL_DATA_C_H +#define SSOL_DATA_C_H + +#include "ssol.h" +#include "ssol_spectrum_c.h" + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +struct ssol_data* +ssol_data_set_real(struct ssol_data* data, const double real) +{ + const double r = real; + ASSERT(data); + ssol_data_clear(data); + data->type = SSOL_DATA_REAL; + data->value.real = r; + return data; +} + +struct ssol_data* +ssol_data_set_spectrum(struct ssol_data* data, struct ssol_spectrum* spectrum) +{ + ASSERT(data && spectrum); + if(data->type == SSOL_DATA_SPECTRUM && data->value.spectrum == spectrum) + return data; + ssol_data_clear(data); + data->type = SSOL_DATA_SPECTRUM; + data->value.spectrum = spectrum; + SSOL(spectrum_ref_get(spectrum)); + return data; +} + +struct ssol_data* +ssol_data_clear(struct ssol_data* data) +{ + ASSERT(data); + if(data->type != SSOL_DATA_SPECTRUM) return data; + ASSERT(data->value.spectrum); + SSOL(spectrum_ref_put(data->value.spectrum)); + *data = SSOL_DATA_NULL; + return data; +} + +struct ssol_data* +ssol_data_copy(struct ssol_data* dst, const struct ssol_data* src) +{ + ASSERT(dst && src); + if(dst == src) return dst; + ssol_data_clear(dst); + switch(src->type) { + case SSOL_DATA_REAL: + dst->value.real = src->value.real; + break; + case SSOL_DATA_SPECTRUM: + ssol_data_set_spectrum(dst, src->value.spectrum); + break; + default: FATAL("Unreachable code.\n"); break; + } + dst->type = src->type; + return dst; +} + +double +ssol_data_get_value(const struct ssol_data* data, const double wavelength) +{ + double val; + ASSERT(data); + + switch(data->type) { + case SSOL_DATA_REAL: + val = data->value.real; + break; + case SSOL_DATA_SPECTRUM: + ASSERT(wavelength >= 0); + val = spectrum_interpolate(data->value.spectrum, wavelength); + break; + default: FATAL("Unreachable code\n"); break; + } + return val; +} + +#endif /* SSOL_DATA_C_H */ + diff --git a/src/ssol_device.c b/src/ssol_device.c @@ -91,7 +91,7 @@ ssol_device_create res = darray_tile_resize(&dev->tiles, dev->nthreads); if(res != RES_OK) goto error; - res = s3d_device_create(logger, mem_allocator, verbose, &dev->s3d); + res = s3d_device_create(logger, mem_allocator, 0, &dev->s3d); if(res != RES_OK) goto error; res = scpr_mesh_create(mem_allocator, &dev->scpr_mesh); diff --git a/src/ssol_draw_draft.c b/src/ssol_draw_draft.c @@ -17,6 +17,7 @@ #include "ssol_camera.h" #include "ssol_device_c.h" #include "ssol_draw.h" +#include "ssol_material_c.h" #include "ssol_object_c.h" #include "ssol_scene_c.h" #include "ssol_shape_c.h" @@ -41,6 +42,7 @@ Li double val[3]) { const float range[2] = {0, FLT_MAX}; + struct ssol_surface_fragment frag; struct ray_data ray_data = RAY_DATA_NULL; struct s3d_hit hit; ASSERT(scn && view && org && dir && val); @@ -52,9 +54,12 @@ Li d3_splat(val, 0); } else { struct ssol_instance* inst; + struct ssol_material* mtl; const struct shaded_shape* sshape; size_t isshape; - float N[3]={0}; + double o[3], wi[3]; + double N[3]={0}; + double cos_N_wi; /* Retrieve the hit shaded shape */ inst = *htable_instance_find(&scn->instances_rt, &hit.prim.inst_id); @@ -65,12 +70,28 @@ Li /* Retrieve and normalized the hit normal */ switch(sshape->shape->type) { - case SHAPE_MESH: f3_normalize(N, hit.normal); break; - case SHAPE_PUNCHED: f3_normalize(N, f3_set_d3(N, ray_data.N)); break; + case SHAPE_MESH: d3_normalize(N, d3_set_f3(N, hit.normal)); break; + case SHAPE_PUNCHED: d3_normalize(N, ray_data.N); break; + break; default: FATAL("Unreachable code"); break; } - ASSERT(f3_is_normalized(N)); - d3_splat(val, fabs(f3_dot(N, dir))); + + d3_set_f3(o, org); + d3_set_f3(wi, dir); + d3_normalize(wi, wi); + if(d3_dot(N, wi) < 0) { + mtl = sshape->mtl_front; + } else { + mtl = sshape->mtl_back; + d3_minus(N, N); + } + + surface_fragment_setup(&frag, o, wi, N, &hit.prim, hit.uv); + material_shade_normal(mtl, &frag, 1/*TODO wavelength*/, N); + + ASSERT(d3_is_normalized(N)); + cos_N_wi = d3_dot(N, d3_minus(wi, wi)); + d3_splat(val, MMAX(cos_N_wi, 0)); } } @@ -84,14 +105,14 @@ draw_pixel const float pix_sz[2], /* Normalized pixel size */ const size_t nsamples, double pixel[3], - void* ctx) + void* data) { - struct darray_float* samples = ctx; + struct darray_float* samples = data; float samp[2]; float ray_org[3], ray_dir[3]; double sum[3] = {0, 0, 0}; size_t i; - ASSERT(scn && cam && view && pix_coords && pix_sz && nsamples && pixel && ctx); + ASSERT(scn && cam && view && pix_coords && pix_sz && nsamples && pixel && data); (void)ithread; FOR_EACH(i, 0, nsamples) { @@ -135,6 +156,10 @@ ssol_draw_draft if(!scn || !spp) return RES_BAD_ARG; darray_float_init(scn->dev->allocator, &samples); + + res = scene_check(scn, FUNC_NAME); + if(res != RES_OK) goto error; + res = darray_float_reserve(&samples, spp * 2/*#dimensions*/); if(res != RES_OK) goto error; diff --git a/src/ssol_draw_pt.c b/src/ssol_draw_pt.c @@ -14,11 +14,13 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "ssol_c.h" +#include "ssol_atmosphere_c.h" #include "ssol_camera.h" #include "ssol_device_c.h" #include "ssol_draw.h" #include "ssol_material_c.h" #include "ssol_object_c.h" +#include "ssol_ranst_sun_wl.h" #include "ssol_scene_c.h" #include "ssol_shape_c.h" #include "ssol_sun_c.h" @@ -37,6 +39,8 @@ struct thread_context { struct ssp_rng* rng; struct ssf_bsdf* bsdf; + struct ranst_sun_wl* ran_wl; + float up[3]; }; static void @@ -45,6 +49,7 @@ thread_context_release(struct thread_context* ctx) ASSERT(ctx); if(ctx->rng) SSP(rng_ref_put(ctx->rng)); if(ctx->bsdf) SSF(bsdf_ref_put(ctx->bsdf)); + if(ctx->ran_wl) ranst_sun_wl_ref_put(ctx->ran_wl); } static res_T @@ -67,12 +72,17 @@ error: static void thread_context_setup (struct thread_context* ctx, - struct ssp_rng* rng) + struct ssp_rng* rng, + struct ranst_sun_wl* ran_wl, + const double up[3]) { - ASSERT(ctx && rng); + ASSERT(ctx && rng && ran_wl && up); if(ctx->rng) SSP(rng_ref_put(ctx->rng)); SSP(rng_ref_get(rng)); + ranst_sun_wl_ref_get(ran_wl); ctx->rng = rng; + ctx->ran_wl = ran_wl; + f3_set_d3(ctx->up, up); } /* Declare the container of the per thread contexts */ @@ -108,13 +118,12 @@ sun_lighting d3_minus(wi, sun->direction); /* The point look backward the sun */ - if(d3_dot(wi, N) < 0) return 0.0; + cos_wi_N = d3_dot(wi, N); + if(cos_wi_N < 0 || eq_eps(cos_wi_N, 0, 1.e-6)) return 0.0; R = ssf_bsdf_eval(bsdf, wo, N, wi); if(R <= 0) return 0.0; - cos_wi_N = d3_dot(wi, N); - f3_set_d3(ray_dir, wi); S3D(scene_view_trace_ray(view, ray_org, ray_dir, ray_range, ray_data, &hit)); if(S3D_HIT_NONE(&hit)) return R * cos_wi_N; @@ -129,13 +138,13 @@ Li(struct ssol_scene* scn, const float dir[3], double val[3]) { - struct ssol_medium medium; + struct ssol_medium medium = SSOL_MEDIUM_VACUUM; struct s3d_hit hit; struct ray_data ray_data = RAY_DATA_NULL; struct ssol_instance* inst; struct ssol_material* mtl; const struct shaded_shape* sshape; - struct surface_fragment frag; + struct ssol_surface_fragment frag; size_t isshape; double throughput = 1.0; double wi[3], o[3], uv[3]; @@ -144,6 +153,8 @@ Li(struct ssol_scene* scn, double L = 0; double R; double pdf; + double cos_wi_Ng; + double wl; const float ray_range[2] = {0, FLT_MAX}; float ray_org[3]; float ray_dir[3]; @@ -159,22 +170,28 @@ Li(struct ssol_scene* scn, f3_set(ray_org, org); f3_set(ray_dir, dir); - /* Assume that the path starts from vacuum */ - medium = SSOL_MEDIUM_VACUUM; + wl = ranst_sun_wl_get(ctx->ran_wl, ctx->rng); + + if(scn->atmosphere) { + ssol_data_copy(&medium.absorption, &scn->atmosphere->absorption); + } for(;;) { + double absorption; S3D(scene_view_trace_ray (view, ray_org, ray_dir, ray_range, &ray_data, &hit)); - if(medium.absorptivity > 0) { - throughput *= exp(-medium.absorptivity * hit.distance); - } - if(S3D_HIT_NONE(&hit)) { /* Background lighting */ - if(ray_dir[2] > 0) L += throughput * 1.e-1; + if(f3_dot(ray_dir, ctx->up) > 0) L += throughput * 1.e-1; break; } + absorption = ssol_data_get_value(&medium.absorption, wl); + if(absorption > 0) { + throughput *= exp(-absorption * hit.distance); + if(throughput <= 0) break; + } + /* Retrieve the hit shaded shape */ inst = *htable_instance_find(&scn->instances_rt, &hit.prim.inst_id); isshape = *htable_shaded_shape_find @@ -188,8 +205,12 @@ Li(struct ssol_scene* scn, /* Retrieve and normalized the hit normal */ switch(sshape->shape->type) { - case SHAPE_MESH: d3_normalize(N, d3_set_f3(N, hit.normal)); break; - case SHAPE_PUNCHED: d3_normalize(N, ray_data.N); break; + case SHAPE_MESH: + d3_normalize(N, d3_set_f3(N, hit.normal)); + break; + case SHAPE_PUNCHED: + d3_normalize(N, ray_data.N); + break; default: FATAL("Unreachable code"); break; } @@ -203,16 +224,25 @@ Li(struct ssol_scene* scn, } surface_fragment_setup(&frag, o, wo, N, &hit.prim, hit.uv); + material_shade_normal(mtl, &frag, wl, N); + + /* Shaded normal may look backward the outgoing direction */ + if(d3_dot(N, wo) > 0) break; + SSF(bsdf_clear(ctx->bsdf)); - res = material_shade_rendering - (mtl, &frag, 1/*TODO wavelength*/, &medium, ctx->bsdf); + res = material_setup_bsdf + (mtl, &frag, wl, &medium, 1/*Rendering*/, ctx->bsdf); if(res != RES_OK) goto error; /* Update the ray */ ray_data.prim_from = hit.prim; ray_data.inst_from = inst; ray_data.side_from = side; - f3_mulf(ray_dir, ray_dir, hit.distance); + switch(sshape->shape->type) { + case SHAPE_MESH: f3_mulf(ray_dir, ray_dir, hit.distance); break; + case SHAPE_PUNCHED: f3_mulf(ray_dir, ray_dir, (float)ray_data.dst); break; + default: FATAL("Unreachable code"); break; + } f3_add(ray_org, ray_org, ray_dir); d3_minus(wo, wo); @@ -221,8 +251,18 @@ Li(struct ssol_scene* scn, (scn->sun, view, &ray_data, ctx->bsdf, wo, N, ray_org); } - R = ssf_bsdf_sample(ctx->bsdf, ctx->rng, wo, frag.Ns, wi, &type, &pdf); + /* Sampling a bounce direction */ + R = ssf_bsdf_sample(ctx->bsdf, ctx->rng, wo, N, wi, &type, &pdf); ASSERT(0 <= R && R <= 1); + + /* Due to the shading normal, the sampled direction may point in the wrong + * direction wrt the sampled BSDF component. */ + cos_wi_Ng = d3_dot(frag.Ng, wi); + if((cos_wi_Ng > 0 && (type & SSF_TRANSMISSION)) + || (cos_wi_Ng < 0 && (type & SSF_REFLECTION))) { + R = 0; + } + f3_set_d3(ray_dir, wi); if(type & SSF_TRANSMISSION) material_get_next_medium(mtl, &medium, &medium); @@ -230,7 +270,7 @@ Li(struct ssol_scene* scn, throughput *= fabs(d3_dot(wi, N)) * R; } else { if(ssp_rng_canonical(ctx->rng) >= R) break; - throughput *= d3_dot(wi, N); + throughput *= fabs(d3_dot(wi, N)); } if(throughput <= 0) break; @@ -242,6 +282,7 @@ Li(struct ssol_scene* scn, d3_splat(val, L); exit: + ssol_medium_clear(&medium); return res; error: d3(val, 1, 1, 0); @@ -271,7 +312,7 @@ draw_pixel ctx = darray_thread_context_data_get(thread_ctxs) + ithread; FOR_EACH(isample, 0, nsamples) { - const int MAX_NFAILURES = 10; + const int MAX_NFAILURES = 100; double weight[3]; float samp[2]; /* Pixel sample */ float ray_org[3], ray_dir[3]; @@ -313,24 +354,31 @@ ssol_draw_pt const size_t width, const size_t height, const size_t spp, + const double up[3], ssol_write_pixels_T writer, void* data) { struct darray_thread_context thread_ctxs; struct ssp_rng_proxy* rng_proxy = NULL; + struct ranst_sun_wl* ran_sun_wl = NULL; size_t i; res_T res = RES_OK; - if(!scn) - return RES_BAD_ARG; + if(!scn || !up) return RES_BAD_ARG; darray_thread_context_init(scn->dev->allocator, &thread_ctxs); + res = scene_check(scn, FUNC_NAME); + if(res != RES_OK) goto error; + /* Create a RNG proxy */ res = ssp_rng_proxy_create (scn->dev->allocator, &ssp_rng_threefry, scn->dev->nthreads, &rng_proxy); if(res != RES_OK) goto error; + res = sun_create_wavelength_distribution(scn->sun, &ran_sun_wl); + if(res != RES_OK) goto error; + /* Create the thread contexts */ res = darray_thread_context_resize(&thread_ctxs, scn->dev->nthreads); if(res != RES_OK) goto error; @@ -343,7 +391,7 @@ ssol_draw_pt res = ssp_rng_proxy_create_rng(rng_proxy, i, &rng); if(res != RES_OK) goto error; - thread_context_setup(ctx, rng); + thread_context_setup(ctx, rng, ran_sun_wl, up); SSP(rng_ref_put(rng)); } @@ -354,6 +402,7 @@ ssol_draw_pt exit: darray_thread_context_release(&thread_ctxs); if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); + if(ran_sun_wl) ranst_sun_wl_ref_put(ran_sun_wl); return (res_T)res; error: goto exit; diff --git a/src/ssol_estimator.c b/src/ssol_estimator.c @@ -118,9 +118,12 @@ ssol_estimator_get_mc_global global->Name.V = data->sqr_weight / N - global->Name.E*global->Name.E; \ global->Name.SE = global->Name.V > 0 ? sqrt(global->Name.V / N) : 0; \ } (void)0 - SETUP_MC_RESULT(cos_loss); + SETUP_MC_RESULT(cos_factor); + SETUP_MC_RESULT(absorbed); SETUP_MC_RESULT(shadowed); SETUP_MC_RESULT(missing); + SETUP_MC_RESULT(atmosphere); + SETUP_MC_RESULT(reflectivity); #undef SETUP_MC_RESULT return RES_OK; } @@ -145,9 +148,8 @@ ssol_estimator_get_mc_sampled_x_receiver memset(rcv, 0, sizeof(rcv[0])); - mc_samp = htable_sampled_find(&estimator->mc_sampled, &samp_instance); - if(!mc_samp || !mc_samp->nb_samples) { + if(!mc_samp) { /* The sampled instance has no MC estimation */ return RES_BAD_ARG; } @@ -161,16 +163,16 @@ ssol_estimator_get_mc_sampled_x_receiver mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back; #define SETUP_MC_RESULT(Name) { \ - const double N = (double)mc_samp->nb_samples; \ + const double N = (double)estimator->realisation_count; \ const struct mc_data* data = &mc_rcv1->Name; \ rcv->Name.E = data->weight / N; \ rcv->Name.V = data->sqr_weight / N - rcv->Name.E*rcv->Name.E; \ rcv->Name.SE = rcv->Name.V > 0 ? sqrt(rcv->Name.V / N) : 0; \ } (void)0 SETUP_MC_RESULT(integrated_irradiance); + SETUP_MC_RESULT(integrated_absorbed_irradiance); SETUP_MC_RESULT(absorptivity_loss); SETUP_MC_RESULT(reflectivity_loss); - SETUP_MC_RESULT(cos_loss); #undef SETUP_MC_RESULT rcv->mc__ = mc_rcv1; rcv->N__ = mc_samp->nb_samples; @@ -178,7 +180,7 @@ ssol_estimator_get_mc_sampled_x_receiver } res_T -ssol_estimator_get_count +ssol_estimator_get_realisation_count (const struct ssol_estimator* estimator, size_t* count) { if (!estimator || !count) return RES_BAD_ARG; @@ -205,6 +207,39 @@ ssol_estimator_get_sampled_area } res_T +ssol_estimator_get_sampled_count + (const struct ssol_estimator* estimator, size_t* count) +{ + if (!estimator || !count) return RES_BAD_ARG; + *count = htable_sampled_size_get(&estimator->mc_sampled); + return RES_OK; +} + +res_T +ssol_estimator_get_mc_sampled + (struct ssol_estimator* estimator, + const struct ssol_instance* samp_instance, + struct ssol_mc_sampled* sampled) +{ + struct mc_sampled* mc = NULL; + if (!estimator || !samp_instance || !sampled) return RES_BAD_ARG; + mc = htable_sampled_find(&estimator->mc_sampled, &samp_instance); + if(!mc) return RES_BAD_ARG; + sampled->nb_samples = mc->nb_samples; + #define SETUP_MC_RESULT(Name) { \ + const double N = (double)estimator->realisation_count; \ + const struct mc_data* data = &mc->Name; \ + sampled->Name.E = data->weight / N; \ + sampled->Name.V = data->sqr_weight/N - sampled->Name.E*sampled->Name.E; \ + sampled->Name.SE = sampled->Name.V > 0 ? sqrt(sampled->Name.V / N) : 0; \ + } (void)0 + SETUP_MC_RESULT(cos_factor); + SETUP_MC_RESULT(shadowed); + #undef SETUP_MC_RESULT + return RES_OK; +} + +res_T ssol_estimator_get_tracked_paths_count (const struct ssol_estimator* estimator, size_t* npaths) { diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h @@ -37,9 +37,9 @@ static const struct mc_data MC_DATA_NULL = MC_DATA_NULL__; #define MC_RECEIVER_DATA \ struct mc_data integrated_irradiance; /* In W */ \ + struct mc_data integrated_absorbed_irradiance; /* In W */ \ struct mc_data absorptivity_loss; /* In W */ \ - struct mc_data reflectivity_loss; /* In W */ \ - struct mc_data cos_loss; /* In W */ + struct mc_data reflectivity_loss; /* In W */ #define MC_RECEIVER_DATA_NULL__ \ MC_DATA_NULL__, \ @@ -148,9 +148,9 @@ mc_receiver_1side_init { ASSERT(mc); mc->integrated_irradiance = MC_DATA_NULL; + mc->integrated_absorbed_irradiance = MC_DATA_NULL; mc->absorptivity_loss = MC_DATA_NULL; mc->reflectivity_loss = MC_DATA_NULL; - mc->cos_loss = MC_DATA_NULL; htable_shape2mc_init(allocator, &mc->shape2mc); } @@ -167,9 +167,9 @@ mc_receiver_1side_copy { ASSERT(dst && src); dst->integrated_irradiance = src->integrated_irradiance; + dst->integrated_absorbed_irradiance = src->integrated_absorbed_irradiance; dst->absorptivity_loss = src->absorptivity_loss; dst->reflectivity_loss = src->reflectivity_loss; - dst->cos_loss = src->cos_loss; return htable_shape2mc_copy(&dst->shape2mc, &src->shape2mc); } @@ -179,9 +179,9 @@ mc_receiver_1side_copy_and_release { ASSERT(dst && src); dst->integrated_irradiance = src->integrated_irradiance; + dst->integrated_absorbed_irradiance = src->integrated_absorbed_irradiance; dst->absorptivity_loss = src->absorptivity_loss; dst->reflectivity_loss = src->reflectivity_loss; - dst->cos_loss = src->cos_loss; return htable_shape2mc_copy_and_release(&dst->shape2mc, &src->shape2mc); } @@ -278,10 +278,8 @@ mc_receiver_copy_and_release ******************************************************************************/ struct mc_sampled { /* Global data for this entity */ - struct mc_data cos_loss; + struct mc_data cos_factor; struct mc_data shadowed; - double area; - double sun_cos; size_t nb_samples; /* By-receptor data for this entity */ @@ -294,10 +292,8 @@ mc_sampled_init struct mc_sampled* samp) { ASSERT(samp); - samp->cos_loss = MC_DATA_NULL; + samp->cos_factor = MC_DATA_NULL; samp->shadowed = MC_DATA_NULL; - samp->area = 0; - samp->sun_cos = 0; samp->nb_samples = 0; htable_receiver_init(allocator, &samp->mc_rcvs); } @@ -313,10 +309,8 @@ static INLINE res_T mc_sampled_copy(struct mc_sampled* dst, const struct mc_sampled* src) { ASSERT(dst && src); - dst->cos_loss = src->cos_loss; + dst->cos_factor = src->cos_factor; dst->shadowed = src->shadowed; - dst->area = src->area; - dst->sun_cos = src->sun_cos; dst->nb_samples = src->nb_samples; return htable_receiver_copy(&dst->mc_rcvs, &src->mc_rcvs); } @@ -325,10 +319,8 @@ static INLINE res_T mc_sampled_copy_and_release(struct mc_sampled* dst, struct mc_sampled* src) { ASSERT(dst && src); - dst->cos_loss = src->cos_loss; + dst->cos_factor = src->cos_factor; dst->shadowed = src->shadowed; - dst->area = src->area; - dst->sun_cos = src->sun_cos; dst->nb_samples = src->nb_samples; return htable_receiver_copy_and_release(&dst->mc_rcvs, &src->mc_rcvs); } @@ -455,9 +447,12 @@ struct ssol_estimator { size_t failed_count; /* Implicit MC computations */ + struct mc_data cos_factor; + struct mc_data absorbed; struct mc_data shadowed; struct mc_data missing; - struct mc_data cos_loss; /* TODO compute it */ + struct mc_data atmosphere; + struct mc_data reflectivity; struct htable_receiver mc_receivers; /* Per receiver MC */ struct htable_sampled mc_sampled; /* Per sampled instance MC */ @@ -518,21 +513,20 @@ get_mc_sampled struct mc_sampled** out_mc_samp) { struct mc_sampled* mc_samp = NULL; - struct mc_sampled mc_samp_null; res_T res = RES_OK; ASSERT(sampled && inst && out_mc_samp); - mc_sampled_init(inst->dev->allocator, &mc_samp_null); - mc_samp = htable_sampled_find(sampled, &inst); if(!mc_samp) { + struct mc_sampled mc_samp_null; + mc_sampled_init(inst->dev->allocator, &mc_samp_null); res = htable_sampled_set(sampled, &inst, &mc_samp_null); + mc_sampled_release(&mc_samp_null); if(res != RES_OK) goto error; mc_samp = htable_sampled_find(sampled, &inst); } exit: - mc_sampled_release(&mc_samp_null); *out_mc_samp = mc_samp; return res; error: diff --git a/src/ssol_image.c b/src/ssol_image.c @@ -13,6 +13,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +#define _POSIX_C_SOURCE 200112L /* nextafter support */ + #include "ssol.h" #include "ssol_image_c.h" #include "ssol_device_c.h" @@ -21,11 +23,35 @@ #include <rsys/ref_count.h> #include <rsys/rsys.h> +#include <math.h> #include <string.h> /******************************************************************************* * Helper functions ******************************************************************************/ +static INLINE double +map_address(const double address, const enum ssol_address_mode mode) +{ + double dbl; + double i; + switch(mode) { + case SSOL_ADDRESS_CLAMP: dbl = CLAMP(address, 0, nextafter(1,0)); break; + case SSOL_ADDRESS_REPEAT: + dbl = modf(address, &i); + if(dbl < 0) dbl = 1.0+dbl; + break; + default: FATAL("Unreachable code.\n"); break; + } + return dbl; +} + +static INLINE const char* +get_pixel(const struct ssol_image* img, const size_t x, const size_t y) +{ + ASSERT(img && x < img->size[0] && y < img->size[1]); + return img->mem + y*img->pitch + x*ssol_sizeof_pixel_format(img->format); +} + static void image_release(ref_T* ref) { @@ -136,7 +162,7 @@ ssol_image_get_layout } res_T -ssol_image_map(const struct ssol_image* img, void** mem) +ssol_image_map(const struct ssol_image* img, char** mem) { if(!img || !mem) return RES_BAD_ARG; *mem = img->mem; @@ -151,6 +177,57 @@ res_T ssol_image_unmap(const struct ssol_image* img) } res_T +ssol_image_sample + (const struct ssol_image* img, + const enum ssol_filter_mode filter, + const enum ssol_address_mode address_u, + const enum ssol_address_mode address_v, + const double uv[2], + void* val) +{ + double* z00, *z01, *z10, *z11; + double x0, y0, x1, y1; + double texsz[2]; + double s, t; + double* pix = val; + double integer; + + if(!img || !uv || !val) return RES_BAD_ARG; + + /* Only double3 pixel format is currently supported */ + if(img->format != SSOL_PIXEL_DOUBLE3) return RES_BAD_ARG; + + x0 = map_address(uv[0], address_u) * (double)img->size[0]; + y0 = map_address(uv[1], address_v) * (double)img->size[1]; + + switch(filter) { + case SSOL_FILTER_NEAREST: + z00 = (double*)get_pixel(img, (size_t)x0, (size_t)y0); + pix[0] = z00[0]; + pix[1] = z00[1]; + pix[2] = z00[2]; + break; + case SSOL_FILTER_LINEAR: + texsz[0] = 1.0/(double)img->size[0]; + texsz[1] = 1.0/(double)img->size[1]; + x1 = map_address(uv[0] + texsz[0], address_u) * (double)img->size[0]; + y1 = map_address(uv[1] + texsz[1], address_v) * (double)img->size[1]; + z00 = (double*)get_pixel(img, (size_t)x0, (size_t)y0); + z01 = (double*)get_pixel(img, (size_t)x0, (size_t)y1); + z10 = (double*)get_pixel(img, (size_t)x1, (size_t)y0); + z11 = (double*)get_pixel(img, (size_t)x1, (size_t)y1); + s = modf(x0, &integer); + t = modf(y0, &integer); + pix[0] = (1-s)*((1-t)*z00[0] + t*z01[0]) + s*((1-t)*z10[0] + t*z11[0]); + pix[1] = (1-s)*((1-t)*z00[1] + t*z01[1]) + s*((1-t)*z10[1] + t*z11[1]); + pix[2] = (1-s)*((1-t)*z00[2] + t*z01[2]) + s*((1-t)*z10[2] + t*z11[2]); + break; + default: FATAL("Unreachable code.\n"); break; + } + return RES_OK; +} + +res_T ssol_image_write (void* image, const size_t origin[2], diff --git a/src/ssol_instance.c b/src/ssol_instance.c @@ -203,7 +203,7 @@ ssol_instance_get_area (const struct ssol_instance* instance, double* area) { - if (!instance || !area) return RES_BAD_ARG;; + if (!instance || !area) return RES_BAD_ARG; /* the area of the 3D surface */ *area = instance->shape_rt_area; return RES_OK; diff --git a/src/ssol_material.c b/src/ssol_material.c @@ -15,11 +15,14 @@ #include "ssol.h" #include "ssol_c.h" -#include "ssol_material_c.h" #include "ssol_device_c.h" +#include "ssol_material_c.h" +#include "ssol_spectrum_c.h" -#include <rsys/double3.h> #include <rsys/double2.h> +#include <rsys/double3.h> +#include <rsys/double33.h> +#include <rsys/float2.h> #include <rsys/float3.h> #include <rsys/float33.h> #include <rsys/ref_count.h> @@ -33,10 +36,61 @@ /******************************************************************************* * Helper functions ******************************************************************************/ +/* Define if the submitted ssol_data are *certainly* equals or not. Note that it + * does not check explicitly the spectrum data since it would be too expensive; + * it compares their checksum and that's why one cannot certify that the data + * are strictly equals. Anyway, since this function is used to detect medium + * inconsistencies, it is actually really sufficient to use this strategy. */ +static INLINE int +ssol_data_ceq(const struct ssol_data* a, const struct ssol_data* b) +{ + int i; + ASSERT(a && b); + + if(a->type != b->type) { + i = 0; + } else { + switch(a->type) { + case SSOL_DATA_REAL: + i = a->value.real == b->value.real; + break; + case SSOL_DATA_SPECTRUM: + i = a->value.spectrum->checksum[0] == b->value.spectrum->checksum[0] + && a->value.spectrum->checksum[1] == b->value.spectrum->checksum[1]; + break; + default: FATAL("Unreachable code\n"); break; + } + } + return i; +} + +/* Define if the submitted media are *certainly* equals. Refer to the + * check_ssol_data for more details. */ +static INLINE int +media_ceq(const struct ssol_medium* a, const struct ssol_medium* b) +{ + ASSERT(a && b); + return ssol_data_ceq(&a->refractive_index, &b->refractive_index) + && ssol_data_ceq(&a->absorption, &b->absorption); +} + +static void +shade_normal_default + (struct ssol_device* dev, + struct ssol_param_buffer* buf, + const double wlen, + const struct ssol_surface_fragment* frag, + double* val) /* Returned value */ +{ + ASSERT(frag && val); + (void)dev, (void)buf, (void)wlen; + d3_set(val, frag->Ns); +} + static res_T -dielectric_shade +setup_dielectric_bsdf (const struct ssol_material* mtl, - const struct surface_fragment* fragment, + const struct ssol_surface_fragment* fragment, const double wavelength, /* In nanometer */ const struct ssol_medium* medium, struct ssf_bsdf* bsdf) @@ -44,27 +98,20 @@ dielectric_shade struct ssf_bxdf* brdf = NULL; struct ssf_bxdf* btdf = NULL; struct ssf_fresnel* fresnel = NULL; - const struct ssol_dielectric_shader* shader; double eta_i, eta_t; - double N[3]; res_T res = RES_OK; ASSERT(mtl && fragment && mtl->type == SSOL_MATERIAL_DIELECTRIC); ASSERT(medium && bsdf); + (void)wavelength, (void)fragment; - shader = &mtl->data.dielectric; - - /* Fetch material attribs */ - shader->normal(mtl->dev, mtl->buf, wavelength, fragment->pos, fragment->Ng, - fragment->Ns, fragment->uv, fragment->dir, N); - - if(!MEDIA_EQ(medium, &mtl->out_medium)) { + if(!media_ceq(medium, &mtl->out_medium)) { log_error(mtl->dev, "Inconsistent medium description.\n"); res = RES_BAD_OP; goto error; } - eta_i = mtl->out_medium.refractive_index; - eta_t = mtl->in_medium.refractive_index; + eta_i = ssol_data_get_value(&mtl->out_medium.refractive_index, wavelength); + eta_t = ssol_data_get_value(&mtl->in_medium.refractive_index, wavelength); #define CALL(Func) { res = Func; if(res != RES_OK) goto error; } (void)0 /* Setup the reflective part */ @@ -91,27 +138,21 @@ error: } static res_T -matte_shade +setup_matte_bsdf (const struct ssol_material* mtl, - const struct surface_fragment* fragment, + const struct ssol_surface_fragment* fragment, const double wavelength, /* In nanometer */ struct ssf_bsdf* bsdf) { struct ssf_bxdf* brdf = NULL; - const struct ssol_matte_shader* shader; - double normal[3]; double reflectivity; res_T res; ASSERT(mtl && fragment && mtl->type == SSOL_MATERIAL_MATTE); ASSERT(bsdf); - shader = &mtl->data.matte; - /* Fetch material attribs */ - shader->normal(mtl->dev, mtl->buf, wavelength, fragment->pos, - fragment->Ng, fragment->Ns, fragment->uv, fragment->dir, normal); - shader->reflectivity(mtl->dev, mtl->buf, wavelength, fragment->pos, - fragment->Ng, fragment->Ns, fragment->uv, fragment->dir, &reflectivity); + mtl->data.matte.reflectivity + (mtl->dev, mtl->buf, wavelength, fragment, &reflectivity); /* Setup the BRDF */ res = ssf_bxdf_create(mtl->dev->allocator, &ssf_lambertian_reflection, &brdf); @@ -131,9 +172,9 @@ error: } static res_T -mirror_shade +setup_mirror_bsdf (const struct ssol_material* mtl, - const struct surface_fragment* fragment, + const struct ssol_surface_fragment* fragment, const double wavelength, /* In nanometer */ const int rendering, struct ssf_bsdf* bsdf) @@ -141,23 +182,17 @@ mirror_shade struct ssf_bxdf* brdf = NULL; struct ssf_fresnel* fresnel = NULL; struct ssf_microfacet_distribution* distrib = NULL; - const struct ssol_mirror_shader* shader; - double normal[3]; double roughness; double reflectivity; res_T res; ASSERT(mtl && fragment && mtl->type == SSOL_MATERIAL_MIRROR); ASSERT(bsdf); - shader = &mtl->data.mirror; - /* Fetch material attribs */ - shader->normal(mtl->dev, mtl->buf, wavelength, fragment->pos, - fragment->Ng, fragment->Ns, fragment->uv, fragment->dir, normal); - shader->reflectivity(mtl->dev, mtl->buf, wavelength, fragment->pos, - fragment->Ng, fragment->Ns, fragment->uv, fragment->dir, &reflectivity); - shader->roughness(mtl->dev, mtl->buf, wavelength, fragment->pos, - fragment->Ng, fragment->Ns, fragment->uv, fragment->dir, &roughness); + mtl->data.mirror.reflectivity + (mtl->dev, mtl->buf, wavelength, fragment, &reflectivity); + mtl->data.mirror.roughness + (mtl->dev, mtl->buf, wavelength, fragment, &roughness); /* Setup the fresnel term */ res = ssf_fresnel_create(mtl->dev->allocator, &ssf_fresnel_constant, &fresnel); @@ -207,31 +242,27 @@ error: } static res_T -thin_dielectric_shade +setup_thin_dielectric_bsdf (const struct ssol_material* mtl, - const struct surface_fragment* fragment, + const struct ssol_surface_fragment* fragment, const double wavelength, /* In nanometer */ struct ssf_bsdf* bsdf) { struct ssf_bxdf* bxdf = NULL; - const struct ssol_thin_dielectric_shader* shader; - double N[3]; double thickness; - double absorptivity; + double absorption; double eta_i; double eta_t; res_T res = RES_OK; ASSERT(mtl && fragment && mtl->type == SSOL_MATERIAL_THIN_DIELECTRIC); ASSERT(bsdf); + (void)wavelength, (void)fragment; - shader = &mtl->data.thin_dielectric.shader; - - /* Fetch material attribs */ - shader->normal(mtl->dev, mtl->buf, wavelength, fragment->pos, - fragment->Ng, fragment->Ns, fragment->uv, fragment->dir, N); - eta_i = mtl->out_medium.refractive_index; - eta_t = mtl->data.thin_dielectric.slab_medium.refractive_index; - absorptivity = mtl->data.thin_dielectric.slab_medium.absorptivity; + eta_i = ssol_data_get_value(&mtl->out_medium.refractive_index, wavelength); + eta_t = ssol_data_get_value + (&mtl->data.thin_dielectric.slab_medium.refractive_index, wavelength); + absorption = ssol_data_get_value + (&mtl->data.thin_dielectric.slab_medium.absorption, wavelength); thickness = mtl->data.thin_dielectric.thickness; /* Setup the BxDF */ @@ -239,7 +270,7 @@ thin_dielectric_shade (mtl->dev->allocator, &ssf_thin_specular_dielectric, &bxdf); if(res != RES_OK) goto error; res = ssf_thin_specular_dielectric_setup - (bxdf, absorptivity, eta_i, eta_t, thickness); + (bxdf, absorption, eta_i, eta_t, thickness); if(res != RES_OK) goto error; /* Setup the BSDF */ @@ -253,39 +284,6 @@ error: goto exit; } -static INLINE res_T -shade - (const struct ssol_material* mtl, - const struct surface_fragment* fragment, - const double wavelength, /* In nanometer */ - const int rendering, /* Is material used for rendering */ - const struct ssol_medium* medium, - struct ssf_bsdf* bsdf) -{ - res_T res = RES_OK; - ASSERT(mtl); - - /* Specific material shading */ - switch(mtl->type) { - case SSOL_MATERIAL_DIELECTRIC: - res = dielectric_shade - (mtl, fragment, wavelength, medium, bsdf); - break; - case SSOL_MATERIAL_MATTE: - res = matte_shade(mtl, fragment, wavelength, bsdf); - break; - case SSOL_MATERIAL_MIRROR: - res = mirror_shade(mtl, fragment, wavelength, rendering, bsdf); - break; - case SSOL_MATERIAL_THIN_DIELECTRIC: - res = thin_dielectric_shade(mtl, fragment, wavelength, bsdf); - break; - case SSOL_MATERIAL_VIRTUAL: /* Nothing to shade */ break; - default: FATAL("Unreachable code\n"); break; - } - return res; -} - static INLINE int check_shader_dielectric(const struct ssol_dielectric_shader* shader) { @@ -318,9 +316,38 @@ check_shader_thin_differential(const struct ssol_thin_dielectric_shader* shader) static INLINE int check_medium(const struct ssol_medium* medium) { - return medium - && medium->refractive_index > 0 - && medium->absorptivity >= 0; + if(!medium) return 0; + + /* Check absorption in [0, INF) */ + switch(medium->absorption.type) { + case SSOL_DATA_REAL: + if(medium->absorption.value.real < 0) + return 0; + break; + case SSOL_DATA_SPECTRUM: + if(!medium->absorption.value.spectrum + || !spectrum_check_data(medium->absorption.value.spectrum, 0, DBL_MAX)) + return 0; + break; + default: FATAL("Unreachable code\n"); break; + } + + /* Check refractive index in ]0, INF) */ + switch(medium->refractive_index.type) { + case SSOL_DATA_REAL: + if(medium->refractive_index.value.real <= 0) + return 0; + break; + case SSOL_DATA_SPECTRUM: + if(!medium->refractive_index.value.spectrum + || !spectrum_check_data + (medium->refractive_index.value.spectrum, DBL_EPSILON, DBL_MAX)) + return 0; + break; + default: FATAL("Unreachable code\n"); break; + } + + return 1; } static void @@ -331,6 +358,11 @@ material_release(ref_T* ref) ASSERT(ref); dev = material->dev; if(material->buf) SSOL(param_buffer_ref_put(material->buf)); + if(material->type == SSOL_MATERIAL_THIN_DIELECTRIC) { + ssol_medium_clear(&material->data.thin_dielectric.slab_medium); + } + ssol_medium_clear(&material->in_medium); + ssol_medium_clear(&material->out_medium); ASSERT(dev && dev->allocator); MEM_RM(dev->allocator, material); SSOL(device_ref_put(dev)); @@ -366,6 +398,7 @@ ssol_material_create material->type = type; material->in_medium = SSOL_MEDIUM_VACUUM; material->out_medium = SSOL_MEDIUM_VACUUM; + material->normal = shade_normal_default; exit: if (out_material) *out_material = material; @@ -461,9 +494,10 @@ ssol_dielectric_setup || !check_medium(outside_medium) || !check_medium(inside_medium)) return RES_BAD_ARG; - material->data.dielectric = *shader; - material->out_medium = *outside_medium; - material->in_medium = *inside_medium; + material->data.dielectric.dummy = 1; + ssol_medium_copy(&material->out_medium, outside_medium); + ssol_medium_copy(&material->in_medium, inside_medium); + material->normal = shader->normal; return RES_OK; } @@ -475,7 +509,9 @@ ssol_mirror_setup || material->type != SSOL_MATERIAL_MIRROR || !check_shader_mirror(shader)) return RES_BAD_ARG; - material->data.mirror = *shader; + material->normal = shader->normal; + material->data.mirror.reflectivity = shader->reflectivity; + material->data.mirror.roughness = shader->roughness; return RES_OK; } @@ -487,7 +523,8 @@ ssol_matte_setup || material->type != SSOL_MATERIAL_MATTE || !check_shader_matte(shader)) return RES_BAD_ARG; - material->data.matte = *shader; + material->normal = shader->normal; + material->data.matte.reflectivity = shader->reflectivity; return RES_OK; } @@ -506,11 +543,11 @@ ssol_thin_dielectric_setup || !check_medium(slab_medium) || thickness < 0) return RES_BAD_ARG; - material->data.thin_dielectric.shader = *shader; - material->data.thin_dielectric.slab_medium = *slab_medium; + ssol_medium_copy(&material->data.thin_dielectric.slab_medium, slab_medium); material->data.thin_dielectric.thickness = thickness; - material->out_medium = *outside_medium; - material->in_medium = *outside_medium; + ssol_medium_copy(&material->out_medium, outside_medium); + ssol_medium_copy(&material->in_medium, outside_medium); + material->normal = shader->normal; return RES_OK; } @@ -526,7 +563,7 @@ ssol_material_create_virtual ******************************************************************************/ void surface_fragment_setup - (struct surface_fragment* fragment, + (struct ssol_surface_fragment* fragment, const double pos[3], const double dir[3], const double normal[3], @@ -535,31 +572,68 @@ surface_fragment_setup { struct s3d_attrib attr; char has_texcoord, has_normal; + struct s3d_attrib uvs[3]; + struct s3d_attrib P[3]; + double duv1[2], duv2[2]; + double dP1[3], dP2[3]; + double det; ASSERT(fragment && pos && dir && primitive && uv); /* Assume that the submitted normal look forward the incoming dir */ ASSERT(d3_dot(normal, dir) <= 0); - /* Setup the incoming direction */ - d3_set(fragment->dir, dir); - - /* Setup the surface position */ - d3_set(fragment->pos, pos); + d3_set(fragment->dir, dir); /* Setup the incoming direction */ + d3_set(fragment->P, pos); /* Setup the surface position */ + d3_normalize(fragment->Ng, normal); /* Normalize the geometry normal */ - /* Normalize the geometry normal */ - d3_set(fragment->Ng, normal); - d3_normalize(fragment->Ng, fragment->Ng); + /* Retrieve the position of the triangle vertices */ + S3D(triangle_get_vertex_attrib(primitive, 0, S3D_POSITION, &P[0])); + S3D(triangle_get_vertex_attrib(primitive, 1, S3D_POSITION, &P[1])); + S3D(triangle_get_vertex_attrib(primitive, 2, S3D_POSITION, &P[2])); /* Retrieve the tex coord */ S3D(primitive_has_attrib(primitive, SSOL_TO_S3D_TEXCOORD, &has_texcoord)); if (!has_texcoord) { d2_set_f2(fragment->uv, uv); + uvs[0].type = uvs[1].type = uvs[2].type = S3D_FLOAT2; + uvs[0].usage = uvs[1].usage = uvs[2].usage = SSOL_TO_S3D_TEXCOORD; + f2(uvs[0].value, 1, 0); + f2(uvs[1].value, 0, 1); + f2(uvs[2].value, 0, 0); } else { S3D(primitive_get_attrib(primitive, SSOL_TO_S3D_TEXCOORD, uv, &attr)); + S3D(triangle_get_vertex_attrib(primitive, 0, SSOL_TO_S3D_TEXCOORD, &uvs[0])); + S3D(triangle_get_vertex_attrib(primitive, 1, SSOL_TO_S3D_TEXCOORD, &uvs[1])); + S3D(triangle_get_vertex_attrib(primitive, 2, SSOL_TO_S3D_TEXCOORD, &uvs[2])); ASSERT(attr.type == S3D_FLOAT2); d2_set_f2(fragment->uv, attr.value); } + /* Compute the partial derivatives. */ + duv1[0] = uvs[1].value[0] - uvs[0].value[0]; + duv1[1] = uvs[1].value[1] - uvs[0].value[1]; + duv2[0] = uvs[2].value[0] - uvs[0].value[0]; + duv2[1] = uvs[2].value[1] - uvs[0].value[1]; + dP1[0] = P[1].value[0] - P[0].value[0]; + dP1[1] = P[1].value[1] - P[0].value[1]; + dP1[2] = P[1].value[2] - P[0].value[2]; + dP2[0] = P[2].value[0] - P[0].value[0]; + dP2[1] = P[2].value[1] - P[0].value[1]; + dP2[2] = P[2].value[2] - P[0].value[2]; + det = duv1[0]*duv2[1] - duv1[1]*duv2[0]; + if(det == 0) { /* Handle zero determinant */ + double basis[9]; + d33_basis(basis, fragment->Ng); + d3_set(fragment->dPdu, basis + 0); + d3_set(fragment->dPdv, basis + 3); + } else { + double a[3], b[3]; + d3_sub(fragment->dPdu, d3_muld(a, dP1, duv2[1]), d3_muld(b, dP2, duv1[1])); + d3_sub(fragment->dPdv, d3_muld(a, dP2, duv1[0]), d3_muld(b, dP1, duv2[0])); + d3_divd(fragment->dPdu, fragment->dPdu, det); + d3_divd(fragment->dPdv, fragment->dPdv, det); + } + /* Retrieve and normalize the shading normal in world space */ S3D(primitive_has_attrib(primitive, SSOL_TO_S3D_NORMAL, &has_normal)); if (!has_normal) { @@ -593,26 +667,47 @@ surface_fragment_setup } } -res_T -material_shade +void +material_shade_normal (const struct ssol_material* mtl, - const struct surface_fragment* fragment, - const double wavelength, /* In nanometer */ - const struct ssol_medium* medium, - struct ssf_bsdf* bsdf) + const struct ssol_surface_fragment* frag, + const double wavelength, + double N[3]) { - return shade(mtl, fragment, wavelength, 0, medium, bsdf); + ASSERT(mtl && frag && N); + mtl->normal(mtl->dev, mtl->buf, wavelength, frag, N); } res_T -material_shade_rendering +material_setup_bsdf (const struct ssol_material* mtl, - const struct surface_fragment* fragment, + const struct ssol_surface_fragment* fragment, const double wavelength, /* In nanometer */ const struct ssol_medium* medium, + const int rendering, /* Is BSDF used for rendering */ struct ssf_bsdf* bsdf) { - return shade(mtl, fragment, wavelength, 1, medium, bsdf); + res_T res = RES_OK; + ASSERT(mtl); + + switch(mtl->type) { + case SSOL_MATERIAL_DIELECTRIC: + res = setup_dielectric_bsdf + (mtl, fragment, wavelength, medium, bsdf); + break; + case SSOL_MATERIAL_MATTE: + res = setup_matte_bsdf(mtl, fragment, wavelength, bsdf); + break; + case SSOL_MATERIAL_MIRROR: + res = setup_mirror_bsdf(mtl, fragment, wavelength, rendering, bsdf); + break; + case SSOL_MATERIAL_THIN_DIELECTRIC: + res = setup_thin_dielectric_bsdf(mtl, fragment, wavelength, bsdf); + break; + case SSOL_MATERIAL_VIRTUAL: /* Nothing to shade */ break; + default: FATAL("Unreachable code\n"); break; + } + return res; } res_T @@ -625,17 +720,18 @@ material_get_next_medium switch(mtl->type) { /* The material is an interface between 2 media */ case SSOL_MATERIAL_DIELECTRIC: - if(MEDIA_EQ(&mtl->out_medium, medium)) { - *next_medium = mtl->in_medium; + if(media_ceq(&mtl->out_medium, medium)) { + ssol_medium_copy(next_medium, &mtl->in_medium); } else { - *next_medium = mtl->out_medium; + ASSERT(media_ceq(&mtl->in_medium, medium)); + ssol_medium_copy(next_medium, &mtl->out_medium); } break; /* The material is not an interface between 2 media */ case SSOL_MATERIAL_MATTE: case SSOL_MATERIAL_MIRROR: case SSOL_MATERIAL_THIN_DIELECTRIC: - *next_medium = *medium; + ssol_medium_copy(next_medium, medium); break; default: FATAL("Unreachable code\n"); break; } diff --git a/src/ssol_material_c.h b/src/ssol_material_c.h @@ -23,23 +23,20 @@ struct s3d_primitive; struct ssf_bsdf; struct ssol_device; -#define MEDIA_EQ(A, B) \ - ( ((A)->refractive_index == (B)->refractive_index) \ - && ((A)->absorptivity == (B)->absorptivity)) - -struct surface_fragment { - double dir[3]; /* World space incoming direction */ - double pos[3]; /* World space position */ - double Ng[3]; /* Normalized world space primitive normal */ - double Ns[3]; /* Normalized world space shading normal */ - double uv[2]; /* Texture coordinates */ +struct dielectric { + int dummy; +}; + +struct matte { + ssol_shader_getter_T reflectivity; +}; + +struct mirror { + ssol_shader_getter_T reflectivity; + ssol_shader_getter_T roughness; }; -#define SURFACE_FRAGMENT_NULL__ {{0,0,0}, {0,0,0}, {0,0,0}, {0,0,0}, {0,0}} -static const struct surface_fragment SURFACE_FRAGMENT_NULL = - SURFACE_FRAGMENT_NULL__; struct thin_dielectric { - struct ssol_thin_dielectric_shader shader; struct ssol_medium slab_medium; double thickness; }; @@ -47,10 +44,12 @@ struct thin_dielectric { struct ssol_material { enum ssol_material_type type; + ssol_shader_getter_T normal; + union { - struct ssol_dielectric_shader dielectric; - struct ssol_matte_shader matte; - struct ssol_mirror_shader mirror; + struct dielectric dielectric; + struct matte matte; + struct mirror mirror; struct thin_dielectric thin_dielectric; } data; @@ -64,28 +63,27 @@ struct ssol_material { extern LOCAL_SYM void surface_fragment_setup - (struct surface_fragment* fragment, + (struct ssol_surface_fragment* fragment, const double pos[3], const double dir[3], const double normal[3], const struct s3d_primitive* primitive, const float uv[2]); -extern LOCAL_SYM res_T -material_shade +extern LOCAL_SYM void +material_shade_normal (const struct ssol_material* mtl, - const struct surface_fragment* fragment, - const double wavelength, /* In nanometer */ - const struct ssol_medium* medium, /* Current medium */ - struct ssf_bsdf* bsdf); /* Bidirectional Scattering Distribution Function */ + const struct ssol_surface_fragment* fragment, + const double wavelength, + double N[3]); -/* Material shading for rendering purposes */ extern LOCAL_SYM res_T -material_shade_rendering +material_setup_bsdf (const struct ssol_material* mtl, - const struct surface_fragment* fragment, + const struct ssol_surface_fragment* fragment, const double wavelength, /* In nanometer */ - const struct ssol_medium* medium, + const struct ssol_medium* medium, /* Current medium */ + const int rendering, /* Is material used for rendering purposes */ struct ssf_bsdf* bsdf); /* Bidirectional Scattering Distribution Function */ extern LOCAL_SYM res_T @@ -95,3 +93,4 @@ material_get_next_medium struct ssol_medium* next_medium); #endif /* SSOL_MATERIAL_C_H */ + diff --git a/src/ssol_mc_receiver.c b/src/ssol_mc_receiver.c @@ -17,6 +17,9 @@ #include "ssol_estimator_c.h" #include "ssol_object_c.h" +#include <rsys/double3.h> +#include <star/s3d.h> + #ifdef COMPILER_CL #pragma warning(push) #pragma warning(disable:4706) /* Assignment within a condition */ @@ -56,9 +59,9 @@ ssol_estimator_get_mc_receiver rcv->Name.SE = rcv->Name.V > 0 ? sqrt(rcv->Name.V / N) : 0; \ } (void)0 SETUP_MC_RESULT(integrated_irradiance); + SETUP_MC_RESULT(integrated_absorbed_irradiance); SETUP_MC_RESULT(absorptivity_loss); SETUP_MC_RESULT(reflectivity_loss); - SETUP_MC_RESULT(cos_loss); #undef SETUP_MC_RESULT rcv->mc__ = mc_rcv1; rcv->N__ = estimator->realisation_count; @@ -106,24 +109,55 @@ ssol_mc_shape_get_mc_primitive prim->Name.SE = 0; \ } (void)0 SETUP_MC_RESULT(integrated_irradiance); + SETUP_MC_RESULT(integrated_absorbed_irradiance); SETUP_MC_RESULT(absorptivity_loss); SETUP_MC_RESULT(reflectivity_loss); - SETUP_MC_RESULT(cos_loss); #undef SETUP_MC_RESULT } else { + struct s3d_attrib attr; + struct s3d_shape* s3d_shape; + double v0[3], v1[3], v2[3], E0[3], E1[3], normal[3]; + double area; + unsigned ids[3]; + res_T res = RES_OK; + + s3d_shape = shape->shape__->shape_rt; + + /* Retrieve the primitive indices */ + res = s3d_mesh_get_triangle_indices(s3d_shape, i, ids); + if(res != RES_OK) return res; + + /* Fetch the primitive vertices */ + S3D(mesh_get_vertex_attrib(s3d_shape, ids[0], S3D_POSITION, &attr)); + d3_set_f3(v0, attr.value); + S3D(mesh_get_vertex_attrib(s3d_shape, ids[1], S3D_POSITION, &attr)); + d3_set_f3(v1, attr.value); + S3D(mesh_get_vertex_attrib(s3d_shape, ids[2], S3D_POSITION, &attr)); + d3_set_f3(v2, attr.value); + + /* Compute the primitive area */ + d3_sub(E0, v1, v0); + d3_sub(E1, v2, v0); + d3_cross(normal, E0, E1); + area = d3_len(normal) * 0.5; + #define SETUP_MC_RESULT(Name) { \ const double N = (double)shape->N__; \ const struct mc_data* data = &mc_prim1->Name; \ prim->Name.E = data->weight / N; \ prim->Name.V = data->sqr_weight/N - prim->Name.E*prim->Name.E; \ prim->Name.SE = prim->Name.V > 0 ? sqrt(prim->Name.V / N) : 0; \ + prim->Name.E /= area; \ + prim->Name.V /= area*area; \ + prim->Name.SE /= area; \ } (void)0 SETUP_MC_RESULT(integrated_irradiance); + SETUP_MC_RESULT(integrated_absorbed_irradiance); SETUP_MC_RESULT(absorptivity_loss); SETUP_MC_RESULT(reflectivity_loss); - SETUP_MC_RESULT(cos_loss); #undef SETUP_MC_RESULT } + return RES_OK; } diff --git a/src/ssol_object.c b/src/ssol_object.c @@ -23,6 +23,7 @@ #include <rsys/ref_count.h> #include <rsys/rsys.h> #include <rsys/mem_allocator.h> +#include <rsys/double3.h> /******************************************************************************* * Helper functions @@ -139,13 +140,11 @@ ssol_object_add_shaded_shape res = s3d_scene_attach_shape(object->scn_rt, shape->shape_rt); if(res != RES_OK) goto error; mask |= BIT(ATTACH_S3D_RT); - object->scn_rt_area += shape->shape_rt_area; /* Add the shape samp to the sampling scene of the object */ res = s3d_scene_attach_shape(object->scn_samp, shape->shape_samp); if(res != RES_OK) goto error; mask |= BIT(ATTACH_S3D_SAMP); - object->scn_samp_area += shape->shape_samp_area; /* Ask for a shaded shape identifier */ i = darray_shaded_shape_size_get(&object->shaded_shapes); @@ -164,6 +163,8 @@ ssol_object_add_shaded_shape mask |= BIT(REGISTER_SAMP); /* Setup the object shaded shape */ + object->scn_rt_area += shape->shape_rt_area; + object->scn_samp_area += shape->shape_samp_area; SSOL(shape_ref_get(shape)); SSOL(material_ref_get(front)); SSOL(material_ref_get(back)); @@ -213,6 +214,7 @@ ssol_object_clear(struct ssol_object* obj) htable_shaded_shape_clear(&obj->shaded_shapes_samp); obj->scn_rt_area = 0; + obj->scn_samp_area = 0; S3D(scene_clear(obj->scn_rt)); S3D(scene_clear(obj->scn_samp)); @@ -240,4 +242,3 @@ object_has_shape(struct ssol_object* obj, const struct ssol_shape* shape) S3D(shape_get_id(shape->shape_rt, &id)); return htable_shaded_shape_find(&obj->shaded_shapes_rt, &id) != NULL; } - diff --git a/src/ssol_param_buffer.c b/src/ssol_param_buffer.c @@ -22,15 +22,19 @@ #define DEFAULT_ALIGNMENT 64 struct param { - size_t size; /* In Bytes */ - size_t offset; /* In Bytes */ + void* mem; + void (*release)(void*); }; +#define DARRAY_NAME param +#define DARRAY_DATA struct param +#include <rsys/dynamic_array.h> + struct ssol_param_buffer { char* pool; size_t capacity; size_t size; - + struct darray_param params; ref_T ref; struct ssol_device* dev; }; @@ -45,8 +49,10 @@ param_buffer_release(ref_T* ref) struct ssol_device* dev; ASSERT(ref); buf = CONTAINER_OF(ref, struct ssol_param_buffer, ref); + SSOL(param_buffer_clear(buf)); dev = buf->dev; if(buf->pool) MEM_RM(dev->allocator, buf->pool); + darray_param_release(&buf->params); MEM_RM(dev->allocator, buf); SSOL(device_ref_put(dev)); } @@ -78,6 +84,8 @@ ssol_param_buffer_create ref_init(&buf->ref); buf->capacity = capacity; buf->size = 0; + darray_param_init(dev->allocator, &buf->params); + buf->pool = MEM_ALLOC_ALIGNED(dev->allocator, capacity, DEFAULT_ALIGNMENT); if(!buf->pool) { res = RES_MEM_ERR; @@ -115,25 +123,32 @@ void* ssol_param_buffer_allocate (struct ssol_param_buffer* buf, const size_t size, /* In Bytes */ - const size_t align) /* Must be a power of 2 in [1, 64] */ + const size_t align, /* Must be a power of 2 in [1, 64] */ + void (*release)(void*)) /* May be NULL */ { + struct param param = { NULL, NULL }; size_t offset; - void* mem = NULL; + res_T res = RES_OK; if(!buf || !size || !IS_POW2(align) || align > DEFAULT_ALIGNMENT) goto error; offset = ALIGN_SIZE(buf->size, align); if(offset + size > buf->capacity) goto error; - - mem = buf->pool + offset; - ASSERT(IS_ALIGNED(mem, align)); + + param.mem = buf->pool + offset; + param.release = release; + ASSERT(IS_ALIGNED(param.mem, align)); + + res = darray_param_push_back(&buf->params, &param); + if(res != RES_OK) goto error; + buf->size = offset + size; exit: - return mem; + return param.mem; error: - mem = NULL; + param.mem = NULL; goto exit; } @@ -147,7 +162,15 @@ ssol_param_buffer_get(struct ssol_param_buffer* buf) res_T ssol_param_buffer_clear(struct ssol_param_buffer* buf) { + size_t i; if(!buf) return RES_BAD_ARG; + + /* Release the parameter */ + FOR_EACH(i, 0, darray_param_size_get(&buf->params)) { + struct param* param = darray_param_data_get(&buf->params)+i; + if(param->release) param->release(param->mem); + } + darray_param_clear(&buf->params); buf->size = 0; return RES_OK; } diff --git a/src/ssol_ranst_sun_wl.c b/src/ssol_ranst_sun_wl.c @@ -158,22 +158,24 @@ ranst_sun_wl_setup const size_t sz) { res_T res = RES_OK; - if (!ran || !wavelengths || !intensities || !sz) + if(!ran || !wavelengths || !intensities || !sz) return RES_BAD_ARG; - if (sz > 1) { + + if(sz <= 1) { + ran->type = WL_DIRAC; + ran->get = &ran_dirac_get; + ran->state.dirac.wavelength = sz ? wavelengths[0] : -1; + } else { ran->type = WL_PIECEWISE; ran->get = &ran_piecewise_get; res = ssp_ranst_piecewise_linear_create (ran->allocator, &ran->state.piecewise.spectrum); - if (res != RES_OK) goto error; + if(res != RES_OK) goto error; res = ssp_ranst_piecewise_linear_setup (ran->state.piecewise.spectrum, wavelengths, intensities, sz); - if (res != RES_OK) goto error; - } else { - ran->type = WL_DIRAC; - ran->get = &ran_dirac_get; - ran->state.dirac.wavelength = wavelengths[0]; + if(res != RES_OK) goto error; } + exit: return res; error: diff --git a/src/ssol_scene.c b/src/ssol_scene.c @@ -14,22 +14,23 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "ssol.h" -#include "ssol_c.h" #include "ssol_atmosphere_c.h" -#include "ssol_scene_c.h" -#include "ssol_sun_c.h" +#include "ssol_c.h" #include "ssol_device_c.h" +#include "ssol_instance_c.h" #include "ssol_material_c.h" -#include "ssol_shape_c.h" #include "ssol_object_c.h" -#include "ssol_instance_c.h" +#include "ssol_scene_c.h" +#include "ssol_shape_c.h" +#include "ssol_spectrum_c.h" +#include "ssol_sun_c.h" +#include <rsys/double3.h> +#include <rsys/float2.h> +#include <rsys/float3.h> #include <rsys/list.h> #include <rsys/mem_allocator.h> #include <rsys/rsys.h> -#include <rsys/float2.h> -#include <rsys/float3.h> -#include <rsys/double3.h> /******************************************************************************* * Helper functions @@ -478,10 +479,10 @@ hit_filter_function /* Project the hit position into the punched shape */ d3_set_f3(dir, dirf); d3_set_f3(org, orgf); - dst = punched_shape_trace_ray(sshape->shape, inst->transform, org, dir, - hit->distance, N); + dst = shape_trace_ray(sshape->shape, inst->transform, org, dir, + hit->distance, N, punched_shape_intersect_local); if(dst >= FLT_MAX) { - /* No projection is found => the ray does not intersect the quadric */ + /* No projection found => the ray does not intersect the quadric */ return 1; } if((float)dst <= rdata->range_min) { @@ -514,10 +515,33 @@ hit_filter_function } /* Save the nearest intersected quadric point */ - if(sshape->shape->type == SHAPE_PUNCHED && rdata->dst >= dst) { + if(sshape->shape->type != SHAPE_MESH && rdata->dst >= dst) { d3_set(rdata->N, N); rdata->dst = dst; } return 0; } + +res_T +scene_check(const struct ssol_scene* scene, const char* caller) +{ + ASSERT(scene && caller); + + if(!scene->sun) { + log_error(scene->dev, "%s: no sun attached.\n", caller); + return RES_BAD_ARG; + } + + if(!scene->sun->spectrum) { + log_error(scene->dev, "%s: sun's spectrum undefined.\n", caller); + return RES_BAD_ARG; + } + + if(scene->sun->dni <= 0) { + log_error(scene->dev, "%s: sun's DNI undefined.\n", caller); + return RES_BAD_ARG; + } + return RES_OK; +} + diff --git a/src/ssol_scene_c.h b/src/ssol_scene_c.h @@ -60,5 +60,10 @@ scene_create_s3d_views double* sampled_area, /* Area of the instance set as "samplable" */ double* sampled_area_proxy); /* Area of the sampled geometries */ +extern LOCAL_SYM res_T +scene_check + (const struct ssol_scene* scene, + const char* caller); + #endif /* SSOL_SCENE_C_H */ diff --git a/src/ssol_shape.c b/src/ssol_shape.c @@ -32,6 +32,7 @@ #include <rsys/math.h> #include <star/scpr.h> +#include <star/s3dut.h> #include <limits.h> /* UINT_MAX constant */ #include <math.h> /* copysign function */ @@ -44,38 +45,53 @@ struct mesh_context { struct quadric_mesh_context { const double* coords; const size_t* ids; - const union priv_quadric_data* quadric; + const union private_data* data; + enum ssol_quadric_type quadric_type; const double* transform; /* 3x4 column major matrix */ + double lower[2]; + double upper[2]; +}; + +struct get_ctx { + size_t nbvert; + double two_pi_over_nbvert; + double radius; }; /******************************************************************************* * Helper functions ******************************************************************************/ -static INLINE int +static FINLINE int check_plane(const struct ssol_quadric_plane* plane) { return plane != NULL; } -static INLINE int +static FINLINE int check_parabol(const struct ssol_quadric_parabol* parabol) { return parabol && parabol->focal > 0; } -static INLINE int +static FINLINE int check_hyperbol(const struct ssol_quadric_hyperbol* hyperbol) { return hyperbol && hyperbol->img_focal > 0 && hyperbol->real_focal > 0; } -static INLINE int +static FINLINE int check_parabolic_cylinder (const struct ssol_quadric_parabolic_cylinder* parabolic_cylinder) { return parabolic_cylinder && parabolic_cylinder->focal > 0; } +static FINLINE int +check_hemisphere(const struct ssol_quadric_hemisphere* hemisphere) +{ + return hemisphere && hemisphere->radius > 0; +} + static INLINE int check_quadric(const struct ssol_quadric* quadric) { @@ -90,6 +106,8 @@ check_quadric(const struct ssol_quadric* quadric) return check_hyperbol(&quadric->data.hyperbol); case SSOL_QUADRIC_PARABOLIC_CYLINDER: return check_parabolic_cylinder(&quadric->data.parabolic_cylinder); + case SSOL_QUADRIC_HEMISPHERE: + return check_hemisphere(&quadric->data.hemisphere); default: return 0; } } @@ -108,8 +126,9 @@ check_punched_surface(const struct ssol_punched_surface* punched_surface) size_t i; if(!punched_surface - || punched_surface->nb_carvings == 0 - || !punched_surface->carvings + || (punched_surface->nb_carvings == 0 + && punched_surface->quadric->type != SSOL_QUADRIC_HEMISPHERE) + || (punched_surface->nb_carvings && !punched_surface->carvings) || !punched_surface->quadric || !check_quadric(punched_surface->quadric)) return 0; @@ -163,6 +182,20 @@ mesh_get_pos(const size_t ivert, double pos[2], void* ctx) } static void +quadric_mesh_get_uv(const unsigned ivert, float uv[2], void* ctx) +{ + const size_t i = ivert*2/*#coords per vertex*/; + const struct quadric_mesh_context* msh = ctx; + double tmp[2]; + ASSERT(uv && ctx); + tmp[0] = (msh->coords[i+0] - msh->lower[0]) / (msh->upper[0] - msh->lower[0]); + tmp[1] = (msh->coords[i+1] - msh->lower[1]) / (msh->upper[1] - msh->lower[1]); + + uv[0] = (float)tmp[0]; + uv[1] = (float)tmp[1]; +} + +static void quadric_mesh_get_ids(const unsigned itri, unsigned ids[3], void* ctx) { const size_t i = itri*3/*#ids per triangle*/; @@ -196,9 +229,10 @@ hyperbol_z (const double p[2], const struct priv_hyperbol_data* hyperbol) { - const double z0 = hyperbol->g_2 + hyperbol->abs_b; + const double z0 = hyperbol->g_square + hyperbol->abs_b; const double r2 = p[0] * p[0] + p[1] * p[1]; - return hyperbol->abs_b * sqrt(1 + r2 * hyperbol->_1_a2) + hyperbol->g_2 - z0; + return hyperbol->abs_b * sqrt(1 + r2 * hyperbol->one_over_a_square) + + hyperbol->g_square - z0; } static FINLINE double @@ -207,7 +241,7 @@ parabol_z const struct priv_parabol_data* parabol) { const double r2 = p[0] * p[0] + p[1] * p[1]; - return r2 * parabol->_1_4f; + return r2 * parabol->one_over_4focal; } static FINLINE double @@ -215,7 +249,19 @@ parabolic_cylinder_z (const double p[2], const struct priv_pcylinder_data* pcyl) { - return (p[1] * p[1]) * pcyl->_1_4f; + return (p[1] * p[1]) * pcyl->one_over_4focal; +} + +static FINLINE double +hemisphere_z + (const double p[2], + const struct priv_hemisphere_data* hemisphere) +{ + const double r2 = p[0] * p[0] + p[1] * p[1]; + const double z2 = hemisphere->sqr_radius - r2; + /* manage numerical unaccuracy */ + ASSERT(z2 >= -hemisphere->sqr_radius * FLT_EPSILON); + return (z2 > 0) ? -sqrt(z2) + hemisphere->radius : hemisphere->radius; } static void @@ -227,7 +273,7 @@ quadric_mesh_parabol_get_pos(const unsigned ivert, float pos[3], void* ctx) ASSERT(pos && ctx); p[0] = msh->coords[i+0]; p[1] = msh->coords[i+1]; - p[2] = parabol_z(p, &msh->quadric->parabol); + p[2] = parabol_z(p, &msh->data->parabol); /* Transform the position in object space */ d33_muld3(p, msh->transform, p); @@ -245,7 +291,7 @@ quadric_mesh_hyperbol_get_pos(const unsigned ivert, float pos[3], void* ctx) ASSERT(pos && ctx); p[0] = msh->coords[i+0]; p[1] = msh->coords[i+1]; - p[2] = hyperbol_z(p, &msh->quadric->hyperbol); + p[2] = hyperbol_z(p, &msh->data->hyperbol); /* Transform the position in object space */ d33_muld3(p, msh->transform, p); @@ -264,7 +310,7 @@ quadric_mesh_parabolic_cylinder_get_pos ASSERT(pos && ctx); p[0] = msh->coords[i+0]; p[1] = msh->coords[i+1]; - p[2] = parabolic_cylinder_z(p, &msh->quadric->pcylinder); + p[2] = parabolic_cylinder_z(p, &msh->data->pcylinder); /* Transform the position in object space */ d33_muld3(p, msh->transform, p); @@ -273,6 +319,25 @@ quadric_mesh_parabolic_cylinder_get_pos f3_set_d3(pos, p); } +static void +quadric_mesh_hemisphere_get_pos + (const unsigned ivert, float pos[3], void* ctx) +{ + const size_t i = ivert * 2/*#coords per vertex*/; + const struct quadric_mesh_context* msh = ctx; + double p[3]; /* Temporary quadric space position */ + ASSERT(pos && ctx); + p[0] = msh->coords[i + 0]; + p[1] = msh->coords[i + 1]; + p[2] = hemisphere_z(p, &msh->data->hemisphere); + + /* Transform the position in object space */ + d33_muld3(p, msh->transform, p); + d3_add(p, p, msh->transform + 9); + + f3_set_d3(pos, p); +} + static FINLINE int aabb_is_degenerated(const double lower[2], const double upper[2]) { @@ -307,6 +372,86 @@ carvings_compute_aabb } } +static double +carvings_compute_radius + (const struct ssol_carving* carvings, const size_t ncarvings) +{ + size_t icarving; + double r2 = -DBL_MAX; + ASSERT(carvings); + + if(!ncarvings) return DBL_MAX; + + FOR_EACH(icarving, 0, ncarvings) { + size_t ivert; + FOR_EACH(ivert, 0, carvings[icarving].nb_vertices) { + double pos[2]; + /* Discard the polygons to subtract */ + if (carvings[icarving].operation == SSOL_SUB) continue; + + carvings[icarving].get(ivert, pos, carvings[icarving].context); + r2 = MMAX(r2, d2_dot(pos, pos)); + } + } + return r2 >= 0 ? sqrt(r2) : DBL_MAX; +} + +static res_T +build_triangulated_disk + (struct darray_double* coords, + struct darray_size_t* ids, + const double radius, + const size_t nsteps) +{ + struct s3dut_mesh_data data; + struct s3dut_mesh* mesh = NULL; + double *c_ptr = NULL; + size_t* i_ptr = NULL; + size_t i; + res_T res = RES_OK; + ASSERT(coords && ids && nsteps && radius > 0); + ASSERT(nsteps < UINT_MAX); + + s3dut_create_hemisphere + (coords->allocator, radius, (unsigned)nsteps, (unsigned)nsteps, &mesh); + if (res != RES_OK) { + fprintf(stderr, "Could not create the hemisphere 3D data.\n"); + goto error; + } + + S3DUT(mesh_get_data(mesh, &data)); + if (!data.nprimitives || !data.nvertices) { + res = RES_BAD_ARG; + goto error; + } + + darray_double_clear(coords); + darray_size_t_clear(ids); + + /* Reserve the memory space for the plane vertices */ + res = darray_double_resize(coords, data.nvertices * 2/*#coords per vertex*/); + if (res != RES_OK) goto error; + + /* Reserve the memory space for the plane indices */ + res = darray_size_t_resize(ids, data.nprimitives * 3/*#ids per triangle*/); + if (res != RES_OK) goto error; + + c_ptr = darray_double_data_get(coords); + FOR_EACH(i, 0, data.nvertices) { + d2_set(c_ptr + i * 2, data.positions + i * 3); /* don't get z */ + } + i_ptr = darray_size_t_data_get(ids); + FOR_EACH(i, 0, data.nprimitives * 3) i_ptr[i] = data.indices[i]; + +exit: + if(mesh) S3DUT(mesh_ref_put(mesh)); + return res; +error: + darray_double_clear(coords); + darray_size_t_clear(ids); + goto exit; +} + static res_T build_triangulated_plane (struct darray_double* coords, @@ -322,7 +467,7 @@ build_triangulated_plane double size_min; double delta; res_T res = RES_OK; - ASSERT(coords && lower && upper && nsteps); + ASSERT(coords && ids && lower && upper && nsteps); ASSERT(!aabb_is_degenerated(lower, upper)); darray_double_clear(coords); @@ -394,7 +539,7 @@ error: } static res_T -clip_triangulated_plane +clip_triangulated_sheet (struct darray_double* coords, struct darray_size_t* ids, struct scpr_mesh* mesh, @@ -513,53 +658,66 @@ quadric_setup_s3d_shape_rt (const struct ssol_shape* shape, const struct darray_double* coords, const struct darray_size_t* ids, + const double lower[2], + const double upper[2], struct s3d_shape* s3dshape, double* rt_area) { struct quadric_mesh_context ctx; - struct s3d_vertex_data vdata; + struct s3d_vertex_data vdata[2]; unsigned nverts; unsigned ntris; res_T res; - ASSERT(shape && coords && ids && s3dshape && rt_area); + ASSERT(shape && coords && ids && lower && upper && s3dshape && rt_area); ASSERT(darray_double_size_get(coords)%2 == 0); ASSERT(darray_size_t_size_get(ids)%3 == 0); ASSERT(darray_double_size_get(coords)/2 <= UINT_MAX); ASSERT(darray_size_t_size_get(ids)/3 <= UINT_MAX); + ASSERT(!aabb_is_degenerated(lower, upper)); nverts = (unsigned)darray_double_size_get(coords) / 2/*#coords per vertex*/; ntris = (unsigned)darray_size_t_size_get(ids) / 3/*#ids per triangle*/; ctx.coords = darray_double_cdata_get(coords); ctx.ids = darray_size_t_cdata_get(ids); - ctx.transform = shape->quadric.transform; + ctx.transform = shape->transform; + d2_set(ctx.lower, lower); + d2_set(ctx.upper, upper); + + vdata[0].usage = S3D_POSITION; + vdata[0].type = S3D_FLOAT3; + vdata[0].get = NULL; + + vdata[1].usage = SSOL_TO_S3D_TEXCOORD; + vdata[1].type = S3D_FLOAT2; + vdata[1].get = quadric_mesh_get_uv; - vdata.usage = S3D_POSITION; - vdata.type = S3D_FLOAT3; - vdata.get = NULL; - ctx.quadric = &shape->priv_quadric; - switch (shape->quadric.type) { + ctx.data = &shape->private_data; + ctx.quadric_type = shape->quadric_type; + switch (shape->quadric_type) { case SSOL_QUADRIC_PARABOL: - vdata.get = quadric_mesh_parabol_get_pos; + vdata[0].get = quadric_mesh_parabol_get_pos; break; case SSOL_QUADRIC_HYPERBOL: - vdata.get = quadric_mesh_hyperbol_get_pos; + vdata[0].get = quadric_mesh_hyperbol_get_pos; break; case SSOL_QUADRIC_PARABOLIC_CYLINDER: - vdata.get = quadric_mesh_parabolic_cylinder_get_pos; + vdata[0].get = quadric_mesh_parabolic_cylinder_get_pos; break; case SSOL_QUADRIC_PLANE: - vdata.get = quadric_mesh_plane_get_pos; + vdata[0].get = quadric_mesh_plane_get_pos; + break; + case SSOL_QUADRIC_HEMISPHERE: + vdata[0].get = quadric_mesh_hemisphere_get_pos; break; default: FATAL("Unreachable code.\n"); break; } res = s3d_mesh_setup_indexed_vertices - (s3dshape, ntris, quadric_mesh_get_ids, nverts, &vdata, 1, &ctx); + (s3dshape, ntris, quadric_mesh_get_ids, nverts, vdata, 2, &ctx); if(res != RES_OK) return res; - ASSERT(vdata.get); *rt_area = mesh_compute_area - (ntris, quadric_mesh_get_ids, nverts, vdata.get, &ctx); + (ntris, quadric_mesh_get_ids, nverts, vdata[0].get, &ctx); return RES_OK; } @@ -570,31 +728,41 @@ quadric_setup_s3d_shape_samp (const struct ssol_quadric* quadric, const struct darray_double* coords, const struct darray_size_t* ids, + const double lower[2], + const double upper[2], struct s3d_shape* shape, double *samp_area) { struct quadric_mesh_context ctx; - struct s3d_vertex_data vdata; + struct s3d_vertex_data vdata[2]; unsigned nverts; unsigned ntris; res_T res; - ASSERT(coords && ids && shape); + ASSERT(coords && ids && shape && ids && lower && samp_area); ASSERT(darray_double_size_get(coords)%2 == 0); ASSERT(darray_size_t_size_get(ids)%3 == 0); ASSERT(darray_double_size_get(coords)/2 <= UINT_MAX); ASSERT(darray_size_t_size_get(ids)/3 <= UINT_MAX); + ASSERT(!aabb_is_degenerated(lower, upper)); nverts = (unsigned)darray_double_size_get(coords) / 2/*#coords per vertex*/; ntris = (unsigned)darray_size_t_size_get(ids) / 3/*#ids per triangle*/; ctx.coords = darray_double_cdata_get(coords); ctx.ids = darray_size_t_cdata_get(ids); ctx.transform = quadric->transform; + d2_set(ctx.lower, lower); + d2_set(ctx.upper, upper); + + vdata[0].usage = S3D_POSITION; + vdata[0].type = S3D_FLOAT3; + vdata[0].get = quadric_mesh_plane_get_pos; + + vdata[1].usage = SSOL_TO_S3D_TEXCOORD; + vdata[1].type = S3D_FLOAT2; + vdata[1].get = quadric_mesh_get_uv; - vdata.usage = S3D_POSITION; - vdata.type = S3D_FLOAT3; - vdata.get = quadric_mesh_plane_get_pos; res = s3d_mesh_setup_indexed_vertices - (shape, ntris, quadric_mesh_get_ids, nverts, &vdata, 1, &ctx); + (shape, ntris, quadric_mesh_get_ids, nverts, vdata, 2, &ctx); if(res != RES_OK) return res; *samp_area = mesh_compute_area (ntris, quadric_mesh_get_ids, nverts, quadric_mesh_plane_get_pos, &ctx); @@ -648,36 +816,42 @@ error: } /* Solve a 2nd degree equation. "hint" is used to select among the 2 solutions - * (if applies) the selected solution is then the closest to hint */ + * (if applies) the selected solution is then the closest to hint ans is + * returned in dist[0]. + * If there is a second solution, it is returned in dist[1]. + * Returns the number of roots. */ static int -quadric_solve_second +solve_second (const double a, const double b, const double c, const double hint, - double* dist) + double dist[2]) { - double t = -1; ASSERT(dist && hint >= 0); - if(a == 0) { - if(b != 0) t = -c / b; /* Degenerated case: 1st degree only */ + if(b != 0) { + dist[0] = -c / b; /* Degenerated case: 1st degree only */ + return 1; + } + return 0; /* 0 distance determined */ } else { /* Standard case: 2nd degree */ const double delta = b*b - 4*a*c; if(delta == 0) { - t = -b / (2*a); + dist[0] = -b / (2*a); + return 1; } else { const double sqrt_delta = sqrt(delta); /* Precise formula */ const double t1 = (-b - copysign(sqrt_delta, b)) / (2*a); const double t2 = c / (a*t1); - /* Choose the closest value to hint */ - t = fabs(t1 - hint) < fabs(t2 - hint) ? t1 : t2; + /* dist[0] is the closest value to hint */ + dist[0] = fabs(t1 - hint) < fabs(t2 - hint) ? t1 : t2; + dist[1] = fabs(t1 - hint) < fabs(t2 - hint) ? t2 : t1; + return 2; } } - *dist = t; - return t >= 0; } static FINLINE void @@ -709,10 +883,10 @@ quadric_hyperbol_gradient_local { ASSERT(quad && pt && grad); { - const double z0 = quad->g_2 + quad->abs_b; + const double z0 = quad->g_square + quad->abs_b; grad[0] = pt[0]; grad[1] = pt[1]; - grad[2] = -(pt[2] + z0 - quad->g_2) * quad->_a2_b2; + grad[2] = -(pt[2] + z0 - quad->g_square) * quad->a_square_over_b_square; } } @@ -728,6 +902,18 @@ quadric_parabolic_cylinder_gradient_local grad[2] = 2 * quad->focal; } +static FINLINE void +quadric_hemisphere_gradient_local + (const struct priv_hemisphere_data* quad, + const double pt[3], + double grad[3]) +{ + ASSERT(pt && grad); + grad[0] = -pt[0]; + grad[1] = -pt[1]; + grad[2] = quad->radius - pt[2]; +} + static FINLINE int quadric_plane_intersect_local (const double org[3], @@ -741,13 +927,13 @@ quadric_plane_intersect_local const double a = 0; const double b = dir[2]; const double c = org[2]; - double dst; - int sol = quadric_solve_second(a, b, c, hint, &dst); + double dst[2]; + const int n = solve_second(a, b, c, hint, dst); - if(!sol) return 0; - d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst)); + if(!n) return 0; + d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); quadric_plane_gradient_local(grad); - *dist = dst; + *dist = *dst; return 1; } @@ -762,17 +948,17 @@ quadric_parabol_intersect_local double* dist) /* in/out: */ { /* Define x^2 + y^2 - 4*focal*z = 0 */ - double dst; + double dst[2]; const double a = dir[0] * dir[0] + dir[1] * dir[1]; const double b = 2 * org[0] * dir[0] + 2 * org[1] * dir[1] - 4 * quad->focal * dir[2]; const double c = org[0] * org[0] + org[1] * org[1] - 4 * quad->focal * org[2]; - const int sol = quadric_solve_second(a, b, c, hint, &dst); + const int n = solve_second(a, b, c, hint, dst); - if(!sol) return 0; - d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst)); + if(!n) return 0; + d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); quadric_parabol_gradient_local(quad, hit_pt, grad); - *dist = dst; + *dist = *dst; return 1; } @@ -786,22 +972,49 @@ quadric_hyperbol_intersect_local double grad[3], double* dist) { - double dst; + double dst[2]; const double b2 = quad->abs_b * quad->abs_b; - const double b2_a2 = b2 * quad->_1_a2; - const double z0 = quad->g_2 + quad->abs_b; + const double b2_a2 = b2 * quad->one_over_a_square; + const double z0 = quad->g_square + quad->abs_b; const double a = b2_a2 * (dir[0] * dir[0] + dir[1] * dir[1]) - dir[2] * dir[2]; const double b = - 2 * (b2_a2 * (org[0] * dir[0] + org[1] * dir[1]) - (org[2] + z0 - quad->g_2) * dir[2]); + 2 * (b2_a2 * (org[0] * dir[0] + org[1] * dir[1]) + - (org[2] + z0 - quad->g_square) * dir[2]); const double c = b2_a2 * (org[0] * org[0] + org[1] * org[1]) + b2 - - (org[2] + z0 - quad->g_2) * (org[2] + z0 - quad->g_2); - const int sol = quadric_solve_second(a, b, c, hint, &dst); + - (org[2] + z0 - quad->g_square) * (org[2] + z0 - quad->g_square); + const int n = solve_second(a, b, c, hint, dst); - if(!sol) return 0; - d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst)); + if(!n) return 0; + d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); quadric_hyperbol_gradient_local(quad, hit_pt, grad); - *dist = dst; + *dist = *dst; + return 1; +} + +static FINLINE int +quadric_hemisphere_intersect_local + (const struct priv_hemisphere_data* quad, + const double org[3], + const double dir[3], + const double hint, + double hit_pt[3], + double grad[3], + double* dist) +{ + double dst[2]; + double z0 = -quad->radius; + const double a = dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]; + const double b = 2 * (org[0] * dir[0] + org[1] * dir[1] + org[2] * dir[2] + z0 * dir[2]); + const double c = + org[0] * org[0] + org[1] * org[1] + org[2] * org[2] - quad->sqr_radius + + 2 * z0 * org[2] + z0 * z0; + const int n = solve_second(a, b, c, hint, dst); + + if(!n) return 0; + d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); + quadric_hemisphere_gradient_local(quad, hit_pt, grad); + *dist = *dst; return 1; } @@ -815,15 +1028,16 @@ quadric_parabolic_cylinder_intersect_local double grad[3], double* dist) { - /* Define y^2 - 4 focal z = 0 */ + double dst[2]; const double a = dir[1] * dir[1]; const double b = 2 * org[1] * dir[1] - 4 * quad->focal * dir[2]; const double c = org[1] * org[1] - 4 * quad->focal * org[2]; - const int sol = quadric_solve_second(a, b, c, hint, dist); + const int n = solve_second(a, b, c, hint, dst); - if(!sol) return 0; - d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dist)); + if(!n) return 0; + d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); quadric_parabolic_cylinder_gradient_local(quad, hit_pt, grad); + *dist = *dst; return 1; } @@ -832,18 +1046,21 @@ punched_shape_set_z_local(const struct ssol_shape* shape, double pt[3]) { ASSERT(shape && pt); ASSERT(shape->type == SHAPE_PUNCHED); - switch (shape->quadric.type) { + switch (shape->quadric_type) { case SSOL_QUADRIC_PLANE: pt[2] = 0; break; case SSOL_QUADRIC_PARABOLIC_CYLINDER: - pt[2] = parabolic_cylinder_z(pt, &shape->priv_quadric.pcylinder); + pt[2] = parabolic_cylinder_z(pt, &shape->private_data.pcylinder); break; case SSOL_QUADRIC_PARABOL: - pt[2] = parabol_z(pt, &shape->priv_quadric.parabol); + pt[2] = parabol_z(pt, &shape->private_data.parabol); break; case SSOL_QUADRIC_HYPERBOL: - pt[2] = hyperbol_z(pt, &shape->priv_quadric.hyperbol); + pt[2] = hyperbol_z(pt, &shape->private_data.hyperbol); + break; + case SSOL_QUADRIC_HEMISPHERE: + pt[2] = hemisphere_z(pt, &shape->private_data.hemisphere); break; default: FATAL("Unreachable code\n"); break; } @@ -857,28 +1074,31 @@ punched_shape_set_normal_local { ASSERT(shape && pt); ASSERT(shape->type == SHAPE_PUNCHED); - switch (shape->quadric.type) { + switch (shape->quadric_type) { case SSOL_QUADRIC_PLANE: quadric_plane_gradient_local(normal); break; case SSOL_QUADRIC_PARABOLIC_CYLINDER: quadric_parabolic_cylinder_gradient_local - (&shape->priv_quadric.pcylinder, pt, normal); + (&shape->private_data.pcylinder, pt, normal); break; - case SSOL_QUADRIC_PARABOL: { + case SSOL_QUADRIC_PARABOL: quadric_parabol_gradient_local - (&shape->priv_quadric.parabol, pt, normal); + (&shape->private_data.parabol, pt, normal); break; case SSOL_QUADRIC_HYPERBOL: quadric_hyperbol_gradient_local - (&shape->priv_quadric.hyperbol, pt, normal); + (&shape->private_data.hyperbol, pt, normal); + break; + case SSOL_QUADRIC_HEMISPHERE: + quadric_hemisphere_gradient_local + (&shape->private_data.hemisphere, pt, normal); break; - } default: FATAL("Unreachable code\n"); break; } } -static FINLINE int +int punched_shape_intersect_local (const struct ssol_shape* shape, const double org[3], @@ -894,21 +1114,25 @@ punched_shape_intersect_local ASSERT(dir[0] || dir[1] || dir[2]); /* Hits on quadrics must be recomputed more accurately */ - switch (shape->quadric.type) { + switch (shape->quadric_type) { case SSOL_QUADRIC_PLANE: hit = quadric_plane_intersect_local(org, dir, hint, pt, N, dist); break; case SSOL_QUADRIC_PARABOLIC_CYLINDER: hit = quadric_parabolic_cylinder_intersect_local - (&shape->priv_quadric.pcylinder, org, dir, hint, pt, N, dist); + (&shape->private_data.pcylinder, org, dir, hint, pt, N, dist); break; case SSOL_QUADRIC_PARABOL: hit = quadric_parabol_intersect_local - (&shape->priv_quadric.parabol, org, dir, hint, pt, N, dist); + (&shape->private_data.parabol, org, dir, hint, pt, N, dist); break; case SSOL_QUADRIC_HYPERBOL: hit = quadric_hyperbol_intersect_local - (&shape->priv_quadric.hyperbol, org, dir, hint, pt, N, dist); + (&shape->private_data.hyperbol, org, dir, hint, pt, N, dist); + break; + case SSOL_QUADRIC_HEMISPHERE: + hit = quadric_hemisphere_intersect_local + (&shape->private_data.hemisphere, org, dir, hint, pt, N, dist); break; default: FATAL("Unreachable code\n"); break; } @@ -937,7 +1161,7 @@ priv_parabol_data_setup { ASSERT(data && parabol); data->focal = parabol->focal; - data->_1_4f = 1 / (4.0 * parabol->focal); + data->one_over_4focal = 1 / (4.0 * parabol->focal); } static FINLINE void @@ -953,10 +1177,10 @@ priv_hyperbol_data_setup f = hyperbol->real_focal / g; a2 = g * g * (f - f * f); - data->g_2 = g * 0.5; + data->g_square = g * 0.5; data->abs_b = g * fabs(f - 0.5); - data->_a2_b2 = a2 / (data->abs_b * data->abs_b); - data->_1_a2 = 1 / a2; + data->a_square_over_b_square = a2 / (data->abs_b * data->abs_b); + data->one_over_a_square = 1 / a2; } static FINLINE void @@ -966,12 +1190,22 @@ priv_parabolic_cylinder_data_setup { ASSERT(data && parabolic_cylinder); data->focal = parabolic_cylinder->focal; - data->_1_4f = 1 / (4.0 * parabolic_cylinder->focal); + data->one_over_4focal = 1 / (4.0 * parabolic_cylinder->focal); +} + +static FINLINE void +priv_hemisphere_data_setup + (struct priv_hemisphere_data* data, + const struct ssol_quadric_hemisphere* hemisphere) +{ + ASSERT(data && hemisphere); + data->radius = hemisphere->radius; + data->sqr_radius = hemisphere->radius * hemisphere->radius; } static INLINE void priv_quadric_data_setup - (union priv_quadric_data* priv_data, + (union private_data* priv_data, const struct ssol_quadric* quadric) { ASSERT(priv_data && quadric); @@ -989,6 +1223,10 @@ priv_quadric_data_setup priv_parabolic_cylinder_data_setup (&priv_data->pcylinder, &quadric->data.parabolic_cylinder); break; + case SSOL_QUADRIC_HEMISPHERE: + priv_hemisphere_data_setup + (&priv_data->hemisphere, &quadric->data.hemisphere); + break; default: FATAL("Unreachable code\n"); break; } } @@ -996,9 +1234,9 @@ priv_quadric_data_setup static INLINE size_t priv_quadric_data_compute_slices_count (const enum ssol_quadric_type type, - const union priv_quadric_data* priv_data, - const double lower[3], - const double upper[3]) + const union private_data* priv_data, + const double lower[2], + const double upper[2]) { size_t nslices; double max_z; @@ -1029,6 +1267,17 @@ priv_quadric_data_compute_slices_count return nslices; } +static INLINE size_t +hemisphere_compute_slices_count + (const struct priv_hemisphere_data* hemisphere, const double radius) +{ + size_t nslices; + ASSERT(hemisphere && radius > 0 && hemisphere->radius >= radius); + /* default ranging from 5 to 16 */ + nslices = (size_t)(5.5 + 11 * radius / hemisphere->radius); + return nslices; +} + /******************************************************************************* * Local functions ******************************************************************************/ @@ -1050,8 +1299,8 @@ punched_shape_project_point ASSERT(shape->type == SHAPE_PUNCHED); /* Compute world<->quadric space transformations */ - d33_muld33(R, transform, shape->quadric.transform); - d33_muld3(T, transform, shape->quadric.transform+9); + d33_muld33(R, transform, shape->transform); + d33_muld3(T, transform, shape->transform+9); d3_add(T, T, transform + 9); d33_invtrans(R_invtrans, R); d3_minus(T_inv, T); @@ -1074,13 +1323,14 @@ punched_shape_project_point } double -punched_shape_trace_ray +shape_trace_ray (struct ssol_shape* shape, const double transform[12], /* Shape to world space transformation */ const double org[3], /* World space position near of the ray origin */ const double dir[3], /* World space ray direction */ const double hint_dst, /* Hint on the hit distance */ - double N_quadric[3]) /* World space normal onto the quadric */ + double N_shape[3], /* World space normal onto the shape */ + intersect_local_fn local) /* the intersection function for this shape */ { double R[9]; /* Quadric to world rotation matrix */ double R_invtrans[9]; /* Inverse transpose of R */ @@ -1092,12 +1342,11 @@ punched_shape_trace_ray double N_local[3]; double dst; /* Hit distance */ int valid; - ASSERT(shape && transform && org && N_quadric); - ASSERT(shape->type == SHAPE_PUNCHED); + ASSERT(shape && transform && org && N_shape); /* Compute world<->quadric space transformations */ - d33_muld33(R, transform, shape->quadric.transform); - d33_muld3(T, transform, shape->quadric.transform+9); + d33_muld33(R, transform, shape->transform); + d33_muld3(T, transform, shape->transform+9); d3_add(T, T, transform + 9); d33_invtrans(R_invtrans, R); d3_minus(T_inv, T); @@ -1109,14 +1358,14 @@ punched_shape_trace_ray /* Transform dir in quadric space */ d3_muld33(dir_local, dir, R_invtrans); - /* Project pos_local onto the quadric and compute its associated normal */ - valid = punched_shape_intersect_local + /* Project pos_local onto the shape and compute its associated normal */ + valid = local (shape, org_local, dir_local, hint_dst, hit_local, N_local, &dst); if(!valid) return INF; - /* Transform the quadric normal in world space */ - d33_muld3(N_quadric, R_invtrans, N_local); - d3_normalize(N_quadric, N_quadric); + /* Transform the shape normal in world space */ + d33_muld3(N_shape, R_invtrans, N_local); + d3_normalize(N_shape, N_shape); return dst; } @@ -1209,11 +1458,11 @@ ssol_shape_get_vertex_attrib /* Transform the fetch attrib */ if(shape->type == SHAPE_PUNCHED) { if(usage == SSOL_POSITION) { - d33_muld3(value, shape->quadric.transform, value); - d3_add(value, shape->quadric.transform + 9, value); + d33_muld3(value, shape->transform, value); + d3_add(value, shape->transform + 9, value); } else if(usage == SSOL_NORMAL) { double R_invtrans[9]; - d33_invtrans(R_invtrans, shape->quadric.transform); + d33_invtrans(R_invtrans, shape->transform); d33_muld3(value, R_invtrans, value); } } @@ -1241,6 +1490,7 @@ ssol_punched_surface_setup const struct ssol_punched_surface* psurf) { double lower[2], upper[2]; /* Carvings AABB */ + double radius = -1; struct darray_double coords; struct darray_size_t ids; size_t nslices; @@ -1257,43 +1507,63 @@ ssol_punched_surface_setup } /* Save quadric for further object instancing */ - shape->quadric = *psurf->quadric; - - carvings_compute_aabb(psurf->carvings, psurf->nb_carvings, lower, upper); - if(aabb_is_degenerated(lower, upper)) { - log_error(shape->dev, - "%s: infinite or null punched surface.\n", - FUNC_NAME); - res = RES_BAD_ARG; - goto error; + d33_set(shape->transform, psurf->quadric->transform); + d3_set(shape->transform+9, psurf->quadric->transform+9); + shape->quadric_type = psurf->quadric->type; + + if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) { + radius = carvings_compute_radius(psurf->carvings, psurf->nb_carvings); + radius = MMIN(radius, psurf->quadric->data.hemisphere.radius); + lower[0] = lower[1] = -radius; + upper[0] = upper[1] = +radius; + } else { + carvings_compute_aabb(psurf->carvings, psurf->nb_carvings, lower, upper); + if(aabb_is_degenerated(lower, upper)) { + log_error(shape->dev, + "%s: infinite or null punched surface.\n", + FUNC_NAME); + res = RES_BAD_ARG; + goto error; + } } /* Setup internal data */ - priv_quadric_data_setup(&shape->priv_quadric, psurf->quadric); + priv_quadric_data_setup(&shape->private_data, psurf->quadric); - /* Define the #slices of the discretized quadric */ + /* Define the #slices of the discreet quadric */ if(psurf->quadric->slices_count_hint != SIZE_MAX) { nslices = psurf->quadric->slices_count_hint; + } else if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) { + nslices = hemisphere_compute_slices_count + (&shape->private_data.hemisphere, radius); } else { nslices = priv_quadric_data_compute_slices_count - (shape->quadric.type, &shape->priv_quadric, lower, upper); + (shape->quadric_type, &shape->private_data, lower, upper); } - res = build_triangulated_plane(&coords, &ids, lower, upper, nslices); + /* Build the quadric mesh */ + if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) { + res = build_triangulated_disk(&coords, &ids, radius, nslices); + } else { + res = build_triangulated_plane(&coords, &ids, lower, upper, nslices); + } if(res != RES_OK) goto error; - res = clip_triangulated_plane - (&coords, &ids, shape->dev->scpr_mesh, psurf->carvings, psurf->nb_carvings); - if(res != RES_OK) goto error; + /* Clip the quadric mesh if necessary */ + if(psurf->nb_carvings) { + res = clip_triangulated_sheet + (&coords, &ids, shape->dev->scpr_mesh, psurf->carvings, psurf->nb_carvings); + if(res != RES_OK) goto error; + } /* Setup the Star-3D shape to ray-trace */ - res = quadric_setup_s3d_shape_rt - (shape, &coords, &ids, shape->shape_rt, &shape->shape_rt_area); + res = quadric_setup_s3d_shape_rt(shape, &coords, &ids, lower, upper, + shape->shape_rt, &shape->shape_rt_area); if(res != RES_OK) goto error; /* Setup the Star-3D shape to sample */ - res = quadric_setup_s3d_shape_samp - (psurf->quadric, &coords, &ids, shape->shape_samp, &shape->shape_samp_area); + res = quadric_setup_s3d_shape_samp(psurf->quadric, &coords, &ids, lower, + upper, shape->shape_samp, &shape->shape_samp_area); if(res != RES_OK) goto error; exit: @@ -1363,7 +1633,7 @@ ssol_mesh_setup shape->shape_rt_area = mesh_compute_area(ntris, get_indices, nverts, get_position, data); - /* The Star-3D shape to sample is the same of the one to ray-traced */ + /* The Star-3D shape to sample is the same that the one to ray-trace */ res = s3d_mesh_copy(shape->shape_rt, shape->shape_samp); if(res != RES_OK) goto error; shape->shape_samp_area = shape->shape_rt_area; @@ -1373,4 +1643,3 @@ exit: error: goto exit; } - diff --git a/src/ssol_shape_c.h b/src/ssol_shape_c.h @@ -28,40 +28,56 @@ enum shape_type { struct priv_parabol_data { double focal; - double _1_4f; + double one_over_4focal; }; struct priv_hyperbol_data { - double g_2; - double _a2_b2; - double _1_a2; + double g_square; + double a_square_over_b_square; + double one_over_a_square; double abs_b; }; struct priv_pcylinder_data { double focal; - double _1_4f; + double one_over_4focal; }; -union priv_quadric_data { +struct priv_hemisphere_data { + double radius; + double sqr_radius; +}; + +union private_data { struct priv_hyperbol_data hyperbol; struct priv_parabol_data parabol; struct priv_pcylinder_data pcylinder; + struct priv_hemisphere_data hemisphere; }; struct ssol_shape { enum shape_type type; + enum ssol_quadric_type quadric_type; /* Defined if type is SHAPE_PUNCHED */ struct s3d_shape* shape_rt; /* Star-3D shape to ray-trace */ struct s3d_shape* shape_samp; /* Star-3D shape to sample */ - union priv_quadric_data priv_quadric; - struct ssol_quadric quadric; + union private_data private_data; + double transform[12]; double shape_rt_area, shape_samp_area; struct ssol_device* dev; ref_T ref; }; +typedef int(*intersect_local_fn) + (const struct ssol_shape* shape, + const double org[3], + const double dir[3], + const double hint, + double pt[3], + double N[3], + double* dist); + /* Project pos onto the punched surface and retrieve its associated normal */ extern LOCAL_SYM void punched_shape_project_point @@ -91,5 +107,26 @@ shape_fetched_raw_vertex_attrib const enum ssol_attrib_usage usage, double value[]); +/* Compute ray/punched shape intersection */ +extern LOCAL_SYM int punched_shape_intersect_local + (const struct ssol_shape* shape, + const double org[3], + const double dir[3], + const double hint, + double pt[3], + double N[3], + double* dist); + +/* Compute ray/shape intersection */ +extern LOCAL_SYM double +shape_trace_ray + (struct ssol_shape* shape, + const double transform[12], /* Shape to world space transformation */ + const double org[3], /* World space position near of the ray origin */ + const double dir[3], /* World space ray direction */ + const double hint_dst, /* Hint on the hit distance */ + double N_quadric[3], /* World space normal onto the quadric */ + intersect_local_fn local); + #endif /* SSOL_SHAPE_C_H */ diff --git a/src/ssol_solver.c b/src/ssol_solver.c @@ -25,7 +25,6 @@ #include "ssol_object_c.h" #include "ssol_sun_c.h" #include "ssol_material_c.h" -#include "ssol_spectrum_c.h" #include "ssol_instance_c.h" #include "ssol_ranst_sun_dir.h" #include "ssol_ranst_sun_wl.h" @@ -44,18 +43,25 @@ #include <limits.h> #include <omp.h> +/* How many percent of random walk realisations may fail before an error occurs */ +#define MAX_PERCENT_FAILURES 0.01 + /******************************************************************************* * Thread context ******************************************************************************/ struct thread_context { struct ssp_rng* rng; struct ssf_bsdf* bsdf; + struct mc_data cos_factor; + struct mc_data absorbed; struct mc_data shadowed; struct mc_data missing; - struct mc_data cos_loss; + struct mc_data atmosphere; + struct mc_data reflectivity; struct htable_receiver mc_rcvs; struct htable_sampled mc_samps; struct darray_path paths; /* paths */ + size_t realisation_count; }; static void @@ -100,9 +106,12 @@ thread_context_copy ASSERT(dst && src); dst->rng = src->rng; dst->bsdf = src->bsdf; + dst->cos_factor = src->cos_factor; + dst->absorbed = src->absorbed; dst->shadowed = src->shadowed; dst->missing = src->missing; - dst->cos_loss = src->cos_loss; + dst->atmosphere = src->atmosphere; + dst->reflectivity = src->reflectivity; res = htable_receiver_copy(&dst->mc_rcvs, &src->mc_rcvs); if(res != RES_OK) return res; res = htable_sampled_copy(&dst->mc_samps, &src->mc_samps); @@ -161,10 +170,12 @@ struct point { double dir[3]; float uv[2]; double wl; /* Sampled wavelength */ - double weight; /* actual weight */ - double absorptivity_loss; - double reflectivity_loss; - double cos_loss; + /* MC weights, before and after hit */ + double incoming_weight, weight; + double cos_factor; /* local cos at the starting point */ + double absorbed_irradiance; /* current hit only */ + double absorptivity_loss_before, absorptivity_loss; + double reflectivity_loss_before, reflectivity_loss; enum ssol_side_flag side; }; @@ -178,11 +189,17 @@ struct point { {0, 0, 0}, /* Direction */ \ {0, 0}, /* UV */ \ 0, /* Wavelength */ \ - 0, 0, 0, 0, /* MC weights */ \ + 0, 0, 0, 0, 0, 0, 0, 0, /* MC weights */ \ SSOL_FRONT /* Side */ \ } static const struct point POINT_NULL = POINT_NULL__; +static FINLINE struct ssol_material* +point_get_material(const struct point* pt) +{ + return pt->side == SSOL_FRONT ? pt->sshape->mtl_front : pt->sshape->mtl_back; +} + static res_T point_init (struct point* pt, @@ -196,9 +213,12 @@ point_init struct ssp_rng* rng, int* is_lit) { + struct ssol_surface_fragment frag; struct s3d_attrib attr; struct s3d_hit hit; struct ray_data ray_data = RAY_DATA_NULL; + struct ssol_material* mtl; + double N[3], Np[3]; float dir[3], pos[3], range[2] = { 0, FLT_MAX }; size_t id; res_T res = RES_OK; @@ -232,42 +252,50 @@ point_init /* Sample a sun direction */ ranst_sun_dir_get(ran_sun_dir, rng, pt->dir); - /* Initialise the Monte Carlo weight */ if(pt->sshape->shape->type != SHAPE_PUNCHED) { - double surface_sun_cos = fabs(d3_dot(pt->N, pt->dir)); - pt->weight = scn->sun->dni * sampled_area_proxy * surface_sun_cos; - pt->cos_loss = scn->sun->dni * sampled_area_proxy * (1 - surface_sun_cos); + d3_set(N, pt->N); } else { - double proxy_sun_cos = fabs(d3_dot(pt->N, pt->dir)); - double cos_ratio, surface_proxy_cos, surface_sun_cos, tmp_n[3]; /* For punched surface, retrieve the sampled position and normal onto the * quadric surface */ punched_shape_project_point - (pt->sshape->shape, pt->inst->transform, pt->pos, pt->pos, tmp_n); - surface_proxy_cos = d3_dot(pt->N, tmp_n); - surface_sun_cos = d3_dot(tmp_n, pt->dir); - cos_ratio = fabs(surface_sun_cos / surface_proxy_cos); - d3_set(pt->N, tmp_n); + (pt->sshape->shape, pt->inst->transform, pt->pos, pt->pos, N); + } + + /* Define the primitive side on which the point lies */ + if(d3_dot(N, pt->dir) < 0) { + pt->side = SSOL_FRONT; + } else { + pt->side = SSOL_BACK; + d3_minus(N, N); /* Force the normal to look forward dir */ + } + + /* Perturb the normal */ + surface_fragment_setup(&frag, pt->pos, pt->dir, N, &pt->prim, pt->uv); + mtl = point_get_material(pt); + material_shade_normal(mtl, &frag, pt->wl, Np); + + /* Initialise the Monte Carlo weight */ + if(pt->sshape->shape->type != SHAPE_PUNCHED) { + double surface_sun_cos = fabs(d3_dot(Np, pt->dir)); + pt->weight = scn->sun->dni * sampled_area_proxy * surface_sun_cos; + pt->cos_factor = surface_sun_cos; + } else { + double cos_ratio, surface_proxy_cos, surface_sun_cos; + surface_proxy_cos = fabs(d3_dot(pt->N, Np)); + surface_sun_cos = fabs(d3_dot(Np, pt->dir)); + cos_ratio = surface_sun_cos / surface_proxy_cos; pt->weight = scn->sun->dni * sampled_area_proxy * cos_ratio; - pt->cos_loss = scn->sun->dni * sampled_area_proxy * (1 - proxy_sun_cos); + pt->cos_factor = surface_sun_cos; } + d3_set(pt->N, N); pt->absorptivity_loss = pt->reflectivity_loss = 0; + ASSERT(d3_dot(pt->N, pt->dir) <= 0); /* Store sampled entity related weights */ res = get_mc_sampled(sampled, pt->inst, &pt->mc_samp); if(res != RES_OK) goto error; - pt->mc_samp->cos_loss.weight += pt->cos_loss; - pt->mc_samp->cos_loss.sqr_weight += pt->cos_loss * pt->cos_loss; pt->mc_samp->nb_samples++; - /* Define the primitive side on which the point lies */ - if(d3_dot(pt->N, pt->dir) < 0) { - pt->side = SSOL_FRONT; - } else { - pt->side = SSOL_BACK; - d3_minus(pt->N, pt->N); /* Force the normal to look forward dir */ - } - /* Initialise the ray data to avoid self intersection */ ray_data.scn = scn; ray_data.prim_from = pt->prim; @@ -340,12 +368,6 @@ point_update_from_hit } } -static FINLINE struct ssol_material* -point_get_material(const struct point* pt) -{ - return pt->side == SSOL_FRONT ? pt->sshape->mtl_front : pt->sshape->mtl_back; -} - static FINLINE res_T point_shade (struct point* pt, @@ -355,9 +377,9 @@ point_shade double dir[3]) { struct ssol_material* mtl; - struct surface_fragment frag; + struct ssol_surface_fragment frag; double r = 1; - double wi[3], pdf; + double wi[3], N[3], pdf; int type; res_T res; ASSERT(pt && bsdf && medium && rng && dir); @@ -378,26 +400,51 @@ point_shade /* Shade the surface fragment */ mtl = point_get_material(pt); SSF(bsdf_clear(bsdf)); - res = material_shade(mtl, &frag, pt->wl, medium, bsdf); + res = material_setup_bsdf(mtl, &frag, pt->wl, medium, 0, bsdf); if(res != RES_OK) return res; + /* Perturbe the normal */ + material_shade_normal(mtl, &frag, pt->wl, N); + /* By convention, Star-SF assumes that incoming and reflected * directions point outward the surface => negate incoming dir */ d3_minus(wi, pt->dir); - if(d3_dot(wi, frag.Ns) <= 0) { + if(d3_dot(wi, N) <= 0) { r = 0; } else { - r = ssf_bsdf_sample(bsdf, rng, wi, frag.Ns, dir, &type, &pdf); + double cos_dir_Ng; + r = ssf_bsdf_sample(bsdf, rng, wi, N, dir, &type, &pdf); ASSERT(0 <= r && r <= 1); + + /* Due to the perturbed normal, the sampled direction may point in the + * wrong direction wrt the sampled BSDF component. */ + cos_dir_Ng = d3_dot(frag.Ng, dir); + if((cos_dir_Ng > 0 && (type & SSF_TRANSMISSION)) + || (cos_dir_Ng < 0 && (type & SSF_REFLECTION))) { + r = 0; + } } - pt->reflectivity_loss += (1 - r) * pt->weight; - pt->weight *= r; + pt->incoming_weight = pt->weight; + pt->absorptivity_loss_before = pt->absorptivity_loss; + pt->reflectivity_loss_before = pt->reflectivity_loss; + pt->absorbed_irradiance = (1 - r) * pt->weight; + pt->reflectivity_loss += pt->absorbed_irradiance; + pt->weight = pt->incoming_weight - pt->absorbed_irradiance; if(type & SSF_TRANSMISSION) material_get_next_medium(mtl, medium, medium); return RES_OK; } +static FINLINE void +point_hit_virtual(struct point* pt) +{ + pt->absorbed_irradiance = 0; + pt->incoming_weight = pt->weight; + pt->absorptivity_loss_before = pt->absorptivity_loss; + pt->reflectivity_loss_before = pt->reflectivity_loss; +} + static FINLINE int point_is_receiver(const struct point* pt) { @@ -425,7 +472,6 @@ point_dump if(!stream) return RES_OK; out.realization_id = irealisation; - out.date = 0; /* TODO */ out.segment_id = (uint32_t)isegment; out.receiver_id = point_get_id(pt); out.wavelength = (float)pt->wl; @@ -441,39 +487,6 @@ point_dump /******************************************************************************* * Helper functions ******************************************************************************/ -static FINLINE res_T -check_scene(const struct ssol_scene* scene, const char* caller) -{ - ASSERT(scene && caller); - - if(!scene->sun) { - log_error(scene->dev, "%s: no sun attached.\n", caller); - return RES_BAD_ARG; - } - - if(!scene->sun->spectrum) { - log_error(scene->dev, "%s: sun's spectrum undefined.\n", caller); - return RES_BAD_ARG; - } - - if(scene->sun->dni <= 0) { - log_error(scene->dev, "%s: sun's DNI undefined.\n", caller); - return RES_BAD_ARG; - } - - if(scene->atmosphere) { - int i; - ASSERT(scene->atmosphere->type == ATMOS_UNIFORM); - i = spectrum_includes - (scene->atmosphere->data.uniform.spectrum, scene->sun->spectrum); - if(!i) { - log_error(scene->dev, "%s: sun/atmosphere spectra mismatch.\n", caller); - return RES_BAD_ARG; - } - } - return RES_OK; -} - /* Compute an empirical length of the path segment coming from/going to the * infinite, wrt the scene bounding box */ static INLINE double @@ -518,9 +531,9 @@ accum_mc_receivers_1side dst->Name.sqr_weight += src->Name.sqr_weight; \ } (void)0 ACCUM_WEIGHT(integrated_irradiance); + ACCUM_WEIGHT(integrated_absorbed_irradiance); ACCUM_WEIGHT(absorptivity_loss); ACCUM_WEIGHT(reflectivity_loss); - ACCUM_WEIGHT(cos_loss); #undef ACCUM_WEIGHT /* Merge the per shape MC */ @@ -555,9 +568,9 @@ accum_mc_receivers_1side mc_prim1_dst->Name.sqr_weight += mc_prim1_src->Name.sqr_weight; \ } (void)0 ACCUM_WEIGHT(integrated_irradiance); + ACCUM_WEIGHT(integrated_absorbed_irradiance); ACCUM_WEIGHT(absorptivity_loss); ACCUM_WEIGHT(reflectivity_loss); - ACCUM_WEIGHT(cos_loss); #undef ACCUM_WEIGHT htable_prim2mc_iterator_next(&it_prim); @@ -585,7 +598,6 @@ accum_mc_sampled(struct mc_sampled* dst, struct mc_sampled* src) dst->Name.weight += src->Name.weight; \ dst->Name.sqr_weight += src->Name.sqr_weight; \ } (void)0 - ACCUM_WEIGHT(cos_loss); ACCUM_WEIGHT(shadowed); #undef ACCUM_WEIGHT @@ -627,43 +639,52 @@ update_mc (const struct point* pt, const size_t irealisation, const size_t ibounce, - struct htable_receiver* mc_rcvs, + struct thread_context* thread_ctx, FILE* output) { struct mc_receiver_1side* mc_rcv1 = NULL; struct mc_receiver_1side* mc_samp_x_rcv1 = NULL; res_T res = RES_OK; - ASSERT(pt && mc_rcvs && point_is_receiver(pt)); + ASSERT(pt && thread_ctx && point_is_receiver(pt)); res = point_dump(pt, irealisation, ibounce, output); if(res != RES_OK) goto error; + /* Global MC accumulation */ + #define ACCUM_WEIGHT(Res, W) { \ + Res.weight += (W); \ + Res.sqr_weight += (W)*(W); \ + } (void)0 + ACCUM_WEIGHT(thread_ctx->absorbed, pt->absorbed_irradiance); + #undef ACCUM_WEIGHT + /* Per receiver MC accumulation */ - res = get_mc_receiver_1side(mc_rcvs, pt->inst, pt->side, &mc_rcv1); + res = get_mc_receiver_1side(&thread_ctx->mc_rcvs, pt->inst, pt->side, &mc_rcv1); if(res != RES_OK) goto error; #define ACCUM_WEIGHT(Name, W) { \ mc_rcv1->Name.weight += (W); \ mc_rcv1->Name.sqr_weight += (W)*(W); \ } (void)0 - ACCUM_WEIGHT(integrated_irradiance, pt->weight); - ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss); - ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss); - ACCUM_WEIGHT(cos_loss, pt->cos_loss); + ACCUM_WEIGHT(integrated_irradiance, pt->incoming_weight); + ACCUM_WEIGHT(integrated_absorbed_irradiance, pt->absorbed_irradiance); + ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss_before); + ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss_before); #undef ACCUM_WEIGHT /* Per-sampled/receiver MC accumulation */ res = mc_sampled_get_mc_receiver_1side (pt->mc_samp, pt->inst, pt->side, &mc_samp_x_rcv1); if(res != RES_OK) goto error; + #define ACCUM_WEIGHT(Name, W) { \ mc_samp_x_rcv1->Name.weight += (W); \ mc_samp_x_rcv1->Name.sqr_weight += (W)*(W); \ } (void)0 - ACCUM_WEIGHT(integrated_irradiance, pt->weight); - ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss); - ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss); - ACCUM_WEIGHT(cos_loss, pt->cos_loss); + ACCUM_WEIGHT(integrated_irradiance, pt->incoming_weight); + ACCUM_WEIGHT(integrated_absorbed_irradiance, pt->absorbed_irradiance); + ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss_before); + ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss_before); #undef ACCUM_WEIGHT /* Per primitive receiver MC accumulation */ @@ -677,14 +698,14 @@ update_mc res = mc_shape_1side_get_mc_primitive(mc_shape1, pt->prim.prim_id, &mc_prim1); if(res != RES_OK) goto error; - #define ACCUM_WEIGHT(Name, W) { \ - mc_prim1->Name.weight += (W); \ - mc_prim1->Name.sqr_weight += (W)*(W); \ + #define ACCUM_WEIGHT(Name, W) { \ + mc_prim1->Name.weight += (W); \ + mc_prim1->Name.sqr_weight += (W)*(W); \ } (void)0 - ACCUM_WEIGHT(integrated_irradiance, pt->weight); - ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss); - ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss); - ACCUM_WEIGHT(cos_loss, pt->cos_loss); + ACCUM_WEIGHT(integrated_irradiance, pt->incoming_weight); + ACCUM_WEIGHT(integrated_absorbed_irradiance, pt->absorbed_irradiance); + ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss_before); + ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss_before); #undef ACCUM_WEIGHT } @@ -707,13 +728,14 @@ trace_radiative_path const struct ssol_path_tracker* tracker, /* May be NULL */ FILE* output) /* May be NULL */ { - struct ssol_medium medium = SSOL_MEDIUM_VACUUM; struct path path; + struct ssol_medium medium = SSOL_MEDIUM_VACUUM; struct s3d_hit hit = S3D_HIT_NULL; struct point pt; float org[3], dir[3], range[2] = { 0, FLT_MAX }; size_t depth = 0; int is_lit = 0; + int hit_a_receiver = 0; res_T res = RES_OK; ASSERT(thread_ctx && scn && view_samp && view_rt && ran_sun_dir && ran_sun_wl); @@ -724,13 +746,7 @@ trace_radiative_path view_samp, view_rt, ran_sun_dir, ran_sun_wl, thread_ctx->rng, &is_lit); if(res != RES_OK) goto error; - if(scn->atmosphere) { - /* Assume that the path starts from an uniform atmosphere */ - medium.absorptivity = atmosphere_uniform_get_absorption - (scn->atmosphere, pt.wl); - } - - if(tracker) { +if(tracker) { /* Add the first point of the starting segment */ if(tracker->sun_ray_length > 0) { double pos[3], wi[3]; @@ -746,14 +762,31 @@ trace_radiative_path if(res != RES_OK) goto error; } + #define ACCUM_WEIGHT(Res, W) { \ + Res.weight += (W); \ + Res.sqr_weight += (W)*(W); \ + } (void)0 + ACCUM_WEIGHT(thread_ctx->cos_factor, pt.cos_factor); + ACCUM_WEIGHT(pt.mc_samp->cos_factor, pt.cos_factor); + if(!is_lit) { /* The starting point is not lit */ - pt.mc_samp->shadowed.weight += pt.weight; - pt.mc_samp->shadowed.sqr_weight += pt.weight; - thread_ctx->shadowed.weight += pt.weight; - thread_ctx->shadowed.sqr_weight += pt.weight * pt.weight; + ACCUM_WEIGHT(pt.mc_samp->shadowed, pt.weight); + ACCUM_WEIGHT(thread_ctx->shadowed, pt.weight); if(tracker) path.type = SSOL_PATH_SHADOW; } else { - int hit_a_receiver = 0; + struct ssol_material* mtl = point_get_material(&pt); + + /* Define the medium in which the sampled point lies */ + switch(mtl->type) { + case SSOL_MATERIAL_DIELECTRIC: + case SSOL_MATERIAL_THIN_DIELECTRIC: + ssol_medium_copy(&medium, &mtl->out_medium); + break; + default: + if(scn->atmosphere) + ssol_data_copy(&medium.absorption, &scn->atmosphere->absorption); + break; + } /* Setup the ray as if it starts from the current point position in order * to handle the points that start from a virtual material */ @@ -763,15 +796,29 @@ trace_radiative_path for(;;) { /* Here we go for the radiative random walk */ struct ray_data ray_data = RAY_DATA_NULL; - struct ssol_material* mtl; + double absorption; + + /* Compute interaction with material */ + mtl = point_get_material(&pt); + if(mtl->type == SSOL_MATERIAL_VIRTUAL) { + point_hit_virtual(&pt); + } else { + /* Modulate the point weight wrt its scattering functions and generate + * an outgoing direction */ + res = point_shade(&pt, thread_ctx->bsdf, &medium, thread_ctx->rng, pt.dir); + if(res != RES_OK) goto error; + } if(point_is_receiver(&pt)) { hit_a_receiver = 1; - res = update_mc(&pt, path_id, depth, &thread_ctx->mc_rcvs, output); + res = update_mc(&pt, path_id, depth, thread_ctx, output); if(res != RES_OK) goto error; } - mtl = point_get_material(&pt); + /* Stop the radiative random walk */ + if(pt.weight == 0) break; + + /* Setup new ray parameters */ if(mtl->type == SSOL_MATERIAL_VIRTUAL) { /* Note that for Virtual materials, the ray parameters 'org' & 'dir' * are not updated to ensure that it pursues its traversal without any @@ -779,15 +826,6 @@ trace_radiative_path range[0] = nextafterf(hit.distance, FLT_MAX); range[1] = FLT_MAX; } else { - /* Modulate the point weight wrt to its scattering functions and - * generate an outgoing direction */ - res = point_shade(&pt, thread_ctx->bsdf, &medium, thread_ctx->rng, pt.dir); - if(res != RES_OK) goto error; - - /* Stop the radiative random walk */ - if(pt.weight == 0) break; - - /* Setup new ray parameters */ f2(range, 0, FLT_MAX); f3_set_d3(org, pt.pos); f3_set_d3(dir, pt.dir); @@ -805,12 +843,12 @@ trace_radiative_path S3D(scene_view_trace_ray(view_rt, org, dir, range, &ray_data, &hit)); if(S3D_HIT_NONE(&hit)) { /* Add the point of the last path segment going to the infinite */ - if(tracker && tracker->infinite_ray_length > 0) { + if (tracker && tracker->infinite_ray_length > 0) { double pos[3], wi[3]; d3_set_f3(wi, dir); d3_add(pos, pt.pos, d3_muld(wi, wi, tracker->infinite_ray_length)); res = path_add_vertex(&path, pos, pt.weight); - if(res != RES_OK) goto error; + if (res != RES_OK) goto error; } break; } @@ -818,8 +856,9 @@ trace_radiative_path depth += mtl->type != SSOL_MATERIAL_VIRTUAL; /* Take into account the medium attenuation */ - if(medium.absorptivity > 0 && hit.distance > 0) { - const double transmissivity = exp(-medium.absorptivity * hit.distance); + absorption = ssol_data_get_value(&medium.absorption, pt.wl); + if(absorption > 0 && hit.distance > 0) { + const double transmissivity = exp(-absorption * hit.distance); ASSERT(0 < transmissivity && transmissivity <= 1); pt.absorptivity_loss += (1 - transmissivity) * pt.weight; pt.weight *= transmissivity; @@ -830,13 +869,14 @@ trace_radiative_path if(tracker) { res = path_add_vertex(&path, pt.pos, pt.weight); - if(res != RES_OK) goto error; + if (res != RES_OK) goto error; } } - if(!hit_a_receiver) { - thread_ctx->missing.weight += pt.weight; - thread_ctx->missing.sqr_weight += pt.weight*pt.weight; - } + ACCUM_WEIGHT(thread_ctx->atmosphere, pt.absorptivity_loss); + /* all the remaining weight is lost */ + ACCUM_WEIGHT(thread_ctx->missing, pt.weight); + #undef ACCUM_WEIGHT + if(tracker) { path.type = hit_a_receiver ? SSOL_PATH_SUCCESS : SSOL_PATH_MISSING; } @@ -846,8 +886,8 @@ trace_radiative_path res = path_register_and_clear(&thread_ctx->paths, &path); if(res != RES_OK) goto error; } - exit: + ssol_medium_clear(&medium); if(tracker) path_release(&path); return res; error: @@ -878,10 +918,13 @@ ssol_solve struct ssp_rng_proxy* rng_proxy = NULL; double sampled_area; double sampled_area_proxy; - int nthreads = 0; int64_t nrealisations = 0; + int64_t max_failures = 0; + int nthreads = 0; int i = 0; - ATOMIC res = RES_OK; + ATOMIC mt_res = RES_OK; + ATOMIC nfailures = 0; + res_T res; ASSERT(nrealisations <= INT_MAX); if(!scn || !rng_state || !realisations_count || !out_estimator) @@ -896,17 +939,20 @@ ssol_solve res = RES_BAD_ARG; goto error; } - nrealisations = (int)realisations_count; - nthreads = (int) scn->dev->nthreads; + nrealisations = (int64_t)realisations_count; + max_failures = (int64_t)((double)nrealisations * MAX_PERCENT_FAILURES); + nthreads = (int)scn->dev->nthreads; - res = check_scene(scn, FUNC_NAME); + res = scene_check(scn, FUNC_NAME); if(res != RES_OK) goto error; /* Create data structures shared by all threads */ res = scene_create_s3d_views(scn, &view_rt, &view_samp, &sampled_area, &sampled_area_proxy); if(res != RES_OK) goto error; - res = sun_create_distributions(scn->sun, &ran_sun_dir, &ran_sun_wl); + res = sun_create_direction_distribution(scn->sun, &ran_sun_dir); + if(res != RES_OK) goto error; + res = sun_create_wavelength_distribution(scn->sun, &ran_sun_wl); if(res != RES_OK) goto error; /* Create a RNG proxy from the submitted RNG state */ @@ -945,7 +991,7 @@ ssol_solve const int ithread = omp_get_thread_num(); res_T res_local; - if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurs */ + if(ATOMIC_GET(&mt_res) != RES_OK) continue; /* An error occured */ /* Fetch per thread data */ thread_ctx = darray_thread_ctx_data_get(&thread_ctxs) + ithread; @@ -953,11 +999,19 @@ ssol_solve /* Execute a MC experiment */ res_local = trace_radiative_path((size_t)i, sampled_area_proxy, thread_ctx, scn, view_samp, view_rt, ran_sun_dir, ran_sun_wl, path_tracker, output); - if(res_local != RES_OK) { - ATOMIC_SET(&res, res_local); - continue; + + if(res_local == RES_BAD_OP) { + if(ATOMIC_INCR(&nfailures) >= max_failures) { + log_error(scn->dev, "Too many unexpected radiative paths.\n"); + ATOMIC_SET(&mt_res, res_local); + } + } else if(res_local != RES_OK) { + ATOMIC_SET(&mt_res, res_local); } + if(res_local != RES_OK) continue; + thread_ctx->realisation_count++; } + estimator->failed_count = (size_t)nfailures; /* Merge per thread global MC estimations */ FOR_EACH(i, 0, nthreads) { @@ -967,9 +1021,13 @@ ssol_solve estimator->Name.weight += thread_ctx->Name.weight; \ estimator->Name.sqr_weight += thread_ctx->Name.sqr_weight; \ } (void)0 + ACCUM_WEIGHT(cos_factor); + ACCUM_WEIGHT(absorbed); ACCUM_WEIGHT(shadowed); ACCUM_WEIGHT(missing); - ACCUM_WEIGHT(cos_loss); + ACCUM_WEIGHT(atmosphere); + ACCUM_WEIGHT(reflectivity); + estimator->realisation_count += thread_ctx->realisation_count; #undef ACCUM_WEIGHT } @@ -1038,8 +1096,8 @@ ssol_solve } } - estimator->realisation_count = realisations_count; estimator->sampled_area = sampled_area; + if(mt_res != RES_OK) res = (res_T)mt_res; exit: darray_thread_ctx_release(&thread_ctxs); @@ -1049,7 +1107,7 @@ exit: if(ran_sun_wl) ranst_sun_wl_ref_put(ran_sun_wl); if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); if(out_estimator) *out_estimator = estimator; - return (res_T)res; + return res; error: if(estimator) { SSOL(estimator_ref_put(estimator)); diff --git a/src/ssol_spectrum.c b/src/ssol_spectrum.c @@ -17,11 +17,11 @@ #include "ssol_spectrum_c.h" #include "ssol_device_c.h" -#include <rsys/rsys.h> +#include <rsys/algorithm.h> +#include <rsys/hash.h> +#include <rsys/math.h> #include <rsys/mem_allocator.h> #include <rsys/ref_count.h> -#include <rsys/math.h> -#include <rsys/algorithm.h> /******************************************************************************* * Helper functions @@ -41,20 +41,6 @@ spectrum_release(ref_T* ref) } static int -spectrum_includes_point - (const struct ssol_spectrum* spectrum, - const double wavelength) -{ - const double* data; - size_t sz; - ASSERT(spectrum); - sz = darray_double_size_get(&spectrum->wavelengths); - ASSERT(sz && sz == darray_double_size_get(&spectrum->intensities)); - data = darray_double_cdata_get(&spectrum->wavelengths); - return data[0] <= wavelength && wavelength <= data[sz - 1]; -} - -static int eq_dbl(const void* key, const void* base) { const double k = *(const double*) key; @@ -67,25 +53,9 @@ eq_dbl(const void* key, const void* base) /******************************************************************************* * Local ssol_spectrum functions ******************************************************************************/ -int -spectrum_includes - (const struct ssol_spectrum* reference, - const struct ssol_spectrum* tested) -{ - const double* test_data; - size_t test_sz; - ASSERT(reference && tested); - - test_sz = darray_double_size_get(&tested->wavelengths); - test_data = darray_double_cdata_get(&tested->wavelengths); - return spectrum_includes_point(reference, test_data[0]) - && spectrum_includes_point(reference, test_data[test_sz - 1]); -} - double spectrum_interpolate - (const struct ssol_spectrum* spectrum, - const double wavelength) + (const struct ssol_spectrum* spectrum, const double wavelength) { const double* wls; const double* ints; @@ -93,17 +63,22 @@ spectrum_interpolate double slope; double intensity; size_t id_next, sz; - ASSERT(spectrum && spectrum_includes_point(spectrum, wavelength)); + ASSERT(spectrum); sz = darray_double_size_get(&spectrum->wavelengths); wls = darray_double_cdata_get(&spectrum->wavelengths); ints = darray_double_cdata_get(&spectrum->intensities); next = search_lower_bound(&wavelength, wls, sz, sizeof(double), &eq_dbl); - ASSERT(next); /* because spectrum_includes_point */ + if(!next) { /* Clamp to upper bound */ + return ints[sz-1]; + } id_next = (size_t)(next - wls); + if(!id_next) { /* Clamp to lower bound */ + return ints[0]; + } + ASSERT(id_next); /* because spectrum_includes_point */ - ASSERT(ints[id_next] >= ints[id_next - 1]); ASSERT(wls[id_next] >= wls[id_next - 1]); slope = (ints[id_next] - ints[id_next-1]) / (wls[id_next] - wls[id_next-1]); @@ -112,6 +87,27 @@ spectrum_interpolate return intensity; } +int +spectrum_check_data + (const struct ssol_spectrum* spectrum, const double lower, const double upper) +{ + size_t sz, i; + double current_wl = 0; + ASSERT(spectrum && lower <= upper); + sz = darray_double_size_get(&spectrum->intensities); + if(!sz) return 0; + if(sz != darray_double_size_get(&spectrum->wavelengths)) return 0; + FOR_EACH(i, 0, sz) { + const double wl = darray_double_cdata_get(&spectrum->wavelengths)[i]; + const double data = darray_double_cdata_get(&spectrum->intensities)[i]; + if(data < lower || data > upper) return 0; + if(wl <= 0) return 0; + if(wl <= current_wl) return 0; + current_wl = wl; + } + return 1; +} + /******************************************************************************* * Exported ssol_spectrum functions ******************************************************************************/ @@ -175,6 +171,7 @@ ssol_spectrum_setup { double* wavelengths; double* intensities; + double current_wl = 0; size_t i; res_T res = RES_OK; @@ -190,7 +187,17 @@ ssol_spectrum_setup wavelengths = darray_double_data_get(&spectrum->wavelengths); intensities = darray_double_data_get(&spectrum->intensities); - FOR_EACH(i, 0, nwlens) get(i, wavelengths+i, intensities+i, ctx); + FOR_EACH(i, 0, nwlens) { + get(i, wavelengths + i, intensities + i, ctx); + if(wavelengths[i] <= current_wl || intensities[i] < 0) { + res = RES_BAD_ARG; + goto error; + } + current_wl = *(wavelengths + i); + } + + spectrum->checksum[0] = hash_fnv64(wavelengths, nwlens*sizeof(double)); + spectrum->checksum[1] = hash_fnv64(intensities, nwlens*sizeof(double)); exit: return res; @@ -198,6 +205,8 @@ error: if(spectrum) { darray_double_clear(&spectrum->wavelengths); darray_double_clear(&spectrum->intensities); + spectrum->checksum[0] = 0; + spectrum->checksum[1] = 0; } goto exit; } diff --git a/src/ssol_spectrum_c.h b/src/ssol_spectrum_c.h @@ -22,17 +22,12 @@ struct ssol_spectrum { struct darray_double wavelengths; struct darray_double intensities; + uint64_t checksum[2]; struct ssol_device* dev; ref_T ref; }; -/* Check that the `tested' spectrum is included into `reference' */ -extern LOCAL_SYM int -spectrum_includes - (const struct ssol_spectrum* reference, - const struct ssol_spectrum* tested); - /* Retrieve the linearly interpolated spectrum intensity for the commited * wavelength */ extern LOCAL_SYM double @@ -40,4 +35,10 @@ spectrum_interpolate (const struct ssol_spectrum* spectrum, const double wavelength); +extern LOCAL_SYM int +spectrum_check_data + (const struct ssol_spectrum* spectrum, + const double lower, /* Inclusive lower bound */ + const double upper); /* Inclusive upper bound */ + #endif /* SSOL_SPECTRUM_C_H */ diff --git a/src/ssol_sun.c b/src/ssol_sun.c @@ -39,7 +39,7 @@ sun_release(ref_T* ref) ASSERT(ref); dev = sun->dev; ASSERT(dev && dev->allocator); - if (sun->spectrum) SSOL(spectrum_ref_put(sun->spectrum)); + if(sun->spectrum) SSOL(spectrum_ref_put(sun->spectrum)); MEM_RM(dev->allocator, sun); SSOL(device_ref_put(dev)); } @@ -50,13 +50,13 @@ sun_create { struct ssol_sun* sun = NULL; res_T res = RES_OK; - if (!dev || !out_sun || type >= SUN_TYPES_COUNT__) { + if(!dev || !out_sun || type >= SUN_TYPES_COUNT__) { return RES_BAD_ARG; } sun = (struct ssol_sun*)MEM_CALLOC (dev->allocator, 1, sizeof(struct ssol_sun)); - if (!sun) { + if(!sun) { res = RES_MEM_ERR; goto error; } @@ -67,10 +67,10 @@ sun_create ref_init(&sun->ref); exit: - if (out_sun) *out_sun = sun; + if(out_sun) *out_sun = sun; return res; error: - if (sun) { + if(sun) { SSOL(sun_ref_put(sun)); sun = NULL; } @@ -102,7 +102,7 @@ ssol_sun_create_buie res_T ssol_sun_ref_get(struct ssol_sun* sun) { - if (!sun) + if(!sun) return RES_BAD_ARG; ref_get(&sun->ref); return RES_OK; @@ -111,7 +111,7 @@ ssol_sun_ref_get(struct ssol_sun* sun) res_T ssol_sun_ref_put(struct ssol_sun* sun) { - if (!sun) + if(!sun) return RES_BAD_ARG; ref_put(&sun->ref, sun_release); return RES_OK; @@ -120,9 +120,9 @@ ssol_sun_ref_put(struct ssol_sun* sun) res_T ssol_sun_set_direction(struct ssol_sun* sun, const double direction[3]) { - if (!sun || !direction) + if(!sun || !direction) return RES_BAD_ARG; - if (0 == d3_normalize(sun->direction, direction)) + if(0 == d3_normalize(sun->direction, direction)) /* zero vector */ return RES_BAD_ARG; return RES_OK; @@ -131,7 +131,7 @@ ssol_sun_set_direction(struct ssol_sun* sun, const double direction[3]) res_T ssol_sun_get_direction(const struct ssol_sun* sun, double direction[3]) { - if (!sun || !direction) + if(!sun || !direction) return RES_BAD_ARG; d3_set(direction, sun->direction); return RES_OK; @@ -140,7 +140,7 @@ ssol_sun_get_direction(const struct ssol_sun* sun, double direction[3]) res_T ssol_sun_set_dni(struct ssol_sun* sun, const double dni) { - if (!sun || dni <= 0) + if(!sun || dni <= 0) return RES_BAD_ARG; sun->dni = dni; return RES_OK; @@ -149,7 +149,7 @@ ssol_sun_set_dni(struct ssol_sun* sun, const double dni) res_T ssol_sun_get_dni(const struct ssol_sun* sun, double* dni) { - if (!sun || !dni) + if(!sun || !dni) return RES_BAD_ARG; *dni = sun->dni; return RES_OK; @@ -158,11 +158,11 @@ ssol_sun_get_dni(const struct ssol_sun* sun, double* dni) res_T ssol_sun_set_spectrum(struct ssol_sun* sun, struct ssol_spectrum* spectrum) { - if (!sun || !spectrum) + if(!sun || !spectrum) return RES_BAD_ARG; - if (spectrum == sun->spectrum) /* no change */ + if(spectrum == sun->spectrum) /* no change */ return RES_OK; - if (sun->spectrum) + if(sun->spectrum) SSOL(spectrum_ref_put(sun->spectrum)); SSOL(spectrum_ref_get(spectrum)); sun->spectrum = spectrum; @@ -170,9 +170,7 @@ ssol_sun_set_spectrum(struct ssol_sun* sun, struct ssol_spectrum* spectrum) } res_T -ssol_sun_set_pillbox_aperture - (struct ssol_sun* sun, - const double angle) +ssol_sun_set_pillbox_aperture(struct ssol_sun* sun, const double angle) { if(!sun || angle <= 0 @@ -201,26 +199,13 @@ ssol_sun_set_buie_param * Local function ******************************************************************************/ res_T -sun_create_distributions - (struct ssol_sun* sun, - struct ranst_sun_dir** out_ran_dir, - struct ranst_sun_wl** out_ran_wl) +sun_create_direction_distribution + (struct ssol_sun* sun, struct ranst_sun_dir** out_ran_dir) { struct ranst_sun_dir* ran_dir = NULL; - struct ranst_sun_wl* ran_wl = NULL; res_T res = RES_OK; - ASSERT(sun && out_ran_dir && out_ran_wl); - - /* Create and setup the spectrum distribution */ - res = ranst_sun_wl_create(sun->dev->allocator, &ran_wl); - if(res != RES_OK) goto error; - res = ranst_sun_wl_setup(ran_wl, - darray_double_cdata_get(&sun->spectrum->wavelengths), - darray_double_cdata_get(&sun->spectrum->intensities), - darray_double_size_get(&sun->spectrum->wavelengths)); - if(res != RES_OK) goto error; + ASSERT(sun && out_ran_dir); - /* Create and setup the the direction distribution */ res = ranst_sun_dir_create(sun->dev->allocator, &ran_dir); if(res != RES_OK) goto error; switch(sun->type) { @@ -237,16 +222,10 @@ sun_create_distributions break; default: FATAL("Unreachable code\n"); break; } - exit: *out_ran_dir = ran_dir; - *out_ran_wl = ran_wl; return res; error: - if(ran_wl) { - CHECK(ranst_sun_wl_ref_put(ran_wl), RES_OK); - ran_wl = NULL; - } if(ran_dir) { CHECK(ranst_sun_dir_ref_put(ran_dir), RES_OK); ran_dir = NULL; @@ -254,4 +233,30 @@ error: goto exit; } +res_T +sun_create_wavelength_distribution + (struct ssol_sun* sun, struct ranst_sun_wl** out_ran_wl) +{ + struct ranst_sun_wl* ran_wl = NULL; + res_T res = RES_OK; + ASSERT(sun && out_ran_wl); + res = ranst_sun_wl_create(sun->dev->allocator, &ran_wl); + if(res != RES_OK) goto error; + + res = ranst_sun_wl_setup(ran_wl, + darray_double_cdata_get(&sun->spectrum->wavelengths), + darray_double_cdata_get(&sun->spectrum->intensities), + darray_double_size_get(&sun->spectrum->wavelengths)); + if(res != RES_OK) goto error; + +exit: + *out_ran_wl = ran_wl; + return res; +error: + if(ran_wl) { + CHECK(ranst_sun_wl_ref_put(ran_wl), RES_OK); + ran_wl = NULL; + } + goto exit; +} diff --git a/src/ssol_sun_c.h b/src/ssol_sun_c.h @@ -54,9 +54,13 @@ struct ssol_sun { }; extern LOCAL_SYM res_T -sun_create_distributions +sun_create_direction_distribution + (struct ssol_sun* sun, + struct ranst_sun_dir** out_ran_dir); + +extern LOCAL_SYM res_T +sun_create_wavelength_distribution (struct ssol_sun* sun, - struct ranst_sun_dir** out_ran_dir, struct ranst_sun_wl** out_ran_wl); #endif /* SSOL_SUN_C_H */ diff --git a/src/test_ssol_atmosphere.c b/src/test_ssol_atmosphere.c @@ -16,14 +16,25 @@ #include "ssol.h" #include "test_ssol_utils.h" +static double wl[3] = { 1, 2, 3 }; +static double ab[3] = { 1, 0.8, 1.1 }; + +static void +get_wlen(const size_t i, double* wlen, double* data, void* ctx) +{ + (void)ctx; + *wlen = wl[i]; + *data = ab[i]; +} + int main(int argc, char** argv) { struct mem_allocator allocator; struct ssol_device* dev; - struct ssol_spectrum* spectrum; - struct ssol_spectrum* spectrum2; + struct ssol_data absorption, absorption2; struct ssol_atmosphere* atm; + struct ssol_spectrum* spectrum; (void) argc, (void) argv; mem_init_proxy_allocator(&allocator, &mem_default_allocator); @@ -31,12 +42,15 @@ main(int argc, char** argv) CHECK(ssol_device_create (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev), RES_OK); + absorption2.type = SSOL_DATA_SPECTRUM; CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK); - CHECK(ssol_spectrum_create(dev, &spectrum2), RES_OK); + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 2, NULL), RES_OK); + absorption2.type = SSOL_DATA_SPECTRUM; + absorption2.value.spectrum = spectrum; - CHECK(ssol_atmosphere_create_uniform(NULL, &atm), RES_BAD_ARG); - CHECK(ssol_atmosphere_create_uniform(dev, NULL), RES_BAD_ARG); - CHECK(ssol_atmosphere_create_uniform(dev, &atm), RES_OK); + CHECK(ssol_atmosphere_create(NULL, &atm), RES_BAD_ARG); + CHECK(ssol_atmosphere_create(dev, NULL), RES_BAD_ARG); + CHECK(ssol_atmosphere_create(dev, &atm), RES_OK); CHECK(ssol_atmosphere_ref_get(NULL), RES_BAD_ARG); CHECK(ssol_atmosphere_ref_get(atm), RES_OK); @@ -44,17 +58,23 @@ main(int argc, char** argv) CHECK(ssol_atmosphere_ref_put(NULL), RES_BAD_ARG); CHECK(ssol_atmosphere_ref_put(atm), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(NULL, spectrum), RES_BAD_ARG); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, NULL), RES_BAD_ARG); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, spectrum), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, spectrum2), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, spectrum2), RES_OK); + absorption.type = SSOL_DATA_REAL; + absorption.value.real = 0.1; + CHECK(ssol_atmosphere_set_absorption(NULL, &absorption), RES_BAD_ARG); + CHECK(ssol_atmosphere_set_absorption(atm, NULL), RES_BAD_ARG); + CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_OK); + CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_OK); + CHECK(ssol_atmosphere_set_absorption(atm, &absorption2), RES_OK); - CHECK(ssol_atmosphere_ref_put(atm), RES_OK); + /* absorption values out of range */ + absorption.value.real = 2; + CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_BAD_ARG); + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK); + CHECK(ssol_atmosphere_set_absorption(atm, &absorption2), RES_BAD_ARG); - CHECK(ssol_spectrum_ref_put(spectrum), RES_OK); - CHECK(ssol_spectrum_ref_put(spectrum2), RES_OK); + CHECK(ssol_spectrum_ref_put(absorption2.value.spectrum), RES_OK); CHECK(ssol_device_ref_put(dev), RES_OK); + CHECK(ssol_atmosphere_ref_put(atm), RES_OK); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_ssol_data.c b/src/test_ssol_data.c @@ -0,0 +1,85 @@ +/* Copyright (C) CNRS 2016-2017 + * + * 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 "ssol.h" +#include "test_ssol_utils.h" + +#include <rsys/math.h> + +static void +get_wlen(const size_t i, double* wlen, double* data, void* ctx) +{ + NCHECK(wlen, NULL); + NCHECK(data, NULL); + *wlen = (double)(i+1); + *data = (double)(i+2); + (void)ctx; +} + +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct ssol_device* dev; + struct ssol_spectrum* spectrum; + struct ssol_data data = SSOL_DATA_NULL; + size_t i; + (void)argc, (void)argv; + + CHECK(mem_init_proxy_allocator(&allocator, &mem_default_allocator), RES_OK); + + CHECK(ssol_device_create + (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev), RES_OK); + CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK); + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 10, NULL), RES_OK); + + CHECK(ssol_data_get_value(&data, 1), ssol_data_get_value(&SSOL_DATA_NULL, 1)); + CHECK(ssol_data_get_value(&data, 4), ssol_data_get_value(&SSOL_DATA_NULL, 4)); + CHECK(ssol_data_get_value(&data, 2), ssol_data_get_value(&SSOL_DATA_NULL, 2)); + CHECK(ssol_data_get_value(&data, 7), ssol_data_get_value(&SSOL_DATA_NULL, 7)); + + CHECK(ssol_data_set_real(&data, 1.25), &data); + CHECK(ssol_data_get_value(&data, 1), 1.25); + CHECK(ssol_data_get_value(&data, 1.234), 1.25); + CHECK(data.type, SSOL_DATA_REAL); + CHECK(data.value.real, 1.25); + + CHECK(ssol_data_set_spectrum(&data, spectrum), &data); + CHECK(ssol_data_set_spectrum(&data, spectrum), &data); + + CHECK(data.type, SSOL_DATA_SPECTRUM); + CHECK(data.value.spectrum, spectrum); + CHECK(ssol_spectrum_ref_put(spectrum), RES_OK); + + FOR_EACH(i, 0, 10) { + CHECK(ssol_data_get_value(&data, (double)(i+1)), (double)(i+2)); + } + + CHECK(eq_eps(ssol_data_get_value(&data, 1.5), 2.5, 1.e-6), 1); + CHECK(eq_eps(ssol_data_get_value(&data, 1.25), 2.25, 1.e-6), 1); + CHECK(ssol_data_get_value(&data, 0.5), 2); + CHECK(ssol_data_get_value(&data, 0.1), 2); + CHECK(ssol_data_get_value(&data, 10), 11); + CHECK(ssol_data_get_value(&data, 10.1), 11); + + CHECK(ssol_device_ref_put(dev), RES_OK); + ssol_data_clear(&data); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHECK(mem_allocated_size(), 0); + return 0; +} + diff --git a/src/test_ssol_draw.c b/src/test_ssol_draw.c @@ -30,6 +30,17 @@ #define PITCH (WIDTH*sizeof(unsigned char[3])) #define PROJ_RATIO ((double)WIDTH/(double)HEIGHT) +static void +get_wlen(const size_t i, double* wlen, double* data, void* ctx) +{ + double wavelengths[3] = { 1, 2, 3 }; + double intensities[3] = { 1, 0.8, 1 }; + CHECK(i < 3, 1); + (void)ctx; + *wlen = wavelengths[i]; + *data = intensities[i]; +} + static res_T write_RGB8 (void* data, @@ -175,6 +186,22 @@ setup_cornell_box(struct ssol_device* dev, struct ssol_scene* scn) CHECK(f3_eq_eps(upper, f3(tmp, 552.f, 559.f, 548.f), 1.e-6f), 1); } + +/* Wrap the ssol_draw_pt function to match the ssol_draw_draft profile */ +static INLINE res_T +draw_pt + (struct ssol_scene* scn, + struct ssol_camera* cam, + const size_t width, + const size_t height, + const size_t spp, + ssol_write_pixels_T writer, + void* data) +{ + const double up[3] = {0, 0, 1}; + return ssol_draw_pt(scn, cam, width, height, spp, up, writer, data); +} + int main(int argc, char** argv) { @@ -182,12 +209,15 @@ main(int argc, char** argv) struct ssol_device* dev; struct ssol_camera* cam; struct ssol_scene* scn; + struct ssol_spectrum* spectrum; struct ssol_sun* sun; - unsigned char* pixels = NULL; + struct image img; + uint8_t* pixels; const double pos[3] = {278.0, -1000.0, 273.0}; const double tgt[3] = {278.0, 0.0, 273.0}; const double up[3] = {0.0, 0.0, 1.0}; double dir[3]; + size_t pitch; res_T (*draw_func) (struct ssol_scene* scn, struct ssol_camera* cam, @@ -206,7 +236,7 @@ main(int argc, char** argv) if(!strcmp(argv[1], "draft")) { draw_func = ssol_draw_draft; } else if(!strcmp(argv[1], "pt")) { - draw_func = ssol_draw_pt; + draw_func = draw_pt; } else { fprintf(stderr, "Usage: %s <draft|pt>\n", argv[0]); return -1; @@ -233,8 +263,10 @@ main(int argc, char** argv) CHECK(ssol_sun_set_dni(sun, 1000), RES_OK); CHECK(ssol_scene_attach_sun(scn, sun), RES_OK); - pixels = MEM_CALLOC(&allocator, HEIGHT, PITCH); - NCHECK(pixels, NULL); + pitch = WIDTH * sizeof_image_format(IMAGE_RGB8); + image_init(&allocator, &img); + image_setup(&img, WIDTH, HEIGHT, pitch, IMAGE_RGB8, NULL); + pixels = (uint8_t*)img.pixels; CHECK(draw_func(NULL, NULL, 0, 0, 0, NULL, pixels), RES_BAD_ARG); CHECK(draw_func(scn, NULL, 0, 0, 0, NULL, pixels), RES_BAD_ARG); @@ -299,14 +331,27 @@ main(int argc, char** argv) CHECK(draw_func(NULL, NULL, WIDTH, HEIGHT, 4, write_RGB8, pixels), RES_BAD_ARG); CHECK(draw_func(scn, NULL, WIDTH, HEIGHT, 4, write_RGB8, pixels), RES_BAD_ARG); CHECK(draw_func(NULL, cam, WIDTH, WIDTH, 4, write_RGB8, pixels), RES_BAD_ARG); + + /* No sun spectrum */ + CHECK(draw_func(scn, cam, WIDTH, HEIGHT, 4, write_RGB8, pixels), RES_BAD_ARG); + + CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK); + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK); + CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK); CHECK(draw_func(scn, cam, WIDTH, HEIGHT, 4, write_RGB8, pixels), RES_OK); - CHECK(image_ppm_write_stream(stdout, WIDTH, HEIGHT, 3, pixels), RES_OK); + CHECK(image_write_ppm_stream(&img, 0, stdout), RES_OK); + + if(draw_func == draw_pt) { + CHECK(ssol_draw_pt + (scn, cam, WIDTH, HEIGHT, 4, NULL, write_RGB8, pixels), RES_BAD_ARG); + } - MEM_RM(&allocator, pixels); + CHECK(image_release(&img), RES_OK); CHECK(ssol_device_ref_put(dev), RES_OK); CHECK(ssol_camera_ref_put(cam), RES_OK); CHECK(ssol_scene_ref_put(scn), RES_OK); + CHECK(ssol_spectrum_ref_put(spectrum), RES_OK); CHECK(ssol_sun_ref_put(sun), RES_OK); check_memory_allocator(&allocator); diff --git a/src/test_ssol_image.c b/src/test_ssol_image.c @@ -16,6 +16,94 @@ #include "ssol.h" #include "test_ssol_utils.h" + +#include <rsys/double2.h> +#include <rsys/double3.h> + +static void +check_sampling(struct ssol_device* dev) +{ + struct ssol_image_layout layout; + struct ssol_image* img; + size_t pixsz; + char* mem; + double uv[2]; + double pix[3]; + double tmp[3]; + + CHECK(ssol_image_create(dev, &img), RES_OK); + CHECK(ssol_image_setup(img, 4, 2, SSOL_PIXEL_DOUBLE3), RES_OK); + CHECK(ssol_image_get_layout(img, &layout), RES_OK); + + pixsz = ssol_sizeof_pixel_format(layout.pixel_format); + + CHECK(ssol_image_map(img, &mem), RES_OK); + d3((double*)(mem + layout.offset + 0*pixsz + 0*layout.row_pitch), 1, 0, 0); + d3((double*)(mem + layout.offset + 1*pixsz + 0*layout.row_pitch), 1, 0, 0); + d3((double*)(mem + layout.offset + 2*pixsz + 0*layout.row_pitch), 1, 1, 0); + d3((double*)(mem + layout.offset + 3*pixsz + 0*layout.row_pitch), 1, 1, 0); + d3((double*)(mem + layout.offset + 0*pixsz + 1*layout.row_pitch), 1, 0, 1); + d3((double*)(mem + layout.offset + 1*pixsz + 1*layout.row_pitch), 1, 0, 1); + d3((double*)(mem + layout.offset + 2*pixsz + 1*layout.row_pitch), 0, 1, 1); + d3((double*)(mem + layout.offset + 3*pixsz + 1*layout.row_pitch), 0, 1, 1); + CHECK(ssol_image_unmap(img), RES_OK); + + #define CLAMPED SSOL_ADDRESS_CLAMP + #define REPEAT SSOL_ADDRESS_REPEAT + #define NEAREST SSOL_FILTER_NEAREST + #define LINEAR SSOL_FILTER_LINEAR + + d2_splat(uv, 0); + CHECK(ssol_image_sample(NULL, NEAREST, CLAMPED, CLAMPED, NULL, NULL), RES_BAD_ARG); + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, NULL, NULL), RES_BAD_ARG); + CHECK(ssol_image_sample(NULL, NEAREST, CLAMPED, CLAMPED, uv, NULL), RES_BAD_ARG); + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, NULL), RES_BAD_ARG); + CHECK(ssol_image_sample(NULL, NEAREST, CLAMPED, CLAMPED, NULL, pix), RES_BAD_ARG); + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, NULL, pix), RES_BAD_ARG); + CHECK(ssol_image_sample(NULL, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_BAD_ARG); + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,0,0)), 1); + + uv[0] = 1.0/4.0; + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,0,0)), 1); + uv[0] = 2.0/4.0; + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,1,0)), 1); + uv[0] = 3.0/4.0; + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,1,0)), 1); + + uv[0] = -1.0/4.0; + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,0,0)), 1); + CHECK(ssol_image_sample(img, NEAREST, REPEAT, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,1,0)), 1); + uv[0] = 4.0/4.0; + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,1,0)), 1); + CHECK(ssol_image_sample(img, NEAREST, REPEAT, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,0,0)), 1); + + uv[1] = 1.0/2.0; + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,0,1,1)), 1); + uv[1] = 2.0/2.0; + CHECK(ssol_image_sample(img, NEAREST, REPEAT, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,0,1)), 1); + CHECK(ssol_image_sample(img, NEAREST, REPEAT, REPEAT, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,0,0)), 1); + CHECK(ssol_image_sample(img, NEAREST, CLAMPED, REPEAT, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,1,1,0)), 1); + + uv[0] = 1.0/4.0 + 1.0/8.0; + uv[1] = 0.0/2.0 + 1.0/4.0; + CHECK(ssol_image_sample(img, LINEAR, CLAMPED, CLAMPED, uv, pix), RES_OK); + CHECK(d3_eq(pix, d3(tmp,0.75,0.5,0.5)), 1); + + CHECK(ssol_image_ref_put(img), RES_OK); +} + int main(int argc, char** argv) { @@ -26,7 +114,7 @@ main(int argc, char** argv) struct ssol_image_layout layout = SSOL_IMAGE_LAYOUT_NULL; size_t org[2]; size_t sz[2]; - void* mem; + char* mem; size_t i, x, y; (void) argc, (void) argv; @@ -103,7 +191,7 @@ main(int argc, char** argv) FOR_EACH(y, 0, layout.height) { const double* row = (const double*) - (((char*)mem + layout.offset) + y * layout.row_pitch); + ((mem + layout.offset) + y * layout.row_pitch); FOR_EACH(x, 0, layout.width) { const double* pixel = row + x*3; if(y < 8) { @@ -117,7 +205,7 @@ main(int argc, char** argv) CHECK(pixel[2], 2); } } else { - if(x < 8) { + if(x < 8) { CHECK(pixel[0], 3); CHECK(pixel[1], 3); CHECK(pixel[2], 3); @@ -129,8 +217,12 @@ main(int argc, char** argv) } } } + CHECK(ssol_image_unmap(NULL), RES_BAD_ARG); + CHECK(ssol_image_unmap(img), RES_OK); CHECK(ssol_image_ref_put(img), RES_OK); + + check_sampling(dev); CHECK(ssol_device_ref_put(dev), RES_OK); check_memory_allocator(&allocator); diff --git a/src/test_ssol_instance.c b/src/test_ssol_instance.c @@ -35,7 +35,7 @@ main(int argc, char** argv) struct ssol_vertex_data attrib = SSOL_VERTEX_DATA_NULL; struct ssol_instantiated_shaded_shape sshape; double transform[12] = {1, 0, 0, 0, 1, 0, 0, 0, 1, 10, 0, 0}; - double val[3]; + double val[3], area; size_t n; unsigned i, count; uint32_t id, id1; @@ -79,6 +79,10 @@ main(int argc, char** argv) CHECK(ssol_instance_set_transform(instance, transform), RES_OK); CHECK(ssol_instance_set_transform(instance, transform), RES_OK); + CHECK(ssol_instance_get_area(instance, NULL), RES_BAD_ARG); + CHECK(ssol_instance_get_area(NULL, &area), RES_BAD_ARG); + CHECK(ssol_instance_get_area(instance, &area), RES_OK); + CHECK(ssol_instance_set_receiver(NULL, 0, 0), RES_BAD_ARG); CHECK(ssol_instance_set_receiver(instance, 0, 0), RES_OK); diff --git a/src/test_ssol_material.c b/src/test_ssol_material.c @@ -162,21 +162,21 @@ test_thin_dielectric(struct ssol_device* dev) CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); shader.normal = get_shader_normal; - mdm0.absorptivity = -1; + ssol_data_set_real(&mdm0.absorption, -1); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); - mdm0.absorptivity = SSOL_MEDIUM_VACUUM.absorptivity; + ssol_data_copy(&mdm0.absorption, &SSOL_MEDIUM_VACUUM.absorption); - mdm0.refractive_index = 0; + ssol_data_set_real(&mdm0.refractive_index, 0); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); - mdm0.refractive_index = SSOL_MEDIUM_VACUUM.refractive_index; + ssol_data_copy(&mdm0.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); - mdm1.absorptivity = -1; + ssol_data_set_real(&mdm1.absorption, -1); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); - mdm1.absorptivity = SSOL_MEDIUM_VACUUM.absorptivity; + ssol_data_copy(&mdm1.absorption, &SSOL_MEDIUM_VACUUM.absorption); - mdm1.refractive_index = 0; + ssol_data_set_real(&mdm1.refractive_index, 0); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); - mdm1.refractive_index = SSOL_MEDIUM_VACUUM.refractive_index; + ssol_data_copy(&mdm1.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); CHECK(ssol_material_ref_put(mtl), RES_OK); } @@ -221,21 +221,21 @@ test_dielectric(struct ssol_device* dev) CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); dielectric.normal = get_shader_normal; - mdm0.refractive_index = 0; + ssol_data_set_real(&mdm0.refractive_index, 0); CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); - mdm0.refractive_index = SSOL_MEDIUM_VACUUM.refractive_index; + ssol_data_copy(&mdm0.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); - mdm1.refractive_index = 0; + ssol_data_set_real(&mdm1.refractive_index, 0); CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); - mdm1.refractive_index = SSOL_MEDIUM_VACUUM.refractive_index; + ssol_data_copy(&mdm1.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); - mdm0.absorptivity = -1; + ssol_data_set_real(&mdm0.absorption, -1); CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); - mdm0.absorptivity = SSOL_MEDIUM_VACUUM.refractive_index; + ssol_data_copy(&mdm0.absorption, &SSOL_MEDIUM_VACUUM.refractive_index); - mdm1.absorptivity = -1; + ssol_data_set_real(&mdm1.absorption, -1); CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); - mdm1.refractive_index = SSOL_MEDIUM_VACUUM.refractive_index; + ssol_data_copy(&mdm1.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); CHECK(ssol_material_ref_put(material), RES_OK); } diff --git a/src/test_ssol_materials.h b/src/test_ssol_materials.h @@ -25,16 +25,12 @@ get_shader_normal (struct ssol_device* dev, struct ssol_param_buffer* buf, const double wavelength, - const double P[3], - const double Ng[3], - const double Ns[3], - const double uv[2], - const double w[3], + const struct ssol_surface_fragment* frag, double* val) { int i; - (void)dev, (void)buf, (void)wavelength, (void)P, (void)Ng, (void)uv, (void)w; - FOR_EACH(i, 0, 3) val[i] = Ns[i]; + (void)dev, (void)buf, (void)wavelength; + FOR_EACH(i, 0, 3) val[i] = frag->Ns[i]; } static INLINE void @@ -42,15 +38,10 @@ get_shader_reflectivity (struct ssol_device* dev, struct ssol_param_buffer* buf, const double wavelength, - const double P[3], - const double Ng[3], - const double Ns[3], - const double uv[2], - const double w[3], + const struct ssol_surface_fragment* frag, double* val) { - (void)dev, (void)buf, (void)wavelength; - (void)P, (void)Ng, (void)Ns, (void)uv, (void) w; + (void)dev, (void)buf, (void)wavelength, (void)frag; *val = 1; } @@ -60,15 +51,10 @@ get_shader_reflectivity_2 (struct ssol_device* dev, struct ssol_param_buffer* buf, const double wavelength, - const double P[3], - const double Ng[3], - const double Ns[3], - const double uv[2], - const double w[3], + const struct ssol_surface_fragment* frag, double* val) { - (void) dev, (void) buf, (void) wavelength; - (void) P, (void) Ng, (void) Ns, (void) uv, (void) w; + (void)dev, (void)buf, (void)wavelength, (void)frag; *val = REFLECTIVITY; } #endif @@ -78,15 +64,10 @@ get_shader_roughness (struct ssol_device* dev, struct ssol_param_buffer* buf, const double wavelength, - const double P[3], - const double Ng[3], - const double Ns[3], - const double uv[2], - const double w[3], + const struct ssol_surface_fragment* frag, double* val) { - (void)dev, (void)buf, (void)wavelength; - (void)P, (void)Ng, (void)Ns, (void)uv, (void) w; + (void)dev, (void)buf, (void)wavelength, (void)frag; *val = 0; } @@ -95,15 +76,10 @@ get_shader_absorption (struct ssol_device* dev, struct ssol_param_buffer* buf, const double wavelength, - const double P[3], - const double Ng[3], - const double Ns[3], - const double uv[2], - const double w[3], + const struct ssol_surface_fragment* frag, double* val) { - (void)dev, (void)buf, (void)wavelength; - (void)P, (void)Ng, (void)Ns, (void)uv, (void) w; + (void)dev, (void)buf, (void)wavelength, (void)frag; *val = 0; } @@ -112,15 +88,10 @@ get_shader_thickness (struct ssol_device* dev, struct ssol_param_buffer* buf, const double wavelength, - const double P[3], - const double Ng[3], - const double Ns[3], - const double uv[2], - const double w[3], + const struct ssol_surface_fragment* frag, double* val) { - (void)dev, (void)buf, (void)wavelength; - (void)P, (void)Ng, (void)Ns, (void)uv, (void) w; + (void)dev, (void)buf, (void)wavelength, (void)frag; *val = 1; } @@ -129,15 +100,10 @@ get_shader_refractive_index (struct ssol_device* dev, struct ssol_param_buffer* buf, const double wavelength, - const double P[3], - const double Ng[3], - const double Ns[3], - const double uv[2], - const double w[3], + const struct ssol_surface_fragment* frag, double* val) { - (void)dev, (void)buf, (void)wavelength; - (void)P, (void)Ng, (void)Ns, (void)uv, (void) w; + (void)dev, (void)buf, (void)wavelength, (void)frag; *val = 1.5; } diff --git a/src/test_ssol_object.c b/src/test_ssol_object.c @@ -26,6 +26,7 @@ main(int argc, char** argv) struct ssol_material* mtl; struct ssol_material* mtl2; struct ssol_object* object; + double a; (void) argc, (void) argv; mem_init_proxy_allocator(&allocator, &mem_default_allocator); @@ -73,6 +74,10 @@ main(int argc, char** argv) CHECK(ssol_object_clear(NULL), RES_BAD_ARG); CHECK(ssol_object_clear(object), RES_OK); + CHECK(ssol_object_get_area(object, NULL), RES_BAD_ARG); + CHECK(ssol_object_get_area(NULL, &a), RES_BAD_ARG); + CHECK(ssol_object_get_area(object, &a), RES_OK); + CHECK(ssol_object_ref_put(object), RES_OK); CHECK(ssol_shape_ref_put(shape), RES_OK); CHECK(ssol_shape_ref_put(shape2), RES_OK); diff --git a/src/test_ssol_param_buffer.c b/src/test_ssol_param_buffer.c @@ -19,18 +19,30 @@ #include <rsys/logger.h> #include <limits.h> +struct param { + char* name; + double d; + int i; + void* ptr; + struct ssol_image* img; +}; + +static void +param_release(void* mem) +{ + struct param* param = mem; + ASSERT(param); + if(param->img) SSOL(image_ref_put(param->img)); +} + int main(int argc, char** argv) { - struct param { - char* name; - double d; - int i; - void* ptr; - }* param; + struct param* param; struct mem_allocator allocator; struct ssol_device* dev; struct ssol_param_buffer* pbuf; + struct ssol_image* img; size_t sz, al; void* mem; (void)argc, (void)argv; @@ -53,14 +65,14 @@ main(int argc, char** argv) sz = sizeof(intptr_t); al = ALIGNOF(intptr_t); - CHECK(mem = ssol_param_buffer_allocate(NULL, 0, 0), NULL); - CHECK(mem = ssol_param_buffer_allocate(pbuf, 0, 0), NULL); - CHECK(mem = ssol_param_buffer_allocate(NULL, sz, 0), NULL); - CHECK(mem = ssol_param_buffer_allocate(pbuf, sz, 0), NULL); - CHECK(mem = ssol_param_buffer_allocate(NULL, 0, al), NULL); - CHECK(mem = ssol_param_buffer_allocate(pbuf, 0, al), NULL); - CHECK(mem = ssol_param_buffer_allocate(NULL, sz, al), NULL); - NCHECK(mem = ssol_param_buffer_allocate(pbuf, sz, al), NULL); + CHECK(mem = ssol_param_buffer_allocate(NULL, 0, 0, NULL), NULL); + CHECK(mem = ssol_param_buffer_allocate(pbuf, 0, 0, NULL), NULL); + CHECK(mem = ssol_param_buffer_allocate(NULL, sz, 0, NULL), NULL); + CHECK(mem = ssol_param_buffer_allocate(pbuf, sz, 0, NULL), NULL); + CHECK(mem = ssol_param_buffer_allocate(NULL, 0, al, NULL), NULL); + CHECK(mem = ssol_param_buffer_allocate(pbuf, 0, al, NULL), NULL); + CHECK(mem = ssol_param_buffer_allocate(NULL, sz, al, NULL), NULL); + NCHECK(mem = ssol_param_buffer_allocate(pbuf, sz, al, NULL), NULL); *(intptr_t*)mem = 0xDECAFBAD; CHECK(*(intptr_t*)ssol_param_buffer_get(pbuf), 0xDECAFBAD); @@ -74,7 +86,7 @@ main(int argc, char** argv) sz = strlen("Foo") + 1; al = 4; - NCHECK(mem = ssol_param_buffer_allocate(pbuf, sz, al), NULL); + NCHECK(mem = ssol_param_buffer_allocate(pbuf, sz, al, NULL), NULL); strcpy(mem, "Foo"); CHECK(strcmp(ssol_param_buffer_get(pbuf), "Foo"), 0); strcpy(mem, "Bar"); @@ -85,12 +97,13 @@ main(int argc, char** argv) sz = sizeof(struct param); al = ALIGNOF(struct param); - NCHECK(param = ssol_param_buffer_allocate(pbuf, sz, al), NULL); - NCHECK(param->name = ssol_param_buffer_allocate(pbuf, 7, 64), NULL); + NCHECK(param = ssol_param_buffer_allocate(pbuf, sz, al, NULL), NULL); + NCHECK(param->name = ssol_param_buffer_allocate(pbuf, 7, 64, NULL), NULL); strcpy(param->name, "0123456"); - NCHECK(param->ptr = ssol_param_buffer_allocate(pbuf, 4, 16), NULL); + NCHECK(param->ptr = ssol_param_buffer_allocate(pbuf, 4, 16, NULL), NULL); param->d = PI; param->i = -123; + param->img = NULL; strcpy(param->ptr, "abc"); NCHECK(param = ssol_param_buffer_get(pbuf), NULL); @@ -102,6 +115,46 @@ main(int argc, char** argv) CHECK(strcmp(param->name, "0123456"), 0); CHECK(strcmp(param->ptr, "abc"), 0); + CHECK(ssol_param_buffer_clear(pbuf), RES_OK); + + sz = sizeof(struct param); + al = ALIGNOF(struct param); + CHECK(ssol_image_create(dev, &img), RES_OK); + CHECK(ssol_image_setup(img, 1280, 720, SSOL_PIXEL_DOUBLE3), RES_OK); + NCHECK(param = ssol_param_buffer_allocate(pbuf, sz, al, &param_release), NULL); + param->d = PI; + param->i = -123; + param->name = NULL; + param->ptr = NULL; + param->img = img; + CHECK(ssol_image_ref_get(img), RES_OK); + + NCHECK(param = ssol_param_buffer_allocate(pbuf, sz, al, &param_release), NULL); + param->d = 123.456; + param->i = -1; + param->name = NULL; + param->ptr = NULL; + param->img = img; + CHECK(ssol_image_ref_get(img), RES_OK); + + NCHECK(param = ssol_param_buffer_allocate(pbuf, sz, al, &param_release), NULL); + param->d = 0.1; + param->i = 789; + param->name = NULL; + param->ptr = NULL; + param->img = img; + CHECK(ssol_image_ref_get(img), RES_OK); + + CHECK(ssol_param_buffer_clear(pbuf), RES_OK); + + NCHECK(param = ssol_param_buffer_allocate(pbuf, sz, al, &param_release), NULL); + param->d = 0.1; + param->i = 789; + param->name = NULL; + param->ptr = NULL; + param->img = img; + CHECK(ssol_image_ref_get(img), RES_OK); + CHECK(ssol_param_buffer_ref_get(NULL), RES_BAD_ARG); CHECK(ssol_param_buffer_ref_get(pbuf), RES_OK); CHECK(ssol_param_buffer_ref_put(NULL), RES_BAD_ARG); @@ -109,9 +162,10 @@ main(int argc, char** argv) CHECK(ssol_param_buffer_ref_put(pbuf), RES_OK); CHECK(ssol_param_buffer_create(dev, 8, &pbuf), RES_OK); - NCHECK(mem = ssol_param_buffer_allocate(pbuf, 2, 1), NULL); - CHECK(mem = ssol_param_buffer_allocate(pbuf, 1, 16), NULL); + NCHECK(mem = ssol_param_buffer_allocate(pbuf, 2, 1, NULL), NULL); + CHECK(mem = ssol_param_buffer_allocate(pbuf, 1, 16, NULL), NULL); + CHECK(ssol_image_ref_put(img), RES_OK); CHECK(ssol_param_buffer_ref_put(pbuf), RES_OK); CHECK(ssol_device_ref_put(dev), RES_OK); diff --git a/src/test_ssol_scene.c b/src/test_ssol_scene.c @@ -45,17 +45,6 @@ instance_func(struct ssol_instance* inst, void* context) return RES_OK; } -static void -get_wlen(const size_t i, double* wlen, double* data, void* ctx) -{ - double wavelengths[3] = { 10, 20, 30 }; - double intensities[3] = { 1, 2.1, 1.5 }; - CHECK(i < 3, 1); - (void)ctx; - *wlen = wavelengths[i]; - *data = intensities[i]; -} - int main(int argc, char** argv) { @@ -87,9 +76,9 @@ main(int argc, char** argv) struct ssol_sun* sun2; struct ssol_scene* scene; struct ssol_scene* scene2; - struct ssol_spectrum* spectrum; struct ssol_atmosphere* atm; struct ssol_atmosphere* atm2; + struct ssol_data absorption; struct ssol_vertex_data vdata; struct scene_ctx ctx; struct desc desc; @@ -179,12 +168,12 @@ main(int argc, char** argv) CHECK(ssol_scene_detach_sun(scene, sun), RES_BAD_ARG); CHECK(ssol_scene_detach_sun(scene2, sun), RES_OK); - CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK); - CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK); - CHECK(ssol_atmosphere_create_uniform(dev, &atm), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, spectrum), RES_OK); - CHECK(ssol_atmosphere_create_uniform(dev, &atm2), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(atm2, spectrum), RES_OK); + CHECK(ssol_atmosphere_create(dev, &atm), RES_OK); + absorption.type = SSOL_DATA_REAL; + absorption.value.real = 0.1; + CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_OK); + CHECK(ssol_atmosphere_create(dev, &atm2), RES_OK); + CHECK(ssol_atmosphere_set_absorption(atm2, &absorption), RES_OK); CHECK(ssol_scene_attach_atmosphere(NULL, atm), RES_BAD_ARG); CHECK(ssol_scene_attach_atmosphere(scene, NULL), RES_BAD_ARG); @@ -266,7 +255,6 @@ main(int argc, char** argv) CHECK(ssol_object_ref_put(object), RES_OK); CHECK(ssol_shape_ref_put(shape), RES_OK); CHECK(ssol_sun_ref_put(sun), RES_OK); - CHECK(ssol_spectrum_ref_put(spectrum), RES_OK); CHECK(ssol_atmosphere_ref_put(atm), RES_OK); CHECK(ssol_atmosphere_ref_put(atm2), RES_OK); CHECK(ssol_material_ref_put(material), RES_OK); diff --git a/src/test_ssol_shape.c b/src/test_ssol_shape.c @@ -151,13 +151,7 @@ main(int argc, char** argv) CHECK(ssol_shape_create_punched_surface(dev, NULL), RES_BAD_ARG); CHECK(ssol_shape_create_punched_surface(NULL, &shape), RES_BAD_ARG); CHECK(ssol_shape_create_punched_surface(dev, &shape), RES_OK); - - CHECK(ssol_shape_ref_get(NULL), RES_BAD_ARG); - CHECK(ssol_shape_ref_get(shape), RES_OK); - - CHECK(ssol_shape_ref_put(NULL), RES_BAD_ARG); - CHECK(ssol_shape_ref_put(shape), RES_OK); - + carving.get = get_polygon_vertices; carving.operation = SSOL_AND; carving.nb_vertices = npolygon_verts; @@ -196,6 +190,10 @@ main(int argc, char** argv) quadric.data.parabol.focal = 1; CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK); + punched_surface.nb_carvings = 0; + CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG); + punched_surface.nb_carvings = 1; + quadric.data.parabol.focal = 0; CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG); quadric.data.parabol.focal = 1; @@ -205,6 +203,10 @@ main(int argc, char** argv) quadric.data.hyperbol.img_focal = 1; CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK); + punched_surface.nb_carvings = 0; + CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG); + punched_surface.nb_carvings = 1; + quadric.data.hyperbol.real_focal = 0; CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG); quadric.data.hyperbol.real_focal = 1; @@ -217,10 +219,26 @@ main(int argc, char** argv) quadric.data.parabolic_cylinder.focal = 1; CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK); + punched_surface.nb_carvings = 0; + CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG); + punched_surface.nb_carvings = 1; + quadric.data.parabolic_cylinder.focal = 0; CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG); quadric.data.parabolic_cylinder.focal = 1; + quadric.type = SSOL_QUADRIC_HEMISPHERE; + quadric.data.hemisphere.radius = 10; + CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK); + + punched_surface.nb_carvings = 0; + CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK); + punched_surface.nb_carvings = 1; + + quadric.data.hemisphere.radius = 0; + CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG); + quadric.data.hemisphere.radius = 10; + CHECK(ssol_shape_ref_put(shape), RES_OK); CHECK(ssol_device_ref_put(dev), RES_OK); diff --git a/src/test_ssol_solver1.c b/src/test_ssol_solver1.c @@ -66,9 +66,11 @@ main(int argc, char** argv) struct ssol_sun* sun; struct ssol_sun* sun_mono; struct ssol_spectrum* spectrum; - struct ssol_spectrum* abs; + struct ssol_spectrum* abs_spectrum; + struct ssol_data abs_data; struct ssol_atmosphere* atm; struct ssol_estimator* estimator; + struct ssol_mc_sampled sampled; struct ssol_mc_global mc_global; struct ssol_mc_receiver mc_rcv; struct ssol_mc_shape mc_shape; @@ -78,19 +80,17 @@ main(int argc, char** argv) double dir[3]; double wavelengths[3] = { 1, 2, 3 }; double intensities[3] = { 1, 0.8, 1 }; - double mismatch[3] = { 1.5, 3.5, 0 }; double ka[3] = { 0, 0, 0 }; double mono = 1.21; double transform1[12]; /* 3x4 column major matrix */ double transform2[12]; /* 3x4 column major matrix */ double dbl; - size_t i, count, fcount; + size_t i, count, fcount, scount; FILE* tmp = NULL; double m, std; double a_m, a_std; uint32_t r_id; unsigned ntris; - (void) argc, (void) argv; d33_splat(transform1, 0); @@ -120,7 +120,6 @@ main(int argc, char** argv) CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK); CHECK(ssol_sun_set_dni(sun, 1000), RES_OK); CHECK(ssol_scene_create(dev, &scene), RES_OK); - CHECK(ssol_scene_attach_sun(scene, sun), RES_OK); CHECK(ssol_solve(NULL, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); CHECK(ssol_solve(scene, NULL, 10, 0, NULL, &estimator), RES_BAD_ARG); @@ -164,6 +163,13 @@ main(int argc, char** argv) CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 0), RES_OK); CHECK(ssol_scene_attach_instance(scene, target), RES_OK); + /* No sun */ + CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); + + CHECK(ssol_scene_attach_sun(scene, sun), RES_OK); + CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_OK); + CHECK(ssol_estimator_ref_put(estimator), RES_OK); + CHECK(ssol_solve (scene, rng, 1, &SSOL_PATH_TRACKER_DEFAULT, NULL, &estimator), RES_OK); @@ -203,12 +209,12 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_sampled_area(estimator, NULL), RES_BAD_ARG); CHECK(ssol_estimator_get_sampled_area(NULL, &dbl), RES_BAD_ARG); CHECK(ssol_estimator_get_sampled_area(estimator, &dbl), RES_OK); - CHECK(eq_eps(dbl, 12, 1.e-6), 1); + CHECK(eq_eps(dbl, 12, DBL_EPSILON), 1); - CHECK(ssol_estimator_get_count(NULL, NULL), RES_BAD_ARG); - CHECK(ssol_estimator_get_count(estimator, NULL), RES_BAD_ARG); - CHECK(ssol_estimator_get_count(NULL, &count), RES_BAD_ARG); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(NULL, NULL), RES_BAD_ARG); + CHECK(ssol_estimator_get_realisation_count(estimator, NULL), RES_BAD_ARG); + CHECK(ssol_estimator_get_realisation_count(NULL, &count), RES_BAD_ARG); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, 1); CHECK(ssol_estimator_get_failed_count(NULL, NULL), RES_BAD_ARG); @@ -217,6 +223,17 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_failed_count(estimator, &fcount), RES_OK); CHECK(fcount, 0); + CHECK(ssol_estimator_get_sampled_count(NULL, NULL), RES_BAD_ARG); + CHECK(ssol_estimator_get_sampled_count(estimator, NULL), RES_BAD_ARG); + CHECK(ssol_estimator_get_sampled_count(NULL, &scount), RES_BAD_ARG); + CHECK(ssol_estimator_get_sampled_count(estimator, &scount), RES_OK); + CHECK(scount, 3); + + CHECK(ssol_estimator_get_mc_sampled(NULL, heliostat, &sampled), RES_BAD_ARG); + CHECK(ssol_estimator_get_mc_sampled(estimator, NULL, &sampled), RES_BAD_ARG); + CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat, NULL), RES_BAD_ARG); + CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat, &sampled), RES_OK); + CHECK(ssol_estimator_get_mc_global(NULL, &mc_global), RES_BAD_ARG); CHECK(ssol_estimator_get_mc_global(estimator, NULL), RES_BAD_ARG); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); @@ -232,6 +249,7 @@ main(int argc, char** argv) CHECK(ssol_instance_sample(secondary, 0), RES_OK); CHECK(ssol_instance_sample(heliostat, 0), RES_OK); CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); + CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat, &sampled), RES_BAD_ARG); CHECK(ssol_instance_sample(target, 1), RES_OK); CHECK(ssol_instance_sample(secondary, 1), RES_OK); @@ -269,20 +287,6 @@ main(int argc, char** argv) CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 0), RES_OK); CHECK(ssol_estimator_ref_put(estimator), RES_OK); - /* Spectra mismatch */ - desc.wavelengths = mismatch; - desc.wavelengths = ka; - desc.count = 2; - CHECK(ssol_spectrum_create(dev, &abs), RES_OK); - CHECK(ssol_spectrum_setup(abs, get_wlen, 2, &desc), RES_OK); - CHECK(ssol_atmosphere_create_uniform(dev, &atm), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, abs), RES_OK); - CHECK(ssol_scene_attach_atmosphere(scene, atm), RES_OK); - CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); - CHECK(ssol_scene_detach_atmosphere(scene, atm), RES_OK); - CHECK(ssol_spectrum_ref_put(abs), RES_OK); - CHECK(ssol_atmosphere_ref_put(atm), RES_OK); - /* Can sample any geometry; variance is high */ NCHECK(tmp = tmpfile(), 0); #define N__ 10000 @@ -290,7 +294,7 @@ main(int argc, char** argv) #define GET_MC_GLOBAL ssol_estimator_get_mc_global CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target, &r_id), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); @@ -308,8 +312,9 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, m, 2 * dbl), 1); - CHECK(eq_eps(mc_global.missing.E, m, 2*mc_global.missing.SE), 1); + CHECK(eq_eps(mc_global.missing.E, 2*m, 2*mc_global.missing.SE), 1); CHECK(GET_MC_RCV(NULL, NULL, SSOL_BACK, NULL), RES_BAD_ARG); CHECK(GET_MC_RCV(estimator, NULL, SSOL_BACK, NULL), RES_BAD_ARG); CHECK(GET_MC_RCV(NULL, target, SSOL_BACK, NULL), RES_BAD_ARG); @@ -338,7 +343,7 @@ main(int argc, char** argv) NCHECK(tmp = tmpfile(), 0); CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); @@ -348,8 +353,10 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, m, 1e-4), 1); + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); @@ -358,18 +365,15 @@ main(int argc, char** argv) CHECK(ssol_estimator_ref_put(estimator), RES_OK); /* Check atmosphere model; with no absorption result is unchanged */ - desc.wavelengths = wavelengths; - desc.intensities = ka; - desc.count = 3; - CHECK(ssol_spectrum_create(dev, &abs), RES_OK); - CHECK(ssol_spectrum_setup(abs, get_wlen, 3, &desc), RES_OK); - CHECK(ssol_atmosphere_create_uniform(dev, &atm), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, abs), RES_OK); + CHECK(ssol_atmosphere_create(dev, &atm), RES_OK); + abs_data.type = SSOL_DATA_REAL; + abs_data.value.real = 0; + CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK); CHECK(ssol_scene_attach_atmosphere(scene, atm), RES_OK); NCHECK(tmp = tmpfile(), 0); CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); @@ -377,18 +381,20 @@ main(int argc, char** argv) CHECK(eq_eps(m, 4 * DNI_cos, MMAX(4 * DNI_cos * 1e-2, std)), 1); CHECK(eq_eps(std, 0, 1e-4), 1); CHECK(ssol_scene_detach_atmosphere(scene, atm), RES_OK); - CHECK(ssol_spectrum_ref_put(abs), RES_OK); CHECK(ssol_atmosphere_ref_put(atm), RES_OK); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, m, 1e-4), 1); + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1); CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1); + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(ssol_estimator_ref_put(estimator), RES_OK); /* Check atmosphere model and imperfect mirror: there are losses */ @@ -407,17 +413,17 @@ main(int argc, char** argv) CHECK(ssol_scene_attach_instance(scene, heliostat2), RES_OK); #define KA 0.03 - ka[0] = ka[1] = ka[2] = KA; - CHECK(ssol_spectrum_create(dev, &abs), RES_OK); - CHECK(ssol_spectrum_setup(abs, get_wlen, 3, &desc), RES_OK); - CHECK(ssol_atmosphere_create_uniform(dev, &atm), RES_OK); - CHECK(ssol_atmosphere_set_uniform_absorption(atm, abs), RES_OK); + abs_data.value.real = KA; + CHECK(ssol_spectrum_create(dev, &abs_spectrum), RES_OK); + CHECK(ssol_spectrum_setup(abs_spectrum, get_wlen, 3, &desc), RES_OK); + CHECK(ssol_atmosphere_create(dev, &atm), RES_OK); + CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK); CHECK(ssol_scene_attach_atmosphere(scene, atm), RES_OK); CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 1), RES_OK); NCHECK(tmp = tmpfile(), 0); CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &a_m, &a_std), RES_OK); CHECK(fclose(tmp), 0); @@ -428,9 +434,13 @@ main(int argc, char** argv) CHECK(eq_eps(a_std, 0, 1e-4), 1); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE); - printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); + printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE); + printf("Atmosphere = %g +/- %g; ", mc_global.atmosphere.E, mc_global.atmosphere.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E + mc_global.atmosphere.E + mc_global.absorbed.E, + m, 1e-4), 1); + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf ("\tIr(target) = %g +/- %g (%.2g %%)\n", @@ -447,11 +457,11 @@ main(int argc, char** argv) mc_rcv.reflectivity_loss.E, mc_rcv.reflectivity_loss.SE, 100 * mc_rcv.reflectivity_loss.E / m); - printf - ("\tCos Loss(target) = %g +/- %g (%.2g %%)\n", - mc_rcv.cos_loss.E, - mc_rcv.cos_loss.SE, - 100 * mc_rcv.cos_loss.E / m); + + CHECK(ssol_estimator_get_sampled_count(estimator, &scount), RES_OK); + CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat, &sampled), RES_BAD_ARG); + CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat2, &sampled), RES_OK); + CHECK(eq_eps(mc_rcv.integrated_irradiance.E, a_m, 1e-8), 1); CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, a_std, 1e-4), 1); CHECK(eq_eps @@ -462,9 +472,7 @@ main(int argc, char** argv) ( mc_rcv.integrated_irradiance.E + mc_rcv.absorptivity_loss.E + mc_rcv.reflectivity_loss.E - + mc_rcv.cos_loss.E, 4 * DNI, 1e-8), 1); - CHECK(eq_eps(mc_rcv.cos_loss.E / (4 * DNI), 1 - COS, 1e-8), 1); - + + (1 - mc_global.cos_factor.E) * 4 * DNI, 4 * DNI, 1e-8), 1); CHECK(ssol_mc_receiver_get_mc_shape(NULL, NULL, NULL), RES_BAD_ARG); CHECK(ssol_mc_receiver_get_mc_shape(&mc_rcv, NULL, NULL), RES_BAD_ARG); CHECK(ssol_mc_receiver_get_mc_shape(NULL, square, NULL), RES_BAD_ARG); @@ -488,8 +496,16 @@ main(int argc, char** argv) dbl = 0; FOR_EACH(i, 0, ntris) { + double v0[3], v1[3], v2[3], E1[3], E2[3], N[3], area; + unsigned ids[3]; + CHECK(ssol_shape_get_triangle_indices(square, (unsigned)i, ids), RES_OK); + CHECK(ssol_shape_get_vertex_attrib(square, ids[0], SSOL_POSITION, v0), RES_OK); + CHECK(ssol_shape_get_vertex_attrib(square, ids[1], SSOL_POSITION, v1), RES_OK); + CHECK(ssol_shape_get_vertex_attrib(square, ids[2], SSOL_POSITION, v2), RES_OK); + area = d3_len(d3_cross(N, d3_sub(E1, v1, v0), d3_sub(E2, v2, v0))) * 0.5; + CHECK(ssol_mc_shape_get_mc_primitive(&mc_shape, (unsigned)i, &mc_prim), RES_OK); - dbl += mc_prim.integrated_irradiance.E; + dbl += mc_prim.integrated_irradiance.E * area; } CHECK(eq_eps(dbl, a_m, 1.e-6), 1); @@ -513,11 +529,16 @@ main(int argc, char** argv) ka[1] = 0.2; ka[0] = ka[2] = 0.1; desc.wavelengths = wavelengths; desc.intensities = ka; - desc.count = 2; - CHECK(ssol_spectrum_setup(abs, get_wlen, 2, &desc), RES_OK); + desc.count = 3; + CHECK(ssol_spectrum_setup(abs_spectrum, get_wlen, 3, &desc), RES_OK); + CHECK(ssol_spectrum_setup(abs_spectrum, get_wlen, 2, &desc), RES_OK); + abs_data.type = SSOL_DATA_SPECTRUM; + abs_data.value.spectrum = abs_spectrum; + CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK); + NCHECK(tmp = tmpfile(), 0); CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); @@ -528,8 +549,10 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, m, 1e-4), 1); + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); @@ -551,10 +574,10 @@ main(int argc, char** argv) CHECK(ssol_material_ref_put(v_mtl), RES_OK); CHECK(ssol_device_ref_put(dev), RES_OK); CHECK(ssol_scene_ref_put(scene), RES_OK); - CHECK(ssol_spectrum_ref_put(abs), RES_OK); CHECK(ssol_atmosphere_ref_put(atm), RES_OK); CHECK(ssol_estimator_ref_put(estimator), RES_OK); CHECK(ssol_spectrum_ref_put(spectrum), RES_OK); + CHECK(ssol_spectrum_ref_put(abs_spectrum), RES_OK); CHECK(ssol_sun_ref_put(sun), RES_OK); CHECK(ssol_sun_ref_put(sun_mono), RES_OK); CHECK(ssp_rng_ref_put(rng), RES_OK); diff --git a/src/test_ssol_solver2.c b/src/test_ssol_solver2.c @@ -182,20 +182,23 @@ main(int argc, char** argv) #define GET_MC_RCV ssol_estimator_get_mc_receiver CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target, &r_id), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); printf("Ir = %g +/- %g\n", m, std); -#define DNI_cos (1000 * cos(PI / 4)) +#define COS cos(PI / 4) +#define DNI_cos (1000 * COS) CHECK(eq_eps(m, 4 * DNI_cos, 4 * DNI_cos * 1e-4), 1); #define SQR(x) ((x)*(x)) CHECK(eq_eps(std, 0, 1e-4), 1); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); - printf("Missing = %g +/- %g", mc_global.missing.E, mc_global.missing.SE); + printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, 4 * DNI_cos, 1e-4), 1); /* nothing absorbed */ + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(GET_MC_RCV(estimator, heliostat1, SSOL_BACK, &mc_rcv), RES_BAD_ARG); CHECK(GET_MC_RCV(estimator, secondary, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(secondary) = %g +/- %g\n", diff --git a/src/test_ssol_solver2b.c b/src/test_ssol_solver2b.c @@ -185,12 +185,13 @@ main(int argc, char** argv) #define N__ 50000 CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target, &r_id), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); printf("Ir = %g +/- %g\n", m, std); -#define DNI_cos (1000 * cos(PI / 4)) +#define COS cos(PI / 4) +#define DNI_cos (1000 * COS) CHECK(eq_eps(m, 2 * DNI_cos, MMAX(2 * DNI_cos * 1e-2, std)), 1); #define SQR(x) ((x)*(x)) CHECK(eq_eps(std, @@ -198,8 +199,10 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, 4 * DNI_cos, 1e-4), 1); /* nothing absorbed */ + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(ssol_estimator_get_mc_receiver (estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", diff --git a/src/test_ssol_solver3.c b/src/test_ssol_solver3.c @@ -70,6 +70,7 @@ main(int argc, char** argv) struct ssol_mc_receiver mc_rcv; double dir[3]; double transform[12]; /* 3x4 column major matrix */ + double area; size_t count; FILE* tmp; double m, std; @@ -139,12 +140,13 @@ main(int argc, char** argv) #define N__ 20000 CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target, &r_id), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); printf("Ir = %g +/- %g\n", m, std); -#define DNI_cos (1000 * cos(PI / 4)) +#define COS cos(PI / 4) +#define DNI_cos (1000 * COS) CHECK(eq_eps(m, 4 * DNI_cos, 4 * DNI_cos * 2e-1), 1); #define SQR(x) ((x)*(x)) CHECK(eq_eps(std, @@ -152,8 +154,10 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, 400 * DNI_cos, 1e-4), 1); /* nothing absorbed */ + CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); CHECK(ssol_estimator_get_mc_receiver (estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", @@ -162,6 +166,8 @@ main(int argc, char** argv) CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1); CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); CHECK(count, 0); + CHECK(ssol_instance_get_area(heliostat, &area), RES_OK); + CHECK(eq_eps(area, 400, DBL_EPSILON), 1); /* Free data */ CHECK(ssol_instance_ref_put(heliostat), RES_OK); diff --git a/src/test_ssol_solver4.c b/src/test_ssol_solver4.c @@ -112,6 +112,7 @@ main(int argc, char** argv) carving.context = &POLY_EDGES__; quadric.type = SSOL_QUADRIC_PARABOL; quadric.data.parabol.focal = FOCAL; + quadric.slices_count_hint = 100; punched.nb_carvings = 1; punched.quadric = &quadric; punched.carvings = &carving; @@ -144,18 +145,19 @@ main(int argc, char** argv) CHECK(ssol_scene_attach_instance(scene, target2), RES_OK); NCHECK(tmp = tmpfile(), 0); -#define N__ 10000 +#define N__ 100000 #define GET_MC_RCV ssol_estimator_get_mc_receiver CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target1, &r_id1), RES_OK); CHECK(ssol_instance_get_id(target2, &r_id2), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id1, count, &m1, &std1), RES_OK); CHECK(pp_sum(tmp, (int32_t)r_id2, count, &m2, &std2), RES_OK); CHECK(fclose(tmp), 0); printf("Ir = %g +/- %g\n", m1, std1); -#define DNI_cos (1000 * cos(0)) +#define COS cos(0) +#define DNI_cos (1000 * COS) CHECK(eq_eps(m1, 400 * DNI_cos, 400 * DNI_cos * 1e-4), 1); CHECK(eq_eps(std1, 0, 1), 1); CHECK(m1, m2); @@ -163,18 +165,19 @@ main(int argc, char** argv) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, 400 * DNI_cos, 1e-2), 1); /* nothing absorbed */ CHECK(GET_MC_RCV(estimator, target1, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target1) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); - CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m1, 1e-8), 1); - CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std1, 1e-4), 1); + CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m1, 1e-2), 1); + CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std1, 1e-2), 1); CHECK(GET_MC_RCV(estimator, target2, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target2) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); - CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m2, 1e-8), 1); - CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std2, 1e-4), 1); + CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m2, 1e-2), 1); + CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std2, 1e-2), 1); CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); CHECK(count, 0); diff --git a/src/test_ssol_solver5.c b/src/test_ssol_solver5.c @@ -140,19 +140,21 @@ main(int argc, char** argv) #define N__ 10000 CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target, &r_id), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK); CHECK(fclose(tmp), 0); printf("Ir = %g +/- %g\n", m, std); -#define DNI_cos (1000 * cos(0)) +#define COS cos(0) +#define DNI_cos (1000 * COS) CHECK(eq_eps(m, 400 * DNI_cos, 20), 1); CHECK(eq_eps(std, 0, 1), 1); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, 400 * DNI_cos, 1e-4), 1); /* nothing absorbed */ CHECK(ssol_estimator_get_mc_receiver (estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", diff --git a/src/test_ssol_solver6.c b/src/test_ssol_solver6.c @@ -181,8 +181,10 @@ main(int argc, char** argv) CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(fclose(tmp), 0); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); + printf("Absorbed = %g +/- %g\n", mc_global.absorbed.E, mc_global.absorbed.SE); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); CHECK(eq_eps(mc_global.shadowed.E, 100000, 2 * 100000/sqrt(N__)), 1); CHECK(eq_eps(mc_global.missing.E, 0, 0), 1); @@ -190,11 +192,13 @@ main(int argc, char** argv) printf("Ir(target1) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); CHECK(eq_eps(mc_rcv.integrated_irradiance.E, 100000, 2*100000/sqrt(N__)), 1); + CHECK(mc_rcv.integrated_irradiance.E, mc_rcv.integrated_absorbed_irradiance.E); CHECK(GET_MC_RCV(estimator, target2, SSOL_BACK, &mc_rcv), RES_OK); printf("Ir(target2) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); CHECK(eq_eps(mc_rcv.integrated_irradiance.E, 0, 1), 1); + CHECK(mc_rcv.integrated_irradiance.E, mc_rcv.integrated_absorbed_irradiance.E); /* Free data */ CHECK(ssol_instance_ref_put(heliostat1), RES_OK); diff --git a/src/test_ssol_solver7.c b/src/test_ssol_solver7.c @@ -19,12 +19,7 @@ #include <rsys/mem_allocator.h> #include <rsys/image.h> -#define SCREEN_GAMMA 2.2 -#define WIDTH 640 -#define HEIGHT 480 -#define PROJ_RATIO (double)WIDTH/(double)HEIGHT - -#define REFLECTIVITY 0 +#define REFLECTIVITY 0.1 #include "test_ssol_materials.h" #define TARGET_SZ 2 @@ -85,7 +80,7 @@ main(int argc, char** argv) struct ssol_estimator* estimator; struct ssol_mc_global mc_global; struct ssol_mc_receiver mc_rcv; - double dir[3]; + double dir[3], area; double transform[12]; /* 3x4 column major matrix */ FILE* tmp; /* primary is a parabol */ @@ -201,20 +196,26 @@ main(int argc, char** argv) CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(fclose(tmp), 0); - printf("Total = %g\n", TOTAL); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); - printf("Shadows = %g +/- %g\n", - mc_global.shadowed.E, mc_global.shadowed.SE); + CHECK(ssol_estimator_get_sampled_area(estimator, &area), RES_OK); + printf("Total = %g\n", area * DNI_cos); + CHECK(eq_eps(area * DNI_cos, TOTAL, TOTAL * 1e-4), 1); + printf("Absorbed = %g +/- %g\n", mc_global.absorbed.E, mc_global.absorbed.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); + printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - printf("Missing = %g +/- %g\n", - mc_global.missing.E, mc_global.missing.SE); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); + printf("Abs(target1) = %g +/- %g\n", + mc_rcv.integrated_absorbed_irradiance.E, + mc_rcv.integrated_absorbed_irradiance.SE); printf("Ir(target1) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); CHECK(eq_eps(mc_rcv.integrated_irradiance.E, TOTAL, TOTAL * 1e-4), 1); - CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, 0, 1e-4), 1); + CHECK(eq_eps(mc_rcv.integrated_absorbed_irradiance.E, + (1 - REFLECTIVITY) * TOTAL, (1 - REFLECTIVITY) *TOTAL * 1e-4), 1); + CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, 0, 1e-2), 1); /* Free data */ CHECK(ssol_instance_ref_put(heliostat), RES_OK); diff --git a/src/test_ssol_solver8.c b/src/test_ssol_solver8.c @@ -144,19 +144,19 @@ main(int argc, char** argv) #define GET_MC_RCV ssol_estimator_get_mc_receiver CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target, &r_id), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(fclose(tmp), 0); CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); CHECK(count, 0); #define S (sqrt(2) * X_SZ * Y_SZ) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); - printf("Cos = %g +/- %g\n", mc_global.cos_loss.E, mc_global.cos_loss.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); - CHECK(eq_eps(mc_global.cos_loss.E, 0, 1e-4), 1); CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1); - CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1); + CHECK(eq_eps(mc_global.missing.E, S * DNI, + 3 * mc_global.missing.SE), 1); /* nothing absorbed */ CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target1) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); diff --git a/src/test_ssol_solver9.c b/src/test_ssol_solver9.c @@ -138,7 +138,7 @@ main(int argc, char** argv) #define GET_MC_RCV ssol_estimator_get_mc_receiver CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK); CHECK(ssol_instance_get_id(target, &r_id), RES_OK); - CHECK(ssol_estimator_get_count(estimator, &count), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(fclose(tmp), 0); CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); @@ -146,13 +146,11 @@ main(int argc, char** argv) #define DNI_TGT_S (DNI * TGT_X * TGT_Y) #define DNI_S (DNI * SZ * SZ) CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); - printf("Cos = %g +/- %g\n", mc_global.cos_loss.E, mc_global.cos_loss.SE); + printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); - /* CHECK(eq_eps(mc_global.cos_loss.E, 4 * DNI_S, DNI_S * 1e-4), 1); */ CHECK(eq_eps(mc_global.shadowed.E, DNI_S, 3 * mc_global.shadowed.SE), 1); - CHECK(eq_eps(mc_global.missing.E, - MMAX(DNI_S, DNI_TGT_S) - MMIN(DNI_S, DNI_TGT_S), + CHECK(eq_eps(mc_global.missing.E, MMAX(DNI_S, DNI_TGT_S), 3 * mc_global.missing.SE), 1); CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target1) = %g +/- %g\n", diff --git a/src/test_ssol_spectrum.c b/src/test_ssol_spectrum.c @@ -77,6 +77,14 @@ main(int argc, char** argv) CHECK(wlens_count, 0); CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, &desc), RES_OK); CHECK(wlens_count, 3); + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, &desc), RES_OK); + + desc.wavelengths[1] = 30; + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, &desc), RES_BAD_ARG); + + desc.wavelengths[1] = 20; + desc.data[1] = -2.1; + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, &desc), RES_BAD_ARG); CHECK(ssol_spectrum_ref_put(spectrum), RES_OK);