Skip to content

Commit

Permalink
Merge pull request #1208 from KrisThielemans/fixTOFalgebra
Browse files Browse the repository at this point in the history 
fixes STIR TOF AcquisitionData algebra
  • Loading branch information
@KrisThielemans
KrisThielemans authored Aug 22, 2023
2 parents 1bb1b11 + abaa910 commit 60bda8a
Show file tree
Hide file tree
Showing 3 changed files with 118 additions and 115 deletions.
87 changes: 31 additions & 56 deletions src/common/Utilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -597,93 +597,91 @@ def test_data_container_algebra(test, x, eps=1e-5):

s = x.norm()
t = numpy.linalg.norm(ax)
test.check_if_equal(1, abs(t - s) <= eps * abs(t))
# needs increased tolerance for large data size
test.check_if_equal_within_tolerance(t, s, 0, eps * 10);

s = x.max()
t = numpy.max(ax)
test.check_if_equal(1, abs(t - s) <= eps * abs(t))
test.check_if_equal_within_tolerance(t, s, 0, eps);

s = x.sum()
t = numpy.sum(ax)
r = numpy.sum(abs(ax))
test.check_if_equal(1, abs(t - s) <= eps * r)
test.check_if_equal_within_tolerance(t, s, 0, eps);

s = x.dot(y)
t = numpy.vdot(ay, ax)
test.check_if_equal(1, abs(t - s) <= eps * abs(t))
# needs increased tolerance for large data size
test.check_if_equal_within_tolerance(t, s, 0, eps * 10);

x2 = x.multiply(2)
ax2 = x2.as_array()
s = numpy.linalg.norm(ax2 - 2*ax)
t = numpy.linalg.norm(ax2)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

x2 *= 0
x.multiply(2, out=x2)
ax2 = x2.as_array()
s = numpy.linalg.norm(ax2 - 2*ax)
t = numpy.linalg.norm(ax2)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

t = x2.norm()
x2 -= x*2
s = x2.norm()
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

y = x.multiply(x)
ax = x.as_array()
ay = y.as_array()
s = numpy.linalg.norm(ay - ax * ax)
t = numpy.linalg.norm(ay)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

y *= 0
x.multiply(x, out=y)
ax = x.as_array()
ay = y.as_array()
s = numpy.linalg.norm(ay - ax * ax)
t = numpy.linalg.norm(ay)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

z = x*y
az = z.as_array()
s = numpy.linalg.norm(az - ax * ay)
t = numpy.linalg.norm(az)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

y = x + 1
ax = x.as_array()
ay = y.as_array()
s = numpy.linalg.norm(ay - (ax + 1))
t = numpy.linalg.norm(ay)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

y *= 0
x.add(1, out=y)
ax = x.as_array()
ay = y.as_array()
s = numpy.linalg.norm(ay - (ax + 1))
t = numpy.linalg.norm(ay)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

z = x/y
az = z.as_array()
s = numpy.linalg.norm(az - ax/ay)
t = numpy.linalg.norm(az)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

z = x/2
az = z.as_array()
s = numpy.linalg.norm(az - ax/2)
t = numpy.linalg.norm(az)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

z *= 0
x.divide(y, out=z)
az = z.as_array()
s = numpy.linalg.norm(az - ax/ay)
t = numpy.linalg.norm(az)
test.check_if_equal(1, abs(s) <= eps * abs(t))
test.check_if_zero_within_tolerance(s, eps * t)

y = x.sapyb(1, x, -1)
s = y.norm()
Expand Down Expand Up @@ -723,7 +721,6 @@ def test_data_container_algebra(test, x, eps=1e-5):
test.check_if_equal(0, s)

y = x.maximum(z)
ax = x.as_array()
ay = y.as_array()
az = z.as_array()
ay -= numpy.maximum(ax, az)
Expand All @@ -732,30 +729,26 @@ def test_data_container_algebra(test, x, eps=1e-5):

y *= 0
x.maximum(z, out=y)
ax = x.as_array()
ay = y.as_array()
az = z.as_array()
ay -= numpy.maximum(ax, az)
s = numpy.linalg.norm(ay)
test.check_if_equal(0, s)

y = x.maximum(0)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.maximum(ax, 0)
s = numpy.linalg.norm(ay)
test.check_if_equal(0, s)

y *= 0
x.maximum(0, out=y)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.maximum(ax, 0)
s = numpy.linalg.norm(ay)
test.check_if_equal(0, s)

y = x.minimum(z)
ax = x.as_array()
ay = y.as_array()
az = z.as_array()
ay -= numpy.minimum(ax, az)
Expand All @@ -764,166 +757,148 @@ def test_data_container_algebra(test, x, eps=1e-5):

y *= 0
x.minimum(z, out=y)
ax = x.as_array()
ay = y.as_array()
az = z.as_array()
ay -= numpy.minimum(ax, az)
s = numpy.linalg.norm(ay)
test.check_if_equal(0, s)

