Skip to content

Commit

Permalink
add TOF loop in computation of Hessian for projdata
Browse files Browse the repository at this point in the history 
Fixes #1305
  • Loading branch information
@KrisThielemans
KrisThielemans committed May 15, 2024
1 parent 5e4fd8c commit 7460c36
Show file tree
Hide file tree
Showing 2 changed files with 137 additions and 80 deletions.
136 changes: 83 additions & 53 deletions src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -912,6 +912,30 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::get_exam_info_up
return exam_info_uptr;
}

static std::vector<ViewgramIndices>
find_basic_viewgram_indices_in_subset(const ProjDataInfo& proj_data_info,
const DataSymmetriesForViewSegmentNumbers& symmetries,
const int min_segment_num,
const int max_segment_num,
const int subset_num,
const int num_subsets)
{
const std::vector<ViewSegmentNumbers> vs_nums_to_process = detail::find_basic_vs_nums_in_subset(
proj_data_info, symmetries, min_segment_num, max_segment_num, subset_num, num_subsets);

std::vector<ViewgramIndices> vg_idx_to_process;
for (auto vs_num : vs_nums_to_process)
{
for (int k = proj_data_info.get_min_tof_pos_num(); k <= proj_data_info.get_max_tof_pos_num(); ++k)
{
ViewgramIndices viewgram_idx = vs_num;
viewgram_idx.timing_pos_num() = k;
vg_idx_to_process.push_back(viewgram_idx);
}
}
return vg_idx_to_process;
}

template <typename TargetT>
Succeeded
PoissonLogLikelihoodWithLinearModelForMeanAndProjData<
Expand All @@ -938,35 +962,36 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<
this->get_projector_pair().get_forward_projector_sptr()->set_input(input);
this->get_projector_pair().get_back_projector_sptr()->start_accumulating_in_new_target();

const std::vector<ViewSegmentNumbers> vs_nums_to_process
= detail::find_basic_vs_nums_in_subset(*this->get_proj_data().get_proj_data_info_sptr(),
*symmetries_sptr,
-this->get_max_segment_num_to_process(),
this->get_max_segment_num_to_process(),
subset_num,
this->get_num_subsets());
const std::vector<ViewgramIndices> vg_idx_to_process
= find_basic_viewgram_indices_in_subset(*this->get_proj_data().get_proj_data_info_sptr(),
*symmetries_sptr,
-this->get_max_segment_num_to_process(),
this->get_max_segment_num_to_process(),
subset_num,
this->get_num_subsets());

info("Forward projecting input image.", 2);
#ifdef STIR_OPENMP
# pragma omp parallel for schedule(dynamic)
#endif
// note: older versions of openmp need an int as loop
for (int i = 0; i < static_cast<int>(vs_nums_to_process.size()); ++i)
for (int i = 0; i < static_cast<int>(vg_idx_to_process.size()); ++i)
{
const auto viewgram_idx = vg_idx_to_process[i];
{
#ifdef STIR_OPENMP
const int thread_num = omp_get_thread_num();
info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d") % thread_num % omp_get_num_threads()
% vs_nums_to_process[i].segment_num() % vs_nums_to_process[i].view_num(),
2);
const int thread_num = omp_get_thread_num();
const int num_threads = omp_get_num_threads();
#else
info(boost::format("calculating segment_num: %d, view_num: %d") % vs_nums_to_process[i].segment_num()
% vs_nums_to_process[i].view_num(),
2);
const int thread_num = 0;
const int num_threads = 1;
#endif
const ViewSegmentNumbers view_segment_num = vs_nums_to_process[i];

info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d, TOF: %d") % thread_num % num_threads
% viewgram_idx.segment_num() % viewgram_idx.view_num() % viewgram_idx.timing_pos_num(),
2);
}
// first compute data-term: y*norm^2
RelatedViewgrams<float> viewgrams = this->get_proj_data().get_related_viewgrams(view_segment_num, symmetries_sptr);
RelatedViewgrams<float> viewgrams = this->get_proj_data().get_related_viewgrams(viewgram_idx, symmetries_sptr);
// TODO add 1 for 1/(y+1) approximation

