commit 4f403b04ecbac79a0fae222e0051812ed0b6ca3b
parent afb9054ca9e495bd48cf8481f3e868ec190f29b7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 27 Jul 2018 16:51:30 +0200
Compute the SW CDFs wrt the CIE XYZ tristimulus
Diffstat:
3 files changed, 229 insertions(+), 3 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -62,6 +62,10 @@ rcmake_prepend_path(HTGOP_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_library(htgop SHARED ${HTGOP_FILES_SRC} ${HTGOP_FILES_INC} ${HTGOP_FILES_INC_API})
target_link_libraries(htgop RSys)
+if(CMAKE_COMPILER_IS_GNUCC)
+ target_link_libraries(htgop m)
+endif()
+
set_target_properties(htgop PROPERTIES
DEFINE_SYMBOL HTGOP_SHARED_BUILD
VERSION ${VERSION}
diff --git a/src/htgop.c b/src/htgop.c
@@ -21,6 +21,16 @@
#include <rsys/logger.h>
#include <rsys/mem_allocator.h>
+#include <math.h>
+
+#define WNUM_TO_WLEN(Num) (1.0e7 / ((double)Num))
+#define WLEN_TO_WNUM(Len) WNUM_TO_WLEN(Len)
+
+#define CIE_XYZ_WAVELENGTH_MIN 380.0 /* In nanometer */
+#define CIE_XYZ_WAVELENGTH_MAX 780.0 /* In nanometer */
+#define CIE_XYZ_WAVENUMBER_MIN WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX) /*In cm^-1*/
+#define CIE_XYZ_WAVENUMBER_MAX WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN) /*In cm^-1*/
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -39,6 +49,178 @@ log_msg
}
}
+static FINLINE double
+trapezoidal_integration
+ (const double lambda_lo, /* Integral lower bound. In nanometer */
+ const double lambda_hi, /* Integral upper bound. In nanometer */
+ double (*f_bar)(const double lambda)) /* Function to integrate */
+{
+ double dlambda;
+ size_t i, n;
+ double integral;
+ ASSERT(lambda_lo <= lambda_hi);
+ ASSERT(lambda_lo > 0);
+
+ n = (size_t)(lambda_hi - lambda_lo) + 1;
+ dlambda = (lambda_hi - lambda_lo) / (double)n;
+
+ FOR_EACH(i, 0, n) {
+ const double lambda1 = lambda_lo + dlambda*(double)(i+0);
+ const double lambda2 = lambda_hi + dlambda*(double)(i+1);
+ const double f1 = f_bar(lambda1);
+ const double f2 = f_bar(lambda2);
+ integral += (f1 + f2)*dlambda*0.5;
+ }
+ return integral;
+}
+
+/* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved
+ * has defined by the 1931 standard observers. These analytical fits are
+ * propsed by C. Wyman, P. P. Sloan & P. Shirley in "Simple Analytic
+ * Approximations to the CIE XYZ Color Matching Functions" - JCGT 2013. */
+static FINLINE double
+fit_x_bar_1931(const double lambda)
+{
+ const double a = (lambda - 442.0) * (lambda < 442.0 ? 0.0624 : 0.0374);
+ const double b = (lambda - 599.8) * (lambda < 599.8 ? 0.0264 : 0.0323);
+ const double c = (lambda - 501.1) * (lambda < 501.1 ? 0.0490 : 0.0382);
+ return 0.362*exp(-0.5*a*a) + 1.056*exp(-0.5f*b*b) - 0.065*exp(-0.5*c*c);
+}
+
+static FINLINE double
+fit_y_bar_1931(const double lambda)
+{
+ const double a = (lambda - 568.8) * (lambda < 568.8 ? 0.0213 : 0.0247);
+ const double b = (lambda - 530.9) * (lambda < 530.9 ? 0.0613 : 0.0322);
+ return 0.821*exp(-0.5*a*a) + 0.286*exp(-0.5*b*b);
+}
+
+static FINLINE double
+fit_z_bar_1931(const double lambda)
+{
+ const double a = (lambda - 437.0) * (lambda < 437.0 ? 0.0845 : 0.0278);
+ const double b = (lambda - 459.0) * (lambda < 459.0 ? 0.0385 : 0.0725);
+ return 1.217*exp(-0.5*a*a) + 0.681*exp(-0.5*b*b);
+}
+
+static res_T
+compute_CIE_XYZ_cdf(struct htgop* htgop)
+{
+ enum { X, Y, Z }; /* Helper constant */
+ size_t iband_min, iband_max;
+ const double* wnums = NULL;
+ double sum[3] = {0,0,0};
+ double* pdf[3] = {NULL, NULL, NULL};
+ double* cdf[3] = {NULL, NULL, NULL};
+ size_t i, n;
+ res_T res = RES_OK;
+ ASSERT(htgop);
+
+ HTGOP(get_sw_spectral_intervals_wave_numbers(htgop, &wnums));
+ HTGOP(get_sw_spectral_intervals_count(htgop, &n));
+ ASSERT(n);
+
+ /* Allocate and reset the memory space for the tristimulus CDF */
+ #define SETUP(Stimulus) { \
+ res = darray_double_resize(&htgop->sw_ ## Stimulus ## _cdf, n); \
+ if(res != RES_OK) { \
+ log_err(htgop, \
+ "Could not reserve the memory space for the CDF " \
+ "of the "STR(X)" stimulus.\n"); \
+ goto error; \
+ } \
+ memset(darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf), 0, \
+ n*sizeof(double)); \
+ cdf[Stimulus] = darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf); \
+ pdf[Stimulus] = cdf[Stimulus]; \
+ } (void)0
+ SETUP(X);
+ SETUP(Y);
+ SETUP(Z);
+ #undef RESERVE
+
+ ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MAX, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MIN, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MAX, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MIN), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MIN, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MAX), 1.e-6));
+
+ /* The sw spectral intervals does not intersect the CIE XYZ interval */
+ if(wnums[0] > CIE_XYZ_WAVENUMBER_MAX || wnums[n] < CIE_XYZ_WAVENUMBER_MIN) {
+ log_warn(htgop,
+ "The loaded SW spectral intervals does not intersect "
+ "the CIE XYZ interval.\n");
+ goto exit;
+ }
+
+ /* Search the lower and upper bound of the bands that overlap the XYZ
+ * spectral interval. Use a simple linear search since the overall number of
+ * bands should be small */
+ iband_min = 0, iband_max = n-1;
+ FOR_EACH(i, 0, n) {
+ if(wnums[i]<=CIE_XYZ_WAVENUMBER_MIN && CIE_XYZ_WAVENUMBER_MIN<=wnums[i+1])
+ iband_min = i;
+ if(wnums[i]<=CIE_XYZ_WAVENUMBER_MAX && CIE_XYZ_WAVENUMBER_MAX<=wnums[i+1])
+ iband_max = i;
+ }
+ ASSERT(iband_min <= iband_max);
+
+ /* Compute the *unormalized* pdf of the tristimulus */
+ FOR_EACH(i, iband_min, iband_max+1) {
+ const double lambda_lo = WNUM_TO_WLEN(wnums[i+1]);
+ const double lambda_hi = WNUM_TO_WLEN(wnums[i+0]);
+ ASSERT(lambda_lo <= lambda_hi);
+ pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
+ pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
+ pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
+ sum[X] += pdf[X][i];
+ sum[Y] += pdf[Y][i];
+ sum[Z] += pdf[Z][i];
+ }
+
+ /* Compute the cdf of the tristimulus. Note that the upper bound is 'n' and not
+ * 'iband_max+1' in order to ensure that the CDF is strictly increasing. */
+ FOR_EACH(i, iband_min, n) {
+ /* Normalize the pdf */
+ pdf[X][i] /= sum[X];
+ pdf[Y][i] /= sum[Y];
+ pdf[Z][i] /= sum[Z];
+ /* Setup the cumulative */
+ if(i == 0) {
+ cdf[X][i] = pdf[X][i];
+ cdf[Y][i] = pdf[Y][i];
+ cdf[Z][i] = pdf[Z][i];
+ } else {
+ cdf[X][i] = pdf[X][i] + cdf[X][i-1];
+ cdf[Y][i] = pdf[Y][i] + cdf[Y][i-1];
+ cdf[Z][i] = pdf[Z][i] + cdf[Z][i-1];
+ }
+ }
+ ASSERT(eq_eps(cdf[X][n-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Y][n-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Z][n-1], 1, 1.e-6));
+
+#ifndef NDEBUG
+ FOR_EACH(i, 1, n) {
+ ASSERT(cdf[X][i] >= cdf[X][i-1]);
+ ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
+ ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
+ }
+#endif
+
+ /* Handle numerical issue */
+ cdf[X][n-1] = 1.0;
+ cdf[Y][n-1] = 1.0;
+ cdf[Z][n-1] = 1.0;
+
+exit:
+ return res;
+error:
+ darray_double_clear(&htgop->sw_X_cdf);
+ darray_double_clear(&htgop->sw_Y_cdf);
+ darray_double_clear(&htgop->sw_Z_cdf);
+ goto exit;
+}
+
static res_T
parse_spectral_intervals
(struct htgop* htgop,
@@ -171,7 +353,12 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
ASSERT(htgop && stream && stream_name);
reader_init(&rdr, stream, stream_name);
- #define CALL(Func) { if((res = Func) != RES_OK) goto error; }
+ #define CALL(Func) { \
+ if((res = Func) != RES_OK) { \
+ log_err(htgop, "%s:%lu: Parsing error.\n", rdr.name, rdr.iline); \
+ goto error; \
+ } \
+ } (void)0
/* Parse the number of levels/layers */
CALL(cstr_to_ulong(read_line(&rdr), &nlvls));
@@ -223,7 +410,7 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
CALL(cstr_to_double(read_line(&rdr), &htgop->ground.lw_emissivity));
CALL(cstr_to_double(read_line(&rdr), &htgop->ground.sw_emissivity));
- /* Parse the numver of long/short wave spectral intervals */
+ /* Parse the number of long/short wave spectral intervals */
CALL(cstr_to_ulong(read_line(&rdr), &nspecints_lw));
CALL(cstr_to_ulong(read_line(&rdr), &nspecints_sw));
@@ -236,11 +423,14 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
#undef CALL
+ /* Compute the cdf for the SW CIE XYZ spectral interval */
+ res = compute_CIE_XYZ_cdf(htgop);
+ if(res != RES_OK) goto error;
+
exit:
reader_release(&rdr);
return res;
error:
- log_err(htgop, "%s:%lu: Parsing error.\n", rdr.name, rdr.iline);
goto exit;
}
@@ -252,6 +442,9 @@ release_htgop(ref_T* ref)
htgop = CONTAINER_OF(ref, struct htgop, ref);
spectral_intervals_release(&htgop->lw_specints);
spectral_intervals_release(&htgop->sw_specints);
+ darray_double_release(&htgop->sw_X_cdf);
+ darray_double_release(&htgop->sw_Y_cdf);
+ darray_double_release(&htgop->sw_Z_cdf);
darray_level_release(&htgop->levels);
darray_layer_release(&htgop->layers);
MEM_RM(htgop->allocator, htgop);
@@ -296,6 +489,9 @@ htgop_create
htgop->verbose = verbose;
spectral_intervals_init(allocator, &htgop->lw_specints);
spectral_intervals_init(allocator, &htgop->sw_specints);
+ darray_double_init(allocator, &htgop->sw_X_cdf);
+ darray_double_init(allocator, &htgop->sw_Y_cdf);
+ darray_double_init(allocator, &htgop->sw_Z_cdf);
darray_level_init(allocator, &htgop->levels);
darray_layer_init(allocator, &htgop->layers);
@@ -598,3 +794,14 @@ log_err(const struct htgop* htgop, const char* msg, ...)
va_end(vargs_list);
}
+void
+log_warn(const struct htgop* htgop, const char* msg, ...)
+{
+ va_list vargs_list;
+ ASSERT(htgop && msg);
+
+ va_start(vargs_list, msg);
+ log_msg(htgop, LOG_WARNING, msg, vargs_list);
+ va_end(vargs_list);
+}
+
diff --git a/src/htgop_c.h b/src/htgop_c.h
@@ -34,6 +34,11 @@ struct htgop {
struct spectral_intervals lw_specints;
struct spectral_intervals sw_specints;
+ /* Cumulative distribution function for the Short Wave CIE XYZ tristimulus */
+ struct darray_double sw_X_cdf;
+ struct darray_double sw_Y_cdf;
+ struct darray_double sw_Z_cdf;
+
struct darray_layer layers; /* Par layer data */
struct darray_level levels; /* Per level data (#level = #layer + 1 ) */
@@ -53,5 +58,15 @@ log_err
#endif
;
+extern LOCAL_SYM void
+log_warn
+ (const struct htgop* htgop,
+ const char* fmt,
+ ...)
+#ifdef COMPILER_GCC
+ __attribute((format(printf, 2, 3)))
+#endif
+ ;
+
#endif /* HTGOP_C_H */
