Skip to content

Commit

Permalink
fixes STIR TOF AcquisitionData algebra
Browse files Browse the repository at this point in the history 
inserted TOF loops. However, this needs the STIR_TOF preprocessor
variable to be defined, which is still to be merged on the STIR TOF branch.
  • Loading branch information
@KrisThielemans
KrisThielemans committed Aug 21, 2023
1 parent 1bb1b11 commit 70f1eba
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 42 deletions.
25 changes: 21 additions & 4 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 @@ -405,10 +409,23 @@ namespace sirf {
return data()->get_max_segment_num();
}
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);
}
#ifdef STIR_TOF
return data()->get_segment_by_sinogram(segment_num, timing_pos_num);
#else
return data()->get_segment_by_sinogram(segment_num);
#endif
}
stir::SegmentBySinogram<float>
get_empty_segment_by_sinogram(const int segment_num, const int timing_pos_num ) const
{
#ifdef STIR_TOF
return data()->get_empty_segment_by_sinogram(segment_num, false, timing_pos_num);
#else
return data()->get_empty_segment_by_sinogram(segment_num);
#endif
}
stir::SegmentBySinogram<float>
get_empty_segment_by_sinogram(const int segment_num) const
{
Expand Down
93 changes: 55 additions & 38 deletions src/xSTIR/cSTIR/stir_data_containers.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
/*
SyneRBI Synergistic Image Reconstruction Framework (SIRF)
Copyright 2017 - 2019 Rutherford Appleton Laboratory STFC
Copyright 2018 - 2020 University College London
Copyright 2017 - 2023 Rutherford Appleton Laboratory STFC
Copyright 2018 - 2023 University College London
This is software developed for the Collaborative Computational
Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
Expand Down Expand Up @@ -32,20 +32,29 @@ using namespace sirf;
std::string STIRAcquisitionData::_storage_scheme;
std::shared_ptr<STIRAcquisitionData> STIRAcquisitionData::_template;

#ifdef STIR_TOF
#define TOF_LOOP for (int k=data()->get_min_tof_pos_num(); k<=data()->get_max_tof_pos_num(); ++k)
#define TOF_ARG , k
#else
#define TOF_LOOP
#define TOF_ARG
#endif

float
STIRAcquisitionData::norm() const
{
double t = 0.0;
TOF_LOOP
for (int s = 0; s <= get_max_segment_num(); ++s)
{
SegmentBySinogram<float> seg = get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) {
double r = *seg_iter++;
t += r*r;
}
if (s != 0) {
seg = get_segment_by_sinogram(-s);
seg = get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) {
double r = *seg_iter++;
t += r*r;
Expand All @@ -60,14 +69,15 @@ STIRAcquisitionData::sum(void* ptr) const
{
int n = get_max_segment_num();
double t = 0;
TOF_LOOP
for (int s = 0; s <= n; ++s)
{
SegmentBySinogram<float> seg = get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
t += *seg_iter++;
if (s != 0) {
seg = get_segment_by_sinogram(-s);
seg = get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
t += *seg_iter++;
}
Expand All @@ -81,14 +91,15 @@ STIRAcquisitionData::max(void* ptr) const
{
int n = get_max_segment_num();
float t = 0;
TOF_LOOP
for (int s = 0; s <= n; ++s)
{
SegmentBySinogram<float> seg = get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
t = std::max(t, *seg_iter++);
if (s != 0) {
seg = get_segment_by_sinogram(-s);
seg = get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
t = std::max(t, *seg_iter++);
}
Expand All @@ -104,10 +115,11 @@ STIRAcquisitionData::dot(const DataContainer& a_x, void* ptr) const
int n = get_max_segment_num();
int nx = x.get_max_segment_num();
double t = 0;
TOF_LOOP
for (int s = 0; s <= n && s <= nx; ++s)
{
SegmentBySinogram<float> seg = get_segment_by_sinogram(s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
SegmentBySinogram<float>::full_iterator sx_iter;
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
Expand All @@ -116,8 +128,8 @@ STIRAcquisitionData::dot(const DataContainer& a_x, void* ptr) const
t += (*seg_iter++) * double(*sx_iter++);
}
if (s != 0) {
seg = get_segment_by_sinogram(-s);
sx = x.get_segment_by_sinogram(-s);
seg = get_segment_by_sinogram(-s TOF_ARG);
sx = x.get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
seg_iter != seg.end_all() && sx_iter != sx.end_all();
/*empty*/)
Expand Down Expand Up @@ -188,11 +200,12 @@ STIRAcquisitionData::inv(float amin, const DataContainer& a_x)
SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x);
int n = get_max_segment_num();
int nx = x.get_max_segment_num();
TOF_LOOP
for (int s = 0; s <= n && s <= nx; ++s)
{
//std::cout << "processing segment " << s << std::endl;
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
SegmentBySinogram<float>::full_iterator sx_iter;
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
Expand All @@ -202,8 +215,8 @@ STIRAcquisitionData::inv(float amin, const DataContainer& a_x)
set_segment(seg);
if (s != 0) {
//std::cout << "processing segment " << -s << std::endl;
seg = get_empty_segment_by_sinogram(-s);
sx = x.get_segment_by_sinogram(-s);
seg = get_empty_segment_by_sinogram(-s TOF_ARG);
sx = x.get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
seg_iter != seg.end_all() && sx_iter != sx.end_all();
/*empty*/) {
Expand All @@ -223,18 +236,19 @@ STIRAcquisitionData::unary_op(
SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x);
int n = get_max_segment_num();
int nx = x.get_max_segment_num();
TOF_LOOP
for (int s = 0; s <= n && s <= nx; ++s) {
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
SegmentBySinogram<float>::full_iterator sx_iter;
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
seg_iter != seg.end_all() && sx_iter != sx.end_all(); /*empty*/)
*seg_iter++ = f(*sx_iter++);
set_segment(seg);
if (s > 0) {
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(-s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(-s);
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(-s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(-s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
SegmentBySinogram<float>::full_iterator sx_iter;
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
Expand All @@ -256,18 +270,19 @@ STIRAcquisitionData::semibinary_op(
SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x);
int n = get_max_segment_num();
int nx = x.get_max_segment_num();
TOF_LOOP
for (int s = 0; s <= n && s <= nx; ++s) {
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
SegmentBySinogram<float>::full_iterator sx_iter;
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
seg_iter != seg.end_all() && sx_iter != sx.end_all(); /*empty*/)
*seg_iter++ = f(*sx_iter++, y);
set_segment(seg);
if (s > 0) {
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(-s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(-s);
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(-s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(-s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
SegmentBySinogram<float>::full_iterator sx_iter;
for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all();
Expand Down Expand Up @@ -296,11 +311,12 @@ STIRAcquisitionData::binary_op(
SegmentBySinogram<float>::full_iterator seg_iter;
SegmentBySinogram<float>::full_iterator sx_iter;
SegmentBySinogram<float>::full_iterator sy_iter;
TOF_LOOP
for (int s = 0; s <= n; ++s)
{
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s);
SegmentBySinogram<float> sy = y.get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sy = y.get_segment_by_sinogram(s TOF_ARG);
if (seg.size_all() != sx.size_all() || seg.size_all() != sy.size_all())
throw std::runtime_error("binary_op error: operands sizes differ");
for (seg_iter = seg.begin_all(),
Expand All @@ -309,9 +325,9 @@ STIRAcquisitionData::binary_op(
*seg_iter++ = f(*sx_iter++, *sy_iter++);
set_segment(seg);
if (s != 0) {
seg = get_empty_segment_by_sinogram(-s);
sx = x.get_segment_by_sinogram(-s);
sy = y.get_segment_by_sinogram(-s);
seg = get_empty_segment_by_sinogram(-s TOF_ARG);
sx = x.get_segment_by_sinogram(-s TOF_ARG);
sy = y.get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(),
sx_iter = sx.begin_all(), sy_iter = sy.begin_all();
seg_iter != seg.end_all(); /*empty*/)
Expand Down Expand Up @@ -340,12 +356,13 @@ STIRAcquisitionData::xapyb(
SegmentBySinogram<float>::full_iterator sx_iter;
SegmentBySinogram<float>::full_iterator sy_iter;
SegmentBySinogram<float>::full_iterator sb_iter;
TOF_LOOP
for (int s = 0; s <= n; ++s)
{
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s);
SegmentBySinogram<float> sy = y.get_segment_by_sinogram(s);
SegmentBySinogram<float> sb = b.get_segment_by_sinogram(s);
SegmentBySinogram<float> seg = get_empty_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sx = x.get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sy = y.get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float> sb = b.get_segment_by_sinogram(s TOF_ARG);
size_t size_seg = seg.size_all();
size_t size_sx = sx.size_all();
size_t size_sy = sy.size_all();
Expand All @@ -358,10 +375,10 @@ STIRAcquisitionData::xapyb(
*seg_iter++ = (*sx_iter++) * a + (*sy_iter++) * (*sb_iter++);
set_segment(seg);
if (s != 0) {
seg = get_empty_segment_by_sinogram(-s);
sx = x.get_segment_by_sinogram(-s);
sy = y.get_segment_by_sinogram(-s);
sb = b.get_segment_by_sinogram(-s);
seg = get_empty_segment_by_sinogram(-s TOF_ARG);
sx = x.get_segment_by_sinogram(-s TOF_ARG);
sy = y.get_segment_by_sinogram(-s TOF_ARG);
sb = b.get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(),
sx_iter = sx.begin_all(), sy_iter = sy.begin_all(), sb_iter = sb.begin_all();
seg_iter != seg.end_all(); /*empty*/)
Expand Down

0 comments on commit 70f1eba

Please sign in to comment.