From e8da29faa0095328d0043fbf0e03d4508aafe2cc Mon Sep 17 00:00:00 2001 From: Tobias Reiter <44025882+tobre1@users.noreply.github.com> Date: Thu, 11 Apr 2024 20:15:10 +0200 Subject: [PATCH] Major code refactor (#53) * Remove rayMetaGeometry and refactor rayBoundary * Refactor rayGeometry * Refactor rayHitCounter * Refactor rayMessage * Remove rayPreCompileMacros header * Refactor raySource * Refactor rayTrace * Refactor rayTraceKernel * Refactor rayUtil * Add option to enable progress bar * Move rayTraceDirection and rayDataLog to rayUtil * Bump version * Add missing header * Enable _USE_MATH_DEFINES * Test _USE_MATH_DEFINES * Test _USE_MATH_DEFINES again * Define M_PI if not defined * Remove old assert --- CMakeLists.txt | 12 +- README.md | 2 +- include/viennaray/rayBoundary.hpp | 275 ++++++++++----------- include/viennaray/rayDataLog.hpp | 24 -- include/viennaray/rayGeometry.hpp | 250 +++++++++---------- include/viennaray/rayHitCounter.hpp | 158 ++++++------ include/viennaray/rayMessage.hpp | 32 ++- include/viennaray/rayMetaGeometry.hpp | 28 --- include/viennaray/rayParticle.hpp | 6 +- include/viennaray/rayPreCompileMacros.hpp | 9 - include/viennaray/rayRNG.hpp | 7 +- include/viennaray/rayReflection.hpp | 21 +- include/viennaray/raySource.hpp | 28 +-- include/viennaray/raySourceGrid.hpp | 68 ++---- include/viennaray/raySourceRandom.hpp | 115 +++++---- include/viennaray/rayTestAsserts.hpp | 3 +- include/viennaray/rayTrace.hpp | 235 +++++++++--------- include/viennaray/rayTraceDirection.hpp | 13 - include/viennaray/rayTraceKernel.hpp | 282 +++++++++++----------- include/viennaray/rayTracingData.hpp | 173 ++++++------- include/viennaray/rayUtil.hpp | 86 +++++-- tests/buildBoundary/buildBoundary.cpp | 18 -- tests/buildBoundary2D/buildBoundary2D.cpp | 14 +- tests/diskAreas/diskAreas.cpp | 9 +- 24 files changed, 875 insertions(+), 993 deletions(-) delete mode 100644 include/viennaray/rayDataLog.hpp delete mode 100644 include/viennaray/rayMetaGeometry.hpp delete mode 100644 include/viennaray/rayPreCompileMacros.hpp delete mode 100644 include/viennaray/rayTraceDirection.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 576d9a0..0c472a2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) project( ViennaRay LANGUAGES CXX - VERSION 2.0.4) + VERSION 2.1.0) # -------------------------------------------------------------------------------------------------------- # Library switches @@ -15,10 +15,12 @@ option(VIENNARAY_DISABLE_COPY "Disable Embree Environment setup" OFF) option(VIENNARAY_LEGACY_OPENMP "Use legacy OpenMP flags" OFF) option(VIENNARAY_USE_WDIST "Enable weighted distribution of ray weights" OFF) +option(VIENNARAY_PRINT_PROGRESS "Enable progress bar" OFF) + option(VIENNARAY_BUILD_EXAMPLES "Build examples" OFF) option(VIENNARAY_BUILD_TESTS "Build tests" OFF) -if (VIENNARAY_LEGACY_OPENMP) +if(VIENNARAY_LEGACY_OPENMP) message(STATUS "[ViennaRay] Using legacy OpenMP flags") find_package(OpenMP REQUIRED) @@ -93,6 +95,10 @@ if(VIENNARAY_USE_WDIST) target_compile_definitions(${PROJECT_NAME} INTERFACE VIENNARAY_USE_WDIST) endif() +if(VIENNARAY_PRINT_PROGRESS) + target_compile_definitions(${PROJECT_NAME} INTERFACE VIENNARAY_PRINT_PROGRESS) +endif() + if(EMBREE_RAY_MASK) message(STATUS "[ViennaRay] Using embree ray masking") target_compile_definitions(${PROJECT_NAME} INTERFACE VIENNARAY_USE_RAY_MASKING) @@ -133,7 +139,7 @@ CPMFindPackage( GIT_REPOSITORY "https://github.com/embree/embree" OPTIONS "EMBREE_TUTORIALS OFF") -if (NOT VIENNARAY_LEGACY_OPENMP) +if(NOT VIENNARAY_LEGACY_OPENMP) find_package(OpenMP REQUIRED) target_link_libraries(${PROJECT_NAME} INTERFACE OpenMP::OpenMP_CXX) endif() diff --git a/README.md b/README.md index b0908da..b156382 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@2.0.1") + CPMAddPackage("gh:viennatools/viennaray@2.1.0") ``` * With a local installation diff --git a/include/viennaray/rayBoundary.hpp b/include/viennaray/rayBoundary.hpp index 4687264..075481f 100644 --- a/include/viennaray/rayBoundary.hpp +++ b/include/viennaray/rayBoundary.hpp @@ -1,10 +1,7 @@ -#ifndef RAY_BOUNDARY_HPP -#define RAY_BOUNDARY_HPP +#pragma once -#include -#include #include -#include +#include enum class rayBoundaryCondition : unsigned { REFLECTIVE = 0, @@ -12,41 +9,40 @@ enum class rayBoundaryCondition : unsigned { IGNORE = 2 }; -template -class rayBoundary : public rayMetaGeometry { - typedef rayPair> boundingBoxType; +template class rayBoundary { + using boundingBoxType = rayPair>; public: - rayBoundary(RTCDevice &pDevice, const boundingBoxType &pBoundingBox, - rayBoundaryCondition pBoundaryConds[D], - std::array &pTraceSettings) - : mbdBox(pBoundingBox), firstDir(pTraceSettings[1]), - secondDir(pTraceSettings[2]), - mBoundaryConds({pBoundaryConds[firstDir], pBoundaryConds[secondDir]}) { - initBoundary(pDevice); + rayBoundary(RTCDevice &device, boundingBoxType const &boundingBox, + rayBoundaryCondition boundaryConds[D], + std::array &traceSettings) + : bdBox_(boundingBox), firstDir_(traceSettings[1]), + secondDir_(traceSettings[2]), + boundaryConds_({boundaryConds[firstDir_], boundaryConds[secondDir_]}) { + initBoundary(device); } - void processHit(RTCRayHit &rayHit, bool &reflect) { + void processHit(RTCRayHit &rayHit, bool &reflect) const { const auto primID = rayHit.hit.primID; if constexpr (D == 2) { assert((primID == 0 || primID == 1 || primID == 2 || primID == 3) && "Assumption"); - if (mBoundaryConds[0] == rayBoundaryCondition::REFLECTIVE) { + if (boundaryConds_[0] == rayBoundaryCondition::REFLECTIVE) { reflectRay(rayHit); reflect = true; return; - } else if (mBoundaryConds[0] == rayBoundaryCondition::PERIODIC) { - auto impactCoords = this->getNewOrigin(rayHit.ray); + } else if (boundaryConds_[0] == rayBoundaryCondition::PERIODIC) { + auto impactCoords = getNewOrigin(rayHit.ray); // periodically move ray origin if (primID == 0 || primID == 1) { // hit at x/y min boundary -> move to max x/y - impactCoords[firstDir] = mbdBox[1][firstDir]; + impactCoords[firstDir_] = bdBox_[1][firstDir_]; } else if (primID == 2 || primID == 3) { // hit at x/y max boundary -> move to min x/y - impactCoords[firstDir] = mbdBox[0][firstDir]; + impactCoords[firstDir_] = bdBox_[0][firstDir_]; } - moveRay(rayHit, impactCoords); + assignRayCoords(rayHit, impactCoords); reflect = true; return; } else { @@ -57,23 +53,23 @@ class rayBoundary : public rayMetaGeometry { assert(false && "Correctness Assumption"); } else { - if (primID == 0 || primID == 1 || primID == 2 || primID == 3) { - if (mBoundaryConds[0] == rayBoundaryCondition::REFLECTIVE) { + if (primID >= 0 && primID <= 3) { + if (boundaryConds_[0] == rayBoundaryCondition::REFLECTIVE) { // use specular reflection reflectRay(rayHit); reflect = true; return; - } else if (mBoundaryConds[0] == rayBoundaryCondition::PERIODIC) { - auto impactCoords = this->getNewOrigin(rayHit.ray); + } else if (boundaryConds_[0] == rayBoundaryCondition::PERIODIC) { + auto impactCoords = getNewOrigin(rayHit.ray); // periodically move ray origin if (primID == 0 || primID == 1) { // hit at firstDir min boundary -> move to max firstDir - impactCoords[firstDir] = mbdBox[1][firstDir]; + impactCoords[firstDir_] = bdBox_[1][firstDir_]; } else if (primID == 2 || primID == 3) { // hit at firstDir max boundary -> move to min fristDir - impactCoords[firstDir] = mbdBox[0][firstDir]; + impactCoords[firstDir_] = bdBox_[0][firstDir_]; } - moveRay(rayHit, impactCoords); + assignRayCoords(rayHit, impactCoords); reflect = true; return; } else { @@ -81,23 +77,23 @@ class rayBoundary : public rayMetaGeometry { reflect = false; return; } - } else if (primID == 4 || primID == 5 || primID == 6 || primID == 7) { - if (mBoundaryConds[1] == rayBoundaryCondition::REFLECTIVE) { + } else if (primID >= 4 && primID <= 7) { + if (boundaryConds_[1] == rayBoundaryCondition::REFLECTIVE) { // use specular reflection reflectRay(rayHit); reflect = true; return; - } else if (mBoundaryConds[1] == rayBoundaryCondition::PERIODIC) { - auto impactCoords = this->getNewOrigin(rayHit.ray); + } else if (boundaryConds_[1] == rayBoundaryCondition::PERIODIC) { + auto impactCoords = getNewOrigin(rayHit.ray); // periodically move ray origin if (primID == 4 || primID == 5) { // hit at secondDir min boundary -> move to max secondDir - impactCoords[secondDir] = mbdBox[1][secondDir]; + impactCoords[secondDir_] = bdBox_[1][secondDir_]; } else if (primID == 6 || primID == 7) { // hit at secondDir max boundary -> move to min secondDir - impactCoords[secondDir] = mbdBox[0][secondDir]; + impactCoords[secondDir_] = bdBox_[0][secondDir_]; } - moveRay(rayHit, impactCoords); + assignRayCoords(rayHit, impactCoords); reflect = true; return; } else { @@ -111,87 +107,90 @@ class rayBoundary : public rayMetaGeometry { } } - RTCGeometry &getRTCGeometry() override final { return mRtcBoundary; } + RTCGeometry const &getRTCGeometry() const { return pRtcBoundary_; } void releaseGeometry() { // Attention: // This function must not be called when the RTCGeometry reference count is // > 1 Doing so leads to leaked memory buffers - if (mTriangleBuffer == nullptr || mVertexBuffer == nullptr || - mRtcBoundary == nullptr) { + if (pTriangleBuffer_ == nullptr || pVertexBuffer_ == nullptr || + pRtcBoundary_ == nullptr) { return; } else { - rtcReleaseGeometry(mRtcBoundary); - mRtcBoundary = nullptr; - mTriangleBuffer = nullptr; - mVertexBuffer = nullptr; + rtcReleaseGeometry(pRtcBoundary_); + pRtcBoundary_ = nullptr; + pTriangleBuffer_ = nullptr; + pVertexBuffer_ = nullptr; } } - rayTriple getPrimNormal(const unsigned int primID) override { - assert(primID < numTriangles && "rayBoundary: primID out of bounds"); - return primNormals[primID]; - } - - boundingBoxType getBoundingBox() const { return mbdBox; } - - rayPair getDirs() const { return {firstDir, secondDir}; } + rayPair getDirs() const { return {firstDir_, secondDir_}; } private: + static rayTriple getNewOrigin(RTCRay &ray) { + assert(rayInternal::IsNormalized( + rayTriple{ray.dir_x, ray.dir_y, ray.dir_z}) && + "MetaGeometry: direction not normalized"); + auto xx = ray.org_x + ray.dir_x * ray.tfar; + auto yy = ray.org_y + ray.dir_y * ray.tfar; + auto zz = ray.org_z + ray.dir_z * ray.tfar; + return {xx, yy, zz}; + } + void initBoundary(RTCDevice &pDevice) { - assert(mVertexBuffer == nullptr && mTriangleBuffer == nullptr && + assert(pVertexBuffer_ == nullptr && pTriangleBuffer_ == nullptr && "Boundary buffer not empty"); - mRtcBoundary = rtcNewGeometry(pDevice, RTC_GEOMETRY_TYPE_TRIANGLE); + pRtcBoundary_ = rtcNewGeometry(pDevice, RTC_GEOMETRY_TYPE_TRIANGLE); - mVertexBuffer = (vertex_f3_t *)rtcSetNewGeometryBuffer( - mRtcBoundary, RTC_BUFFER_TYPE_VERTEX, + pVertexBuffer_ = (vertex_f3_t *)rtcSetNewGeometryBuffer( + pRtcBoundary_, RTC_BUFFER_TYPE_VERTEX, 0, // the slot - RTC_FORMAT_FLOAT3, sizeof(vertex_f3_t), numVertices); + RTC_FORMAT_FLOAT3, sizeof(vertex_f3_t), numVertices_); - auto xmin = mbdBox[0][0]; - auto xmax = mbdBox[1][0]; - auto ymin = mbdBox[0][1]; - auto ymax = mbdBox[1][1]; - auto zmin = mbdBox[0][2]; - auto zmax = mbdBox[1][2]; + auto xmin = bdBox_[0][0]; + auto xmax = bdBox_[1][0]; + auto ymin = bdBox_[0][1]; + auto ymax = bdBox_[1][1]; + auto zmin = bdBox_[0][2]; + auto zmax = bdBox_[1][2]; // Vertices - mVertexBuffer[0].xx = (float)xmin; - mVertexBuffer[0].yy = (float)ymin; - mVertexBuffer[0].zz = (float)zmin; + pVertexBuffer_[0].xx = (float)xmin; + pVertexBuffer_[0].yy = (float)ymin; + pVertexBuffer_[0].zz = (float)zmin; - mVertexBuffer[1].xx = (float)xmax; - mVertexBuffer[1].yy = (float)ymin; - mVertexBuffer[1].zz = (float)zmin; + pVertexBuffer_[1].xx = (float)xmax; + pVertexBuffer_[1].yy = (float)ymin; + pVertexBuffer_[1].zz = (float)zmin; - mVertexBuffer[2].xx = (float)xmax; - mVertexBuffer[2].yy = (float)ymax; - mVertexBuffer[2].zz = (float)zmin; + pVertexBuffer_[2].xx = (float)xmax; + pVertexBuffer_[2].yy = (float)ymax; + pVertexBuffer_[2].zz = (float)zmin; - mVertexBuffer[3].xx = (float)xmin; - mVertexBuffer[3].yy = (float)ymax; - mVertexBuffer[3].zz = (float)zmin; + pVertexBuffer_[3].xx = (float)xmin; + pVertexBuffer_[3].yy = (float)ymax; + pVertexBuffer_[3].zz = (float)zmin; - mVertexBuffer[4].xx = (float)xmin; - mVertexBuffer[4].yy = (float)ymin; - mVertexBuffer[4].zz = (float)zmax; + pVertexBuffer_[4].xx = (float)xmin; + pVertexBuffer_[4].yy = (float)ymin; + pVertexBuffer_[4].zz = (float)zmax; - mVertexBuffer[5].xx = (float)xmax; - mVertexBuffer[5].yy = (float)ymin; - mVertexBuffer[5].zz = (float)zmax; + pVertexBuffer_[5].xx = (float)xmax; + pVertexBuffer_[5].yy = (float)ymin; + pVertexBuffer_[5].zz = (float)zmax; - mVertexBuffer[6].xx = (float)xmax; - mVertexBuffer[6].yy = (float)ymax; - mVertexBuffer[6].zz = (float)zmax; + pVertexBuffer_[6].xx = (float)xmax; + pVertexBuffer_[6].yy = (float)ymax; + pVertexBuffer_[6].zz = (float)zmax; - mVertexBuffer[7].xx = (float)xmin; - mVertexBuffer[7].yy = (float)ymax; - mVertexBuffer[7].zz = (float)zmax; + pVertexBuffer_[7].xx = (float)xmin; + pVertexBuffer_[7].yy = (float)ymax; + pVertexBuffer_[7].zz = (float)zmax; - mTriangleBuffer = (triangle_t *)rtcSetNewGeometryBuffer( - mRtcBoundary, RTC_BUFFER_TYPE_INDEX, + pTriangleBuffer_ = (triangle_t *)rtcSetNewGeometryBuffer( + pRtcBoundary_, RTC_BUFFER_TYPE_INDEX, 0, // slot - RTC_FORMAT_UINT3, sizeof(triangle_t), numTriangles); + RTC_FORMAT_UINT3, sizeof(triangle_t), numTriangles_); constexpr rayQuadruple> xMinMaxPlanes = { 0, 3, 7, 0, 7, 4, 6, 2, 1, 6, 1, 5}; @@ -203,62 +202,53 @@ class rayBoundary : public rayMetaGeometry { xMinMaxPlanes, yMinMaxPlanes, zMinMaxPlanes}; for (size_t idx = 0; idx < 4; ++idx) { - mTriangleBuffer[idx].v0 = Planes[firstDir][idx][0]; - mTriangleBuffer[idx].v1 = Planes[firstDir][idx][1]; - mTriangleBuffer[idx].v2 = Planes[firstDir][idx][2]; + pTriangleBuffer_[idx].v0 = Planes[firstDir_][idx][0]; + pTriangleBuffer_[idx].v1 = Planes[firstDir_][idx][1]; + pTriangleBuffer_[idx].v2 = Planes[firstDir_][idx][2]; - mTriangleBuffer[idx + 4].v0 = Planes[secondDir][idx][0]; - mTriangleBuffer[idx + 4].v1 = Planes[secondDir][idx][1]; - mTriangleBuffer[idx + 4].v2 = Planes[secondDir][idx][2]; - } - - for (size_t idx = 0; idx < numTriangles; ++idx) { - auto triangle = getTriangleCoords(idx); - auto triNorm = rayInternal::ComputeNormal(triangle); - rayInternal::Normalize(triNorm); - primNormals[idx] = triNorm; + pTriangleBuffer_[idx + 4].v0 = Planes[secondDir_][idx][0]; + pTriangleBuffer_[idx + 4].v1 = Planes[secondDir_][idx][1]; + pTriangleBuffer_[idx + 4].v2 = Planes[secondDir_][idx][2]; } #ifdef VIENNARAY_USE_RAY_MASKING - rtcSetGeometryMask(mRtcBoundary, -1); + rtcSetGeometryMask(pRtcBoundary_, -1); #endif - rtcCommitGeometry(mRtcBoundary); + rtcCommitGeometry(pRtcBoundary_); assert(rtcGetDeviceError(pDevice) == RTC_ERROR_NONE && "RTC Error: rtcCommitGeometry"); } rayTriple> getTriangleCoords(const size_t primID) { - assert(primID < numTriangles && "rtBounday: primID out of bounds"); - auto tt = mTriangleBuffer[primID]; - return {(NumericType)mVertexBuffer[tt.v0].xx, - (NumericType)mVertexBuffer[tt.v0].yy, - (NumericType)mVertexBuffer[tt.v0].zz, - (NumericType)mVertexBuffer[tt.v1].xx, - (NumericType)mVertexBuffer[tt.v1].yy, - (NumericType)mVertexBuffer[tt.v1].zz, - (NumericType)mVertexBuffer[tt.v2].xx, - (NumericType)mVertexBuffer[tt.v2].yy, - (NumericType)mVertexBuffer[tt.v2].zz}; + assert(primID < numTriangles_ && "rtBounday: primID out of bounds"); + auto tt = pTriangleBuffer_[primID]; + return {(NumericType)pVertexBuffer_[tt.v0].xx, + (NumericType)pVertexBuffer_[tt.v0].yy, + (NumericType)pVertexBuffer_[tt.v0].zz, + (NumericType)pVertexBuffer_[tt.v1].xx, + (NumericType)pVertexBuffer_[tt.v1].yy, + (NumericType)pVertexBuffer_[tt.v1].zz, + (NumericType)pVertexBuffer_[tt.v2].xx, + (NumericType)pVertexBuffer_[tt.v2].yy, + (NumericType)pVertexBuffer_[tt.v2].zz}; } - void reflectRay(RTCRayHit &rayHit) { - auto dir = - *reinterpret_cast *>(&rayHit.ray.dir_x); - auto normal = - *reinterpret_cast *>(&rayHit.hit.Ng_x); + void reflectRay(RTCRayHit &rayHit) const { + auto dir = *reinterpret_cast *>( + &rayHit.ray.dir_x); + auto normal = *reinterpret_cast *>( + &rayHit.hit.Ng_x); rayInternal::Normalize(dir); rayInternal::Normalize(normal); - dir = rayReflectionSpecular(dir, normal); + dir = rayReflectionSpecular(dir, normal); // normal gets reused for new origin here - normal = this->getNewOrigin(rayHit.ray); + normal = getNewOrigin(rayHit.ray); #ifdef ARCH_X86 reinterpret_cast<__m128 &>(rayHit.ray) = - _mm_set_ps(1e-4f, (rtcNumericType)normal[2], (rtcNumericType)normal[1], - (rtcNumericType)normal[0]); + _mm_set_ps(1e-4f, normal[2], normal[1], normal[0]); reinterpret_cast<__m128 &>(rayHit.ray.dir_x) = - _mm_set_ps(0.0f, (rtcNumericType)dir[2], (rtcNumericType)dir[1], - (rtcNumericType)dir[0]); + _mm_set_ps(0.0f, dir[2], dir[1], dir[0]); #else rayHit.ray.org_x = normal[0]; rayHit.ray.org_y = normal[1]; @@ -272,11 +262,14 @@ class rayBoundary : public rayMetaGeometry { #endif } - void moveRay(RTCRayHit &rayHit, const rayTriple &coords) { + void + assignRayCoords(RTCRayHit &rayHit, + const rayTriple &coords) const { #ifdef ARCH_X86 reinterpret_cast<__m128 &>(rayHit.ray) = - _mm_set_ps(1e-4f, (rtcNumericType)coords[2], (rtcNumericType)coords[1], - (rtcNumericType)coords[0]); + _mm_set_ps(1e-4f, (rayInternal::rtcNumericType)coords[2], + (rayInternal::rtcNumericType)coords[1], + (rayInternal::rtcNumericType)coords[0]); rayHit.ray.time = 0.0f; #else rayHit.ray.org_x = coords[0]; @@ -293,23 +286,21 @@ class rayBoundary : public rayMetaGeometry { // in single precision floating point types. float xx, yy, zz; }; - vertex_f3_t *mVertexBuffer = nullptr; + vertex_f3_t *pVertexBuffer_ = nullptr; struct triangle_t { // The triangle geometry uses an index buffer that contains an array // of three 32-bit indices per triangle. uint32_t v0, v1, v2; }; - triangle_t *mTriangleBuffer = nullptr; - - RTCGeometry mRtcBoundary = nullptr; - const boundingBoxType mbdBox; - const int firstDir = 0; - const int secondDir = 1; - const std::array mBoundaryConds = {}; - static constexpr size_t numTriangles = 8; - static constexpr size_t numVertices = 8; - std::array, numTriangles> primNormals; -}; + triangle_t *pTriangleBuffer_ = nullptr; -#endif // RAY_BOUNDARY_HPP \ No newline at end of file +private: + RTCGeometry pRtcBoundary_ = nullptr; + boundingBoxType const &bdBox_; + const int firstDir_ = 0; + const int secondDir_ = 1; + const std::array boundaryConds_; + static constexpr size_t numTriangles_ = 8; + static constexpr size_t numVertices_ = 8; +}; diff --git a/include/viennaray/rayDataLog.hpp b/include/viennaray/rayDataLog.hpp deleted file mode 100644 index 882b1d6..0000000 --- a/include/viennaray/rayDataLog.hpp +++ /dev/null @@ -1,24 +0,0 @@ -#ifndef RAY_DATA_LOG -#define RAY_DATA_LOG - -#include -#include - -template struct rayDataLog { - - std::vector> data; - - void merge(rayDataLog &pOther) { - assert(pOther.data.size() == data.size() && - "Size mismatch when merging logs"); - for (std::size_t i = 0; i < data.size(); i++) { - assert(pOther.data[i].size() == data[i].size() && - "Size mismatch when merging log data"); - for (std::size_t j = 0; j < data[i].size(); j++) { - data[i][j] += pOther.data[i][j]; - } - } - } -}; - -#endif // RAY_DATA_LOG \ No newline at end of file diff --git a/include/viennaray/rayGeometry.hpp b/include/viennaray/rayGeometry.hpp index 1fb74af..d228fd5 100644 --- a/include/viennaray/rayGeometry.hpp +++ b/include/viennaray/rayGeometry.hpp @@ -1,19 +1,13 @@ -#ifndef RAY_GEOMETRY_HPP -#define RAY_GEOMETRY_HPP +#pragma once -#include #include -template -class rayGeometry : public rayMetaGeometry { -private: - typedef std::vector> pointNeighborhoodType; +template class rayGeometry { + using pointNeighborhoodType = std::vector>; public: - rayGeometry() {} - template - void initGeometry(RTCDevice &pDevice, + void initGeometry(RTCDevice &device, std::vector> &points, std::vector> &normals, NumericType discRadii) { @@ -26,145 +20,143 @@ class rayGeometry : public rayMetaGeometry { // overwriting the geometry without releasing it beforehand causes the old // buffer to leak releaseGeometry(); - mRTCGeometry = - rtcNewGeometry(pDevice, RTC_GEOMETRY_TYPE_ORIENTED_DISC_POINT); - assert(rtcGetDeviceError(pDevice) == RTC_ERROR_NONE && + pRtcGeometry_ = + rtcNewGeometry(device, RTC_GEOMETRY_TYPE_ORIENTED_DISC_POINT); + assert(rtcGetDeviceError(device) == RTC_ERROR_NONE && "RTC Error: rtcNewGeometry"); - mNumPoints = points.size(); + numPoints_ = points.size(); // The buffer data is managed internally (embree) and automatically freed // when the geometry is destroyed. - mPointBuffer = (point_4f_t *)rtcSetNewGeometryBuffer( - mRTCGeometry, RTC_BUFFER_TYPE_VERTEX, + pPointBuffer_ = (point_4f_t *)rtcSetNewGeometryBuffer( + pRtcGeometry_, RTC_BUFFER_TYPE_VERTEX, 0, // slot - RTC_FORMAT_FLOAT4, sizeof(point_4f_t), mNumPoints); - assert(rtcGetDeviceError(pDevice) == RTC_ERROR_NONE && + RTC_FORMAT_FLOAT4, sizeof(point_4f_t), numPoints_); + assert(rtcGetDeviceError(device) == RTC_ERROR_NONE && "RTC Error: rtcSetNewGeometryBuffer points"); - mDiscRadii = discRadii; + discRadii_ = discRadii; for (int i = 0; i < D; i++) { - mMinCoords[i] = nummax; - mMaxCoords[i] = nummin; + minCoords_[i] = std::numeric_limits::max(); + maxCoords_[i] = std::numeric_limits::lowest(); } - for (size_t i = 0; i < mNumPoints; ++i) { - mPointBuffer[i].xx = (float)points[i][0]; - mPointBuffer[i].yy = (float)points[i][1]; - mPointBuffer[i].radius = (float)discRadii; - if (points[i][0] < mMinCoords[0]) - mMinCoords[0] = points[i][0]; - if (points[i][1] < mMinCoords[1]) - mMinCoords[1] = points[i][1]; - if (points[i][0] > mMaxCoords[0]) - mMaxCoords[0] = points[i][0]; - if (points[i][1] > mMaxCoords[1]) - mMaxCoords[1] = points[i][1]; + for (size_t i = 0; i < numPoints_; ++i) { + pPointBuffer_[i].xx = (float)points[i][0]; + pPointBuffer_[i].yy = (float)points[i][1]; + pPointBuffer_[i].radius = (float)discRadii_; + if (points[i][0] < minCoords_[0]) + minCoords_[0] = points[i][0]; + if (points[i][1] < minCoords_[1]) + minCoords_[1] = points[i][1]; + if (points[i][0] > maxCoords_[0]) + maxCoords_[0] = points[i][0]; + if (points[i][1] > maxCoords_[1]) + maxCoords_[1] = points[i][1]; if constexpr (D == 2) { - mPointBuffer[i].zz = 0.f; - mMinCoords[2] = 0.; - mMaxCoords[2] = 0.; + pPointBuffer_[i].zz = 0.f; + minCoords_[2] = 0.; + maxCoords_[2] = 0.; } else { - mPointBuffer[i].zz = (float)points[i][2]; - if (points[i][2] < mMinCoords[2]) - mMinCoords[2] = points[i][2]; - if (points[i][2] > mMaxCoords[2]) - mMaxCoords[2] = points[i][2]; + pPointBuffer_[i].zz = (float)points[i][2]; + if (points[i][2] < minCoords_[2]) + minCoords_[2] = points[i][2]; + if (points[i][2] > maxCoords_[2]) + maxCoords_[2] = points[i][2]; } } - mNormalVecBuffer = (normal_vec_3f_t *)rtcSetNewGeometryBuffer( - mRTCGeometry, RTC_BUFFER_TYPE_NORMAL, + pNormalVecBuffer_ = (normal_vec_3f_t *)rtcSetNewGeometryBuffer( + pRtcGeometry_, RTC_BUFFER_TYPE_NORMAL, 0, // slot - RTC_FORMAT_FLOAT3, sizeof(normal_vec_3f_t), mNumPoints); - assert(rtcGetDeviceError(pDevice) == RTC_ERROR_NONE && + RTC_FORMAT_FLOAT3, sizeof(normal_vec_3f_t), numPoints_); + assert(rtcGetDeviceError(device) == RTC_ERROR_NONE && "RTC Error: rtcSetNewGeometryBuffer normals"); - for (size_t i = 0; i < mNumPoints; ++i) { - mNormalVecBuffer[i].xx = (float)normals[i][0]; - mNormalVecBuffer[i].yy = (float)normals[i][1]; + for (size_t i = 0; i < numPoints_; ++i) { + pNormalVecBuffer_[i].xx = (float)normals[i][0]; + pNormalVecBuffer_[i].yy = (float)normals[i][1]; if constexpr (D == 2) { - mNormalVecBuffer[i].zz = 0.f; + pNormalVecBuffer_[i].zz = 0.f; } else { - mNormalVecBuffer[i].zz = (float)normals[i][2]; + pNormalVecBuffer_[i].zz = (float)normals[i][2]; } } #ifdef VIENNARAY_USE_RAY_MASKING - rtcSetGeometryMask(mRTCGeometry, -1); + rtcSetGeometryMask(pRtcGeometry_, -1); #endif - rtcCommitGeometry(mRTCGeometry); - assert(rtcGetDeviceError(pDevice) == RTC_ERROR_NONE && + rtcCommitGeometry(pRtcGeometry_); + assert(rtcGetDeviceError(device) == RTC_ERROR_NONE && "RTC Error: rtcCommitGeometry"); initPointNeighborhood(points); - if (mMaterialIds.empty()) { - mMaterialIds.resize(mNumPoints, 0); + if (materialIds_.empty()) { + materialIds_.resize(numPoints_, 0); } } template void setMaterialIds(std::vector &pMaterialIds) { - assert(pMaterialIds.size() == mNumPoints && + assert(pMaterialIds.size() == numPoints_ && "rayGeometry: Material IDs size mismatch"); - mMaterialIds.clear(); - mMaterialIds.reserve(mNumPoints); + materialIds_.clear(); + materialIds_.reserve(numPoints_); for (const auto id : pMaterialIds) { - mMaterialIds.push_back(static_cast(id)); + materialIds_.push_back(static_cast(id)); } } rayPair> getBoundingBox() const { - return {mMinCoords, mMaxCoords}; + return {minCoords_, maxCoords_}; } rayTriple getPoint(const unsigned int primID) const { - assert(primID < mNumPoints && "rayGeometry: Prim ID out of bounds"); - auto const &pnt = mPointBuffer[primID]; + assert(primID < numPoints_ && "rayGeometry: Prim ID out of bounds"); + auto const &pnt = pPointBuffer_[primID]; return {(NumericType)pnt.xx, (NumericType)pnt.yy, (NumericType)pnt.zz}; } - std::vector getNeighborIndicies(const unsigned int idx) const { - assert(idx < mNumPoints && "rayGeometry: Index out of bounds"); - return mPointNeighborhood[idx]; + std::vector const & + getNeighborIndicies(const unsigned int idx) const { + assert(idx < numPoints_ && "rayGeometry: Index out of bounds"); + return pointNeighborhood_[idx]; } - size_t getNumPoints() const { return mNumPoints; } + size_t getNumPoints() const { return numPoints_; } - NumericType getDiscRadius() const { return mDiscRadii; } + NumericType getDiscRadius() const { return discRadii_; } - RTCGeometry &getRTCGeometry() override final { return mRTCGeometry; } + RTCGeometry const &getRTCGeometry() const { return pRtcGeometry_; } - rayTriple - getPrimNormal(const unsigned int primID) override final { - assert(primID < mNumPoints && "rayGeometry: Prim ID out of bounds"); - auto const &normal = mNormalVecBuffer[primID]; + rayTriple getPrimNormal(const unsigned int primID) const { + assert(primID < numPoints_ && "rayGeometry: Prim ID out of bounds"); + auto const &normal = pNormalVecBuffer_[primID]; return {(NumericType)normal.xx, (NumericType)normal.yy, (NumericType)normal.zz}; } - rayQuadruple &getPrimRef(unsigned int primID) { - assert(primID < mNumPoints && "rayGeometry: Prim ID out of bounds"); - return *reinterpret_cast *>( - &mPointBuffer[primID]); + rayQuadruple &getPrimRef(unsigned int primID) { + assert(primID < numPoints_ && "rayGeometry: Prim ID out of bounds"); + return *reinterpret_cast *>( + &pPointBuffer_[primID]); } - rayTriple &getNormalRef(unsigned int primID) { - assert(primID < mNumPoints && "rayGeometry: Prim ID out of bounds"); - return *reinterpret_cast *>( - &mNormalVecBuffer[primID]); + rayTriple &getNormalRef(unsigned int primID) { + assert(primID < numPoints_ && "rayGeometry: Prim ID out of bounds"); + return *reinterpret_cast *>( + &pNormalVecBuffer_[primID]); } - std::vector &getMaterialIds() { return mMaterialIds; } - - int getMaterialId(const unsigned int primID) const override final { - assert(primID < mNumPoints && "rayGeometry Prim ID out of bounds"); - return mMaterialIds[primID]; + int getMaterialId(const unsigned int primID) const { + assert(primID < numPoints_ && "rayGeometry Prim ID out of bounds"); + return materialIds_[primID]; } - bool checkGeometryEmpty() { - if (mPointBuffer == nullptr || mNormalVecBuffer == nullptr || - mRTCGeometry == nullptr) { + bool checkGeometryEmpty() const { + if (pPointBuffer_ == nullptr || pNormalVecBuffer_ == nullptr || + pRtcGeometry_ == nullptr) { return true; } return false; @@ -174,14 +166,14 @@ class rayGeometry : public rayMetaGeometry { // Attention: // This function must not be called when the RTCGeometry reference count is // > 1. Doing so leads to leaked memory buffers - if (mPointBuffer == nullptr || mNormalVecBuffer == nullptr || - mRTCGeometry == nullptr) { + if (pPointBuffer_ == nullptr || pNormalVecBuffer_ == nullptr || + pRtcGeometry_ == nullptr) { return; } else { - rtcReleaseGeometry(mRTCGeometry); - mPointBuffer = nullptr; - mNormalVecBuffer = nullptr; - mRTCGeometry = nullptr; + rtcReleaseGeometry(pRtcGeometry_); + pPointBuffer_ = nullptr; + pNormalVecBuffer_ = nullptr; + pRtcGeometry_ = nullptr; } } @@ -189,16 +181,16 @@ class rayGeometry : public rayMetaGeometry { template void initPointNeighborhood(std::vector> &points) { - mPointNeighborhood.clear(); - mPointNeighborhood.resize(mNumPoints, std::vector{}); + pointNeighborhood_.clear(); + pointNeighborhood_.resize(numPoints_, std::vector{}); if constexpr (D == 3) { std::vector side1; std::vector side2; // create copy of bounding box - rayTriple min = mMinCoords; - rayTriple max = mMaxCoords; + rayTriple min = minCoords_; + rayTriple max = maxCoords_; std::vector dirs; for (int i = 0; i < 3; ++i) { @@ -212,7 +204,7 @@ class rayGeometry : public rayMetaGeometry { NumericType pivot = (max[dirs[dirIdx]] + min[dirs[dirIdx]]) / 2; // divide point data - for (unsigned int idx = 0; idx < mNumPoints; ++idx) { + for (unsigned int idx = 0; idx < numPoints_; ++idx) { if (points[idx][dirs[dirIdx]] <= pivot) { side1.push_back(idx); } else { @@ -222,11 +214,11 @@ class rayGeometry : public rayMetaGeometry { createNeighborhood(points, side1, side2, min, max, dirIdx, dirs, pivot); } else { // TODO: 2D divide and conquer algorithm - for (unsigned int idx1 = 0; idx1 < mNumPoints; ++idx1) { - for (unsigned int idx2 = idx1 + 1; idx2 < mNumPoints; ++idx2) { - if (checkDistance(points[idx1], points[idx2], 2 * mDiscRadii)) { - mPointNeighborhood[idx1].push_back(idx2); - mPointNeighborhood[idx2].push_back(idx1); + for (unsigned int idx1 = 0; idx1 < numPoints_; ++idx1) { + for (unsigned int idx2 = idx1 + 1; idx2 < numPoints_; ++idx2) { + if (checkDistance(points[idx1], points[idx2], 2 * discRadii_)) { + pointNeighborhood_[idx1].push_back(idx2); + pointNeighborhood_[idx2].push_back(idx1); } } } @@ -260,8 +252,8 @@ class rayGeometry : public rayMetaGeometry { auto const &pi1 = sides[idx1]; auto const &pi2 = sides[idx2]; assert(pi1 != pi2 && "Assumption"); - mPointNeighborhood[pi1].push_back(pi2); - mPointNeighborhood[pi2].push_back(pi1); + pointNeighborhood_[pi1].push_back(pi2); + pointNeighborhood_[pi2].push_back(pi1); } } return; @@ -288,7 +280,7 @@ class rayGeometry : public rayMetaGeometry { } else { s1r2set.push_back(side1[idx]); } - if (point[dirs[dirIdx]] + 2 * mDiscRadii <= pivot) { + if (point[dirs[dirIdx]] + 2 * discRadii_ <= pivot) { continue; } side1Cand.push_back(side1[idx]); @@ -301,7 +293,7 @@ class rayGeometry : public rayMetaGeometry { } else { s2r2set.push_back(side2[idx]); } - if (point[dirs[dirIdx]] - 2 * mDiscRadii >= pivot) { + if (point[dirs[dirIdx]] - 2 * discRadii_ >= pivot) { continue; } side2Cand.push_back(side2[idx]); @@ -315,11 +307,11 @@ class rayGeometry : public rayMetaGeometry { const auto &point2 = points[side2Cand[ci2]]; assert(std::abs(point1[dirs[dirIdx]] - point2[dirs[dirIdx]]) <= - (4 * mDiscRadii) && + (4 * discRadii_) && "Correctness Assertion"); - if (checkDistance(point1, point2, 2 * mDiscRadii)) { - mPointNeighborhood[side1Cand[ci1]].push_back(side2Cand[ci2]); - mPointNeighborhood[side2Cand[ci2]].push_back(side1Cand[ci1]); + if (checkDistance(point1, point2, 2 * discRadii_)) { + pointNeighborhood_[side1Cand[ci1]].push_back(side2Cand[ci2]); + pointNeighborhood_[side2Cand[ci2]].push_back(side1Cand[ci1]); } } } @@ -341,9 +333,9 @@ class rayGeometry : public rayMetaGeometry { } template - bool checkDistance(const std::array &p1, - const std::array &p2, - const NumericType &dist) { + static bool checkDistance(const std::array &p1, + const std::array &p2, + const NumericType &dist) { for (int i = 0; i < D; ++i) { if (std::abs(p1[i] - p2[i]) >= dist) return false; @@ -354,6 +346,7 @@ class rayGeometry : public rayMetaGeometry { return false; } +private: // "RTC_GEOMETRY_TYPE_POINT: // The vertex buffer stores each control vertex in the form of a single // precision position and radius stored in (x, y, z, r) order in memory @@ -363,7 +356,7 @@ class rayGeometry : public rayMetaGeometry { struct point_4f_t { float xx, yy, zz, radius; }; - point_4f_t *mPointBuffer = nullptr; + point_4f_t *pPointBuffer_ = nullptr; // "RTC_GEOMETRY_TYPE_POINT: // [...] the normal buffer stores a single precision normal per control @@ -372,19 +365,14 @@ class rayGeometry : public rayMetaGeometry { struct normal_vec_3f_t { float xx, yy, zz; }; - normal_vec_3f_t *mNormalVecBuffer = nullptr; - - RTCGeometry mRTCGeometry = nullptr; - - size_t mNumPoints; - NumericType mDiscRadii; - constexpr static NumericType nummax = std::numeric_limits::max(); - constexpr static NumericType nummin = - std::numeric_limits::lowest(); - rayTriple mMinCoords{nummax, nummax, nummax}; - rayTriple mMaxCoords{nummin, nummin, nummin}; - pointNeighborhoodType mPointNeighborhood; - std::vector mMaterialIds; -}; + normal_vec_3f_t *pNormalVecBuffer_ = nullptr; + + RTCGeometry pRtcGeometry_ = nullptr; -#endif // RAY_GEOMETRY_HPP \ No newline at end of file + size_t numPoints_; + NumericType discRadii_; + rayTriple minCoords_; + rayTriple maxCoords_; + pointNeighborhoodType pointNeighborhood_; + std::vector materialIds_; +}; diff --git a/include/viennaray/rayHitCounter.hpp b/include/viennaray/rayHitCounter.hpp index 9bf7a61..30dad25 100644 --- a/include/viennaray/rayHitCounter.hpp +++ b/include/viennaray/rayHitCounter.hpp @@ -1,44 +1,44 @@ -#ifndef RAY_HITCOUNTER_HPP -#define RAY_HITCOUNTER_HPP +#pragma once #include template class rayHitCounter { public: - rayHitCounter() : mTotalCnts(0) {} + rayHitCounter() : totalCounts_(0) {} // elements initialized to 0. rayHitCounter(size_t size) - : mCnts(size, 0), mTotalCnts(0), mDiskAreas(size, 0), mS1s(size, 0), - mS2s(size, 0) {} + : counts_(size, 0), totalCounts_(0), diskAreas_(size, 0), S1s_(size, 0), + S2s_(size, 0) {} // copy construct the vector members - rayHitCounter(rayHitCounter const &pA) - : mCnts(pA.mCnts), mTotalCnts(pA.mTotalCnts), mDiskAreas(pA.mDiskAreas), - mS1s(pA.mS1s), mS2s(pA.mS2s) {} + rayHitCounter(rayHitCounter const &other) + : counts_(other.counts_), totalCounts_(other.totalCounts_), + diskAreas_(other.diskAreas_), S1s_(other.S1s_), S2s_(other.S2s_) {} // move the vector members - rayHitCounter(rayHitCounter const &&pA) - : mCnts(std::move(pA.mCnts)), mTotalCnts(std::move(pA.mTotalCnts)), - mDiskAreas(std::move(pA.mDiskAreas)), mS1s(std::move(pA.mS1s)), - mS2s(std::move(pA.mS2s)) {} + rayHitCounter(rayHitCounter const &&other) + : counts_(std::move(other.counts_)), + totalCounts_(std::move(other.totalCounts_)), + diskAreas_(std::move(other.diskAreas_)), S1s_(std::move(other.S1s_)), + S2s_(std::move(other.S2s_)) {} // A copy constructor which can accumulate values from two instances // Precondition: the size of the accumulators are equal - rayHitCounter(rayHitCounter const &pA1, - rayHitCounter const &pA2) - : rayHitCounter(pA1) { - for (size_t idx = 0; idx < mCnts.size(); ++idx) { - mCnts[idx] += pA2.mCnts[idx]; - mS1s[idx] += pA2.mS1s[idx]; - mS2s[idx] += pA2.mS2s[idx]; + rayHitCounter(rayHitCounter const &A1, + rayHitCounter const &A2) + : rayHitCounter(A1) { + for (size_t idx = 0; idx < counts_.size(); ++idx) { + counts_[idx] += A2.counts_[idx]; + S1s_[idx] += A2.S1s_[idx]; + S2s_[idx] += A2.S2s_[idx]; } - mTotalCnts = pA1.mTotalCnts + pA2.mTotalCnts; - for (size_t idx = 0; idx < pA1.mDiskAreas.size(); ++idx) { - mDiskAreas[idx] = pA1.mDiskAreas[idx] > pA2.mDiskAreas[idx] - ? pA1.mDiskAreas[idx] - : pA2.mDiskAreas[idx]; + totalCounts_ = A1.totalCounts_ + A2.totalCounts_; + for (size_t idx = 0; idx < A1.diskAreas_.size(); ++idx) { + diskAreas_[idx] = A1.diskAreas_[idx] > A2.diskAreas_[idx] + ? A1.diskAreas_[idx] + : A2.diskAreas_[idx]; } } @@ -46,15 +46,15 @@ template class rayHitCounter { operator=(rayHitCounter const &pOther) { if (this != &pOther) { // copy from pOther to this - mCnts.clear(); - mCnts = pOther.mCnts; - mTotalCnts = pOther.mTotalCnts; - mDiskAreas.clear(); - mDiskAreas = pOther.mDiskAreas; - mS1s.clear(); - mS1s = pOther.mS1s; - mS2s.clear(); - mS2s = pOther.mS2s; + counts_.clear(); + counts_ = pOther.counts_; + totalCounts_ = pOther.totalCounts_; + diskAreas_.clear(); + diskAreas_ = pOther.diskAreas_; + S1s_.clear(); + S1s_ = pOther.S1s_; + S2s_.clear(); + S2s_ = pOther.S2s_; } return *this; } @@ -63,77 +63,77 @@ template class rayHitCounter { operator=(rayHitCounter const &&pOther) { if (this != &pOther) { // move from pOther to this - mCnts.clear(); - mCnts = std::move(pOther.mCnts); - mTotalCnts = pOther.mTotalCnts; - mDiskAreas.clear(); - mDiskAreas = std::move(pOther.mDiskAreas); - mS1s.clear(); - mS1s = std::move(pOther.mS1s); - mS2s.clear(); - mS2s = std::move(pOther.mS2s); + counts_.clear(); + counts_ = std::move(pOther.counts_); + totalCounts_ = pOther.totalCounts_; + diskAreas_.clear(); + diskAreas_ = std::move(pOther.diskAreas_); + S1s_.clear(); + S1s_ = std::move(pOther.S1s_); + S2s_.clear(); + S2s_ = std::move(pOther.S2s_); } return *this; } void use(unsigned int primID, NumericType value) { - mCnts[primID] += 1; - mTotalCnts += 1; - mS1s[primID] += value; - mS2s[primID] += value * value; + counts_[primID] += 1; + totalCounts_ += 1; + S1s_[primID] += value; + S2s_[primID] += value * value; } void merge(rayHitCounter const &pOther, const bool calcFlux) { if (calcFlux) { - for (size_t idx = 0; idx < mCnts.size(); ++idx) { - mCnts[idx] += pOther.mCnts[idx]; - mS1s[idx] += pOther.mS1s[idx]; - mS2s[idx] += pOther.mS2s[idx]; + for (size_t idx = 0; idx < counts_.size(); ++idx) { + counts_[idx] += pOther.counts_[idx]; + S1s_[idx] += pOther.S1s_[idx]; + S2s_[idx] += pOther.S2s_[idx]; } } - mTotalCnts += pOther.mTotalCnts; - for (size_t idx = 0; idx < mDiskAreas.size(); ++idx) { - mDiskAreas[idx] = mDiskAreas[idx] > pOther.mDiskAreas[idx] - ? mDiskAreas[idx] - : pOther.mDiskAreas[idx]; + totalCounts_ += pOther.totalCounts_; + for (size_t idx = 0; idx < diskAreas_.size(); ++idx) { + diskAreas_[idx] = diskAreas_[idx] > pOther.diskAreas_[idx] + ? diskAreas_[idx] + : pOther.diskAreas_[idx]; } } void resize(const size_t numPoints, const bool calcFlux) { - mDiskAreas.resize(numPoints); - mTotalCnts = 0; + diskAreas_.resize(numPoints); + totalCounts_ = 0; if (calcFlux) { - mCnts.resize(numPoints); - mS1s.resize(numPoints); - mS2s.resize(numPoints); + counts_.resize(numPoints); + S1s_.resize(numPoints); + S2s_.resize(numPoints); } } void clear() { - mDiskAreas.clear(); - mCnts.clear(); - mS1s.clear(); - mS2s.clear(); + diskAreas_.clear(); + counts_.clear(); + S1s_.clear(); + S2s_.clear(); } - std::vector getValues() const { return mS1s; } + std::vector getValues() const { return S1s_; } - std::vector getCounts() const { return mCnts; } + std::vector getCounts() const { return counts_; } - size_t getTotalCounts() const { return mTotalCnts; } + size_t getTotalCounts() const { return totalCounts_; } - const std::vector &getDiskAreas() const { return mDiskAreas; } + const std::vector &getDiskAreas() const { return diskAreas_; } std::vector getRelativeError() { auto result = std::vector( - mS1s.size(), + S1s_.size(), std::numeric_limits::max()); // size, initial values - if (mTotalCnts == 0) { + if (totalCounts_ == 0) { return result; } for (size_t idx = 0; idx < result.size(); ++idx) { - auto s1square = mS1s[idx] * mS1s[idx]; + auto s1square = S1s_[idx] * S1s_[idx]; if (s1square == 0) { continue; } @@ -141,25 +141,23 @@ template class rayHitCounter { // sqrt(N) For details and an exact formula see the book Exploring Monte // Carlo Methods by Dunn and Shultis page 83 and 84. result[idx] = - (NumericType)(std::sqrt(mS2s[idx] / s1square - 1.0 / mTotalCnts)); + (NumericType)(std::sqrt(S2s_[idx] / s1square - 1.0 / totalCounts_)); } return result; } void setDiskAreas(std::vector &pDiskAreas) { - mDiskAreas = pDiskAreas; + diskAreas_ = pDiskAreas; } private: - std::vector mCnts; - size_t mTotalCnts; - std::vector mDiskAreas; + std::vector counts_; + size_t totalCounts_; + std::vector diskAreas_; // S1 denotes the sum of sample values - std::vector mS1s; + std::vector S1s_; // S2 denotes the sum of squared sample values // these are need to compute the relative error - std::vector mS2s; + std::vector S2s_; }; - -#endif // RAY_HITCOUNTER_HPP \ No newline at end of file diff --git a/include/viennaray/rayMessage.hpp b/include/viennaray/rayMessage.hpp index 5100a32..7f2f22e 100644 --- a/include/viennaray/rayMessage.hpp +++ b/include/viennaray/rayMessage.hpp @@ -1,14 +1,14 @@ -#ifndef RAY_MESSAGE_HPP -#define RAY_MESSAGE_HPP +#pragma once #include +#include /// Singleton class for thread-safe logging. class rayMessage { - std::string message; + std::string message_; - bool error = false; - const unsigned tabWidth = 2; + bool error_ = false; + constexpr static unsigned tabWidth_ = 4; rayMessage() {} @@ -25,8 +25,8 @@ class rayMessage { rayMessage &addWarning(std::string s) { #pragma omp critical { - message += - "\n[ViennaRay]" + std::string(tabWidth, ' ') + "WARNING: " + s + "\n"; + message_ += "\n[ViennaRay]" + std::string(tabWidth_, ' ') + + "WARNING: " + s + "\n"; } return *this; } @@ -34,11 +34,11 @@ class rayMessage { rayMessage &addError(std::string s, bool shouldAbort = true) { #pragma omp critical { - message += - "\n[ViennaRay]" + std::string(tabWidth, ' ') + "ERROR: " + s + "\n"; + message_ += + "\n[ViennaRay]" + std::string(tabWidth_, ' ') + "ERROR: " + s + "\n"; } // always abort once error message should be printed - error = true; + error_ = true; // abort now if asked if (shouldAbort) print(); @@ -48,18 +48,16 @@ class rayMessage { rayMessage &addDebug(std::string s) { #pragma omp critical { - message += - "[ViennaRay]" + std::string(tabWidth, ' ') + "DEBUG: " + s + "\n"; + message_ += + "[ViennaRay]" + std::string(tabWidth_, ' ') + "DEBUG: " + s + "\n"; } return *this; } void print(std::ostream &out = std::cout) { - out << message; - message.clear(); - if (error) + out << message_; + message_.clear(); + if (error_) abort(); } }; - -#endif // RAY_MESSAGE_HPP \ No newline at end of file diff --git a/include/viennaray/rayMetaGeometry.hpp b/include/viennaray/rayMetaGeometry.hpp deleted file mode 100644 index ee07499..0000000 --- a/include/viennaray/rayMetaGeometry.hpp +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef RAY_METAGEOMETRY_HPP -#define RAY_METAGEOMETRY_HPP - -#if VIENNARAY_EMBREE_VERSION < 4 -#include -#else -#include -#endif -#include - -template class rayMetaGeometry { -public: - virtual ~rayMetaGeometry() {} - virtual RTCGeometry &getRTCGeometry() = 0; - virtual rayTriple getPrimNormal(const unsigned int primID) = 0; - virtual rayTriple getNewOrigin(RTCRay &ray) { - assert(rayInternal::IsNormalized( - rayTriple{ray.dir_x, ray.dir_y, ray.dir_z}) && - "MetaGeometry: direction not normalized"); - auto xx = ray.org_x + ray.dir_x * ray.tfar; - auto yy = ray.org_y + ray.dir_y * ray.tfar; - auto zz = ray.org_z + ray.dir_z * ray.tfar; - return {xx, yy, zz}; - } - virtual int getMaterialId(const unsigned int primID) const { return 0; }; -}; - -#endif // RAY_METAGEOMETRY_HPP \ No newline at end of file diff --git a/include/viennaray/rayParticle.hpp b/include/viennaray/rayParticle.hpp index 0b4dea2..4a00617 100644 --- a/include/viennaray/rayParticle.hpp +++ b/include/viennaray/rayParticle.hpp @@ -1,7 +1,5 @@ -#ifndef RAY_PARTICLE_HPP -#define RAY_PARTICLE_HPP +#pragma once -#include #include #include #include @@ -139,5 +137,3 @@ class rayTestParticle void logData(rayDataLog &log) override final {} }; - -#endif // RAY_PARTICLE_HPP \ No newline at end of file diff --git a/include/viennaray/rayPreCompileMacros.hpp b/include/viennaray/rayPreCompileMacros.hpp deleted file mode 100644 index f5f58ba..0000000 --- a/include/viennaray/rayPreCompileMacros.hpp +++ /dev/null @@ -1,9 +0,0 @@ -#ifndef RAY_PRE_COMPILE_MACROS -#define RAY_PRE_COMPILE_MACROS - -#if defined(__x86_64__) || defined(_M_X64) -#define ARCH_X86 -#include -#endif - -#endif \ No newline at end of file diff --git a/include/viennaray/rayRNG.hpp b/include/viennaray/rayRNG.hpp index 9b50f6e..723f2ca 100644 --- a/include/viennaray/rayRNG.hpp +++ b/include/viennaray/rayRNG.hpp @@ -1,11 +1,10 @@ -#ifndef RAY_RNG_HPP -#define RAY_RNG_HPP +#pragma once #include #include /// Use mersenne twister 19937 as random number generator. -typedef std::mt19937_64 rayRNG; +using rayRNG = std::mt19937_64; namespace rayInternal { @@ -25,5 +24,3 @@ static unsigned int tea(unsigned int val0, unsigned int val1) { return v0; } } // namespace rayInternal - -#endif // RAY_RNG_HPP \ No newline at end of file diff --git a/include/viennaray/rayReflection.hpp b/include/viennaray/rayReflection.hpp index 4d6f412..defb7bb 100644 --- a/include/viennaray/rayReflection.hpp +++ b/include/viennaray/rayReflection.hpp @@ -1,12 +1,10 @@ -#ifndef RAY_REFLECTION_HPP -#define RAY_REFLECTION_HPP +#pragma once -#include #include #include // Specular reflection -template +template static rayTriple rayReflectionSpecular(const rayTriple &rayDir, const rayTriple &geomNormal) { @@ -34,7 +32,7 @@ rayReflectionDiffuse(const rayTriple &geomNormal, rayRNG &RNG) { "rayReflectionDiffuse: Surface normal is not normalized"); auto randomDirection = - rayInternal::PickRandomPointOnUnitSphere(RNG); + rayInternal::pickRandomPointOnUnitSphere(RNG); randomDirection[0] += geomNormal[0]; randomDirection[1] += geomNormal[1]; if constexpr (D == 3) @@ -63,9 +61,9 @@ static rayTriple rayReflectionConedCosine( do { // sample phi uniformly in [0, 2pi] - NumericType phi = uniDist(RNG) * 2 * PI; + NumericType phi = uniDist(RNG) * 2 * M_PI; // theta on sphere - assert(maxConeAngle >= 0. && maxConeAngle <= PI / 2. && + assert(maxConeAngle >= 0. && maxConeAngle <= M_PI / 2. && "Cone angle not allowed"); NumericType cosTheta = std::cos(maxConeAngle); // sample z uniformly on [cos(theta),1] @@ -137,7 +135,7 @@ static rayTriple rayReflectionConedCosineOld( sqrt_1m_u = std::sqrt(1. - u); angle = avgReflAngle * sqrt_1m_u; } while (uniDist(RNG) * angle * u > - std::cos(rayInternal::PI / 2. * sqrt_1m_u) * std::sin(angle)); + std::cos(M_PI_2 * sqrt_1m_u) * std::sin(angle)); // Random Azimuthal Rotation NumericType costheta = @@ -222,8 +220,7 @@ rayReflectionConedCosineOld2(const rayTriple &rayDir, const NumericType incAngle = std::acos(std::max(std::min(cosTheta, 1.), 0.)); - NumericType avgReflAngle = - std::max(rayInternal::PI / 2. - incAngle, minAvgConeAngle); + NumericType avgReflAngle = std::max(M_PI_2 - incAngle, minAvgConeAngle); std::uniform_real_distribution uniDist; NumericType u, sqrt_1m_u; @@ -234,7 +231,7 @@ rayReflectionConedCosineOld2(const rayTriple &rayDir, sqrt_1m_u = std::sqrt(1. - u); angle = avgReflAngle * sqrt_1m_u; } while (uniDist(RNG) * angle * u > - std::cos(rayInternal::PI / 2. * sqrt_1m_u) * std::sin(angle)); + std::cos(M_PI_2 * sqrt_1m_u) * std::sin(angle)); cosTheta = std::cos(angle); @@ -289,5 +286,3 @@ rayReflectionConedCosineOld2(const rayTriple &rayDir, return randomDir; } } // namespace rayInternal - -#endif // RAY_REFLECTION_HPP \ No newline at end of file diff --git a/include/viennaray/raySource.hpp b/include/viennaray/raySource.hpp index a3aa6b2..6864a9c 100644 --- a/include/viennaray/raySource.hpp +++ b/include/viennaray/raySource.hpp @@ -1,21 +1,19 @@ -#ifndef RAY_SOURCE_HPP -#define RAY_SOURCE_HPP +#pragma once -#if VIENNARAY_EMBREE_VERSION < 4 -#include -#else -#include -#endif -#include #include #include -template class raySource { +template class raySource { +protected: + raySource() = default; + ~raySource() = default; + public: - virtual ~raySource() {} - virtual void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) {} - virtual size_t getNumPoints() const { return 0; } - virtual void printIndexCounter(){}; -}; + Derived &derived() { return static_cast(*this); } + const Derived &derived() const { return static_cast(*this); } -#endif // RAY_SOURCE_HPP \ No newline at end of file + void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) const { + derived().fillRay(ray, idx, RngState); + } + size_t getNumPoints() const { return derived().getNumPoints(); } +}; diff --git a/include/viennaray/raySourceGrid.hpp b/include/viennaray/raySourceGrid.hpp index b8b4b5e..93f6bf1 100644 --- a/include/viennaray/raySourceGrid.hpp +++ b/include/viennaray/raySourceGrid.hpp @@ -1,26 +1,20 @@ -#ifndef RAY_SOURCEGRID_HPP -#define RAY_SOURCEGRID_HPP +#pragma once -#include #include template -class raySourceGrid : public raySource { +class raySourceGrid : public raySource> { public: raySourceGrid(std::vector> &sourceGrid, - NumericType passedCosinePower, - const std::array &passedTraceSettings) - : mSourceGrid(sourceGrid), mNumPoints(sourceGrid.size()), - cosinePower(passedCosinePower), rayDir(passedTraceSettings[0]), - firstDir(passedTraceSettings[1]), secondDir(passedTraceSettings[2]), - minMax(passedTraceSettings[3]), posNeg(passedTraceSettings[4]), - ee(((NumericType)2) / (passedCosinePower + 1)), - indexCounter(sourceGrid.size(), 0) {} + NumericType cosinePower, + const std::array &traceSettings) + : sourceGrid_(sourceGrid), numPoints_(sourceGrid.size()), + rayDir_(traceSettings[0]), firstDir_(traceSettings[1]), + secondDir_(traceSettings[2]), minMax_(traceSettings[3]), + posNeg_(traceSettings[4]), ee_(((NumericType)2) / (cosinePower + 1)) {} - void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) override final { - auto index = idx % mNumPoints; - indexCounter[index]++; - auto origin = mSourceGrid[idx % mNumPoints]; + void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) const { + auto origin = sourceGrid_[idx % numPoints_]; auto direction = getDirection(RngState); #ifdef ARCH_X86 @@ -42,29 +36,23 @@ class raySourceGrid : public raySource { #endif } - size_t getNumPoints() const override final { return mNumPoints; } - - void printIndexCounter() override final { - for (const auto &idx : indexCounter) { - std::cout << idx << std::endl; - } - } + size_t getNumPoints() const { return numPoints_; } private: - rayTriple getDirection(rayRNG &RngState) { + rayTriple getDirection(rayRNG &RngState) const { rayTriple direction{0., 0., 0.}; std::uniform_real_distribution uniDist; auto r1 = uniDist(RngState); auto r2 = uniDist(RngState); - NumericType tt = pow(r2, ee); - direction[rayDir] = posNeg * sqrtf(tt); - direction[firstDir] = cosf(two_pi * r1) * sqrtf(1 - tt); + NumericType tt = pow(r2, ee_); + direction[rayDir_] = posNeg_ * sqrtf(tt); + direction[firstDir_] = cosf(M_PI * 2.f * r1) * sqrtf(1.f - tt); if constexpr (D == 2) { - direction[secondDir] = 0; + direction[secondDir_] = 0; } else { - direction[secondDir] = sinf(two_pi * r1) * sqrtf(1 - tt); + direction[secondDir_] = sinf(M_PI * 2.f * r1) * sqrtf(1.f - tt); } rayInternal::Normalize(direction); @@ -72,17 +60,13 @@ class raySourceGrid : public raySource { return direction; } - const std::vector> &mSourceGrid; - const size_t mNumPoints; - const NumericType cosinePower; - const int rayDir; - const int firstDir; - const int secondDir; - const int minMax; - const NumericType posNeg; - const NumericType ee; - std::vector indexCounter; - constexpr static NumericType two_pi = rayInternal::PI * 2; +private: + const std::vector> &sourceGrid_; + const size_t numPoints_; + const int rayDir_; + const int firstDir_; + const int secondDir_; + const int minMax_; + const NumericType posNeg_; + const NumericType ee_; }; - -#endif // RT_RAYSOURCE_HPP \ No newline at end of file diff --git a/include/viennaray/raySourceRandom.hpp b/include/viennaray/raySourceRandom.hpp index 0edfb76..62c924c 100644 --- a/include/viennaray/raySourceRandom.hpp +++ b/include/viennaray/raySourceRandom.hpp @@ -1,29 +1,28 @@ -#ifndef RAY_SOURCERANDOM_HPP -#define RAY_SOURCERANDOM_HPP +#pragma once #include template -class raySourceRandom : public raySource { - typedef rayPair> boundingBoxType; +class raySourceRandom : public raySource> { + using boundingBoxType = rayPair>; public: raySourceRandom( - boundingBoxType pBoundingBox, NumericType pCosinePower, - std::array &pTraceSettings, const size_t pNumPoints, - const bool pCustomDirection, - const std::array, 3> &pOrthonormalBasis) - : bdBox(pBoundingBox), rayDir(pTraceSettings[0]), - firstDir(pTraceSettings[1]), secondDir(pTraceSettings[2]), - minMax(pTraceSettings[3]), posNeg(pTraceSettings[4]), - ee(((NumericType)2) / (pCosinePower + 1)), mNumPoints(pNumPoints), - customDirection(pCustomDirection), orthonormalBasis(pOrthonormalBasis) { + boundingBoxType boundingBox, NumericType cosinePower, + std::array &pTraceSettings, const size_t numPoints_, + const bool customDirection, + const std::array, 3> &orthonormalBasis) + : bdBox_(boundingBox), rayDir_(pTraceSettings[0]), + firstDir_(pTraceSettings[1]), secondDir_(pTraceSettings[2]), + minMax_(pTraceSettings[3]), posNeg_(pTraceSettings[4]), + ee_(((NumericType)2) / (cosinePower + 1)), numPoints_(numPoints_), + customDirection_(customDirection), orthonormalBasis_(orthonormalBasis) { } - void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) override final { + void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) const { auto origin = getOrigin(RngState); rayTriple direction; - if (customDirection) { + if (customDirection_) { direction = getCustomDirection(RngState); } else { direction = getDirection(RngState); @@ -48,50 +47,50 @@ class raySourceRandom : public raySource { #endif } - size_t getNumPoints() const override final { return mNumPoints; } + size_t getNumPoints() const { return numPoints_; } private: - rayTriple getOrigin(rayRNG &RngState) { + rayTriple getOrigin(rayRNG &RngState) const { rayTriple origin{0., 0., 0.}; std::uniform_real_distribution uniDist; auto r1 = uniDist(RngState); - origin[rayDir] = bdBox[minMax][rayDir]; - origin[firstDir] = - bdBox[0][firstDir] + (bdBox[1][firstDir] - bdBox[0][firstDir]) * r1; + origin[rayDir_] = bdBox_[minMax_][rayDir_]; + origin[firstDir_] = bdBox_[0][firstDir_] + + (bdBox_[1][firstDir_] - bdBox_[0][firstDir_]) * r1; if constexpr (D == 2) { - origin[secondDir] = 0.; + origin[secondDir_] = 0.; } else { auto r2 = uniDist(RngState); - origin[secondDir] = bdBox[0][secondDir] + - (bdBox[1][secondDir] - bdBox[0][secondDir]) * r2; + origin[secondDir_] = bdBox_[0][secondDir_] + + (bdBox_[1][secondDir_] - bdBox_[0][secondDir_]) * r2; } return origin; } - rayTriple getDirection(rayRNG &RngState) { + rayTriple getDirection(rayRNG &RngState) const { rayTriple direction{0., 0., 0.}; std::uniform_real_distribution uniDist; auto r1 = uniDist(RngState); auto r2 = uniDist(RngState); - const NumericType tt = pow(r2, ee); - direction[rayDir] = posNeg * sqrtf(tt); - direction[firstDir] = cosf(two_pi * r1) * sqrtf(1 - tt); + const NumericType tt = pow(r2, ee_); + direction[rayDir_] = posNeg_ * sqrtf(tt); + direction[firstDir_] = cosf(M_PI * 2.f * r1) * sqrtf(1 - tt); if constexpr (D == 2) { - direction[secondDir] = 0; + direction[secondDir_] = 0; rayInternal::Normalize(direction); } else { - direction[secondDir] = sinf(two_pi * r1) * sqrtf(1 - tt); + direction[secondDir_] = sinf(M_PI * 2.f * r1) * sqrtf(1 - tt); } return direction; } - rayTriple getCustomDirection(rayRNG &RngState) { + rayTriple getCustomDirection(rayRNG &RngState) const { rayTriple direction; std::uniform_real_distribution uniDist; @@ -100,42 +99,40 @@ class raySourceRandom : public raySource { auto r1 = uniDist(RngState); auto r2 = uniDist(RngState); - const NumericType tt = pow(r2, ee); + const NumericType tt = pow(r2, ee_); rndDirection[0] = sqrtf(tt); - rndDirection[1] = cosf(two_pi * r1) * sqrtf(1 - tt); - rndDirection[2] = sinf(two_pi * r1) * sqrtf(1 - tt); - - direction[0] = orthonormalBasis[0][0] * rndDirection[0] + - orthonormalBasis[1][0] * rndDirection[1] + - orthonormalBasis[2][0] * rndDirection[2]; - direction[1] = orthonormalBasis[0][1] * rndDirection[0] + - orthonormalBasis[1][1] * rndDirection[1] + - orthonormalBasis[2][1] * rndDirection[2]; - direction[2] = orthonormalBasis[0][2] * rndDirection[0] + - orthonormalBasis[1][2] * rndDirection[1] + - orthonormalBasis[2][2] * rndDirection[2]; - } while ((posNeg < 0. && direction[rayDir] > 0.) || - (posNeg > 0. && direction[rayDir] < 0.)); + rndDirection[1] = cosf(M_PI * 2.f * r1) * sqrtf(1 - tt); + rndDirection[2] = sinf(M_PI * 2.f * r1) * sqrtf(1 - tt); + + direction[0] = orthonormalBasis_[0][0] * rndDirection[0] + + orthonormalBasis_[1][0] * rndDirection[1] + + orthonormalBasis_[2][0] * rndDirection[2]; + direction[1] = orthonormalBasis_[0][1] * rndDirection[0] + + orthonormalBasis_[1][1] * rndDirection[1] + + orthonormalBasis_[2][1] * rndDirection[2]; + direction[2] = orthonormalBasis_[0][2] * rndDirection[0] + + orthonormalBasis_[1][2] * rndDirection[1] + + orthonormalBasis_[2][2] * rndDirection[2]; + } while ((posNeg_ < 0. && direction[rayDir_] > 0.) || + (posNeg_ > 0. && direction[rayDir_] < 0.)); if constexpr (D == 2) { - direction[secondDir] = 0; + direction[secondDir_] = 0; rayInternal::Normalize(direction); } return direction; } - const boundingBoxType bdBox; - const int rayDir; - const int firstDir; - const int secondDir; - const int minMax; - const NumericType posNeg; - const NumericType ee; - const size_t mNumPoints; - constexpr static double two_pi = rayInternal::PI * 2; - const bool customDirection = false; - const std::array, 3> &orthonormalBasis; +private: + const boundingBoxType bdBox_; + const int rayDir_; + const int firstDir_; + const int secondDir_; + const int minMax_; + const NumericType posNeg_; + const NumericType ee_; + const size_t numPoints_; + const bool customDirection_ = false; + const std::array, 3> &orthonormalBasis_; }; - -#endif // RAY_SOURCERANDOM_HPP \ No newline at end of file diff --git a/include/viennaray/rayTestAsserts.hpp b/include/viennaray/rayTestAsserts.hpp index 3af782f..c8bb857 100644 --- a/include/viennaray/rayTestAsserts.hpp +++ b/include/viennaray/rayTestAsserts.hpp @@ -1,7 +1,6 @@ #pragma once -#include + #include -#include #ifdef _MSC_VER #define __PRETTY_FUNCTION__ __FUNCSIG__ diff --git a/include/viennaray/rayTrace.hpp b/include/viennaray/rayTrace.hpp index afbb709..9fd4ea7 100644 --- a/include/viennaray/rayTrace.hpp +++ b/include/viennaray/rayTrace.hpp @@ -1,64 +1,65 @@ -#ifndef RAY_TRACE_HPP -#define RAY_TRACE_HPP +#pragma once #include #include #include #include #include -#include #include #include +#include template class rayTrace { public: - rayTrace() : mDevice(rtcNewDevice("hugepages=1")) {} + rayTrace() : device_(rtcNewDevice("hugepages=1")) {} + + rayTrace(const rayTrace &) = delete; + rayTrace &operator=(const rayTrace &) = delete; + rayTrace(rayTrace &&) = delete; + rayTrace &operator=(rayTrace &&) = delete; ~rayTrace() { - mGeometry.releaseGeometry(); - rtcReleaseDevice(mDevice); + geometry_.releaseGeometry(); + rtcReleaseDevice(device_); } /// Run the ray tracer void apply() { checkSettings(); initMemoryFlags(); - auto boundingBox = mGeometry.getBoundingBox(); + auto boundingBox = geometry_.getBoundingBox(); rayInternal::adjustBoundingBox( - boundingBox, mSourceDirection, mDiskRadius); - auto traceSettings = rayInternal::getTraceSettings(mSourceDirection); + boundingBox, sourceDirection_, diskRadius_); + auto traceSettings = rayInternal::getTraceSettings(sourceDirection_); auto boundary = rayBoundary( - mDevice, boundingBox, mBoundaryConditions, traceSettings); + device_, boundingBox, boundaryConditions_, traceSettings); - std::array, 3> orthoBasis; - if (usePrimaryDirection) - orthoBasis = rayInternal::getOrthonormalBasis(primaryDirection); + std::array, 3> orthonormalBasis; + if (usePrimaryDirection_) + orthonormalBasis = rayInternal::getOrthonormalBasis(primaryDirection_); auto raySource = raySourceRandom( - boundingBox, mParticle->getSourceDistributionPower(), traceSettings, - mGeometry.getNumPoints(), usePrimaryDirection, orthoBasis); + boundingBox, pParticle_->getSourceDistributionPower(), traceSettings, + geometry_.getNumPoints(), usePrimaryDirection_, orthonormalBasis); - auto localDataLabels = mParticle->getLocalDataLabels(); + auto localDataLabels = pParticle_->getLocalDataLabels(); if (!localDataLabels.empty()) { - mLocalData.setNumberOfVectorData(localDataLabels.size()); - auto numPoints = mGeometry.getNumPoints(); - auto localDataLabels = mParticle->getLocalDataLabels(); + localData_.setNumberOfVectorData(localDataLabels.size()); + auto numPoints = geometry_.getNumPoints(); + auto localDataLabels = pParticle_->getLocalDataLabels(); for (int i = 0; i < localDataLabels.size(); ++i) { - mLocalData.setVectorData(i, numPoints, 0., localDataLabels[i]); + localData_.setVectorData(i, numPoints, 0., localDataLabels[i]); } } - auto tracer = rayTraceKernel( - mDevice, mGeometry, boundary, raySource, mParticle, mDataLog, - mNumberOfRaysPerPoint, mNumberOfRaysFixed, mUseRandomSeeds, mCalcFlux, - mRunNumber++); - - tracer.setTracingData(&mLocalData, mGlobalData); - tracer.setHitCounter(&mHitCounter); - tracer.setRayTraceInfo(&mRTInfo); + rayTraceKernel tracer(device_, geometry_, boundary, raySource, pParticle_, + dataLog_, numberOfRaysPerPoint_, numberOfRaysFixed_, + useRandomSeeds_, calcFlux_, runNumber_++, hitCounter_, + RTInfo_); + tracer.setTracingData(&localData_, pGlobalData_); tracer.apply(); - if (mCheckError) + if (checkError_) checkRelativeError(); boundary.releaseGeometry(); @@ -67,12 +68,12 @@ template class rayTrace { /// Set the particle type used for ray tracing /// The particle is a user defined object that has to interface the /// rayParticle class. - template - void setParticleType(std::unique_ptr &p) { - static_assert(std::is_base_of, - ParticleType>::value && - "Particle object does not interface correct class"); - mParticle = p->clone(); + template , + ParticleType>, + bool> = true> + void setParticleType(std::unique_ptr &particle) { + pParticle_ = particle->clone(); } /// Set the ray tracing geometry @@ -85,9 +86,9 @@ template class rayTrace { static_assert((D != 3 || Dim != 2) && "Setting 2D geometry in 3D trace object"); - mGridDelta = gridDelta; - mDiskRadius = gridDelta * rayInternal::DiskFactor; - mGeometry.initGeometry(mDevice, points, normals, mDiskRadius); + gridDelta_ = gridDelta; + diskRadius_ = gridDelta * rayInternal::DiskFactor; + geometry_.initGeometry(device_, points, normals, diskRadius_); } /// Set the ray tracing geometry @@ -99,45 +100,45 @@ template class rayTrace { static_assert((D != 3 || Dim != 2) && "Setting 2D geometry in 3D trace object"); - mGridDelta = gridDelta; - mDiskRadius = diskRadii; - mGeometry.initGeometry(mDevice, points, normals, mDiskRadius); + gridDelta_ = gridDelta; + diskRadius_ = diskRadii; + geometry_.initGeometry(device_, points, normals, diskRadius_); } /// Set material ID's for each geometry point. /// If not set, all material ID's are default 0. - template void setMaterialIds(std::vector &pMaterialIds) { - mGeometry.setMaterialIds(pMaterialIds); + template void setMaterialIds(std::vector &materialIds) { + geometry_.setMaterialIds(materialIds); } /// Set the boundary conditions. /// There has to be a boundary condition defined for each space dimension, /// however the boundary condition in direction of the tracing direction is /// ignored. - void setBoundaryConditions(rayBoundaryCondition pBoundaryConditions[D]) { + void setBoundaryConditions(rayBoundaryCondition boundaryConditions[D]) { for (size_t i = 0; i < D; ++i) { - mBoundaryConditions[i] = pBoundaryConditions[i]; + boundaryConditions_[i] = boundaryConditions[i]; } } /// Set the number of rays per geometry point. /// The total number of rays, that are traced, is the set number set here /// times the number of points in the geometry. - void setNumberOfRaysPerPoint(const size_t pNum) { - mNumberOfRaysPerPoint = pNum; - mNumberOfRaysFixed = 0; + void setNumberOfRaysPerPoint(const size_t numRaysPerPoint) { + numberOfRaysPerPoint_ = numRaysPerPoint; + numberOfRaysFixed_ = 0; } /// Set the number of total rays traced to a fixed amount, /// independent of the geometry - void setNumberOfRaysFixed(const size_t pNum) { - mNumberOfRaysFixed = pNum; - mNumberOfRaysPerPoint = 0; + void setNumberOfRaysFixed(const size_t numRaysFixed) { + numberOfRaysFixed_ = numRaysFixed; + numberOfRaysPerPoint_ = 0; } /// Set the source direction, where the rays should be traced from. - void setSourceDirection(const rayTraceDirection pDirection) { - mSourceDirection = pDirection; + void setSourceDirection(const rayTraceDirection direction) { + sourceDirection_ = direction; } /// Set the primary direction of the source distribution. This can be used to @@ -145,20 +146,20 @@ template class rayTrace { /// not change the position of the source plane. Therefore, one has the be /// careful that the resulting distribution does not lie completely above the /// source plane. - void setPrimaryDirection(const rayTriple pPrimaryDirection) { - primaryDirection = pPrimaryDirection; - usePrimaryDirection = true; + void setPrimaryDirection(const rayTriple primaryDirection) { + primaryDirection_ = primaryDirection; + usePrimaryDirection_ = true; } /// Set whether random seeds for the internal random number generators /// should be used. - void setUseRandomSeeds(const bool useRand) { mUseRandomSeeds = useRand; } + void setUseRandomSeeds(const bool useRand) { useRandomSeeds_ = useRand; } /// Set whether the flux and hit counts should be recorded. If not needed, /// this should be turned off to increase performance. If set to false, the /// functions getTotalFlux(), getNormalizedFlux(), getHitCounts() and /// getRelativeError() can not be used. - void setCalculateFlux(const bool calcFlux) { mCalcFlux = calcFlux; } + void setCalculateFlux(const bool calcFlux) { calcFlux_ = calcFlux; } /// Set whether to check the relative error after a tracing. If the relative /// error at a surface point is larger than 0.05 a warning is printed. The @@ -166,18 +167,18 @@ template class rayTrace { /// Code, Version 5, Vol. I Overview and Theory, LA - UR - 03 - 1987, Los /// Alamos Nat.Lab., Los Alamos, NM" void setCheckRelativeError(const bool checkError) { - mCheckError = checkError; + checkError_ = checkError; } /// Returns the total flux on each disk. std::vector getTotalFlux() const { - return mHitCounter.getValues(); + return hitCounter_.getValues(); } /// Returns the normalized flux on each disk. std::vector getNormalizedFlux(rayNormalizationType normalization, bool averageNeighborhood = false) { - auto flux = mHitCounter.getValues(); + auto flux = hitCounter_.getValues(); normalizeFlux(flux, normalization); if (averageNeighborhood) { smoothFlux(flux); @@ -190,11 +191,11 @@ template class rayTrace { /// value. void normalizeFlux(std::vector &flux, rayNormalizationType norm = rayNormalizationType::SOURCE) { - assert(flux.size() == mGeometry.getNumPoints() && + assert(flux.size() == geometry_.getNumPoints() && "Unequal number of points in normalizeFlux"); - auto diskArea = mHitCounter.getDiskAreas(); - const auto totalDiskArea = mDiskRadius * mDiskRadius * rayInternal::PI; + auto diskArea = hitCounter_.getDiskAreas(); + const auto totalDiskArea = diskRadius_ * diskRadius_ * M_PI; switch (norm) { case rayNormalizationType::MAX: { @@ -208,9 +209,9 @@ template class rayTrace { case rayNormalizationType::SOURCE: { NumericType sourceArea = getSourceArea(); - auto numTotalRays = mNumberOfRaysFixed == 0 - ? flux.size() * mNumberOfRaysPerPoint - : mNumberOfRaysFixed; + auto numTotalRays = numberOfRaysFixed_ == 0 + ? flux.size() * numberOfRaysPerPoint_ + : numberOfRaysFixed_; NumericType normFactor = sourceArea / numTotalRays; #pragma omp parallel for for (int idx = 0; idx < flux.size(); ++idx) { @@ -227,12 +228,12 @@ template class rayTrace { /// Helper function to smooth the recorded flux by averaging over the /// neighborhood in a post-processing step. void smoothFlux(std::vector &flux) { - assert(flux.size() == mGeometry.getNumPoints() && + assert(flux.size() == geometry_.getNumPoints() && "Unequal number of points in smoothFlux"); auto oldFlux = flux; #pragma omp parallel for - for (int idx = 0; idx < mGeometry.getNumPoints(); idx++) { - auto neighborhood = mGeometry.getNeighborIndicies(idx); + for (int idx = 0; idx < geometry_.getNumPoints(); idx++) { + auto neighborhood = geometry_.getNeighborIndicies(idx); for (auto const &nbi : neighborhood) { flux[idx] += oldFlux[nbi]; } @@ -241,45 +242,47 @@ template class rayTrace { } /// Returns the total number of hits for each geometry point. - std::vector getHitCounts() const { return mHitCounter.getCounts(); } + std::vector getHitCounts() const { return hitCounter_.getCounts(); } /// Returns the relative error of the flux for each geometry point std::vector getRelativeError() { - return mHitCounter.getRelativeError(); + return hitCounter_.getRelativeError(); } /// Returns the disk area for each geometry point - std::vector getDiskAreas() { return mHitCounter.getDiskAreas(); } + std::vector getDiskAreas() { return hitCounter_.getDiskAreas(); } - rayTracingData &getLocalData() { return mLocalData; } + rayTracingData &getLocalData() { return localData_; } - rayTracingData *getGlobalData() { return mGlobalData; } + rayTracingData *getGlobalData() { return pGlobalData_; } - void setGlobalData(rayTracingData &data) { mGlobalData = &data; } + void setGlobalData(rayTracingData &data) { + pGlobalData_ = &data; + } - rayTraceInfo getRayTraceInfo() { return mRTInfo; } + rayTraceInfo getRayTraceInfo() { return RTInfo_; } - rayDataLog &getDataLog() { return mDataLog; } + rayDataLog &getDataLog() { return dataLog_; } private: NumericType getSourceArea() { - const auto boundingBox = mGeometry.getBoundingBox(); + const auto boundingBox = geometry_.getBoundingBox(); NumericType sourceArea = 0; - if (mSourceDirection == rayTraceDirection::NEG_X || - mSourceDirection == rayTraceDirection::POS_X) { + if (sourceDirection_ == rayTraceDirection::NEG_X || + sourceDirection_ == rayTraceDirection::POS_X) { sourceArea = (boundingBox[1][1] - boundingBox[0][1]); if constexpr (D == 3) { sourceArea *= (boundingBox[1][2] - boundingBox[0][2]); } - } else if (mSourceDirection == rayTraceDirection::NEG_Y || - mSourceDirection == rayTraceDirection::POS_Y) { + } else if (sourceDirection_ == rayTraceDirection::NEG_Y || + sourceDirection_ == rayTraceDirection::POS_Y) { sourceArea = (boundingBox[1][0] - boundingBox[0][0]); if constexpr (D == 3) { sourceArea *= (boundingBox[1][2] - boundingBox[0][2]); } - } else if (mSourceDirection == rayTraceDirection::NEG_Z || - mSourceDirection == rayTraceDirection::POS_Z) { + } else if (sourceDirection_ == rayTraceDirection::NEG_Z || + sourceDirection_ == rayTraceDirection::POS_Z) { assert(D == 3 && "Error in flux normalization"); sourceArea = (boundingBox[1][0] - boundingBox[0][0]); sourceArea *= (boundingBox[1][1] - boundingBox[0][1]); @@ -312,7 +315,7 @@ template class rayTrace { } } if (!allPassed) { - mRTInfo.warning = true; + RTInfo_.warning = true; rayMessage::getInstance() .addWarning( "Large relative error detected. Consider using more rays.") @@ -321,24 +324,24 @@ template class rayTrace { } void checkSettings() { - if (mParticle == nullptr) { - mRTInfo.error = true; + if (pParticle_ == nullptr) { + RTInfo_.error = true; rayMessage::getInstance().addError( "No particle was specified in rayTrace. Aborting."); } - if (mGeometry.checkGeometryEmpty()) { - mRTInfo.error = true; + if (geometry_.checkGeometryEmpty()) { + RTInfo_.error = true; rayMessage::getInstance().addError( "No geometry was passed to rayTrace. Aborting."); } - if ((D == 2 && mSourceDirection == rayTraceDirection::POS_Z) || - (D == 2 && mSourceDirection == rayTraceDirection::NEG_Z)) { - mRTInfo.error = true; + if ((D == 2 && sourceDirection_ == rayTraceDirection::POS_Z) || + (D == 2 && sourceDirection_ == rayTraceDirection::NEG_Z)) { + RTInfo_.error = true; rayMessage::getInstance().addError( "Invalid source direction in 2D geometry. Aborting."); } - if (mDiskRadius > mGridDelta) { - mRTInfo.warning = true; + if (diskRadius_ > gridDelta_) { + RTInfo_.warning = true; rayMessage::getInstance() .addWarning("Disk radius should be smaller than grid delta. Hit " "count normalization not correct.") @@ -356,26 +359,24 @@ template class rayTrace { } private: - RTCDevice mDevice; - rayGeometry mGeometry; - std::unique_ptr> mParticle = nullptr; - size_t mNumberOfRaysPerPoint = 1000; - size_t mNumberOfRaysFixed = 0; - NumericType mDiskRadius = 0; - NumericType mGridDelta = 0; - rayBoundaryCondition mBoundaryConditions[D] = {}; - rayTraceDirection mSourceDirection = rayTraceDirection::POS_Z; - rayTriple primaryDirection = {0.}; - bool usePrimaryDirection = false; - bool mUseRandomSeeds = false; - size_t mRunNumber = 0; - bool mCalcFlux = true; - bool mCheckError = true; - rayHitCounter mHitCounter; - rayTracingData mLocalData; - rayTracingData *mGlobalData = nullptr; - rayTraceInfo mRTInfo; - rayDataLog mDataLog; + RTCDevice device_; + rayGeometry geometry_; + std::unique_ptr> pParticle_ = nullptr; + size_t numberOfRaysPerPoint_ = 1000; + size_t numberOfRaysFixed_ = 0; + NumericType diskRadius_ = 0; + NumericType gridDelta_ = 0; + rayBoundaryCondition boundaryConditions_[D] = {}; + rayTraceDirection sourceDirection_ = rayTraceDirection::POS_Z; + rayTriple primaryDirection_ = {0.}; + bool usePrimaryDirection_ = false; + bool useRandomSeeds_ = false; + size_t runNumber_ = 0; + bool calcFlux_ = true; + bool checkError_ = true; + rayHitCounter hitCounter_; + rayTracingData localData_; + rayTracingData *pGlobalData_ = nullptr; + rayTraceInfo RTInfo_; + rayDataLog dataLog_; }; - -#endif // RAY_TRACE_HPP \ No newline at end of file diff --git a/include/viennaray/rayTraceDirection.hpp b/include/viennaray/rayTraceDirection.hpp deleted file mode 100644 index 9e8463f..0000000 --- a/include/viennaray/rayTraceDirection.hpp +++ /dev/null @@ -1,13 +0,0 @@ -#ifndef RAY_TRACEDIRECTION_HPP -#define RAY_TRACEDIRECTION_HPP - -enum struct rayTraceDirection : unsigned { - POS_X = 0, - NEG_X = 1, - POS_Y = 2, - NEG_Y = 3, - POS_Z = 4, - NEG_Z = 5 -}; - -#endif \ No newline at end of file diff --git a/include/viennaray/rayTraceKernel.hpp b/include/viennaray/rayTraceKernel.hpp index 76ae492..ce99a7a 100644 --- a/include/viennaray/rayTraceKernel.hpp +++ b/include/viennaray/rayTraceKernel.hpp @@ -1,5 +1,4 @@ -#ifndef RAY_TRACEKERNEL_HPP -#define RAY_TRACEKERNEL_HPP +#pragma once #include #include @@ -10,34 +9,32 @@ #include #include -#define PRINT_PROGRESS false -#define PRINT_RESULT false - -template class rayTraceKernel { - +template +class rayTraceKernel { public: - rayTraceKernel(RTCDevice &pDevice, rayGeometry &pRTCGeometry, - rayBoundary &pRTCBoundary, - raySource &pSource, - std::unique_ptr> &pParticle, - rayDataLog &pDataLog, - const size_t pNumOfRayPerPoint, const size_t pNumOfRayFixed, - const bool pUseRandomSeed, const bool pCalcFlux, - const size_t pRunNumber) - : mDevice(pDevice), mGeometry(pRTCGeometry), mBoundary(pRTCBoundary), - mSource(pSource), mParticle(pParticle->clone()), - mNumRays(pNumOfRayFixed == 0 - ? pSource.getNumPoints() * pNumOfRayPerPoint - : pNumOfRayFixed), - mUseRandomSeeds(pUseRandomSeed), mRunNumber(pRunNumber), - mCalcFlux(pCalcFlux), dataLog(pDataLog) { - assert(rtcGetDeviceProperty(mDevice, RTC_DEVICE_PROPERTY_VERSION) >= + rayTraceKernel(RTCDevice &device, rayGeometry &rtcGeometry, + rayBoundary &rtcBoundary, + raySource &source, + std::unique_ptr> &particle, + rayDataLog &dataLog, const size_t numRaysPerPoint, + const size_t numRaysFixed, const bool useRandomSeed, + const bool calcFlux, const size_t runNumber, + rayHitCounter &hitCounter, + rayTraceInfo &traceInfo) + : device_(device), geometry_(rtcGeometry), boundary_(rtcBoundary), + source_(source), pParticle_(particle->clone()), + numRays_(numRaysFixed == 0 ? source.getNumPoints() * numRaysPerPoint + : numRaysFixed), + useRandomSeeds_(useRandomSeed), runNumber_(runNumber), + calcFlux_(calcFlux), hitCounter_(hitCounter), traceInfo_(traceInfo), + dataLog_(dataLog) { + assert(rtcGetDeviceProperty(device_, RTC_DEVICE_PROPERTY_VERSION) >= 30601 && "Error: The minimum version of Embree is 3.6.1"); } void apply() { - auto rtcScene = rtcNewScene(mDevice); + auto rtcScene = rtcNewScene(device_); // RTC scene flags rtcSetSceneFlags(rtcScene, RTC_SCENE_FLAG_NONE); @@ -47,44 +44,42 @@ template class rayTraceKernel { // RTC_BUILD_QUALITY_MEDIUM. auto bbquality = RTC_BUILD_QUALITY_HIGH; rtcSetSceneBuildQuality(rtcScene, bbquality); - auto rtcGeometry = mGeometry.getRTCGeometry(); - auto rtcBoundary = mBoundary.getRTCGeometry(); + auto rtcGeometry = geometry_.getRTCGeometry(); + auto rtcBoundary = boundary_.getRTCGeometry(); auto boundaryID = rtcAttachGeometry(rtcScene, rtcBoundary); auto geometryID = rtcAttachGeometry(rtcScene, rtcGeometry); - assert(rtcGetDeviceError(mDevice) == RTC_ERROR_NONE && + assert(rtcGetDeviceError(device_) == RTC_ERROR_NONE && "Embree device error"); size_t geohitc = 0; size_t nongeohitc = 0; size_t totaltraces = 0; - const bool calcFlux = mCalcFlux; // thread local data storage const int numThreads = omp_get_max_threads(); std::vector> threadLocalData(numThreads); for (auto &data : threadLocalData) { - data = *localData; + data = *pLocalData_; } // thread local data log std::vector> threadLocalDataLog(numThreads); for (auto &data : threadLocalDataLog) { - data = dataLog; - assert(data.data.size() == dataLog.data.size()); + data = dataLog_; + assert(data.data.size() == dataLog_.data.size()); for (auto &d : data.data) { std::fill(d.begin(), d.end(), 0.); } } // hit counters - assert(hitCounter != nullptr && "Hit counter is nullptr"); - hitCounter->clear(); - hitCounter->resize(mGeometry.getNumPoints(), calcFlux); std::vector> threadLocalHitCounter(numThreads); - if (calcFlux) { + if (calcFlux_) { + hitCounter_.clear(); + hitCounter_.resize(geometry_.getNumPoints(), calcFlux_); for (auto &hitC : threadLocalHitCounter) { - hitC = *hitCounter; + hitC = hitCounter_; } } @@ -99,14 +94,14 @@ template class rayTraceKernel { RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; const int threadID = omp_get_thread_num(); - unsigned int seed = mRunNumber; - if (mUseRandomSeeds) { + unsigned int seed = runNumber_; + if (useRandomSeeds_) { std::random_device rd; seed = static_cast(rd()); } // thread-local particle object - auto particle = mParticle->clone(); + auto particle = pParticle_->clone(); auto &myLocalData = threadLocalData[threadID]; auto &myHitCounter = threadLocalHitCounter[threadID]; @@ -123,7 +118,7 @@ template class rayTraceKernel { [[maybe_unused]] size_t progressCount = 0; #pragma omp for schedule(dynamic) - for (long long idx = 0; idx < mNumRays; ++idx) { + for (long long idx = 0; idx < numRays_; ++idx) { // particle specific RNG seed auto particleSeed = rayInternal::tea<3>(idx, seed); rayRNG RngState(particleSeed); @@ -132,20 +127,21 @@ template class rayTraceKernel { particle->logData(myDataLog); NumericType rayWeight = initialRayWeight; - mSource.fillRay(rayHit.ray, idx, RngState); // fills also tnear + source_.fillRay(rayHit.ray, idx, RngState); // fills also tnear #ifdef VIENNARAY_USE_RAY_MASKING rayHit.ray.mask = -1; #endif - if constexpr (PRINT_PROGRESS) { - printProgress(progressCount); - } +#ifdef VIENNARAY_PRINT_PROGRESS + printProgress(progressCount); +#endif bool reflect = false; bool hitFromBack = false; do { - rayHit.ray.tfar = std::numeric_limits::max(); + rayHit.ray.tfar = + std::numeric_limits::max(); rayHit.hit.instID[0] = RTC_INVALID_GEOMETRY_ID; rayHit.hit.geomID = RTC_INVALID_GEOMETRY_ID; // rayHit.ray.tnear = 1e-4f; // tnear is also set in the particle @@ -169,20 +165,23 @@ template class rayTraceKernel { /* -------- Boundary hit -------- */ if (rayHit.hit.geomID == boundaryID) { - mBoundary.processHit(rayHit, reflect); + boundary_.processHit(rayHit, reflect); continue; } // Calculate point of impact const auto &ray = rayHit.ray; - const rtcNumericType xx = ray.org_x + ray.dir_x * ray.tfar; - const rtcNumericType yy = ray.org_y + ray.dir_y * ray.tfar; - const rtcNumericType zz = ray.org_z + ray.dir_z * ray.tfar; + const rayInternal::rtcNumericType xx = + ray.org_x + ray.dir_x * ray.tfar; + const rayInternal::rtcNumericType yy = + ray.org_y + ray.dir_y * ray.tfar; + const rayInternal::rtcNumericType zz = + ray.org_z + ray.dir_z * ray.tfar; /* -------- Hit from back -------- */ const auto rayDir = rayTriple{ray.dir_x, ray.dir_y, ray.dir_z}; - const auto geomNormal = mGeometry.getPrimNormal(rayHit.hit.primID); + const auto geomNormal = geometry_.getPrimNormal(rayHit.hit.primID); if (rayInternal::DotProduct(rayDir, geomNormal) > 0) { // If the dot product of the ray direction and the surface normal is // greater than zero, then we hit the back face of the disk. @@ -213,13 +212,13 @@ template class rayTraceKernel { std::vector hitDiskIds(1, rayHit.hit.primID); #ifdef VIENNARAY_USE_WDIST - std::vector + std::vector impactDistances; // distances between point of impact and disk // origins of hit disks { // distance on first disk hit - const auto &disk = mGeometry.getPrimRef(rayHit.hit.primID); - const auto &diskOrigin = - *reinterpret_cast const *>(&disk); + const auto &disk = geometry_.getPrimRef(rayHit.hit.primID); + const auto &diskOrigin = *reinterpret_cast< + rayTriple const *>(&disk); impactDistances.push_back( rayInternal::Distance({xx, yy, zz}, diskOrigin) + 1e-6f); // add eps to avoid division by 0 @@ -227,8 +226,8 @@ template class rayTraceKernel { #endif // check for additional intersections for (const auto &id : - mGeometry.getNeighborIndicies(rayHit.hit.primID)) { - rtcNumericType distance; + geometry_.getNeighborIndicies(rayHit.hit.primID)) { + rayInternal::rtcNumericType distance; if (checkLocalIntersection(ray, id, distance)) { hitDiskIds.push_back(id); #ifdef VIENNARAY_USE_WDIST @@ -239,14 +238,14 @@ template class rayTraceKernel { const size_t numDisksHit = hitDiskIds.size(); #ifdef VIENNARAY_USE_WDIST - rtcNumericType invDistanceWeightSum = 0; + rayInternal::rtcNumericType invDistanceWeightSum = 0; for (const auto &d : impactDistances) invDistanceWeightSum += 1 / d; #endif // for each disk hit for (size_t diskId = 0; diskId < numDisksHit; ++diskId) { - const auto matID = mGeometry.getMaterialId(hitDiskIds[diskId]); - const auto normal = mGeometry.getPrimNormal(hitDiskIds[diskId]); + const auto matID = geometry_.getMaterialId(hitDiskIds[diskId]); + const auto normal = geometry_.getPrimNormal(hitDiskIds[diskId]); #ifdef VIENNARAY_USE_WDIST auto distRayWeight = rayWeight / impactDistances[diskId] / invDistanceWeightSum * numDisksHit; @@ -255,16 +254,17 @@ template class rayTraceKernel { #endif particle->surfaceCollision(distRayWeight, rayDir, normal, hitDiskIds[diskId], matID, myLocalData, - globalData, RngState); + pGlobalData_, RngState); } // get sticking probability and reflection direction - const auto stickingnDirection = particle->surfaceReflection( + const auto stickingDirection = particle->surfaceReflection( rayWeight, rayDir, geomNormal, rayHit.hit.primID, - mGeometry.getMaterialId(rayHit.hit.primID), globalData, RngState); - auto valueToDrop = rayWeight * stickingnDirection.first; + geometry_.getMaterialId(rayHit.hit.primID), pGlobalData_, + RngState); + auto valueToDrop = rayWeight * stickingDirection.first; - if (calcFlux) { + if (calcFlux_) { for (size_t diskId = 0; diskId < numDisksHit; ++diskId) { #ifdef VIENNARAY_USE_WDIST auto distRayWeight = valueToDrop / impactDistances[diskId] / @@ -290,19 +290,22 @@ template class rayTraceKernel { #ifdef ARCH_X86 reinterpret_cast<__m128 &>(rayHit.ray) = _mm_set_ps(1e-4f, zz, yy, xx); - reinterpret_cast<__m128 &>(rayHit.ray.dir_x) = - _mm_set_ps(0.0f, (rtcNumericType)stickingnDirection.second[2], - (rtcNumericType)stickingnDirection.second[1], - (rtcNumericType)stickingnDirection.second[0]); + reinterpret_cast<__m128 &>(rayHit.ray.dir_x) = _mm_set_ps( + 0.0f, (rayInternal::rtcNumericType)stickingDirection.second[2], + (rayInternal::rtcNumericType)stickingDirection.second[1], + (rayInternal::rtcNumericType)stickingDirection.second[0]); #else rayHit.ray.org_x = xx; rayHit.ray.org_y = yy; rayHit.ray.org_z = zz; rayHit.ray.tnear = 1e-4f; - rayHit.ray.dir_x = (rtcNumericType)stickingnDirection.second[0]; - rayHit.ray.dir_y = (rtcNumericType)stickingnDirection.second[1]; - rayHit.ray.dir_z = (rtcNumericType)stickingnDirection.second[2]; + rayHit.ray.dir_x = + (rayInternal::rtcNumericType)stickingDirection.second[0]; + rayHit.ray.dir_y = + (rayInternal::rtcNumericType)stickingDirection.second[1]; + rayHit.ray.dir_z = + (rayInternal::rtcNumericType)stickingDirection.second[2]; rayHit.ray.time = 0.0f; #endif } while (reflect); @@ -316,19 +319,19 @@ template class rayTraceKernel { // merge hit counters and data logs for (int i = 0; i < numThreads; ++i) { - hitCounter->merge(threadLocalHitCounter[i], calcFlux); - dataLog.merge(threadLocalDataLog[i]); + hitCounter_.merge(threadLocalHitCounter[i], calcFlux_); + dataLog_.merge(threadLocalDataLog[i]); } // merge local data - if (!localData->getVectorData().empty()) { + if (!pLocalData_->getVectorData().empty()) { // merge vector data #pragma omp parallel for - for (int i = 0; i < localData->getVectorData().size(); ++i) { - switch (localData->getVectorMergeType(i)) { + for (int i = 0; i < pLocalData_->getVectorData().size(); ++i) { + switch (pLocalData_->getVectorMergeType(i)) { case rayTracingDataMergeEnum::SUM: { - for (size_t j = 0; j < localData->getVectorData(i).size(); ++j) { + for (size_t j = 0; j < pLocalData_->getVectorData(i).size(); ++j) { for (int k = 0; k < numThreads; ++k) { - localData->getVectorData(i)[j] += + pLocalData_->getVectorData(i)[j] += threadLocalData[k].getVectorData(i)[j]; } } @@ -336,9 +339,10 @@ template class rayTraceKernel { } case rayTracingDataMergeEnum::APPEND: { - localData->getVectorData(i).clear(); + pLocalData_->getVectorData(i).clear(); for (int k = 0; k < numThreads; ++k) { - localData->appendVectorData(i, threadLocalData[k].getVectorData(i)); + pLocalData_->appendVectorData(i, + threadLocalData[k].getVectorData(i)); } break; } @@ -353,22 +357,24 @@ template class rayTraceKernel { } } - if (!localData->getScalarData().empty()) { + if (!pLocalData_->getScalarData().empty()) { // merge scalar data - for (int i = 0; i < localData->getScalarData().size(); ++i) { - switch (localData->getScalarMergeType(i)) { + for (int i = 0; i < pLocalData_->getScalarData().size(); ++i) { + switch (pLocalData_->getScalarMergeType(i)) { case rayTracingDataMergeEnum::SUM: { for (int j = 0; j < numThreads; ++j) { - localData->getScalarData(i) += threadLocalData[j].getScalarData(i); + pLocalData_->getScalarData(i) += + threadLocalData[j].getScalarData(i); } break; } case rayTracingDataMergeEnum::AVERAGE: { for (int j = 0; j < numThreads; ++j) { - localData->getScalarData(i) += threadLocalData[j].getScalarData(i); + pLocalData_->getScalarData(i) += + threadLocalData[j].getScalarData(i); } - localData->getScalarData(i) /= (NumericType)numThreads; + pLocalData_->getScalarData(i) /= (NumericType)numThreads; break; } @@ -382,33 +388,25 @@ template class rayTraceKernel { } } - if (mTraceInfo != nullptr) { - mTraceInfo->numRays = mNumRays; - mTraceInfo->totalRaysTraced = totaltraces; - mTraceInfo->totalDiskHits = hitCounter->getTotalCounts(); - mTraceInfo->nonGeometryHits = nongeohitc; - mTraceInfo->geometryHits = geohitc; - mTraceInfo->time = (endTime - time) * 1e-3; - } + traceInfo_.numRays = numRays_; + traceInfo_.totalRaysTraced = totaltraces; + traceInfo_.totalDiskHits = hitCounter_.getTotalCounts(); + traceInfo_.nonGeometryHits = nongeohitc; + traceInfo_.geometryHits = geohitc; + traceInfo_.time = (endTime - time) * 1e-3; rtcReleaseScene(rtcScene); } - void setTracingData(rayTracingData *plocalData, - const rayTracingData *pglobalData) { - localData = plocalData; - globalData = pglobalData; + void setTracingData(rayTracingData *pLocalData, + const rayTracingData *pGlobalData) { + pLocalData_ = pLocalData; + pGlobalData_ = pGlobalData; } - void setHitCounter(rayHitCounter *phitCounter) { - hitCounter = phitCounter; - } - - void setRayTraceInfo(rayTraceInfo *pTraceInfo) { mTraceInfo = pTraceInfo; } - private: - bool rejectionControl(NumericType &rayWeight, const NumericType &initWeight, - rayRNG &RNG) { + static bool rejectionControl(NumericType &rayWeight, + const NumericType &initWeight, rayRNG &RNG) { // Choosing a good value for the weight lower threshold is important NumericType lowerThreshold = 0.1 * initWeight; NumericType renewWeight = 0.3 * initWeight; @@ -434,20 +432,19 @@ template class rayTraceKernel { return true; } - std::vector computeDiskAreas() { + std::vector computeDiskAreas() const { constexpr double eps = 1e-3; - auto bdBox = mGeometry.getBoundingBox(); - const auto numOfPrimitives = mGeometry.getNumPoints(); - const auto boundaryDirs = mBoundary.getDirs(); + auto bdBox = geometry_.getBoundingBox(); + const auto numOfPrimitives = geometry_.getNumPoints(); + const auto boundaryDirs = boundary_.getDirs(); auto areas = std::vector(numOfPrimitives, 0); #pragma omp for for (long idx = 0; idx < numOfPrimitives; ++idx) { - auto const &disk = mGeometry.getPrimRef(idx); + auto const &disk = geometry_.getPrimRef(idx); if constexpr (D == 3) { - areas[idx] = - disk[3] * disk[3] * (NumericType)rayInternal::PI; // full disk area + areas[idx] = disk[3] * disk[3] * M_PI; // full disk area if (std::fabs(disk[boundaryDirs[0]] - bdBox[0][boundaryDirs[0]]) < eps || @@ -466,7 +463,7 @@ template class rayTraceKernel { } } else { areas[idx] = 2 * disk[3]; - auto normal = mGeometry.getNormalRef(idx); + auto normal = geometry_.getNormalRef(idx); // test min boundary if (std::abs(disk[boundaryDirs[0]] - bdBox[0][boundaryDirs[0]]) < @@ -502,7 +499,7 @@ template class rayTraceKernel { return areas; } - void printProgress(size_t &progressCount) { + [[maybe_unused]] void printProgress(size_t &progressCount) const { if (omp_get_thread_num() != 0) { return; } @@ -512,11 +509,11 @@ template class rayTraceKernel { constexpr auto emptySymbol = '-'; constexpr auto barEndSymbol = ']'; constexpr auto percentageStringFormatLength = 3; // 3 digits - if (progressCount % (int)std::ceil((rtcNumericType)mNumRays / + if (progressCount % (int)std::ceil((NumericType)numRays_ / omp_get_num_threads() / barLength) == 0) { auto fillLength = - (int)std::ceil(progressCount / ((rtcNumericType)mNumRays / + (int)std::ceil(progressCount / ((NumericType)numRays_ / omp_get_num_threads() / barLength)); auto percentageString = std::to_string((fillLength * 100) / barLength); percentageString = @@ -536,17 +533,21 @@ template class rayTraceKernel { progressCount += 1; } - bool checkLocalIntersection(RTCRay const &ray, const unsigned int primID, - rtcNumericType &impactDistance) { + bool + checkLocalIntersection(RTCRay const &ray, const unsigned int primID, + rayInternal::rtcNumericType &impactDistance) const { auto const &rayOrigin = - *reinterpret_cast const *>(&ray.org_x); + *reinterpret_cast const *>( + &ray.org_x); auto const &rayDirection = - *reinterpret_cast const *>(&ray.dir_x); + *reinterpret_cast const *>( + &ray.dir_x); - const auto &normal = mGeometry.getNormalRef(primID); - const auto &disk = mGeometry.getPrimRef(primID); + const auto &normal = geometry_.getNormalRef(primID); + const auto &disk = geometry_.getPrimRef(primID); const auto &diskOrigin = - *reinterpret_cast const *>(&disk); + *reinterpret_cast const *>( + &disk); auto prodOfDirections = rayInternal::DotProduct(normal, rayDirection); if (prodOfDirections > 0.f) { @@ -571,7 +572,7 @@ template class rayTraceKernel { } // copy ray direction - auto rayDirectionC = rayTriple{ + auto rayDirectionC = rayTriple{ rayDirection[0], rayDirection[1], rayDirection[2]}; rayInternal::Scale(tt, rayDirectionC); auto hitpoint = rayInternal::Sum(rayOrigin, rayDirectionC); @@ -584,20 +585,19 @@ template class rayTraceKernel { return false; } - RTCDevice &mDevice; - rayGeometry &mGeometry; - rayBoundary &mBoundary; - raySource &mSource; - std::unique_ptr> const mParticle = nullptr; - const long long mNumRays; - const bool mUseRandomSeeds; - const size_t mRunNumber; - const bool mCalcFlux; - rayTracingData *localData = nullptr; - const rayTracingData *globalData = nullptr; - rayHitCounter *hitCounter = nullptr; - rayTraceInfo *mTraceInfo = nullptr; - rayDataLog &dataLog; +private: + RTCDevice &device_; + rayGeometry &geometry_; + rayBoundary const &boundary_; + raySource const &source_; + std::unique_ptr> const pParticle_ = nullptr; + const long long numRays_; + const bool useRandomSeeds_; + const size_t runNumber_; + const bool calcFlux_; + rayTracingData *pLocalData_ = nullptr; + rayTracingData const *pGlobalData_ = nullptr; + rayHitCounter &hitCounter_; + rayTraceInfo &traceInfo_; + rayDataLog &dataLog_; }; - -#endif // RAY_TRACEKERNEL_HPP \ No newline at end of file diff --git a/include/viennaray/rayTracingData.hpp b/include/viennaray/rayTracingData.hpp index dc38ea4..e56e04a 100644 --- a/include/viennaray/rayTracingData.hpp +++ b/include/viennaray/rayTracingData.hpp @@ -1,9 +1,11 @@ -#ifndef RAY_TRACINGDATA_HPP -#define RAY_TRACINGDATA_HPP +#pragma once -#include +#include -enum struct rayTracingDataMergeEnum : unsigned { +#include +#include + +enum class rayTracingDataMergeEnum : unsigned { SUM = 0, APPEND = 1, AVERAGE = 2, @@ -15,121 +17,122 @@ template class rayTracingData { using vectorDataType = std::vector; using mergeType = std::vector; - std::vector scalarData; - std::vector vectorData; - std::vector scalarDataLabels; - std::vector vectorDataLabels; - mergeType scalarDataMerge; - mergeType vectorDataMerge; + std::vector scalarData_; + std::vector vectorData_; + std::vector scalarDataLabels_; + std::vector vectorDataLabels_; + mergeType scalarDataMerge_; + mergeType vectorDataMerge_; public: rayTracingData() {} rayTracingData(const rayTracingData &otherData) - : scalarData(otherData.scalarData), vectorData(otherData.vectorData), - scalarDataLabels(otherData.scalarDataLabels), - vectorDataLabels(otherData.vectorDataLabels), - scalarDataMerge(otherData.scalarDataMerge), - vectorDataMerge(otherData.vectorDataMerge) {} + : scalarData_(otherData.scalarData_), vectorData_(otherData.vectorData_), + scalarDataLabels_(otherData.scalarDataLabels_), + vectorDataLabels_(otherData.vectorDataLabels_), + scalarDataMerge_(otherData.scalarDataMerge_), + vectorDataMerge_(otherData.vectorDataMerge_) {} rayTracingData(rayTracingData &&otherData) - : scalarData(std::move(otherData.scalarData)), - vectorData(std::move(otherData.vectorData)), - scalarDataLabels(std::move(otherData.scalarDataLabels)), - vectorDataLabels(std::move(otherData.vectorDataLabels)), - scalarDataMerge(std::move(otherData.scalarDataMerge)), - vectorDataMerge(std::move(otherData.vectorDataMerge)) {} + : scalarData_(std::move(otherData.scalarData_)), + vectorData_(std::move(otherData.vectorData_)), + scalarDataLabels_(std::move(otherData.scalarDataLabels_)), + vectorDataLabels_(std::move(otherData.vectorDataLabels_)), + scalarDataMerge_(std::move(otherData.scalarDataMerge_)), + vectorDataMerge_(std::move(otherData.vectorDataMerge_)) {} rayTracingData &operator=(const rayTracingData &otherData) { - scalarData = otherData.scalarData; - vectorData = otherData.vectorData; - scalarDataLabels = otherData.scalarDataLabels; - vectorDataLabels = otherData.vectorDataLabels; - scalarDataMerge = otherData.scalarDataMerge; - vectorDataMerge = otherData.vectorDataMerge; + scalarData_ = otherData.scalarData_; + vectorData_ = otherData.vectorData_; + scalarDataLabels_ = otherData.scalarDataLabels_; + vectorDataLabels_ = otherData.vectorDataLabels_; + scalarDataMerge_ = otherData.scalarDataMerge_; + vectorDataMerge_ = otherData.vectorDataMerge_; return *this; } rayTracingData &operator=(rayTracingData &&otherData) { - scalarData = std::move(otherData.scalarData); - vectorData = std::move(otherData.vectorData); - scalarDataLabels = std::move(otherData.scalarDataLabels); - vectorDataLabels = std::move(otherData.vectorDataLabels); - scalarDataMerge = std::move(otherData.scalarDataMerge); - vectorDataMerge = std::move(otherData.vectorDataMerge); + scalarData_ = std::move(otherData.scalarData_); + vectorData_ = std::move(otherData.vectorData_); + scalarDataLabels_ = std::move(otherData.scalarDataLabels_); + vectorDataLabels_ = std::move(otherData.vectorDataLabels_); + scalarDataMerge_ = std::move(otherData.scalarDataMerge_); + vectorDataMerge_ = std::move(otherData.vectorDataMerge_); return *this; } void appendVectorData(int num, const std::vector &vec) { - vectorData[num].insert(vectorData[num].end(), vec.begin(), vec.end()); + vectorData_[num].insert(vectorData_[num].end(), vec.begin(), vec.end()); } void setNumberOfVectorData(int size) { - vectorData.clear(); - vectorData.resize(size); - vectorDataMerge.resize(size, rayTracingDataMergeEnum::SUM); - vectorDataLabels.resize(size, "vectorData"); + vectorData_.clear(); + vectorData_.resize(size); + vectorDataMerge_.resize(size, rayTracingDataMergeEnum::SUM); + vectorDataLabels_.resize(size, "vectorData"); } void setNumberOfScalarData(int size) { - scalarData.clear(); - scalarData.resize(size); - scalarDataMerge.resize(size, rayTracingDataMergeEnum::SUM); - scalarDataLabels.resize(size, "scalarData"); + scalarData_.clear(); + scalarData_.resize(size); + scalarDataMerge_.resize(size, rayTracingDataMergeEnum::SUM); + scalarDataLabels_.resize(size, "scalarData"); } void setScalarData(int num, NumericType value, std::string label = "scalarData") { - if (num >= vectorData.size()) + if (num >= vectorData_.size()) rayMessage::getInstance() .addError("Setting scalar data in rayTracingData out of range.") .print(); - scalarData[num] = value; - scalarDataLabels[num] = label; + scalarData_[num] = value; + scalarDataLabels_[num] = label; } void setVectorData(int num, std::vector &vector, std::string label = "vectorData") { - if (num >= vectorData.size()) + if (num >= vectorData_.size()) rayMessage::getInstance() .addError("Setting vector data in rayTracingData out of range.") .print(); - vectorData[num] = vector; - vectorDataLabels[num] = label; + vectorData_[num] = vector; + vectorDataLabels_[num] = label; } void setVectorData(int num, std::vector &&vector, std::string label = "vectorData") { - if (num >= vectorData.size()) + if (num >= vectorData_.size()) rayMessage::getInstance() .addError("Setting vector data in rayTracingData out of range.") .print(); - vectorData[num] = std::move(vector); - vectorDataLabels[num] = label; + vectorData_[num] = std::move(vector); + vectorDataLabels_[num] = label; } void setVectorData(int num, size_t size, NumericType value, std::string label = "vectorData") { - if (num >= vectorData.size()) + if (num >= vectorData_.size()) rayMessage::getInstance() .addError("Setting vector data in rayTracingData out of range.") .print(); - vectorData[num].resize(size, value); - vectorDataLabels[num] = label; + vectorData_[num].resize(size, value); + vectorDataLabels_[num] = label; } void setVectorData(int num, NumericType value, std::string label = "vectorData") { - if (num >= vectorData.size()) + if (num >= vectorData_.size()) rayMessage::getInstance() .addError("Setting vector data in rayTracingData out of range.") .print(); - vectorData[num].fill(vectorData[num].begin(), vectorData[num].end(), value); - vectorDataLabels[num] = label; + vectorData_[num].fill(vectorData_[num].begin(), vectorData_[num].end(), + value); + vectorDataLabels_[num] = label; } void resizeAllVectorData(size_t size, NumericType val = 0) { - for (auto &vector : vectorData) { + for (auto &vector : vectorData_) { vector.clear(); vector.resize(size, val); } @@ -137,71 +140,71 @@ template class rayTracingData { void setVectorMergeType(const std::vector &mergeType) { - vectorDataMerge = mergeType; + vectorDataMerge_ = mergeType; } void setVectorMergeType(int num, rayTracingDataMergeEnum mergeType) { - vectorDataMerge[num] = mergeType; + vectorDataMerge_[num] = mergeType; } void setScalarMergeType(const std::vector &mergeType) { - scalarDataMerge = mergeType; + scalarDataMerge_ = mergeType; } void setScalarMergeType(int num, rayTracingDataMergeEnum mergeType) { - scalarDataMerge[num] = mergeType; + scalarDataMerge_[num] = mergeType; } - vectorDataType &getVectorData(int i) { return vectorData[i]; } + vectorDataType &getVectorData(int i) { return vectorData_[i]; } - const vectorDataType &getVectorData(int i) const { return vectorData[i]; } + const vectorDataType &getVectorData(int i) const { return vectorData_[i]; } vectorDataType &getVectorData(std::string label) { int idx = getVectorDataIndex(label); - return vectorData[idx]; + return vectorData_[idx]; } - std::vector &getVectorData() { return vectorData; } + std::vector &getVectorData() { return vectorData_; } const std::vector &getVectorData() const { - return vectorData; + return vectorData_; } - scalarDataType &getScalarData(int i) { return scalarData[i]; } + scalarDataType &getScalarData(int i) { return scalarData_[i]; } - const scalarDataType &getScalarData(int i) const { return scalarData[i]; } + const scalarDataType &getScalarData(int i) const { return scalarData_[i]; } scalarDataType &getScalarData(std::string label) { int idx = getScalarDataIndex(label); - return scalarData[idx]; + return scalarData_[idx]; } - std::vector &getScalarData() { return scalarData; } + std::vector &getScalarData() { return scalarData_; } const std::vector &getScalarData() const { - return scalarData; + return scalarData_; } std::string getVectorDataLabel(int i) const { - if (i >= vectorDataLabels.size()) + if (i >= vectorDataLabels_.size()) rayMessage::getInstance() .addError("Getting vector data label in rayTracingData out of range.") .print(); - return vectorDataLabels[i]; + return vectorDataLabels_[i]; } std::string getScalarDataLabel(int i) const { - if (i >= scalarDataLabels.size()) + if (i >= scalarDataLabels_.size()) rayMessage::getInstance() .addError("Getting scalar data label in rayTracingData out of range.") .print(); - return scalarDataLabels[i]; + return scalarDataLabels_[i]; } - int getVectorDataIndex(std::string label) { - for (int i = 0; i < vectorDataLabels.size(); ++i) { - if (vectorDataLabels[i] == label) { + int getVectorDataIndex(std::string label) const { + for (int i = 0; i < vectorDataLabels_.size(); ++i) { + if (vectorDataLabels_[i] == label) { return i; } } @@ -211,9 +214,9 @@ template class rayTracingData { return -1; } - int getScalarDataIndex(std::string label) { - for (int i = 0; i < scalarDataLabels.size(); ++i) { - if (scalarDataLabels[i] == label) { + int getScalarDataIndex(std::string label) const { + for (int i = 0; i < scalarDataLabels_.size(); ++i) { + if (scalarDataLabels_[i] == label) { return i; } } @@ -224,12 +227,10 @@ template class rayTracingData { } const rayTracingDataMergeEnum getVectorMergeType(int num) const { - return vectorDataMerge[num]; + return vectorDataMerge_[num]; } const rayTracingDataMergeEnum getScalarMergeType(int num) const { - return scalarDataMerge[num]; + return scalarDataMerge_[num]; } }; - -#endif \ No newline at end of file diff --git a/include/viennaray/rayUtil.hpp b/include/viennaray/rayUtil.hpp index c26396b..5b76052 100644 --- a/include/viennaray/rayUtil.hpp +++ b/include/viennaray/rayUtil.hpp @@ -1,5 +1,20 @@ -#ifndef RAY_UTIL_HPP -#define RAY_UTIL_HPP +#pragma once + +#include +#include + +#if VIENNARAY_EMBREE_VERSION < 4 +#include +#else +#include +#endif + +#if defined(__x86_64__) || defined(_M_X64) +#define ARCH_X86 +#include +#endif + +#include #include #include @@ -8,22 +23,45 @@ #include #include #include -#include -#include -#include -#include #include +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + template using rayPair = std::array; template using rayTriple = std::array; template using rayQuadruple = std::array; -// embree uses float internally -using rtcNumericType = float; +enum class rayNormalizationType : unsigned { SOURCE = 0, MAX = 1 }; + +enum class rayTraceDirection : unsigned { + POS_X = 0, + NEG_X = 1, + POS_Y = 2, + NEG_Y = 3, + POS_Z = 4, + NEG_Z = 5 +}; -enum struct rayNormalizationType : unsigned { SOURCE = 0, MAX = 1 }; +template struct rayDataLog { + + std::vector> data; + + void merge(rayDataLog &pOther) { + assert(pOther.data.size() == data.size() && + "Size mismatch when merging logs"); + for (std::size_t i = 0; i < data.size(); i++) { + assert(pOther.data[i].size() == data[i].size() && + "Size mismatch when merging log data"); + for (std::size_t j = 0; j < data[i].size(); j++) { + data[i][j] += pOther.data[i][j]; + } + } + } +}; struct rayTraceInfo { size_t numRays; @@ -37,7 +75,9 @@ struct rayTraceInfo { }; namespace rayInternal { -constexpr double PI = 3.14159265358979323846; + +// embree uses float internally +using rtcNumericType = float; template constexpr double DiskFactor = 0.5 * (D == 3 ? 1.7320508 : 1.41421356237) * @@ -260,7 +300,7 @@ std::array getTraceSettings(rayTraceDirection sourceDir) { /* ------------------------------------------------------ */ template -static rayTriple PickRandomPointOnUnitSphere(rayRNG &RNG) { +static rayTriple pickRandomPointOnUnitSphere(rayRNG &RNG) { std::uniform_real_distribution uniDist; NumericType x, y, z, x2, y2, x2py2; do { @@ -277,15 +317,15 @@ static rayTriple PickRandomPointOnUnitSphere(rayRNG &RNG) { return rayTriple{x, y, z}; } -// Returns some orthonormal basis containing a the input vector pVector +// Returns some orthonormal basis containing a the input vector // (possibly scaled) as the first element of the return value. // This function is deterministic, i.e., for one input it will return always // the same result. template rayTriple> -getOrthonormalBasis(const rayTriple &pVector) { +getOrthonormalBasis(const rayTriple &vec) { rayTriple> rr; - rr[0] = pVector; + rr[0] = vec; // Calculate a vector (rr[1]) which is perpendicular to rr[0] // https://math.stackexchange.com/questions/137362/how-to-find-perpendicular-vector-to-another-vector#answer-211195 @@ -314,12 +354,11 @@ getOrthonormalBasis(const rayTriple &pVector) { rayInternal::Normalize(rr[2]); // Sanity check - const double eps = 1e-6; - assert(std::abs(rayInternal::DotProduct(rr[0], rr[1])) < eps && + assert(std::abs(rayInternal::DotProduct(rr[0], rr[1])) < 1e-6 && "Error in orthonormal basis computation"); - assert(std::abs(rayInternal::DotProduct(rr[1], rr[2])) < eps && + assert(std::abs(rayInternal::DotProduct(rr[1], rr[2])) < 1e-6 && "Error in orthonormal basis computation"); - assert(std::abs(rayInternal::DotProduct(rr[2], rr[0])) < eps && + assert(std::abs(rayInternal::DotProduct(rr[2], rr[0])) < 1e-6 && "Error in orthonormal basis computation"); return rr; } @@ -376,7 +415,7 @@ void readGridFromFile(std::string fileName, NumericType &gridDelta, template void writeVTK(std::string filename, const std::vector> &points, - const std::vector &mcestimates) { + const std::vector &flux) { std::ofstream f(filename.c_str()); f << "# vtk DataFile Version 2.0" << std::endl; @@ -401,12 +440,11 @@ void writeVTK(std::string filename, for (unsigned i = 0; i < points.size(); ++i) f << 1 << std::endl; - f << "CELL_DATA " << mcestimates.size() << std::endl; + f << "CELL_DATA " << flux.size() << std::endl; f << "SCALARS flux float" << std::endl; f << "LOOKUP_TABLE default" << std::endl; - for (unsigned j = 0; j < mcestimates.size(); ++j) { - f << ((std::abs(mcestimates[j]) < 1e-6) ? 0.0 : mcestimates[j]) - << std::endl; + for (unsigned j = 0; j < flux.size(); ++j) { + f << ((std::abs(flux[j]) < 1e-6) ? 0.0 : flux[j]) << std::endl; } f.close(); @@ -489,5 +527,3 @@ void printBoundingBox(const rayPair> &bdBox) { } /* ------------------------------------------------------- */ } // namespace rayInternal - -#endif // RAY_UTIL_HPP \ No newline at end of file diff --git a/tests/buildBoundary/buildBoundary.cpp b/tests/buildBoundary/buildBoundary.cpp index 9227259..39afd00 100644 --- a/tests/buildBoundary/buildBoundary.cpp +++ b/tests/buildBoundary/buildBoundary.cpp @@ -36,12 +36,6 @@ int main() { // assert boundary is extended in x direction RAYTEST_ASSERT_ISCLOSE(boundingBox[1][0], (1 + 2 * gridDelta), eps) - // assert boundary normal vectors are perpendicular to x direction - auto xplane = rayTriple{1., 0., 0.}; - for (unsigned int i = 0; i < 8; i++) { - auto normal = boundary.getPrimNormal(i); - RAYTEST_ASSERT_ISNORMAL(normal, xplane, eps) - } boundary.releaseGeometry(); } @@ -62,12 +56,6 @@ int main() { // assert boundary is extended in y direction RAYTEST_ASSERT_ISCLOSE(boundingBox[1][1], (1 + 2 * gridDelta), eps) - // assert boundary normal vectors are perpendicular to y direction - auto yplane = rayTriple{0., 1., 0.}; - for (unsigned int i = 0; i < 8; i++) { - auto normal = boundary.getPrimNormal(i); - RAYTEST_ASSERT_ISNORMAL(normal, yplane, eps) - } boundary.releaseGeometry(); } @@ -88,12 +76,6 @@ int main() { // assert boundary is extended in x direction RAYTEST_ASSERT_ISCLOSE(boundingBox[1][2], (1 + 2 * gridDelta), eps) - // assert boundary normal vectors are perpendicular to x direction - auto zplane = rayTriple{0., 0., 1.}; - for (unsigned int i = 0; i < 8; i++) { - auto normal = boundary.getPrimNormal(i); - RAYTEST_ASSERT_ISNORMAL(normal, zplane, eps) - } boundary.releaseGeometry(); } diff --git a/tests/buildBoundary2D/buildBoundary2D.cpp b/tests/buildBoundary2D/buildBoundary2D.cpp index 103f461..f7460ee 100644 --- a/tests/buildBoundary2D/buildBoundary2D.cpp +++ b/tests/buildBoundary2D/buildBoundary2D.cpp @@ -36,12 +36,7 @@ int main() { // assert boundary is extended in x direction RAYTEST_ASSERT_ISCLOSE(boundingBox[1][0], (1 + 2 * gridDelta), eps) - // assert boundary normal vectors are perpendicular to x direction - auto xplane = rayTriple{1., 0., 0.}; - for (unsigned int i = 0; i < 8; i++) { - auto normal = boundary.getPrimNormal(i); - RAYTEST_ASSERT_ISNORMAL(normal, xplane, eps) - } + boundary.releaseGeometry(); } { @@ -63,12 +58,7 @@ int main() { // assert boundary is extended in y direction RAYTEST_ASSERT_ISCLOSE(boundingBox[1][1], (1 + 2 * gridDelta), eps) - // assert boundary normal vectors are perpendicular to y direction - auto yplane = rayTriple{0., 1., 0.}; - for (unsigned int i = 0; i < 8; i++) { - auto normal = boundary.getPrimNormal(i); - RAYTEST_ASSERT_ISNORMAL(normal, yplane, eps) - } + boundary.releaseGeometry(); } rtcReleaseDevice(device); diff --git a/tests/diskAreas/diskAreas.cpp b/tests/diskAreas/diskAreas.cpp index 0859eb6..53314e1 100644 --- a/tests/diskAreas/diskAreas.cpp +++ b/tests/diskAreas/diskAreas.cpp @@ -1,4 +1,3 @@ -#include #include #include #include @@ -48,15 +47,15 @@ int main() { auto cp = particle.clone(); rayDataLog log; - auto tracer = rayTraceKernel( - device, geometry, boundary, raySource, cp, log, 1, 0, false, true, 0); + rayTraceInfo info; + auto tracer = rayTraceKernel(device, geometry, boundary, raySource, cp, log, + 1, 0, false, true, 0, hitCounter, info); tracer.setTracingData(&localData, &globalData); - tracer.setHitCounter(&hitCounter); tracer.apply(); auto diskAreas = hitCounter.getDiskAreas(); auto boundaryDirs = boundary.getDirs(); - auto wholeDiskArea = diskRadius * diskRadius * rayInternal::PI; + auto wholeDiskArea = diskRadius * diskRadius * M_PI; for (unsigned int idx = 0; idx < geometry.getNumPoints(); ++idx) { auto const &disk = geometry.getPrimRef(idx); if (std::fabs(disk[boundaryDirs[0]] - boundingBox[0][boundaryDirs[0]]) <