this->get_normalisation().apply(viewgrams);
Expand All @@ -978,7 +1003,7 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<
RelatedViewgrams<float> tmp_viewgrams;
// set tmp_viewgrams to geometric forward projection of input
{
tmp_viewgrams = this->get_proj_data().get_empty_related_viewgrams(view_segment_num, symmetries_sptr);
tmp_viewgrams = this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr);
this->get_projector_pair().get_forward_projector_sptr()->forward_project(tmp_viewgrams);
}

Expand Down Expand Up @@ -1046,43 +1071,46 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::actual_accumulat
this->get_projector_pair().get_forward_projector_sptr()->set_input(input);
this->get_projector_pair().get_back_projector_sptr()->start_accumulating_in_new_target();

const std::vector<ViewSegmentNumbers> vs_nums_to_process
= detail::find_basic_vs_nums_in_subset(*this->get_proj_data().get_proj_data_info_sptr(),
*symmetries_sptr,
-this->get_max_segment_num_to_process(),
this->get_max_segment_num_to_process(),
subset_num,
this->get_num_subsets());
const std::vector<ViewgramIndices> vg_idx_to_process
= find_basic_viewgram_indices_in_subset(*this->get_proj_data().get_proj_data_info_sptr(),
*symmetries_sptr,
-this->get_max_segment_num_to_process(),
this->get_max_segment_num_to_process(),
subset_num,
this->get_num_subsets());

// Create and populate the input_viewgrams_vec with empty values.
// This is needed to make the order of the vector correct w.r.t vs_nums_to_process.
// This is needed to make the order of the vector correct w.r.t vg_idx_to_process.
// OMP may mess this up
// Try: std::vector<RelatedViewgrams<float>> input_viewgrams_vec(vs_nums_to_process.size());
// Try: std::vector<RelatedViewgrams<float>> input_viewgrams_vec(vg_idx_to_process.size());
std::vector<RelatedViewgrams<float>> input_viewgrams_vec;
for (int i = 0; i < static_cast<int>(vs_nums_to_process.size()); ++i)
for (int i = 0; i < static_cast<int>(vg_idx_to_process.size()); ++i)
{
const ViewSegmentNumbers view_segment_num = vs_nums_to_process[i];
input_viewgrams_vec.push_back(this->get_proj_data().get_empty_related_viewgrams(view_segment_num, symmetries_sptr));
const auto viewgram_idx = vg_idx_to_process[i];
input_viewgrams_vec.push_back(this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr));
}

// Forward project input image
info("Forward projecting input image.", 2);
#ifdef STIR_OPENMP
# pragma omp parallel for schedule(dynamic)
#endif
for (int i = 0; i < static_cast<int>(vs_nums_to_process.size()); ++i)
{ // Loop over eah of the viewgrams in input_viewgrams_vec, forward projecting input into them
for (int i = 0; i < static_cast<int>(vg_idx_to_process.size()); ++i)
{ // Loop over each of the viewgrams in input_viewgrams_vec, forward projecting input into them
const auto viewgram_idx = vg_idx_to_process[i];
{
#ifdef STIR_OPENMP
const int thread_num = omp_get_thread_num();
info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d") % thread_num % omp_get_num_threads()
% vs_nums_to_process[i].segment_num() % vs_nums_to_process[i].view_num(),
2);
const int thread_num = omp_get_thread_num();
const int num_threads = omp_get_num_threads();
#else
info(boost::format("calculating segment_num: %d, view_num: %d") % vs_nums_to_process[i].segment_num()
% vs_nums_to_process[i].view_num(),
2);
const int thread_num = 0;
const int num_threads = 1;
#endif
input_viewgrams_vec[i] = this->get_proj_data().get_empty_related_viewgrams(vs_nums_to_process[i], symmetries_sptr);
info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d, TOF: %d") % thread_num % num_threads
% viewgram_idx.segment_num() % viewgram_idx.view_num() % viewgram_idx.timing_pos_num(),
2);
}
input_viewgrams_vec[i] = this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr);
this->get_projector_pair().get_forward_projector_sptr()->forward_project(input_viewgrams_vec[i]);
}

