Skip to content

Commit

Permalink
Add source area function (#58)
Browse files Browse the repository at this point in the history
  • Loading branch information
@tobre1
tobre1 authored May 23, 2024
1 parent df6adb9 commit 702cf6a
Show file tree
Hide file tree
Showing 10 changed files with 75 additions and 63 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 1 addition & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum
* Installation with CPM

```cmake
CPMAddPackage("gh:viennatools/[email protected].1")
CPMAddPackage("gh:viennatools/[email protected].3")
```

* With a local installation
Expand Down Expand Up @@ -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: [email protected]

ViennaRay was developed under the aegis of the 'Institute for Microelectronics' at the 'TU Wien'.
Expand Down
6 changes: 4 additions & 2 deletions include/viennaray/rayBoundary.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ template <typename NumericType, int D> 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);
Expand Down Expand Up @@ -120,6 +120,8 @@ template <typename NumericType, int D> class rayBoundary {
}
}

auto getBoundaryConditions() const { return boundaryConds_; }

RTCGeometry const &getRTCGeometry() const { return pRtcBoundary_; }

void releaseGeometry() {
Expand Down Expand Up @@ -234,7 +236,7 @@ template <typename NumericType, int D> class rayBoundary {
}

rayTriple<rayTriple<NumericType>> 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,
Expand Down
1 change: 1 addition & 0 deletions include/viennaray/raySource.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,6 @@ template <typename NumericType> class raySource {
virtual rayPair<rayTriple<NumericType>>
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.; }
};
24 changes: 19 additions & 5 deletions include/viennaray/raySourceGrid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,18 @@

template <typename NumericType, int D>
class raySourceGrid : public raySource<NumericType> {
using boundingBoxType = rayPair<rayTriple<NumericType>>;

public:
raySourceGrid(std::vector<rayTriple<NumericType>> &sourceGrid,
raySourceGrid(const boundingBoxType &boundingBox,
std::vector<rayTriple<NumericType>> &sourceGrid,
NumericType cosinePower,
const std::array<int, 5> &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<rayTriple<NumericType>>
getOriginAndDirection(const size_t idx, rayRNG &RngState) const override {
Expand Down Expand Up @@ -45,7 +49,17 @@ class raySourceGrid : public raySource<NumericType> {
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<rayTriple<NumericType>> &sourceGrid_;
const size_t numPoints_;
const int rayDir_;
Expand Down
11 changes: 10 additions & 1 deletion include/viennaray/raySourceRandom.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ class raySourceRandom : public raySource<NumericType> {

public:
raySourceRandom(
boundingBoxType boundingBox, NumericType cosinePower,
const boundingBoxType &boundingBox, NumericType cosinePower,
std::array<int, 5> &pTraceSettings, const size_t numPoints_,
const bool customDirection,
const std::array<std::array<NumericType, 3>, 3> &orthonormalBasis)
Expand All @@ -34,6 +34,15 @@ class raySourceRandom : public raySource<NumericType> {

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<NumericType> getOrigin(rayRNG &RngState) const {
rayTriple<NumericType> origin{0., 0., 0.};
Expand Down
34 changes: 7 additions & 27 deletions include/viennaray/rayTrace.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,13 @@ template <class NumericType, int D> 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_;
Expand Down Expand Up @@ -285,32 +291,6 @@ template <class NumericType, int D> class rayTrace {
[[nodiscard]] rayDataLog<NumericType> &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();
Expand Down
47 changes: 26 additions & 21 deletions include/viennaray/rayTraceKernel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,17 +171,17 @@ template <typename NumericType, int D> class rayTraceKernel {
std::uniform_real_distribution<NumericType> 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<rayInternal::rtcNumericType, 3> origin = {
static_cast<rayInternal::rtcNumericType>(ray.org_x +
ray.dir_x * rndm),
ray.dir_x * rnd),
static_cast<rayInternal::rtcNumericType>(ray.org_y +
ray.dir_y * rndm),
ray.dir_y * rnd),
static_cast<rayInternal::rtcNumericType>(ray.org_z +
ray.dir_z * rndm)};
ray.dir_z * rnd)};

std::array<rayInternal::rtcNumericType, 3> direction{0, 0, 0};
for (int i = 0; i < D; ++i) {
Expand Down Expand Up @@ -435,9 +435,9 @@ template <typename NumericType, int D> 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;
}
Expand All @@ -450,6 +450,7 @@ template <typename NumericType, int D> class rayTraceKernel {
std::vector<NumericType> 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<NumericType>(numOfPrimitives, 0);
Expand All @@ -461,18 +462,20 @@ template <typename NumericType, int D> 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;
}
Expand All @@ -481,8 +484,9 @@ template <typename NumericType, int D> 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) {
Expand All @@ -496,8 +500,9 @@ template <typename NumericType, int D> 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) {
Expand Down Expand Up @@ -590,8 +595,8 @@ template <typename NumericType, int D> class rayTraceKernel {
auto rayDirectionC = rayTriple<rayInternal::rtcNumericType>{
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;
Expand Down
7 changes: 4 additions & 3 deletions tests/createSourceGrid/createSourceGrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,16 @@ int main() {
auto grid = rayInternal::createSourceGrid<NumericType, D>(
boundingBox, points.size(), gridDelta, traceSettings);

rayRNG rngstate(0);
rayRNG rngState(0);
{
// build source in positive z direction;
auto source = raySourceGrid<NumericType, D>(grid, 1., traceSettings);
auto source =
raySourceGrid<NumericType, D>(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]);

Expand Down
2 changes: 2 additions & 0 deletions tests/traceInterface/traceInterface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@ class MySource : public raySource<NumericType> {
}

size_t getNumPoints() const override { return 0; }

NumericType getSourceArea() const override { return 1; }
};

int main() {
Expand Down

0 comments on commit 702cf6a

Please sign in to comment.