y = x.minimum(0)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.minimum(ax, 0)
s = numpy.linalg.norm(ay)
test.check_if_equal(0, s)

y *= 0
x.minimum(0, out=y)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.minimum(ax, 0)
s = numpy.linalg.norm(ay)
test.check_if_equal(0, s)

y = x.exp()
ax = x.as_array()
ay = y.as_array()
ay -= numpy.exp(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y *= 0
x.exp(out=y)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.exp(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y = x.log()
ax = x.as_array()
ay = y.as_array()
az = numpy.log(ax)
numpy.nan_to_num(ay, copy=False, posinf=0.0, neginf=0.0)
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
ay -= az
s = numpy.linalg.norm(ay)
t = numpy.linalg.norm(az)
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y *= 0
x.log(out=y)
ax = x.as_array()
ay = y.as_array()
az = numpy.log(ax)
numpy.nan_to_num(ay, copy=False, posinf=0.0, neginf=0.0)
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
ay -= az
s = numpy.linalg.norm(ay)
t = numpy.linalg.norm(az)
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y = x.sqrt()
ax = x.as_array()
ay = y.as_array()
ay -= numpy.sqrt(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y *= 0
x.sqrt(out=y)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.sqrt(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y = x.sign()
ax = x.as_array()
ay = y.as_array()
ay -= numpy.sign(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y *= 0
x.sign(out=y)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.sign(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y = x.abs()
ax = x.as_array()
ay = y.as_array()
ay -= numpy.abs(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

y *= 0
x.abs(out=y)
ax = x.as_array()
ay = y.as_array()
ay -= numpy.abs(ax)
s = numpy.linalg.norm(ay)
t = y.norm()
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

p = -0.5
z = x.power(p)
ax = x.as_array()
az = z.as_array()
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
t = numpy.linalg.norm(az)
az -= numpy.power(ax, p)
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
s = numpy.linalg.norm(az)
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

z *= 0
x.power(p, out=z)
ax = x.as_array()
az = z.as_array()
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
t = numpy.linalg.norm(az)
az -= numpy.power(ax, p)
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
s = numpy.linalg.norm(az)
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

ax = x.as_array()
ay = -numpy.ones_like(ax)/2
y.fill(ay)
z = x.power(y)
ax = x.as_array()
ay = y.as_array()
az = z.as_array()
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
t = numpy.linalg.norm(az)
az -= numpy.power(ax, ay)
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
s = numpy.linalg.norm(az)
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)

z *= 0
x.power(y, out=z)
ax = x.as_array()
ay = y.as_array()
az = z.as_array()
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
t = numpy.linalg.norm(az)
az -= numpy.power(ax, ay)
numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0)
s = numpy.linalg.norm(az)
test.check_if_equal(1, s <= eps * t)
test.check_if_zero_within_tolerance(s, eps * t)


class DataContainerAlgebraTests(object):
Expand Down
34 changes: 27 additions & 7 deletions src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,7 +249,11 @@ namespace sirf {
if (_is_empty != -1)
return _is_empty ? 0 : 1;
try {
get_segment_by_sinogram(0);
#ifdef STIR_TOF
get_segment_by_sinogram(0,0);
#else
get_segment_by_sinogram(0);
#endif
}
catch (...) {
_is_empty = 1;
Expand Down Expand Up @@ -404,16 +408,32 @@ namespace sirf {
{
return data()->get_max_segment_num();
}
#ifdef STIR_TOF
stir::SegmentBySinogram<float>
get_segment_by_sinogram(const int segment_num) const
get_segment_by_sinogram(const int segment_num, const int timing_pos_num) const
{
return data()->get_segment_by_sinogram(segment_num);
}
return data()->get_segment_by_sinogram(segment_num, timing_pos_num);
}
#else
stir::SegmentBySinogram<float>
get_empty_segment_by_sinogram(const int segment_num) const
get_segment_by_sinogram(const int segment_num) const
{
return data()->get_empty_segment_by_sinogram(segment_num);
}
return data()->get_segment_by_sinogram(segment_num);
}
#endif
#ifdef STIR_TOF
stir::SegmentBySinogram<float>
get_empty_segment_by_sinogram(const int segment_num, const int timing_pos_num) const
{
return data()->get_empty_segment_by_sinogram(segment_num, false, timing_pos_num);
}
#else
stir::SegmentBySinogram<float>
get_empty_segment_by_sinogram(const int segment_num) const
{
return data()->get_empty_segment_by_sinogram(segment_num);
}
#endif
void set_segment(const stir::SegmentBySinogram<float>& s)
{
if (data()->set_segment(s) != stir::Succeeded::yes)
Expand Down
Loading

0 comments on commit 60bda8a

Please sign in to comment.