Expand All @@ -1091,35 +1119,37 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::actual_accumulat
#ifdef STIR_OPENMP
# pragma omp parallel for schedule(dynamic)
#endif
for (int i = 0; i < static_cast<int>(vs_nums_to_process.size()); ++i)
for (int i = 0; i < static_cast<int>(vg_idx_to_process.size()); ++i)
{
const auto viewgram_idx = vg_idx_to_process[i];
{
#ifdef STIR_OPENMP
const int thread_num = omp_get_thread_num();
info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d") % thread_num % omp_get_num_threads()
% vs_nums_to_process[i].segment_num() % vs_nums_to_process[i].view_num(),
2);
const int thread_num = omp_get_thread_num();
const int num_threads = omp_get_num_threads();
#else
info(boost::format("calculating segment_num: %d, view_num: %d") % vs_nums_to_process[i].segment_num()
% vs_nums_to_process[i].view_num(),
2);
const int thread_num = 0;
const int num_threads = 1;
#endif
info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d, TOF: %d") % thread_num % num_threads
% viewgram_idx.segment_num() % viewgram_idx.view_num() % viewgram_idx.timing_pos_num(),
2);
}
// Compute ybar_sq_viewgram = [ F(current_image_est) + additive ]^2
RelatedViewgrams<float> ybar_sq_viewgram;
{
ybar_sq_viewgram = this->get_proj_data().get_empty_related_viewgrams(vs_nums_to_process[i], symmetries_sptr);
ybar_sq_viewgram = this->get_proj_data().get_empty_related_viewgrams(vg_idx_to_process[i], symmetries_sptr);
this->get_projector_pair().get_forward_projector_sptr()->forward_project(ybar_sq_viewgram);

// add additive sinogram to forward projection
if (!(is_null_ptr(this->get_additive_proj_data_sptr())))
ybar_sq_viewgram += this->get_additive_proj_data().get_related_viewgrams(vs_nums_to_process[i], symmetries_sptr);
ybar_sq_viewgram += this->get_additive_proj_data().get_related_viewgrams(vg_idx_to_process[i], symmetries_sptr);
// square ybar
ybar_sq_viewgram *= ybar_sq_viewgram;
}

// Compute: final_viewgram * F(input) / ybar_sq_viewgram
// final_viewgram starts as measured data
RelatedViewgrams<float> final_viewgram
= this->get_proj_data().get_related_viewgrams(vs_nums_to_process[i], symmetries_sptr);
RelatedViewgrams<float> final_viewgram = this->get_proj_data().get_related_viewgrams(vg_idx_to_process[i], symmetries_sptr);
{
// Mult input_viewgram
final_viewgram *= input_viewgrams_vec[i];
Expand Down
81 changes: 54 additions & 27 deletions src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests
*/
PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests(char const* const proj_data_filename = 0,
char const* const density_filename = 0);
void construct_input_data(shared_ptr<target_type>& density_sptr);
void construct_input_data(shared_ptr<target_type>& density_sptr, const bool TOF_or_not);

void run_tests() override;

Expand Down Expand Up @@ -174,18 +174,37 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::test_approximate_Hes
}

