Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
IP
elsa
Commits
432c010b
Commit
432c010b
authored
Dec 10, 2019
by
Jens Petit
Committed by
Tobias Lasser
Dec 10, 2019
Browse files
#32
fix warnings and clang-tidy issues everywhere
parent
47915505
Pipeline
#191720
passed with stages
in 7 minutes and 59 seconds
Changes
56
Pipelines
1
Expand all
Hide whitespace changes
Inline
Side-by-side
elsa/projectors/JosephsMethod.h
View file @
432c010b
...
...
@@ -119,8 +119,9 @@ namespace elsa
* angle \param[in] from index of the current vector position \param[in] to index of the
* current result position \param[in] mainDirection specifies the main direction of the ray
*/
template
<
bool
adjoint
>
void
linear
(
const
DataContainer
<
data_t
>&
vector
,
DataContainer
<
data_t
>&
result
,
const
RealVector_t
&
fractionals
,
bool
adjoint
,
in
t
domainDim
,
const
RealVector_t
&
fractionals
,
index_
t
domainDim
,
const
IndexVector_t
&
currentVoxel
,
float
intersection
,
index_t
from
,
index_t
to
,
int
mainDirection
)
const
;
...
...
elsa/projectors/SiddonsMethod.cpp
View file @
432c010b
...
...
@@ -11,8 +11,8 @@ namespace elsa
const
DataDescriptor
&
rangeDescriptor
,
const
std
::
vector
<
Geometry
>&
geometryList
)
:
LinearOperator
<
data_t
>
(
domainDescriptor
,
rangeDescriptor
),
_
geometryList
{
geometryList
}
,
_
boundingBox
(
domainDescriptor
.
getNumberOfCoefficientsPerDimension
())
_
boundingBox
(
domainDescriptor
.
getNumberOfCoefficientsPerDimension
())
,
_
geometryList
{
geometryList
}
{
auto
dim
=
_domainDescriptor
->
getNumberOfDimensions
();
if
(
dim
!=
_rangeDescriptor
->
getNumberOfDimensions
())
{
...
...
@@ -71,26 +71,25 @@ namespace elsa
void
SiddonsMethod
<
data_t
>::
traverseVolume
(
const
DataContainer
<
data_t
>&
vector
,
DataContainer
<
data_t
>&
result
)
const
{
index_t
maxIterations
{
0
}
;
if
(
adjoint
)
{
maxIterations
=
vector
.
getSize
();
const
index_t
maxIterations
=
adjoint
?
vector
.
getSize
()
:
result
.
getSize
()
;
if
constexpr
(
adjoint
)
{
result
=
0
;
// initialize volume to 0, because we are not going to hit every voxel!
}
else
maxIterations
=
result
.
getSize
();
}
const
auto
rangeDim
=
_rangeDescriptor
->
getNumberOfDimensions
();
// --> loop either over every voxel that should updated or every detector
// cell that should be calculated
#pragma omp parallel for
for
(
size
_t
rangeIndex
=
0
;
rangeIndex
<
maxIterations
;
++
rangeIndex
)
{
for
(
index
_t
rangeIndex
=
0
;
rangeIndex
<
maxIterations
;
++
rangeIndex
)
{
// --> get the current ray to the detector center
auto
ray
=
computeRayToDetector
(
rangeIndex
,
rangeDim
);
// --> setup traversal algorithm
TraverseAABB
traverse
(
_boundingBox
,
ray
);
if
(
!
adjoint
)
if
constexpr
(
!
adjoint
)
result
[
rangeIndex
]
=
0
;
// --> initial index to access the data vector
...
...
@@ -101,7 +100,7 @@ namespace elsa
auto
weight
=
traverse
.
updateTraverseAndGetDistance
();
// --> update result depending on the operation performed
if
(
adjoint
)
if
constexpr
(
adjoint
)
#pragma omp atomic
result
[
dataIndexForCurrentVoxel
]
+=
vector
[
rangeIndex
]
*
weight
;
else
...
...
elsa/projectors/TraverseAABB.cpp
View file @
432c010b
...
...
@@ -14,7 +14,7 @@ namespace elsa
initStepDirection
(
r
.
direction
());
// select the initial voxel (closest to the entry point)
selectClosestVoxel
(
r
.
direction
()
);
selectClosestVoxel
();
// initialize the step sizes for the step parameter
initDelta
(
r
.
direction
());
...
...
@@ -65,7 +65,7 @@ namespace elsa
void
TraverseAABB
::
calculateAABBIntersections
(
const
TraverseAABB
::
Ray
&
r
)
{
// entry and exit point parameters
real_t
tmin
,
tmax
;
real_t
tmin
;
// --> calculate intersection parameter and if the volume is hit
auto
opt
=
Intersection
::
withRay
(
_aabb
,
r
);
...
...
@@ -91,7 +91,7 @@ namespace elsa
_stepDirection
=
rd
.
array
().
sign
();
}
void
TraverseAABB
::
selectClosestVoxel
(
const
RealVector_t
&
rd
)
void
TraverseAABB
::
selectClosestVoxel
()
{
RealVector_t
lowerCorner
=
_entryPoint
.
array
().
floor
();
lowerCorner
=
...
...
elsa/projectors/TraverseAABB.h
View file @
432c010b
...
...
@@ -86,15 +86,13 @@ namespace elsa
/// constant vector containing the maximum number
const
RealVector_t
_MAX
{
RealVector_t
(
_aabb
.
_dim
).
setConstant
(
std
::
numeric_limits
<
real_t
>::
max
())};
/// constant to decide whether we are in next voxel
const
real_t
_NEXT_VOXEL_THRESHOLD
{
0.01
};
/// compute the entry and exit points of ray r with the volume (aabb)
void
calculateAABBIntersections
(
const
Ray
&
r
);
/// setup the step directions (which is basically the sign of the ray direction rd)
void
initStepDirection
(
const
RealVector_t
&
rd
);
/// select the closest voxel to the entry point
(considering the ray direction rd)
void
selectClosestVoxel
(
const
RealVector_t
&
rd
);
/// select the closest voxel to the entry point
void
selectClosestVoxel
();
/// setup the step sizes considering the ray direction rd
void
initDelta
(
const
RealVector_t
&
rd
);
/// setup the maximum step parameters considering the ray direction rd
...
...
elsa/projectors/TraverseAABBJosephsMethod.cpp
View file @
432c010b
...
...
@@ -16,7 +16,7 @@ namespace elsa
if
(
!
isInBoundingBox
())
return
;
selectClosestVoxel
(
rt
.
direction
(),
_entryPoint
);
selectClosestVoxel
(
_entryPoint
);
// exit direction stored in _exitDirection
(
_exitPoint
-
(
_aabb
.
_max
-
_aabb
.
_min
)
/
2
).
cwiseAbs
().
maxCoeff
(
&
_exitDirection
);
...
...
@@ -44,7 +44,7 @@ namespace elsa
//*0.5 because a change of 1 spacing corresponds to a change of
// fractionals from -0.5 to 0.5 0.5 or -0.5 = voxel border is still
// ok
if
(
f
abs
(
_fractionals
(
i
))
>
0.5
)
{
if
(
std
::
abs
(
_fractionals
(
i
))
>
0.5
)
{
_fractionals
(
i
)
-=
_stepDirection
(
i
);
_currentPos
(
i
)
+=
_stepDirection
(
i
);
// --> is the traverse algorithm still in the bounding box?
...
...
@@ -56,7 +56,7 @@ namespace elsa
if
(
_isInAABB
)
{
// now ignore main direction
_nextStep
.
cwiseAbs
().
maxCoeff
(
&
_ignoreDirection
);
_nextStep
/=
f
abs
(
_nextStep
(
_ignoreDirection
));
_nextStep
/=
std
::
abs
(
_nextStep
(
_ignoreDirection
));
_fractionals
(
_ignoreDirection
)
=
0
;
_intersectionLength
=
_nextStep
.
norm
();
_stage
=
INTERIOR
;
...
...
@@ -76,7 +76,7 @@ namespace elsa
_ignoreDirection
=
_exitDirection
;
// move to last sampling point
RealVector_t
currentPosition
=
_exitPoint
-
_intersectionLength
*
_nextStep
/
2
;
selectClosestVoxel
(
_nextStep
,
currentPosition
);
selectClosestVoxel
(
currentPosition
);
initFractionals
(
currentPosition
);
_isInAABB
=
true
;
_stage
=
LAST
;
...
...
@@ -96,7 +96,8 @@ namespace elsa
void
TraverseAABBJosephsMethod
::
initFractionals
(
const
RealVector_t
&
currentPosition
)
{
for
(
int
i
=
0
;
i
<
_aabb
.
_dim
;
i
++
)
_fractionals
(
i
)
=
(
fabs
(
currentPosition
(
i
))
-
_currentPos
(
i
)
-
0.5
);
_fractionals
(
i
)
=
static_cast
<
real_t
>
(
std
::
abs
(
currentPosition
(
i
))
-
_currentPos
(
i
)
-
0.5
);
}
void
TraverseAABBJosephsMethod
::
moveToFirstSamplingPoint
(
const
RealVector_t
&
rd
,
...
...
@@ -106,11 +107,11 @@ namespace elsa
rd
.
cwiseAbs
().
maxCoeff
(
&
_ignoreDirection
);
// find distance to next plane orthogonal to main direction
int
nextBoundary
=
_currentPos
(
_ignoreDirection
);
auto
nextBoundary
=
static_cast
<
int
>
(
_currentPos
(
_ignoreDirection
)
)
;
if
(
_stepDirection
(
_ignoreDirection
)
>
0
)
nextBoundary
+=
1
;
real_t
distToBoundary
=
(
nextBoundary
-
_entryPoint
[
_ignoreDirection
])
/
rd
[
_ignoreDirection
];
real_t
distToBoundary
=
(
static_cast
<
real_t
>
(
nextBoundary
)
-
_entryPoint
[
_ignoreDirection
])
/
rd
[
_ignoreDirection
];
// intialize _nextStep as the step for interior pixels
_nextStep
(
_ignoreDirection
)
=
_stepDirection
(
_ignoreDirection
);
...
...
@@ -136,7 +137,7 @@ namespace elsa
_nextStep
=
_nextStep
/
2
+
_nextStep
*
(
distToBoundary
/
(
2
*
_nextStep
.
norm
()));
_intersectionLength
=
distToBoundary
;
}
selectClosestVoxel
(
rd
,
currentPosition
);
selectClosestVoxel
(
currentPosition
);
// init fractionals
initFractionals
(
currentPosition
);
...
...
@@ -147,12 +148,11 @@ namespace elsa
return
_currentPos
(
index
)
<
_aabb
.
_max
(
index
)
&&
_currentPos
(
index
)
>=
_aabb
.
_min
(
index
);
}
void
TraverseAABBJosephsMethod
::
selectClosestVoxel
(
const
RealVector_t
&
rd
,
const
RealVector_t
&
currentPosition
)
void
TraverseAABBJosephsMethod
::
selectClosestVoxel
(
const
RealVector_t
&
currentPosition
)
{
RealVector_t
lowerCorner
=
currentPosition
.
array
().
floor
();
lowerCorner
=
((
lowerCorner
.
array
()
==
currentPosition
.
array
())
&&
(
_stepDirection
.
array
()
<
0.0
)
((
(
lowerCorner
.
array
()
==
currentPosition
.
array
())
&&
(
_stepDirection
.
array
()
<
0.0
)
)
||
currentPosition
.
array
()
==
_aabb
.
_max
.
array
())
.
select
(
lowerCorner
.
array
()
-
1
,
lowerCorner
);
...
...
@@ -180,7 +180,7 @@ namespace elsa
real_t
TraverseAABBJosephsMethod
::
calculateAABBIntersections
(
const
Ray
&
ray
)
{
real_t
tmin
,
tmax
;
real_t
tmin
;
// --> calculate intersection parameter and if the volume is hit
auto
opt
=
Intersection
::
withRay
(
_aabb
,
ray
);
...
...
elsa/projectors/TraverseAABBJosephsMethod.h
View file @
432c010b
...
...
@@ -125,9 +125,8 @@ namespace elsa
/// compute the entry and exit points of ray r with the volume (aabb), returns the length of
/// the intersection
real_t
calculateAABBIntersections
(
const
Ray
&
ray
);
/// select the closest voxel to the entry point (considering the ray direction and the
/// current positon)
void
selectClosestVoxel
(
const
RealVector_t
&
rd
,
const
RealVector_t
&
currentPosition
);
/// select the closest voxel to the entry point (considering the current position)
void
selectClosestVoxel
(
const
RealVector_t
&
currentPosition
);
/// check if the current index is still in the bounding box
bool
isCurrentPositionInAABB
(
index_t
index
)
const
;
/// advances the traversal algorithm to the first sampling point
...
...
elsa/projectors/tests/test_BinaryMethod.cpp
View file @
432c010b
...
...
@@ -62,7 +62,7 @@ SCENARIO("Testing BinaryMethod with only one ray")
WHEN
(
"We have a single ray with 180 degrees"
)
{
geom
.
emplace_back
(
100
,
5
,
180
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
180
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
THEN
(
"A^t A x should be close to the original data"
)
...
...
@@ -85,7 +85,7 @@ SCENARIO("Testing BinaryMethod with only one ray")
WHEN
(
"We have a single ray with 90 degrees"
)
{
geom
.
emplace_back
(
100
,
5
,
90
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
90
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
THEN
(
"A^t A x should be close to the original data"
)
...
...
@@ -108,7 +108,7 @@ SCENARIO("Testing BinaryMethod with only one ray")
WHEN
(
"We have a single ray with 90 degrees"
)
{
geom
.
emplace_back
(
100
,
5
,
270
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
270
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
THEN
(
"A^t A x should be close to the original data"
)
...
...
@@ -131,7 +131,7 @@ SCENARIO("Testing BinaryMethod with only one ray")
WHEN
(
"We have a single ray with 90 degrees"
)
{
geom
.
emplace_back
(
100
,
5
,
45
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
45
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
THEN
(
"A^t A x should be close to the original data"
)
...
...
@@ -151,7 +151,7 @@ SCENARIO("Testing BinaryMethod with only one ray")
WHEN
(
"We have a single ray with 90 degrees"
)
{
geom
.
emplace_back
(
100
,
5
,
225
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
225
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
// This test case is a little awkward, but the Problem is inside of Geometry, with tiny
...
...
@@ -204,9 +204,9 @@ SCENARIO("Testing BinaryMethod with only 1 rays for 4 angles")
WHEN
(
"We have a single ray with 0, 90, 180, 270 degrees"
)
{
geom
.
emplace_back
(
100
,
5
,
0
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
90
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
180
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
270
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
90
*
pi
_t
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
180
*
pi
_t
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
270
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
THEN
(
"A^t A x should be close to the original data"
)
...
...
@@ -323,7 +323,7 @@ SCENARIO("Testing BinaryMethod")
GIVEN
(
"A traversal with 5 rays at 180 degrees"
)
{
std
::
vector
<
Geometry
>
geom
;
geom
.
emplace_back
(
100
,
5
,
180
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
180
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
...
...
@@ -356,7 +356,7 @@ SCENARIO("Testing BinaryMethod")
GIVEN
(
"A traversal with 5 rays at 90 degrees"
)
{
std
::
vector
<
Geometry
>
geom
;
geom
.
emplace_back
(
100
,
5
,
90
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
90
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
...
...
@@ -389,7 +389,7 @@ SCENARIO("Testing BinaryMethod")
GIVEN
(
"A traversal with 5 rays at 270 degrees"
)
{
std
::
vector
<
Geometry
>
geom
;
geom
.
emplace_back
(
100
,
5
,
270
*
pi
/
180.
,
domain
,
range
);
geom
.
emplace_back
(
100
,
5
,
270
*
pi
_t
/
180.
,
domain
,
range
);
auto
op
=
BinaryMethod
(
domain
,
range
,
geom
);
...
...
@@ -440,7 +440,7 @@ SCENARIO("Calls to functions of super class")
std
::
vector
<
Geometry
>
geom
;
for
(
index_t
i
=
0
;
i
<
numImgs
;
i
++
)
{
real_t
angle
=
i
*
2
*
pi
/
50
;
auto
angle
=
static_cast
<
real_t
>
(
i
*
2
)
*
pi
_t
/
50
;
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
angle
,
volumeDescriptor
,
sinoDescriptor
);
}
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
...
...
@@ -636,8 +636,8 @@ SCENARIO("Rays not intersecting the bounding box are present")
WHEN
(
"Tracing along a x-axis-aligned ray with a negative y-coordinate of origin"
)
{
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
pi
/
2
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
-
volSize
);
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
pi
_t
/
2
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
-
volSize
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
...
...
@@ -661,8 +661,8 @@ SCENARIO("Rays not intersecting the bounding box are present")
WHEN
(
"Tracing along a x-axis-aligned ray with a y-coordinate of origin beyond the bounding "
"box"
)
{
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
pi
/
2
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
volSize
);
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
pi
_t
/
2
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
volSize
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
...
...
@@ -704,8 +704,9 @@ SCENARIO("Rays not intersecting the bounding box are present")
const
index_t
numCases
=
9
;
real_t
alpha
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
};
real_t
beta
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
pi
/
2
,
pi
/
2
,
pi
/
2
};
real_t
gamma
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
pi
/
2
,
pi
/
2
,
pi
/
2
,
pi
/
2
,
pi
/
2
,
pi
/
2
};
real_t
beta
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
pi_t
/
2
,
pi_t
/
2
,
pi_t
/
2
};
real_t
gamma
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
pi_t
/
2
,
pi_t
/
2
,
pi_t
/
2
,
pi_t
/
2
,
pi_t
/
2
,
pi_t
/
2
};
real_t
offsetx
[
numCases
]
=
{
volSize
,
0.0
,
volSize
,
0.0
,
0.0
,
0.0
,
volSize
,
0.0
,
volSize
};
real_t
offsety
[
numCases
]
=
{
0.0
,
volSize
,
volSize
,
volSize
,
0.0
,
volSize
,
0.0
,
0.0
,
0.0
};
real_t
offsetz
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
volSize
,
volSize
,
0.0
,
volSize
,
volSize
};
...
...
@@ -761,7 +762,7 @@ SCENARIO("Axis-aligned rays are present")
std
::
vector
<
Geometry
>
geom
;
const
index_t
numCases
=
4
;
const
real_t
angles
[
numCases
]
=
{
0.0
,
pi
/
2
,
pi
,
3
*
pi
/
2
};
const
real_t
angles
[
numCases
]
=
{
0.0
,
pi
_t
/
2
,
pi
_t
,
3
*
pi
_t
/
2
};
RealVector_t
backProj
[
2
];
backProj
[
0
].
resize
(
volSize
*
volSize
);
backProj
[
1
].
resize
(
volSize
*
volSize
);
...
...
@@ -905,8 +906,8 @@ SCENARIO("Axis-aligned rays are present")
std
::
vector
<
Geometry
>
geom
;
const
index_t
numCases
=
6
;
real_t
beta
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
pi
/
2
,
3
*
pi
/
2
};
real_t
gamma
[
numCases
]
=
{
0.0
,
pi
,
pi
/
2
,
3
*
pi
/
2
,
pi
/
2
,
3
*
pi
/
2
};
real_t
beta
[
numCases
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
pi
_t
/
2
,
3
*
pi
_t
/
2
};
real_t
gamma
[
numCases
]
=
{
0.0
,
pi
_t
,
pi
_t
/
2
,
3
*
pi
_t
/
2
,
pi
_t
/
2
,
3
*
pi
_t
/
2
};
std
::
string
al
[
numCases
]
=
{
"z"
,
"-z"
,
"x"
,
"-x"
,
"y"
,
"-y"
};
RealVector_t
backProj
[
numCases
];
...
...
@@ -1099,11 +1100,11 @@ SCENARIO("Axis-aligned rays are present")
WHEN
(
"Both x- and y-axis-aligned rays are present"
)
{
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
0
,
volumeDescriptor
,
sinoDescriptor
);
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
90
*
pi
/
180.
,
volumeDescriptor
,
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
90
*
pi
_t
/
180.
,
volumeDescriptor
,
sinoDescriptor
);
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
180
*
pi
/
180.
,
volumeDescriptor
,
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
180
*
pi
_t
/
180.
,
volumeDescriptor
,
sinoDescriptor
);
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
270
*
pi
/
180.
,
volumeDescriptor
,
geom
.
emplace_back
(
20
*
volSize
,
volSize
,
270
*
pi
_t
/
180.
,
volumeDescriptor
,
sinoDescriptor
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
...
...
@@ -1157,8 +1158,8 @@ SCENARIO("Axis-aligned rays are present")
WHEN
(
"x-, y and z-axis-aligned rays are present"
)
{
real_t
beta
[
numImgs
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
pi
/
2
,
3
*
pi
/
2
};
real_t
gamma
[
numImgs
]
=
{
0.0
,
pi
,
pi
/
2
,
3
*
pi
/
2
,
pi
/
2
,
3
*
pi
/
2
};
real_t
beta
[
numImgs
]
=
{
0.0
,
0.0
,
0.0
,
0.0
,
pi
_t
/
2
,
3
*
pi
_t
/
2
};
real_t
gamma
[
numImgs
]
=
{
0.0
,
pi
_t
,
pi
_t
/
2
,
3
*
pi
_t
/
2
,
pi
_t
/
2
,
3
*
pi
_t
/
2
};
for
(
index_t
i
=
0
;
i
<
numImgs
;
i
++
)
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
gamma
[
i
],
...
...
@@ -1227,7 +1228,7 @@ SCENARIO("Projection under an angle")
{
// In this case the ray enters and exits the volume through the borders along the main
// direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
/
6
,
volumeDescriptor
,
sinoDescriptor
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
_t
/
6
,
volumeDescriptor
,
sinoDescriptor
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1299,8 +1300,8 @@ SCENARIO("Projection under an angle")
{
// In this case the ray exits through a border along the main ray direction, but enters
// through a border not along the main direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
/
6
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
sqrt
(
3
));
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
_t
/
6
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
sqrt
(
3
));
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1350,8 +1351,8 @@ SCENARIO("Projection under an angle")
{
// In this case the ray enters through a border along the main ray direction, but exits
// through a border not along the main direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
/
6
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
-
sqrt
(
3
));
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
_t
/
6
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
-
sqrt
(
3
));
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1399,8 +1400,8 @@ SCENARIO("Projection under an angle")
WHEN
(
"Projecting under an angle of 30 degrees and ray only intersects a single pixel"
)
{
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
/
6
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
-
2
-
sqrt
(
3
)
/
2
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
pi
_t
/
6
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
-
2
-
sqrt
(
3
)
/
2
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1438,7 +1439,8 @@ SCENARIO("Projection under an angle")
{
// In this case the ray enters and exits the volume through the borders along the main
// direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi
/
3
,
volumeDescriptor
,
sinoDescriptor
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi_t
/
3
,
volumeDescriptor
,
sinoDescriptor
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1507,8 +1509,8 @@ SCENARIO("Projection under an angle")
{
// In this case the ray exits through a border along the main ray direction, but enters
// through a border not along the main direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi
/
3
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
sqrt
(
3
));
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi
_t
/
3
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
sqrt
(
3
));
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1559,8 +1561,8 @@ SCENARIO("Projection under an angle")
{
// In this case the ray enters through a border along the main ray direction, but exits
// through a border not along the main direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi
/
3
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
-
sqrt
(
3
));
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi
_t
/
3
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
-
sqrt
(
3
));
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1610,8 +1612,8 @@ SCENARIO("Projection under an angle")
WHEN
(
"Projecting under an angle of 120 degrees and ray only intersects a single pixel"
)
{
// This is a special case that is handled separately in both forward and backprojection
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi
/
3
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
-
2
-
sqrt
(
3
)
/
2
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
-
2
*
pi
_t
/
3
,
volumeDescriptor
,
sinoDescriptor
,
0.0
,
0.0
,
-
2
-
sqrt
(
3
)
/
2
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"Ray intersects the correct pixels"
)
...
...
@@ -1665,7 +1667,7 @@ SCENARIO("Projection under an angle")
WHEN
(
"A ray with an angle of 30 degrees goes through the center of the volume"
)
{
// In this case the ray enters and exits the volume along the main direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
/
6
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
_t
/
6
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"The ray intersects the correct voxels"
)
...
...
@@ -1715,8 +1717,8 @@ SCENARIO("Projection under an angle")
WHEN
(
"A ray with an angle of 30 degrees enters through the right border"
)
{
// In this case the ray enters through a border orthogonal to a non-main direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
/
6
,
0.0
,
0.0
,
0.0
,
0.0
,
1
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
_t
/
6
,
0.0
,
0.0
,
0.0
,
0.0
,
1
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"The ray intersects the correct voxels"
)
...
...
@@ -1764,8 +1766,8 @@ SCENARIO("Projection under an angle")
WHEN
(
"A ray with an angle of 30 degrees exits through the left border"
)
{
// In this case the ray exit through a border orthogonal to a non-main direction
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
/
6
,
0.0
,
0.0
,
0.0
,
0.0
,
-
1
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
_t
/
6
,
0.0
,
0.0
,
0.0
,
0.0
,
-
1
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"The ray intersects the correct voxels"
)
...
...
@@ -1813,8 +1815,8 @@ SCENARIO("Projection under an angle")
WHEN
(
"A ray with an angle of 30 degrees only intersects a single voxel"
)
{
// special case - no interior voxels, entry and exit voxels are the same
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
/
6
,
0.0
,
0.0
,
0.0
,
0.0
,
-
2
);
geom
.
emplace_back
(
volSize
*
20
,
volSize
,
volumeDescriptor
,
sinoDescriptor
,
pi
_t
/
6
,
0.0
,
0.0
,
0.0
,
0.0
,
-
2
);
BinaryMethod
op
(
volumeDescriptor
,
sinoDescriptor
,
geom
);
THEN
(
"The ray intersects the correct voxels"
)
...
...
elsa/projectors/tests/test_Geometry.cpp
View file @
432c010b
...
...
@@ -65,7 +65,7 @@ SCENARIO("Testing 2D geometries")
THEN
(
"rays make sense"
)
{
// test outer + central pixels
for
(
real_t
detPixel
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
RealVector_t
pixel
(
1
);
pixel
<<
detPixel
;
auto
[
ro
,
rd
]
=
g
.
computeRayTo
(
pixel
);
...
...
@@ -121,7 +121,7 @@ SCENARIO("Testing 2D geometries")
THEN
(
"rays make sense"
)
{
// test outer + central pixels
for
(
real_t
detPixel
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
RealVector_t
pixel
(
1
);
pixel
<<
detPixel
;
auto
[
ro
,
rd
]
=
g
.
computeRayTo
(
pixel
);
...
...
@@ -139,7 +139,7 @@ SCENARIO("Testing 2D geometries")
WHEN
(
"testing geometry with 90 degree rotation and no offsets"
)
{
real_t
angle
=
pi
/
2
;
// 90 degrees
real_t
angle
=
pi
_t
/
2
;
// 90 degrees
Geometry
g
(
s2c
,
c2d
,
angle
,
ddVol
,
ddDet
);
THEN
(
"a copy is the same"
)
...
...
@@ -178,7 +178,7 @@ SCENARIO("Testing 2D geometries")
THEN
(
"rays make sense"
)
{
// test outer + central pixels
for
(
real_t
detPixel
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
RealVector_t
pixel
(
1
);
pixel
<<
detPixel
;
auto
[
ro
,
rd
]
=
g
.
computeRayTo
(
pixel
);
...
...
@@ -196,7 +196,7 @@ SCENARIO("Testing 2D geometries")
WHEN
(
"testing geometry with 45 degree rotation and offset center of rotation"
)
{
real_t
angle
=
pi
/
4
;
// 45 degrees
real_t
angle
=
pi
_t
/
4
;
// 45 degrees
real_t
cx
=
-
1
;
real_t
cy
=
2
;
Geometry
g
(
s2c
,
c2d
,
angle
,
ddVol
,
ddDet
,
0
,
cx
,
cy
);
...
...
@@ -242,7 +242,7 @@ SCENARIO("Testing 2D geometries")
THEN
(
"rays make sense"
)
{
// test outer + central pixels
for
(
real_t
detPixel
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
RealVector_t
pixel
(
1
);
pixel
<<
detPixel
;
auto
[
ro
,
rd
]
=
g
.
computeRayTo
(
pixel
);
...
...
@@ -321,8 +321,8 @@ SCENARIO("Testing 3D geometries")
THEN
(
"rays make sense"
)
{
// test outer + central pixels
for
(
real_t
detPixel1
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel2
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel1
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel2
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
RealVector_t
pixel
(
2
);
pixel
<<
detPixel1
,
detPixel2
;
auto
[
ro
,
rd
]
=
g
.
computeRayTo
(
pixel
);
...
...
@@ -392,8 +392,8 @@ SCENARIO("Testing 3D geometries")
THEN
(
"rays make sense"
)
{
// test outer + central pixels
for
(
real_t
detPixel1
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel2
:
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel1
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
for
(
real_t
detPixel2
:
std
::
initializer_list
<
real_t
>
{
0.5
,
2.5
,
4.5
})
{
RealVector_t
pixel
(
2
);
pixel
<<
detPixel1
,
detPixel2
;
auto
[
ro
,
rd
]
=
g
.
computeRayTo
(
pixel
);
...
...
@@ -423,7 +423,7 @@ SCENARIO("Testing 3D geometries")
WHEN
(
"testing geometry with 90 degree rotation and no offsets"
)
{