diff --git a/include/mitsuba/core/spectrum.h b/include/mitsuba/core/spectrum.h index 72723cb12..914939616 100644 --- a/include/mitsuba/core/spectrum.h +++ b/include/mitsuba/core/spectrum.h @@ -3,11 +3,8 @@ #include #include -#include -#include #include #include -#include #include #include #include @@ -149,13 +146,15 @@ extern MTS_EXPORT_CORE void cie_alloc(); */ template > Result cie1931_xyz(Float wavelength, mask_t active = true) { - using Int32 = int32_array_t; + using Int32 = int32_array_t; + using ScalarFloat = scalar_t; - Float t = (wavelength - MTS_CIE_MIN) * - ((MTS_CIE_SAMPLES - 1) / (MTS_CIE_MAX - MTS_CIE_MIN)); + Float t = (wavelength - (ScalarFloat) MTS_CIE_MIN) * + ((MTS_CIE_SAMPLES - 1) / + ((ScalarFloat) MTS_CIE_MAX - (ScalarFloat) MTS_CIE_MIN)); - active &= wavelength >= MTS_CIE_MIN && - wavelength <= MTS_CIE_MAX; + active &= wavelength >= (ScalarFloat) MTS_CIE_MIN && + wavelength <= (ScalarFloat) MTS_CIE_MAX; Int32 i0 = clamp(Int32(t), zero(), Int32(MTS_CIE_SAMPLES - 2)), i1 = i0 + 1; @@ -168,7 +167,7 @@ Result cie1931_xyz(Float wavelength, mask_t active = true) { v1_z = gather(cie1931_z_data, i1, active); Float w1 = t - Float(i0), - w0 = 1.f - w1; + w0 = (ScalarFloat) 1.f - w1; return Result(fmadd(w0, v0_x, w1 * v1_x), fmadd(w0, v0_y, w1 * v1_y), @@ -183,12 +182,14 @@ Result cie1931_xyz(Float wavelength, mask_t active = true) { template Float cie1931_y(Float wavelength, mask_t active = true) { using Int32 = int32_array_t; + using ScalarFloat = scalar_t; - Float t = (wavelength - MTS_CIE_MIN) * - ((MTS_CIE_SAMPLES - 1) / (MTS_CIE_MAX - MTS_CIE_MIN)); + Float t = (wavelength - (ScalarFloat) MTS_CIE_MIN) * + ((MTS_CIE_SAMPLES - 1) / + ((ScalarFloat) MTS_CIE_MAX - (ScalarFloat) MTS_CIE_MIN)); - active &= wavelength >= MTS_CIE_MIN && - wavelength <= MTS_CIE_MAX; + active &= wavelength >= (ScalarFloat) MTS_CIE_MIN && + wavelength <= (ScalarFloat) MTS_CIE_MAX; Int32 i0 = clamp(Int32(t), zero(), Int32(MTS_CIE_SAMPLES - 2)), i1 = i0 + 1; @@ -197,7 +198,7 @@ Float cie1931_y(Float wavelength, mask_t active = true) { v1 = gather(cie1931_y_data, i1, active); Float w1 = t - Float(i0), - w0 = 1.f - w1; + w0 = (ScalarFloat) 1.f - w1; return select(active, fmadd(w0, v0, w1 * v1), 0.f); } @@ -311,85 +312,13 @@ std::pair, Spectrum> sample_wavelength(Float sample) { } template -inline void spectrum_from_file(const std::string &filename, - std::vector &wavelengths, - std::vector &values) { - auto fs = Thread::thread()->file_resolver(); - fs::path file_path = fs->resolve(filename); - if (!fs::exists(file_path)) - Log(Error, "\"%s\": file does not exist!", file_path); - - Log(Info, "Loading spectral data file \"%s\" ..", file_path); - ref file = new FileStream(file_path); - - std::string line, rest; - float wav, value; - while (true) { - try { - line = file->read_line(); - if (line.size() == 0 || line[0] == '#') - continue; - - std::istringstream iss(line); - iss >> wav; - iss >> value; - if (iss >> rest) - Log(Error, "\"%s\": excess tokens after wavlengths-value pair in file:\n%s!", file_path, line); - wavelengths.push_back(wav); - values.push_back(value); - } catch (std::exception &) { - break; - } - } -} +MTS_EXPORT_CORE void spectrum_from_file(const std::string &filename, + std::vector &wavelengths, + std::vector &values); template -inline Color spectrum_to_rgb(const std::vector &wavelengths, - const std::vector &values, - bool bounded=true) { - Color color = 0.f; - - const int steps = 1000; - for (int i = 0; i < steps; ++i) { - Scalar x = MTS_WAVELENGTH_MIN + - (i / (Scalar) (steps - 1)) * - (MTS_WAVELENGTH_MAX - MTS_WAVELENGTH_MIN); - - if (x < wavelengths.front() || x > wavelengths.back()) - continue; - - // Find interval containing 'x' - uint32_t index = math::find_interval( - (uint32_t) wavelengths.size(), - [&](uint32_t idx) { - return wavelengths[idx] <= x; - }); - - Scalar x0 = wavelengths[index], - x1 = wavelengths[index + 1], - y0 = values[index], - y1 = values[index + 1]; - - // Linear interpolant at 'x' - Scalar y = (x*y0 - x1*y0 - x*y1 + x0*y1) / (x0 - x1); - - Color xyz = cie1931_xyz(x); - color += xyz * y; - } - - // Last specified value repeats implicitly - color *= (MTS_WAVELENGTH_MAX - MTS_WAVELENGTH_MIN) / (Scalar) steps; - color = xyz_to_srgb(color); - - if (bounded && any(color < 0.f || color > 1.f)) { - Log(Warn, "Spectrum: clamping out-of-gamut color %s", color); - color = clamp(color, 0.f, 1.f); - } else if (!bounded && any(color < 0.f)) { - Log(Warn, "Spectrum: clamping out-of-gamut color %s", color); - color = max(color, 0.f); - } - - return color; -} +MTS_EXPORT_CORE Color spectrum_to_rgb(const std::vector &wavelengths, + const std::vector &values, + bool bounded = true); NAMESPACE_END(mitsuba) diff --git a/src/libcore/spectrum.cpp b/src/libcore/spectrum.cpp index 20e89ce65..e1912041b 100644 --- a/src/libcore/spectrum.cpp +++ b/src/libcore/spectrum.cpp @@ -1,7 +1,107 @@ +#include +#include #include NAMESPACE_BEGIN(mitsuba) + +template +void spectrum_from_file(const std::string &filename, std::vector &wavelengths, + std::vector &values) { + auto fs = Thread::thread()->file_resolver(); + fs::path file_path = fs->resolve(filename); + if (!fs::exists(file_path)) + Log(Error, "\"%s\": file does not exist!", file_path); + + Log(Info, "Loading spectral data file \"%s\" ..", file_path); + ref file = new FileStream(file_path); + + std::string line, rest; + Scalar wav, value; + while (true) { + try { + line = file->read_line(); + if (line.size() == 0 || line[0] == '#') + continue; + + std::istringstream iss(line); + iss >> wav; + iss >> value; + if (iss >> rest) + Log(Error, "\"%s\": excess tokens after wavlengths-value pair in file:\n%s!", file_path, line); + wavelengths.push_back(wav); + values.push_back(value); + } catch (std::exception &) { + break; + } + } +} + +template +Color spectrum_to_rgb(const std::vector &wavelengths, + const std::vector &values, bool bounded) { + Color color = (Scalar) 0.f; + + const int steps = 1000; + for (int i = 0; i < steps; ++i) { + Scalar x = (Scalar) MTS_WAVELENGTH_MIN + + (i / (Scalar)(steps - 1)) * ((Scalar) MTS_WAVELENGTH_MAX - + (Scalar) MTS_WAVELENGTH_MIN); + + if (x < wavelengths.front() || x > wavelengths.back()) + continue; + + // Find interval containing 'x' + uint32_t index = math::find_interval( + (uint32_t) wavelengths.size(), + [&](uint32_t idx) { + return wavelengths[idx] <= x; + }); + + Scalar x0 = wavelengths[index], + x1 = wavelengths[index + 1], + y0 = values[index], + y1 = values[index + 1]; + + // Linear interpolant at 'x' + Scalar y = (x*y0 - x1*y0 - x*y1 + x0*y1) / (x0 - x1); + + Color xyz = cie1931_xyz(x); + color += xyz * y; + } + + // Last specified value repeats implicitly + color *= ((Scalar) MTS_WAVELENGTH_MAX - (Scalar) MTS_WAVELENGTH_MIN) / (Scalar) steps; + color = xyz_to_srgb(color); + + if (bounded && any(color < (Scalar) 0.f || color > (Scalar) 1.f)) { + Log(Warn, "Spectrum: clamping out-of-gamut color %s", color); + color = clamp(color, (Scalar) 0.f, (Scalar) 1.f); + } else if (!bounded && any(color < (Scalar) 0.f)) { + Log(Warn, "Spectrum: clamping out-of-gamut color %s", color); + color = max(color, (Scalar) 0.f); + } + + return color; +} + + + +/// Explicit instantiations +template MTS_EXPORT void spectrum_from_file(const std::string &filename, + std::vector &wavelengths, + std::vector &values); +template MTS_EXPORT void spectrum_from_file(const std::string &filename, + std::vector &wavelengths, + std::vector &values); + +template MTS_EXPORT Color +spectrum_to_rgb(const std::vector &wavelengths, + const std::vector &values, bool bounded); +template MTS_EXPORT Color +spectrum_to_rgb(const std::vector &wavelengths, + const std::vector &values, bool bounded); + // ======================================================================= //! @{ \name CIE 1931 2 degree observer implementation // =======================================================================