void
PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data(shared_ptr<target_type>& density_sptr)
PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data(shared_ptr<target_type>& density_sptr,
const bool TOF_or_not)
{
if (this->proj_data_filename == 0)
{
shared_ptr<ProjDataInfo> proj_data_info_sptr;
// construct a small scanner and sinogram
shared_ptr<Scanner> scanner_sptr(new Scanner(Scanner::E953));
scanner_sptr->set_num_rings(5);
shared_ptr<ProjDataInfo> proj_data_info_sptr(ProjDataInfo::ProjDataInfoCTI(scanner_sptr,
if (TOF_or_not)
{
std::cerr << "------ TOF data ----\n";
shared_ptr<Scanner> scanner_sptr(new Scanner(Scanner::Discovery690));
scanner_sptr->set_num_rings(4);
proj_data_info_sptr = std::move(ProjDataInfo::construct_proj_data_info(scanner_sptr,
/*span=*/3,
/*max_delta=*/4,
/*max_delta=*/2,
/*num_views=*/16,
/*num_tang_poss=*/16));
/*num_tang_poss=*/16,
/* arccorrected=*/false,
/* TOF_mash_factor=*/11));
}
else
{
std::cerr << "------ non-TOF data ----\n";
shared_ptr<Scanner> scanner_sptr(new Scanner(Scanner::E953));
scanner_sptr->set_num_rings(5);
proj_data_info_sptr.reset(ProjDataInfo::ProjDataInfoCTI(scanner_sptr,
/*span=*/3,
/*max_delta=*/4,
/*num_views=*/16,
/*num_tang_poss=*/16));
}
shared_ptr<ExamInfo> exam_info_sptr(new ExamInfo(ImagingModality::PT));
proj_data_sptr.reset(new ProjDataInMemory(exam_info_sptr, proj_data_info_sptr));
for (int seg_num = proj_data_sptr->get_min_segment_num(); seg_num <= proj_data_sptr->get_max_segment_num(); ++seg_num)
Expand Down Expand Up @@ -251,24 +270,23 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data
// multiplicative term
shared_ptr<BinNormalisation> bin_norm_sptr(new TrivialBinNormalisation());
{

mult_proj_data_sptr.reset(new ProjDataInMemory(proj_data_sptr->get_exam_info_sptr(),
proj_data_sptr->get_proj_data_info_sptr()->create_shared_clone()));
proj_data_sptr->get_proj_data_info_sptr()->create_non_tof_clone()));
for (int seg_num = proj_data_sptr->get_min_segment_num(); seg_num <= proj_data_sptr->get_max_segment_num(); ++seg_num)
{
for (int timing_pos_num = proj_data_sptr->get_min_tof_pos_num(); timing_pos_num <= proj_data_sptr->get_max_tof_pos_num();
++timing_pos_num)
{
SegmentByView<float> segment = proj_data_sptr->get_empty_segment_by_view(seg_num, false, timing_pos_num);
// fill in some crazy values
float value = 0;
for (SegmentByView<float>::full_iterator iter = segment.begin_all(); iter != segment.end_all(); ++iter)
{
value = float(fabs(seg_num * value - .2)); // needs to be positive for Poisson
*iter = value;
}
segment /= 0.5F * segment.find_max(); // normalise to 2 (to avoid large numbers)
mult_proj_data_sptr->set_segment(segment);
}
{
auto segment = mult_proj_data_sptr->get_empty_segment_by_view(seg_num);
// fill in some crazy values
float value = 0;
for (auto iter = segment.begin_all(); iter != segment.end_all(); ++iter)
{
value = float(fabs(seg_num * value - .2)); // needs to be positive for Poisson
*iter = value;
}
segment /= 0.5F * segment.find_max(); // normalise to 2 (to avoid large numbers)
mult_proj_data_sptr->set_segment(segment);
}
}
bin_norm_sptr.reset(new BinNormalisationFromProjData(mult_proj_data_sptr));
}
Expand Down Expand Up @@ -305,8 +323,8 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data
shared_ptr<ProjectorByBinPair> proj_pair_sptr(new ProjectorByBinPairUsingProjMatrixByBin(proj_matrix_sptr));
objective_function.set_projector_pair_sptr(proj_pair_sptr);
/*
void set_frame_num(const int);
void set_frame_definitions(const TimeFrameDefinitions&);
void set_frame_num(const int);
void set_frame_definitions(const TimeFrameDefinitions&);
*/
objective_function.set_normalisation_sptr(bin_norm_sptr);
objective_function.set_additive_proj_data_sptr(add_proj_data_sptr);
Expand All @@ -323,9 +341,18 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::run_tests()
Verbosity::set(0);

#if 1
shared_ptr<target_type> density_sptr;
construct_input_data(density_sptr);
this->run_tests_for_objective_function(*this->objective_function_sptr, *density_sptr);
{
shared_ptr<target_type> density_sptr;
construct_input_data(density_sptr, /*TOF_or_not=*/false);
this->run_tests_for_objective_function(*this->objective_function_sptr, *density_sptr);
}
if (this->proj_data_filename == 0)
{
shared_ptr<target_type> density_sptr;
construct_input_data(density_sptr, /*TOF_or_not=*/true);
this->run_tests_for_objective_function(*this->objective_function_sptr, *density_sptr);
}

#else
// alternative that gets the objective function from an OSMAPOSL .par file
// currently disabled
Expand Down

0 comments on commit 7460c36

Please sign in to comment.