diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ef4564..e52eea9 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.1.1) + VERSION 2.1.3) # -------------------------------------------------------------------------------------------------------- # Library switches diff --git a/README.md b/README.md index 3601115..5228c6c 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.1.1") + CPMAddPackage("gh:viennatools/viennaray@2.1.3") ``` * With a local installation @@ -110,8 +110,6 @@ If you want to contribute to ViennaRay, make sure to follow the [LLVM Coding gui ## Authors -Current contributors: Tobias Reiter - Contact us via: viennatools@iue.tuwien.ac.at ViennaRay was developed under the aegis of the 'Institute for Microelectronics' at the 'TU Wien'. diff --git a/include/viennaray/rayBoundary.hpp b/include/viennaray/rayBoundary.hpp index 8ecef37..ef67a2c 100644 --- a/include/viennaray/rayBoundary.hpp +++ b/include/viennaray/rayBoundary.hpp @@ -79,7 +79,7 @@ template class rayBoundary { // hit at firstDir min boundary -> move to max firstDir impactCoords[firstDir_] = bdBox_[1][firstDir_]; } else if (primID == 2 || primID == 3) { - // hit at firstDir max boundary -> move to min fristDir + // hit at firstDir max boundary -> move to min firstDir impactCoords[firstDir_] = bdBox_[0][firstDir_]; } assignRayCoords(rayHit, impactCoords); @@ -120,6 +120,8 @@ template class rayBoundary { } } + auto getBoundaryConditions() const { return boundaryConds_; } + RTCGeometry const &getRTCGeometry() const { return pRtcBoundary_; } void releaseGeometry() { @@ -234,7 +236,7 @@ template class rayBoundary { } rayTriple> getTriangleCoords(const size_t primID) { - assert(primID < numTriangles_ && "rtBounday: primID out of bounds"); + assert(primID < numTriangles_ && "rtBoundary: primID out of bounds"); auto tt = pTriangleBuffer_[primID]; return {(NumericType)pVertexBuffer_[tt.v0].xx, (NumericType)pVertexBuffer_[tt.v0].yy, diff --git a/include/viennaray/raySource.hpp b/include/viennaray/raySource.hpp index 84fc959..a63f8ce 100644 --- a/include/viennaray/raySource.hpp +++ b/include/viennaray/raySource.hpp @@ -10,5 +10,6 @@ template class raySource { virtual rayPair> getOriginAndDirection(const size_t idx, rayRNG &RngState) const = 0; virtual size_t getNumPoints() const = 0; + virtual NumericType getSourceArea() const = 0; virtual NumericType getInitialRayWeight(const size_t idx) const { return 1.; } }; diff --git a/include/viennaray/raySourceGrid.hpp b/include/viennaray/raySourceGrid.hpp index 47bf39d..d6bdbfc 100644 --- a/include/viennaray/raySourceGrid.hpp +++ b/include/viennaray/raySourceGrid.hpp @@ -4,14 +4,18 @@ template class raySourceGrid : public raySource { + using boundingBoxType = rayPair>; + public: - raySourceGrid(std::vector> &sourceGrid, + raySourceGrid(const boundingBoxType &boundingBox, + std::vector> &sourceGrid, 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)) {} + : bdBox_(boundingBox), 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)) {} rayPair> getOriginAndDirection(const size_t idx, rayRNG &RngState) const override { @@ -45,7 +49,17 @@ class raySourceGrid : public raySource { return direction; } + NumericType getSourceArea() const override { + if constexpr (D == 2) { + return (bdBox_[1][firstDir_] - bdBox_[0][firstDir_]); + } else { + return (bdBox_[1][firstDir_] - bdBox_[0][firstDir_]) * + (bdBox_[1][secondDir_] - bdBox_[0][secondDir_]); + } + } + private: + const boundingBoxType bdBox_; const std::vector> &sourceGrid_; const size_t numPoints_; const int rayDir_; diff --git a/include/viennaray/raySourceRandom.hpp b/include/viennaray/raySourceRandom.hpp index 16e297f..b8b7e0b 100644 --- a/include/viennaray/raySourceRandom.hpp +++ b/include/viennaray/raySourceRandom.hpp @@ -8,7 +8,7 @@ class raySourceRandom : public raySource { public: raySourceRandom( - boundingBoxType boundingBox, NumericType cosinePower, + const boundingBoxType &boundingBox, NumericType cosinePower, std::array &pTraceSettings, const size_t numPoints_, const bool customDirection, const std::array, 3> &orthonormalBasis) @@ -34,6 +34,15 @@ class raySourceRandom : public raySource { size_t getNumPoints() const override { return numPoints_; } + NumericType getSourceArea() const override { + if constexpr (D == 2) { + return (bdBox_[1][firstDir_] - bdBox_[0][firstDir_]); + } else { + return (bdBox_[1][firstDir_] - bdBox_[0][firstDir_]) * + (bdBox_[1][secondDir_] - bdBox_[0][secondDir_]); + } + } + private: rayTriple getOrigin(rayRNG &RngState) const { rayTriple origin{0., 0., 0.}; diff --git a/include/viennaray/rayTrace.hpp b/include/viennaray/rayTrace.hpp index 82c5e72..39c49a1 100644 --- a/include/viennaray/rayTrace.hpp +++ b/include/viennaray/rayTrace.hpp @@ -220,7 +220,13 @@ template class rayTrace { } case rayNormalizationType::SOURCE: { - NumericType sourceArea = getSourceArea(); + if (!pSource_) { + rayMessage::getInstance() + .addWarning("No source was specified in for the normalization.") + .print(); + break; + } + NumericType sourceArea = pSource_->getSourceArea(); auto numTotalRays = numberOfRaysFixed_ == 0 ? flux.size() * numberOfRaysPerPoint_ : numberOfRaysFixed_; @@ -285,32 +291,6 @@ template class rayTrace { [[nodiscard]] rayDataLog &getDataLog() { return dataLog_; } private: - NumericType getSourceArea() { - const auto boundingBox = geometry_.getBoundingBox(); - NumericType sourceArea = 0; - - 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 (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 (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]); - } - - return sourceArea; - } - void checkRelativeError() { auto error = getRelativeError(); const int numPoints = error.size(); diff --git a/include/viennaray/rayTraceKernel.hpp b/include/viennaray/rayTraceKernel.hpp index 0515724..bc5f036 100644 --- a/include/viennaray/rayTraceKernel.hpp +++ b/include/viennaray/rayTraceKernel.hpp @@ -171,17 +171,17 @@ template class rayTraceKernel { std::uniform_real_distribution dist(0., 1.); NumericType scatterProbability = 1 - std::exp(-rayHit.ray.tfar / lambda); - auto rndm = dist(RngState); - if (rndm < scatterProbability) { + auto rnd = dist(RngState); + if (rnd < scatterProbability) { const auto &ray = rayHit.ray; std::array origin = { static_cast(ray.org_x + - ray.dir_x * rndm), + ray.dir_x * rnd), static_cast(ray.org_y + - ray.dir_y * rndm), + ray.dir_y * rnd), static_cast(ray.org_z + - ray.dir_z * rndm)}; + ray.dir_z * rnd)}; std::array direction{0, 0, 0}; for (int i = 0; i < D; ++i) { @@ -435,9 +435,9 @@ template class rayTraceKernel { // We want to set the weight of (the reflection of) the ray to the value of // renewWeight. In order to stay unbiased we kill the reflection with a // probability of (1 - rayWeight / renewWeight). - auto rndm = RNG(); + auto rnd = RNG(); auto killProbability = 1.0 - rayWeight / renewWeight; - if (rndm < (killProbability * RNG.max())) { + if (rnd < (killProbability * RNG.max())) { // kill the ray return false; } @@ -450,6 +450,7 @@ template class rayTraceKernel { std::vector computeDiskAreas() const { constexpr double eps = 1e-3; auto bdBox = geometry_.getBoundingBox(); + auto boundaryConds = boundary_.getBoundaryConditions(); const auto numOfPrimitives = geometry_.getNumPoints(); const auto boundaryDirs = boundary_.getDirs(); auto areas = std::vector(numOfPrimitives, 0); @@ -461,18 +462,20 @@ template class rayTraceKernel { if constexpr (D == 3) { areas[idx] = disk[3] * disk[3] * M_PI; // full disk area - if (std::fabs(disk[boundaryDirs[0]] - bdBox[0][boundaryDirs[0]]) < - eps || - std::fabs(disk[boundaryDirs[0]] - bdBox[1][boundaryDirs[0]]) < - eps) { + if ((boundaryConds[boundaryDirs[0]] != rayBoundaryCondition::IGNORE) && + (std::fabs(disk[boundaryDirs[0]] - bdBox[0][boundaryDirs[0]]) < + eps || + std::fabs(disk[boundaryDirs[0]] - bdBox[1][boundaryDirs[0]]) < + eps)) { // disk intersects boundary in first direction areas[idx] /= 2; } - if (std::fabs(disk[boundaryDirs[1]] - bdBox[0][boundaryDirs[1]]) < - eps || - std::fabs(disk[boundaryDirs[1]] - bdBox[1][boundaryDirs[1]]) < - eps) { + if ((boundaryConds[boundaryDirs[1]] != rayBoundaryCondition::IGNORE) && + (std::fabs(disk[boundaryDirs[1]] - bdBox[0][boundaryDirs[1]]) < + eps || + std::fabs(disk[boundaryDirs[1]] - bdBox[1][boundaryDirs[1]]) < + eps)) { // disk intersects boundary in second direction areas[idx] /= 2; } @@ -481,8 +484,9 @@ template class rayTraceKernel { auto normal = geometry_.getNormalRef(idx); // test min boundary - if (std::abs(disk[boundaryDirs[0]] - bdBox[0][boundaryDirs[0]]) < - disk[3]) { + if ((boundaryConds[boundaryDirs[0]] != rayBoundaryCondition::IGNORE) && + (std::abs(disk[boundaryDirs[0]] - bdBox[0][boundaryDirs[0]]) < + disk[3])) { NumericType insideTest = 1 - normal[boundaryDirs[0]] * normal[boundaryDirs[0]]; if (insideTest > 1e-4) { @@ -496,8 +500,9 @@ template class rayTraceKernel { } // test max boundary - if (std::abs(disk[boundaryDirs[0]] - bdBox[1][boundaryDirs[0]]) < - disk[3]) { + if ((boundaryConds[boundaryDirs[0]] != rayBoundaryCondition::IGNORE) && + (std::abs(disk[boundaryDirs[0]] - bdBox[1][boundaryDirs[0]]) < + disk[3])) { NumericType insideTest = 1 - normal[boundaryDirs[0]] * normal[boundaryDirs[0]]; if (insideTest > 1e-4) { @@ -590,8 +595,8 @@ template class rayTraceKernel { auto rayDirectionC = rayTriple{ rayDirection[0], rayDirection[1], rayDirection[2]}; rayInternal::Scale(tt, rayDirectionC); - auto hitpoint = rayInternal::Sum(rayOrigin, rayDirectionC); - auto distance = rayInternal::Distance(hitpoint, diskOrigin); + auto hitPoint = rayInternal::Sum(rayOrigin, rayDirectionC); + auto distance = rayInternal::Distance(hitPoint, diskOrigin); auto const &radius = disk[3]; if (radius > distance) { impactDistance = distance; diff --git a/tests/createSourceGrid/createSourceGrid.cpp b/tests/createSourceGrid/createSourceGrid.cpp index d0598ba..51e0ffa 100644 --- a/tests/createSourceGrid/createSourceGrid.cpp +++ b/tests/createSourceGrid/createSourceGrid.cpp @@ -26,15 +26,16 @@ int main() { auto grid = rayInternal::createSourceGrid( boundingBox, points.size(), gridDelta, traceSettings); - rayRNG rngstate(0); + rayRNG rngState(0); { // build source in positive z direction; - auto source = raySourceGrid(grid, 1., traceSettings); + auto source = + raySourceGrid(boundingBox, grid, 1., traceSettings); auto numGridPoints = source.getNumPoints(); alignas(128) auto rayhit = RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; for (size_t i = 0; i < numGridPoints; ++i) { - auto originAndDirection = source.getOriginAndDirection(i, rngstate); + auto originAndDirection = source.getOriginAndDirection(i, rngState); rayInternal::fillRay(rayhit.ray, originAndDirection[0], originAndDirection[1]); diff --git a/tests/traceInterface/traceInterface.cpp b/tests/traceInterface/traceInterface.cpp index 1ff7d98..a76e9d2 100644 --- a/tests/traceInterface/traceInterface.cpp +++ b/tests/traceInterface/traceInterface.cpp @@ -15,6 +15,8 @@ class MySource : public raySource { } size_t getNumPoints() const override { return 0; } + + NumericType getSourceArea() const override { return 1; } }; int main() {