diff --git a/README.md b/README.md index 79786ff..0e645db 100644 --- a/README.md +++ b/README.md @@ -63,7 +63,7 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum * Installation with CPM ```cmake - CPMAddPackage("gh:viennatools/viennaray@3.1.0") + CPMAddPackage("gh:viennatools/viennaray@3.1.1") ``` * With a local installation @@ -101,7 +101,7 @@ ctest --test-dir build ## Contributing -If you want to contribute to ViennaRay, make sure to follow the [LLVM Coding guidelines](https://llvm.org/docs/CodingStandards.html). Before creating a pull request, make sure ALL files have been formatted by clang-format, which can be done using the `format-project.sh` script in the root directory. +If you want to contribute to ViennaRay, make sure to follow the [LLVM Coding guidelines](https://llvm.org/docs/CodingStandards.html). Before creating a pull request, make sure ALL files have been formatted by clang-format. ## Authors diff --git a/include/viennaray/rayReflection.hpp b/include/viennaray/rayReflection.hpp index 3d08e32..1bf09a5 100644 --- a/include/viennaray/rayReflection.hpp +++ b/include/viennaray/rayReflection.hpp @@ -60,28 +60,31 @@ template // (https://math.stackexchange.com/a/182936) std::uniform_real_distribution uniDist; - Vec3D direction; + assert(maxConeAngle >= 0. && maxConeAngle <= M_PI / 2. && + "Cone angle not allowed"); + Vec3D direction; + // compute specular direction + auto const dirOldInv = Inv(rayDir); + auto const specDirection = + (2 * DotProduct(geomNormal, dirOldInv) * geomNormal) - dirOldInv; + auto const basis = getOrthonormalBasis(specDirection); + + // theta from coned cosine distribution (accept reject sampling) + NumericType u, sqrt_1m_u, theta; + do { // generate a random angle between 0 and specular angle + u = std::sqrt(uniDist(rngState)); // x + sqrt_1m_u = std::sqrt(1. - u); // y + theta = maxConeAngle * sqrt_1m_u; // theta + } while (uniDist(rngState) * theta * u > + std::cos(M_PI_2 * sqrt_1m_u) * std::sin(theta)); + NumericType cosTheta = std::cos(theta); + NumericType sinTheta = std::sin(theta); + + // random rotation do { // sample phi uniformly in [0, 2pi] NumericType phi = uniDist(rngState) * 2 * M_PI; - // theta on sphere - assert(maxConeAngle >= 0. && maxConeAngle <= M_PI / 2. && - "Cone angle not allowed"); - NumericType cosTheta = std::cos(maxConeAngle); - // sample z uniformly on [cos(theta),1] - NumericType z = uniDist(rngState) * (1 - cosTheta) + cosTheta; - - // compute specular direction - auto dirOldInv = Inv(rayDir); - auto specDirection = - (2 * DotProduct(geomNormal, dirOldInv) * geomNormal) - dirOldInv; - - // rotate - auto basis = getOrthonormalBasis(specDirection); - NumericType theta = std::acos(z); - cosTheta = std::cos(theta); - NumericType sinTheta = std::sin(theta); NumericType cosPhi = std::cos(phi); NumericType sinPhi = std::sin(phi); @@ -92,7 +95,7 @@ template direction[2] = sinTheta * (cosPhi * basis[1][2] + sinPhi * basis[2][2]) + cosTheta * basis[0][2]; - } while (DotProduct(direction, geomNormal) < 0.); + } while (DotProduct(direction, geomNormal) <= 0.); if constexpr (D == 2) { direction[2] = 0; @@ -112,9 +115,9 @@ using namespace viennacore; // Coned cosine reflection (deprecated) template [[nodiscard]] Vec3D -ReflectionConedCosineOld(NumericType avgReflAngle, - const Vec3D &rayDir, - const Vec3D &geomNormal, RNG &rngState) { +ReflectionConedCosineOld(const Vec3D &rayDir, + const Vec3D &geomNormal, RNG &rngState, + const NumericType maxConeAngle) { assert(IsNormalized(geomNormal) && "ReflectionSpecular: Surface normal is not normalized"); @@ -123,7 +126,6 @@ ReflectionConedCosineOld(NumericType avgReflAngle, // Calculate specular direction auto dirOldInv = Inv(rayDir); - auto specDirection = (2 * DotProduct(geomNormal, dirOldInv) * geomNormal) - dirOldInv; @@ -135,57 +137,57 @@ ReflectionConedCosineOld(NumericType avgReflAngle, // loop until ray is reflected away from the surface normal // this loop takes care of the case where part of the cone points // into the geometry + + // accept reject sampling + do { // generate a random angle between 0 and specular angle + u = std::sqrt(uniDist(rngState)); // x + sqrt_1m_u = std::sqrt(1. - u); // y + angle = maxConeAngle * sqrt_1m_u; // theta + } while (uniDist(rngState) * angle * u > + std::cos(M_PI_2 * sqrt_1m_u) * std::sin(angle)); + + NumericType cosTheta = std::cos(angle); + assert(cosTheta >= 0. && "ReflectionConedCosine: cosTheta < 0"); + assert(cosTheta <= 1. + 1e-6 && "ReflectionConedCosine: cosTheta > 1"); + + // Random Azimuthal Rotation + NumericType cosPhi, sinPhi; + NumericType r2; do { - do { // generate a random angle between 0 and specular angle - u = std::sqrt(uniDist(rngState)); - sqrt_1m_u = std::sqrt(1. - u); - angle = avgReflAngle * sqrt_1m_u; - } while (uniDist(rngState) * angle * u > - std::cos(M_PI_2 * sqrt_1m_u) * std::sin(angle)); - - // Random Azimuthal Rotation - NumericType cosTheta = - std::max(std::min(std::cos(angle), NumericType(1)), NumericType(0)); - NumericType cosPhi, sinPhi; - NumericType r2; - - do { - cosPhi = uniDist(rngState) - 0.5; - sinPhi = uniDist(rngState) - 0.5; - r2 = cosPhi * cosPhi + sinPhi * sinPhi; - } while (r2 >= 0.25 || r2 <= std::numeric_limits::epsilon()); - - // Rotate - cosTheta = std::min(cosTheta, NumericType(1)); - - NumericType a0; - NumericType a1; - - if (std::fabs(specDirection[0]) <= std::fabs(specDirection[1])) { - a0 = specDirection[0]; - a1 = specDirection[1]; - } else { - a0 = specDirection[1]; - a1 = specDirection[0]; - } - - const NumericType a0_a0_m1 = 1. - a0 * a0; - const NumericType tmp = std::sqrt( - std::max(NumericType(1) - cosTheta * cosTheta, NumericType(0)) / - (r2 * a0_a0_m1)); - const NumericType tmp_sinPhi = tmp * sinPhi; - const NumericType tmp_cosPhi = tmp * cosPhi; - const NumericType cosTheta_p_a0_tmp_sinPhi = cosTheta + a0 * tmp_sinPhi; - - randomDir[0] = a0 * cosTheta - a0_a0_m1 * tmp_sinPhi; - randomDir[1] = - a1 * cosTheta_p_a0_tmp_sinPhi + specDirection[2] * tmp_cosPhi; - randomDir[2] = - specDirection[2] * cosTheta_p_a0_tmp_sinPhi - a1 * tmp_cosPhi; - - if (a0 != specDirection[0]) - std::swap(randomDir[0], randomDir[1]); - } while (DotProduct(randomDir, geomNormal) <= 0.); + cosPhi = uniDist(rngState) - 0.5; + sinPhi = uniDist(rngState) - 0.5; + r2 = cosPhi * cosPhi + sinPhi * sinPhi; + } while (r2 >= 0.25 || r2 <= std::numeric_limits::epsilon()); + + // Rotate + NumericType a0; + NumericType a1; + + if (std::fabs(specDirection[0]) <= std::fabs(specDirection[1])) { + a0 = specDirection[0]; + a1 = specDirection[1]; + } else { + a0 = specDirection[1]; + a1 = specDirection[0]; + } + + const double a0_a0_m1 = 1. - a0 * a0; + const double tmp = + std::sqrt(std::max(1. - cosTheta * cosTheta, 0.) / (r2 * a0_a0_m1)); + const double tmp_sinPhi = tmp * sinPhi; + const double tmp_cosPhi = tmp * cosPhi; + const NumericType cosTheta_p_a0_tmp_sinPhi = cosTheta + a0 * tmp_sinPhi; + + randomDir[0] = a0 * cosTheta - a0_a0_m1 * tmp_sinPhi; + randomDir[1] = a1 * cosTheta_p_a0_tmp_sinPhi + specDirection[2] * tmp_cosPhi; + randomDir[2] = specDirection[2] * cosTheta_p_a0_tmp_sinPhi - a1 * tmp_cosPhi; + + if (a0 != specDirection[0]) + std::swap(randomDir[0], randomDir[1]); + + assert(DotProduct(randomDir, geomNormal) > 0. && "ReflectionConedCosine: " + "New direction points into " + "the geometry"); if constexpr (D == 2) { randomDir[2] = 0; @@ -202,7 +204,7 @@ template [[nodiscard]] Vec3D ReflectionConedCosineOld2(const Vec3D &rayDir, const Vec3D &geomNormal, RNG &rngState, - NumericType &minAvgConeAngle = 0.) { + NumericType minAvgConeAngle = 0.) { assert(IsNormalized(geomNormal) && "ReflectionSpecular: Surface normal is not normalized"); diff --git a/tests/reflection/CMakeLists.txt b/tests/reflection/CMakeLists.txt new file mode 100644 index 0000000..534667e --- /dev/null +++ b/tests/reflection/CMakeLists.txt @@ -0,0 +1,7 @@ +project(reflection LANGUAGES CXX) + +add_executable(${PROJECT_NAME} "${PROJECT_NAME}.cpp") +target_link_libraries(${PROJECT_NAME} PRIVATE ViennaRay) + +add_dependencies(ViennaRay_Tests ${PROJECT_NAME}) +add_test(NAME ${PROJECT_NAME} COMMAND $) diff --git a/tests/reflection/reflection.cpp b/tests/reflection/reflection.cpp new file mode 100644 index 0000000..63620f1 --- /dev/null +++ b/tests/reflection/reflection.cpp @@ -0,0 +1,93 @@ +#include +#include + +#include + +using namespace viennaray; + +int main() { + + using NumericType = double; + constexpr int N = 50000; + RNG rngState(12351263); + Timer timer; + + { + std::array, N> directions; + + // diffuse reflection + Vec3D normal = {0.0, 0.0, 1.0}; + + timer.start(); + for (int i = 0; i < N; ++i) { + directions[i] = ReflectionDiffuse(normal, rngState); + } + timer.finish(); + + std::cout << "Time for " << N + << " diffuse reflections: " << timer.currentDuration * 1e-6 + << " ms" << std::endl; + + std::ofstream file("diffuse_reflection.txt"); + for (auto const &dir : directions) { + file << dir[0] << " " << dir[1] << " " << dir[2] << std::endl; + } + file.close(); + } + + Vec3D normal = {0.0, 0.0, 1.0}; + NumericType const minAngle = 85.0 * M_PI / 180.0; + NumericType incAngle = 75.0 * M_PI / 180.0; + Vec3D rayDir = {0.0, -std::sin(incAngle), -std::cos(incAngle)}; + auto coneAngle = M_PI_2 - std::min(incAngle, minAngle); + + std::cout << "minAngle: " << minAngle << std::endl; + std::cout << "incAngle: " << incAngle << std::endl; + std::cout << "incAngle [deg]: " << incAngle * 180.0 / M_PI << std::endl; + std::cout << "coneAngle: " << coneAngle << std::endl; + + { + std::array, N> directions; + + // coned specular reflection + timer.start(); + for (int i = 0; i < N; ++i) { + directions[i] = ReflectionConedCosine( + rayDir, normal, rngState, coneAngle); + } + timer.finish(); + + std::cout << "Time for " << N + << " coned specular reflections: " << timer.currentDuration * 1e-6 + << " ms" << std::endl; + + std::ofstream file("coned_specular_reflection.txt"); + for (auto const &dir : directions) { + file << dir[0] << " " << dir[1] << " " << dir[2] << std::endl; + } + file.close(); + } + + { + std::array, N> directions; + + // coned specular reflection + timer.start(); + for (int i = 0; i < N; ++i) { + directions[i] = rayInternal::ReflectionConedCosineOld( + rayDir, normal, rngState, coneAngle); + } + timer.finish(); + + std::cout << "Time for " << N << " coned specular reflections (old): " + << timer.currentDuration * 1e-6 << " ms" << std::endl; + + std::ofstream file("coned_specular_reflection_old.txt"); + for (auto const &dir : directions) { + file << dir[0] << " " << dir[1] << " " << dir[2] << std::endl; + } + file.close(); + } + + return 0; +} \ No newline at end of file