merge

ViennaTools · Feb 4, 2025 · cee1fc3 · cee1fc3
2 parents 8c53721 + ca626e4
commit cee1fc3
Show file tree

Hide file tree

Showing 4 changed files with 178 additions and 76 deletions.
diff --git 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/[email protected].0")
+  CPMAddPackage("gh:viennatools/[email protected].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
@@ -60,28 +60,31 @@ template <typename NumericType, int D>
   // (https://math.stackexchange.com/a/182936)
   std::uniform_real_distribution<NumericType> uniDist;
 
-  Vec3D<NumericType> direction;
+  assert(maxConeAngle >= 0. && maxConeAngle <= M_PI / 2. &&
+         "Cone angle not allowed");
 
+  Vec3D<NumericType> 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 <typename NumericType, int D>
     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 <typename NumericType, int D>
 [[nodiscard]] Vec3D<NumericType>
-ReflectionConedCosineOld(NumericType avgReflAngle,
-                         const Vec3D<NumericType> &rayDir,
-                         const Vec3D<NumericType> &geomNormal, RNG &rngState) {
+ReflectionConedCosineOld(const Vec3D<NumericType> &rayDir,
+                         const Vec3D<NumericType> &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<NumericType>::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<NumericType>::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 <typename NumericType, int D>
 [[nodiscard]] Vec3D<NumericType>
 ReflectionConedCosineOld2(const Vec3D<NumericType> &rayDir,
                           const Vec3D<NumericType> &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
@@ -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 $<TARGET_FILE:${PROJECT_NAME}>)
diff --git a/tests/reflection/reflection.cpp b/tests/reflection/reflection.cpp
@@ -0,0 +1,93 @@
+#include <rayReflection.hpp>
+#include <vcTimer.hpp>
+
+#include <fstream>
+
+using namespace viennaray;
+
+int main() {
+
+  using NumericType = double;
+  constexpr int N = 50000;
+  RNG rngState(12351263);
+  Timer timer;
+
+  {
+    std::array<Vec3D<NumericType>, N> directions;
+
+    // diffuse reflection
+    Vec3D<NumericType> normal = {0.0, 0.0, 1.0};
+
+    timer.start();
+    for (int i = 0; i < N; ++i) {
+      directions[i] = ReflectionDiffuse<NumericType, 3>(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<NumericType> 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<NumericType> 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<Vec3D<NumericType>, N> directions;
+
+    // coned specular reflection
+    timer.start();
+    for (int i = 0; i < N; ++i) {
+      directions[i] = ReflectionConedCosine<NumericType, 3>(
+          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<Vec3D<NumericType>, N> directions;
+
+    // coned specular reflection
+    timer.start();
+    for (int i = 0; i < N; ++i) {
+      directions[i] = rayInternal::ReflectionConedCosineOld<NumericType, 3>(
+          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;
+}