From 76f1ba1055ac916a30c3808dc68fbd7c59a93089 Mon Sep 17 00:00:00 2001 From: cgf Date: Wed, 14 Oct 2020 16:04:51 +0800 Subject: [PATCH] add virtual crystal support for find_ML_normfactors3D.cxx and apply_normfactors3D.cxx --- src/buildblock/ML_norm.cxx | 2848 ++++++++++++----------- src/utilities/apply_normfactors3D.cxx | 219 +- src/utilities/find_ML_normfactors3D.cxx | 54 +- 3 files changed, 1588 insertions(+), 1533 deletions(-) diff --git a/src/buildblock/ML_norm.cxx b/src/buildblock/ML_norm.cxx index 2f1d265a3..b9ee6c496 100644 --- a/src/buildblock/ML_norm.cxx +++ b/src/buildblock/ML_norm.cxx @@ -38,464 +38,464 @@ using std::max; START_NAMESPACE_STIR -DetPairData::DetPairData() -{} + DetPairData::DetPairData() + {} -DetPairData::DetPairData(const IndexRange<2>& range) -:base_type(range), num_detectors(range.get_length()) -{ -} + DetPairData::DetPairData(const IndexRange<2>& range) + :base_type(range), num_detectors(range.get_length()) + { + } -DetPairData& -DetPairData::operator=(const DetPairData& other) -{ - base_type::operator=(other); - num_detectors = other.num_detectors; - return *this; -} + DetPairData& + DetPairData::operator=(const DetPairData& other) + { + base_type::operator=(other); + num_detectors = other.num_detectors; + return *this; + } -float & DetPairData::operator()(const int a, const int b) -{ - return (*this)[a][b=get_min_index(a)) - return b<=get_max_index(a); - else - return b+num_detectors<=get_max_index(a); -} + bool + DetPairData:: + is_in_data(const int a, const int b) const + { + if (b>=get_min_index(a)) + return b<=get_max_index(a); + else + return b+num_detectors<=get_max_index(a); + } -void DetPairData::fill(const float d) -{ - base_type::fill(d); -} + void DetPairData::fill(const float d) + { + base_type::fill(d); + } -void DetPairData::grow(const IndexRange<2>& range) -{ - base_type::grow(range); - num_detectors=range.get_length(); -} + void DetPairData::grow(const IndexRange<2>& range) + { + base_type::grow(range); + num_detectors=range.get_length(); + } -int DetPairData::get_min_index() const -{ - return base_type::get_min_index(); -} + int DetPairData::get_min_index() const + { + return base_type::get_min_index(); + } -int DetPairData::get_max_index() const -{ - return base_type::get_max_index(); -} + int DetPairData::get_max_index() const + { + return base_type::get_max_index(); + } -int DetPairData::get_min_index(const int a) const -{ - return (*this)[a].get_min_index(); -} + int DetPairData::get_min_index(const int a) const + { + return (*this)[a].get_min_index(); + } -int DetPairData::get_max_index(const int a) const -{ - return (*this)[a].get_max_index(); -} + int DetPairData::get_max_index(const int a) const + { + return (*this)[a].get_max_index(); + } -float DetPairData::sum() const -{ - return base_type::sum(); -} + float DetPairData::sum() const + { + return base_type::sum(); + } -float DetPairData::sum(const int a) const -{ - return (*this)[a].sum(); -} + float DetPairData::sum(const int a) const + { + return (*this)[a].sum(); + } -float DetPairData::find_max() const -{ - return base_type::find_max(); -} + float DetPairData::find_max() const + { + return base_type::find_max(); + } -float DetPairData::find_min() const -{ - return base_type::find_min(); -} + float DetPairData::find_min() const + { + return base_type::find_min(); + } -int DetPairData::get_num_detectors() const -{ - return num_detectors; -} + int DetPairData::get_num_detectors() const + { + return num_detectors; + } -void make_det_pair_data(DetPairData& det_pair_data, - const ProjDataInfo& proj_data_info_general_type, - const int segment_num, - const int ax_pos_num) -{ - const ProjDataInfoCylindricalNoArcCorr& proj_data_info = - dynamic_cast(proj_data_info_general_type); - - const int num_detectors = - proj_data_info.get_scanner_ptr()->get_num_detectors_per_ring(); - const int fan_size = - 2*max(proj_data_info.get_max_tangential_pos_num(), - -proj_data_info.get_min_tangential_pos_num()) + 1; - // fan will range from -half_fan_size to +half_fan_size (i.e. an odd number of elements) - const int half_fan_size = fan_size/2; - - IndexRange<2> fan_indices; - fan_indices.grow(0,num_detectors-1); - for (int a = 0; a < num_detectors; ++a) - { - fan_indices[a] = - IndexRange<1>(a+num_detectors/2-half_fan_size, - a+num_detectors/2+half_fan_size); - } - det_pair_data.grow(fan_indices); - det_pair_data.fill(0); -} + void make_det_pair_data(DetPairData& det_pair_data, + const ProjDataInfo& proj_data_info_general_type, + const int segment_num, + const int ax_pos_num) + { + const ProjDataInfoCylindricalNoArcCorr& proj_data_info = + dynamic_cast(proj_data_info_general_type); + + const int num_detectors = + proj_data_info.get_scanner_ptr()->get_num_detectors_per_ring(); + const int fan_size = + 2*max(proj_data_info.get_max_tangential_pos_num(), + -proj_data_info.get_min_tangential_pos_num()) + 1; + // fan will range from -half_fan_size to +half_fan_size (i.e. an odd number of elements) + const int half_fan_size = fan_size/2; + + IndexRange<2> fan_indices; + fan_indices.grow(0,num_detectors-1); + for (int a = 0; a < num_detectors; ++a) + { + fan_indices[a] = + IndexRange<1>(a+num_detectors/2-half_fan_size, + a+num_detectors/2+half_fan_size); + } + det_pair_data.grow(fan_indices); + det_pair_data.fill(0); + } -void make_det_pair_data(DetPairData& det_pair_data, - const ProjData& proj_data, - const int segment_num, - const int ax_pos_num) -{ - make_det_pair_data(det_pair_data, - *proj_data.get_proj_data_info_sptr(), - segment_num, - ax_pos_num); - const int num_detectors = - det_pair_data.get_num_detectors(); - const ProjDataInfoCylindricalNoArcCorr& proj_data_info = - dynamic_cast(*proj_data.get_proj_data_info_sptr()); - - shared_ptr > - pos_sino_ptr(new Sinogram(proj_data.get_sinogram(ax_pos_num,segment_num))); - shared_ptr > neg_sino_ptr; - if (segment_num == 0) - neg_sino_ptr = pos_sino_ptr; - else - neg_sino_ptr. - reset(new Sinogram(proj_data.get_sinogram(ax_pos_num,-segment_num))); - - - for (int view_num = 0; view_num < num_detectors/2; view_num++) - for (int tang_pos_num = proj_data.get_min_tangential_pos_num(); - tang_pos_num <= proj_data.get_max_tangential_pos_num(); - ++tang_pos_num) - { - int det_num_a = 0; - int det_num_b = 0; + void make_det_pair_data(DetPairData& det_pair_data, + const ProjData& proj_data, + const int segment_num, + const int ax_pos_num) + { + make_det_pair_data(det_pair_data, + *proj_data.get_proj_data_info_sptr(), + segment_num, + ax_pos_num); + const int num_detectors = + det_pair_data.get_num_detectors(); + const ProjDataInfoCylindricalNoArcCorr& proj_data_info = + dynamic_cast(*proj_data.get_proj_data_info_sptr()); + + shared_ptr > + pos_sino_ptr(new Sinogram(proj_data.get_sinogram(ax_pos_num,segment_num))); + shared_ptr > neg_sino_ptr; + if (segment_num == 0) + neg_sino_ptr = pos_sino_ptr; + else + neg_sino_ptr. + reset(new Sinogram(proj_data.get_sinogram(ax_pos_num,-segment_num))); - proj_data_info.get_det_num_pair_for_view_tangential_pos_num(det_num_a, det_num_b, view_num, tang_pos_num); - det_pair_data(det_num_a,det_num_b) = - (*pos_sino_ptr)[view_num][tang_pos_num]; - det_pair_data(det_num_b,det_num_a) = - (*neg_sino_ptr)[view_num][tang_pos_num]; - } -} + for (int view_num = 0; view_num < num_detectors/2; view_num++) + for (int tang_pos_num = proj_data.get_min_tangential_pos_num(); + tang_pos_num <= proj_data.get_max_tangential_pos_num(); + ++tang_pos_num) + { + int det_num_a = 0; + int det_num_b = 0; -void set_det_pair_data(ProjData& proj_data, - const DetPairData& det_pair_data, - const int segment_num, - const int ax_pos_num) -{ - const shared_ptr proj_data_info_sptr = - proj_data.get_proj_data_info_sptr(); - const ProjDataInfoCylindricalNoArcCorr& proj_data_info = - dynamic_cast(*proj_data_info_sptr); - - const int num_detectors = det_pair_data.get_num_detectors(); - assert(proj_data_info.get_scanner_ptr()->get_num_detectors_per_ring() == num_detectors); - - shared_ptr > - pos_sino_ptr(new Sinogram(proj_data.get_empty_sinogram(ax_pos_num,segment_num))); - shared_ptr > neg_sino_ptr; - if (segment_num != 0) - neg_sino_ptr. - reset(new Sinogram(proj_data.get_empty_sinogram(ax_pos_num,-segment_num))); - - - for (int view_num = 0; view_num < num_detectors/2; view_num++) - for (int tang_pos_num = proj_data.get_min_tangential_pos_num(); - tang_pos_num <= proj_data.get_max_tangential_pos_num(); - ++tang_pos_num) - { - int det_num_a = 0; - int det_num_b = 0; + proj_data_info.get_det_num_pair_for_view_tangential_pos_num(det_num_a, det_num_b, view_num, tang_pos_num); - proj_data_info.get_det_num_pair_for_view_tangential_pos_num(det_num_a, det_num_b, view_num, tang_pos_num); + det_pair_data(det_num_a,det_num_b) = + (*pos_sino_ptr)[view_num][tang_pos_num]; + det_pair_data(det_num_b,det_num_a) = + (*neg_sino_ptr)[view_num][tang_pos_num]; + } + } - (*pos_sino_ptr)[view_num][tang_pos_num] = - det_pair_data(det_num_a,det_num_b); - if (segment_num!=0) - (*neg_sino_ptr)[view_num][tang_pos_num] = - det_pair_data(det_num_b,det_num_a); - } - proj_data.set_sinogram(*pos_sino_ptr); - if (segment_num != 0) - proj_data.set_sinogram(*neg_sino_ptr); -} + void set_det_pair_data(ProjData& proj_data, + const DetPairData& det_pair_data, + const int segment_num, + const int ax_pos_num) + { + const shared_ptr proj_data_info_sptr = + proj_data.get_proj_data_info_sptr(); + const ProjDataInfoCylindricalNoArcCorr& proj_data_info = + dynamic_cast(*proj_data_info_sptr); + + const int num_detectors = det_pair_data.get_num_detectors(); + assert(proj_data_info.get_scanner_ptr()->get_num_detectors_per_ring() == num_detectors); + + shared_ptr > + pos_sino_ptr(new Sinogram(proj_data.get_empty_sinogram(ax_pos_num,segment_num))); + shared_ptr > neg_sino_ptr; + if (segment_num != 0) + neg_sino_ptr. + reset(new Sinogram(proj_data.get_empty_sinogram(ax_pos_num,-segment_num))); + + + for (int view_num = 0; view_num < num_detectors/2; view_num++) + for (int tang_pos_num = proj_data.get_min_tangential_pos_num(); + tang_pos_num <= proj_data.get_max_tangential_pos_num(); + ++tang_pos_num) + { + int det_num_a = 0; + int det_num_b = 0; + proj_data_info.get_det_num_pair_for_view_tangential_pos_num(det_num_a, det_num_b, view_num, tang_pos_num); -void apply_block_norm(DetPairData& det_pair_data, const BlockData& block_data, const bool apply) -{ - const int num_detectors = det_pair_data.get_num_detectors(); - const int num_blocks = block_data.get_length(); - const int num_crystals_per_block = num_detectors/num_blocks; - assert(num_blocks * num_crystals_per_block == num_detectors); - - for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) - for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) - { - // note: add 2*num_detectors to newb to avoid using mod with negative numbers - if (det_pair_data(a,b) == 0) - continue; - if (apply) - det_pair_data(a,b) *= - block_data[a/num_crystals_per_block][(b/num_crystals_per_block)%num_blocks]; - else - det_pair_data(a,b) /= - block_data[a/num_crystals_per_block][(b/num_crystals_per_block)%num_blocks]; - } -} + (*pos_sino_ptr)[view_num][tang_pos_num] = + det_pair_data(det_num_a,det_num_b); + if (segment_num!=0) + (*neg_sino_ptr)[view_num][tang_pos_num] = + det_pair_data(det_num_b,det_num_a); + } + proj_data.set_sinogram(*pos_sino_ptr); + if (segment_num != 0) + proj_data.set_sinogram(*neg_sino_ptr); + } -void apply_geo_norm(DetPairData& det_pair_data, const GeoData& geo_data, const bool apply) -{ - const int num_detectors = det_pair_data.get_num_detectors(); - const int num_crystals_per_block = geo_data.get_length()*2; - for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) - for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) - { - if (det_pair_data(a,b) == 0) - continue; - int newa = a % num_crystals_per_block; - int newb = b - (a - newa); - if (newa > num_crystals_per_block - 1 - newa) - { - newa = num_crystals_per_block - 1 - newa; - newb = - newb + num_crystals_per_block - 1; - } - // note: add 2*num_detectors to newb to avoid using mod with negative numbers - if (apply) - det_pair_data(a,b) *= - geo_data[newa][(2*num_detectors + newb)%num_detectors]; - else - det_pair_data(a,b) /= - geo_data[newa][(2*num_detectors + newb)%num_detectors]; - } -} + void apply_block_norm(DetPairData& det_pair_data, const BlockData& block_data, const bool apply) + { + const int num_detectors = det_pair_data.get_num_detectors(); + const int num_blocks = block_data.get_length(); + const int num_crystals_per_block = num_detectors/num_blocks; + assert(num_blocks * num_crystals_per_block == num_detectors); -void apply_efficiencies(DetPairData& det_pair_data, const Array<1,float>& efficiencies, const bool apply) -{ - const int num_detectors = det_pair_data.get_num_detectors(); - for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) - for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) - { - if (det_pair_data(a,b) == 0) - continue; - if (apply) - det_pair_data(a,b) *= - efficiencies[a]*efficiencies[b%num_detectors]; - else - det_pair_data(a,b) /= - efficiencies[a]*efficiencies[b%num_detectors]; - } -} + for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) + for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) + { + // note: add 2*num_detectors to newb to avoid using mod with negative numbers + if (det_pair_data(a,b) == 0) + continue; + if (apply) + det_pair_data(a,b) *= + block_data[a/num_crystals_per_block][(b/num_crystals_per_block)%num_blocks]; + else + det_pair_data(a,b) /= + block_data[a/num_crystals_per_block][(b/num_crystals_per_block)%num_blocks]; + } + } + void apply_geo_norm(DetPairData& det_pair_data, const GeoData& geo_data, const bool apply) + { + const int num_detectors = det_pair_data.get_num_detectors(); + const int num_crystals_per_block = geo_data.get_length()*2; -void make_fan_sum_data(Array<1,float>& data_fan_sums, const DetPairData& det_pair_data) -{ - for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) - data_fan_sums[a] = det_pair_data.sum(a); -} + for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) + for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) + { + if (det_pair_data(a,b) == 0) + continue; + int newa = a % num_crystals_per_block; + int newb = b - (a - newa); + if (newa > num_crystals_per_block - 1 - newa) + { + newa = num_crystals_per_block - 1 - newa; + newb = - newb + num_crystals_per_block - 1; + } + // note: add 2*num_detectors to newb to avoid using mod with negative numbers + if (apply) + det_pair_data(a,b) *= + geo_data[newa][(2*num_detectors + newb)%num_detectors]; + else + det_pair_data(a,b) /= + geo_data[newa][(2*num_detectors + newb)%num_detectors]; + } + } -void make_geo_data(GeoData& geo_data, const DetPairData& det_pair_data) -{ - const int num_detectors = det_pair_data.get_num_detectors(); - const int num_crystals_per_block = geo_data.get_length()*2; - const int num_blocks = num_detectors / num_crystals_per_block; - assert(num_blocks * num_crystals_per_block == num_detectors); + void apply_efficiencies(DetPairData& det_pair_data, const Array<1,float>& efficiencies, const bool apply) + { + const int num_detectors = det_pair_data.get_num_detectors(); + for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) + for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) + { + if (det_pair_data(a,b) == 0) + continue; + if (apply) + det_pair_data(a,b) *= + efficiencies[a]*efficiencies[b%num_detectors]; + else + det_pair_data(a,b) /= + efficiencies[a]*efficiencies[b%num_detectors]; + } + } - // TODO optimise - DetPairData work = det_pair_data; - work.fill(0); - for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) - for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) - { - // mirror symmetry - work(a,b) = - det_pair_data(a,b) + - det_pair_data(num_detectors-1-a,(2*num_detectors-1-b)%num_detectors); - } + void make_fan_sum_data(Array<1,float>& data_fan_sums, const DetPairData& det_pair_data) + { + for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) + data_fan_sums[a] = det_pair_data.sum(a); + } - geo_data.fill(0); + void make_geo_data(GeoData& geo_data, const DetPairData& det_pair_data) + { + const int num_detectors = det_pair_data.get_num_detectors(); + const int num_crystals_per_block = geo_data.get_length()*2; + const int num_blocks = num_detectors / num_crystals_per_block; + assert(num_blocks * num_crystals_per_block == num_detectors); - for (int crystal_num_a = 0; crystal_num_a < num_crystals_per_block/2; ++crystal_num_a) - for (int det_num_b = det_pair_data.get_min_index(crystal_num_a); det_num_b <= det_pair_data.get_max_index(crystal_num_a); ++det_num_b) - { - for (int block_num = 0; block_num& efficiencies, - const Array<1,float>& data_fan_sums, - const DetPairData& model) -{ - const int num_detectors = efficiencies.get_length(); + // TODO optimise + DetPairData work = det_pair_data; + work.fill(0); + + for (int a = det_pair_data.get_min_index(); a <= det_pair_data.get_max_index(); ++a) + for (int b = det_pair_data.get_min_index(a); b <= det_pair_data.get_max_index(a); ++b) + { + // mirror symmetry + work(a,b) = + det_pair_data(a,b) + + det_pair_data(num_detectors-1-a,(2*num_detectors-1-b)%num_detectors); + } + + geo_data.fill(0); - for (int a = 0; a < num_detectors; ++a) + for (int crystal_num_a = 0; crystal_num_a < num_crystals_per_block/2; ++crystal_num_a) + for (int det_num_b = det_pair_data.get_min_index(crystal_num_a); det_num_b <= det_pair_data.get_max_index(crystal_num_a); ++det_num_b) + { + for (int block_num = 0; block_num=threshold || - measured_geo_data[a][b] < 10000*norm_geo_data[a][b]) - ? measured_geo_data[a][b] / norm_geo_data[a][b] - : 0; - } -} - -void iterate_block_norm(BlockData& norm_block_data, - const BlockData& measured_block_data, - const DetPairData& model) -{ - make_block_data(norm_block_data, model); - //norm_block_data = measured_block_data / norm_block_data; - const int num_blocks = norm_block_data.get_length(); - const float threshold = measured_block_data.find_max()/10000.F; - for (int a = 0; a < num_blocks; ++a) - for (int b = 0; b < num_blocks; ++b) - { - norm_block_data[a][b] = - (measured_block_data[a][b]>=threshold || - measured_block_data[a][b] < 10000*norm_block_data[a][b]) - ? measured_block_data[a][b] / norm_block_data[a][b] - : 0; - } -} + void iterate_efficiencies(Array<1,float>& efficiencies, + const Array<1,float>& data_fan_sums, + const DetPairData& model) + { + const int num_detectors = efficiencies.get_length(); -double KL(const DetPairData& d1, const DetPairData& d2, const double threshold) -{ - double sum=0; - for (int a = d1.get_min_index(); a <= d1.get_max_index(); ++a) + for (int a = 0; a < num_detectors; ++a) + { + if (data_fan_sums[a] == 0) + efficiencies[a] = 0; + else + { + //const float denominator = inner_product(efficiencies,model[a]); + float denominator = 0; + for (int b = model.get_min_index(a); b <= model.get_max_index(a); ++b) + denominator += efficiencies[b%num_detectors]*model(a,b); + efficiencies[a] = data_fan_sums[a] / denominator; + } + } + } + + void iterate_geo_norm(GeoData& norm_geo_data, + const GeoData& measured_geo_data, + const DetPairData& model) { - double bsum=0; - for (int b = d1.get_min_index(a); b <= d1.get_max_index(a); ++b) - bsum += KL(d1(a,b), d2(a,b), threshold); - sum += bsum; + make_geo_data(norm_geo_data, model); + //norm_geo_data = measured_geo_data / norm_geo_data; + const int num_detectors = model.get_num_detectors(); + const int num_crystals_per_block = measured_geo_data.get_length()*2; + const float threshold = measured_geo_data.find_max()/10000.F; + for (int a = 0; a < num_crystals_per_block/2; ++a) + for (int b = 0; b < num_detectors; ++b) + { + norm_geo_data[a][b] = + (measured_geo_data[a][b]>=threshold || + measured_geo_data[a][b] < 10000*norm_geo_data[a][b]) + ? measured_geo_data[a][b] / norm_geo_data[a][b] + : 0; + } } - return static_cast(sum); -} + void iterate_block_norm(BlockData& norm_block_data, + const BlockData& measured_block_data, + const DetPairData& model) + { + make_block_data(norm_block_data, model); + //norm_block_data = measured_block_data / norm_block_data; + const int num_blocks = norm_block_data.get_length(); + const float threshold = measured_block_data.find_max()/10000.F; + for (int a = 0; a < num_blocks; ++a) + for (int b = 0; b < num_blocks; ++b) + { + norm_block_data[a][b] = + (measured_block_data[a][b]>=threshold || + measured_block_data[a][b] < 10000*norm_block_data[a][b]) + ? measured_block_data[a][b] / norm_block_data[a][b] + : 0; + } + } -////////// ******* GeoData3D ******** ////// -GeoData3D::GeoData3D() -{} + double KL(const DetPairData& d1, const DetPairData& d2, const double threshold) + { + double sum=0; + for (int a = d1.get_min_index(); a <= d1.get_max_index(); ++a) + { + double bsum=0; + for (int b = d1.get_min_index(a); b <= d1.get_max_index(a); ++b) + bsum += KL(d1(a,b), d2(a,b), threshold); + sum += bsum; + } + return static_cast(sum); + } -GeoData3D::~GeoData3D() -{} +////////// ******* GeoData3D ******** ////// + GeoData3D::GeoData3D() + {} -GeoData3D:: -GeoData3D(const int num_axial_crystals_per_block, const int half_num_transaxial_crystals_per_block, const int num_rings, const int num_detectors_per_ring) -:num_axial_crystals_per_block(num_axial_crystals_per_block), half_num_transaxial_crystals_per_block(half_num_transaxial_crystals_per_block), -num_rings(num_rings), num_detectors_per_ring(num_detectors_per_ring) -{ - - // assert(max_ring_diff= (num_axial_crystals_per_block -1) + GeoData3D::~GeoData3D() + {} - - IndexRange<4> fan_indices; - fan_indices.grow(0,num_axial_crystals_per_block-1); - for (int ra = 0; ra < num_axial_crystals_per_block; ++ra) + + GeoData3D:: + GeoData3D(const int num_axial_crystals_per_block, const int half_num_transaxial_crystals_per_block, const int num_rings, const int num_detectors_per_ring) + :num_axial_crystals_per_block(num_axial_crystals_per_block), half_num_transaxial_crystals_per_block(half_num_transaxial_crystals_per_block), + num_rings(num_rings), num_detectors_per_ring(num_detectors_per_ring) { - const int min_rb = ra; - const int max_rb = num_rings-1; // I assumed max ring diff is num_ring_1 - fan_indices[ra].grow(0,half_num_transaxial_crystals_per_block-1); - for (int a = 0; a < half_num_transaxial_crystals_per_block; ++a) + + // assert(max_ring_diff= (num_axial_crystals_per_block -1) + + + IndexRange<4> fan_indices; + fan_indices.grow(0,num_axial_crystals_per_block-1); + for (int ra = 0; ra < num_axial_crystals_per_block; ++ra) { - // store only 1 half of data as ra,a,rb,b = rb,b,ra,a - fan_indices[ra][a].grow(min_rb, max_rb); - for (int rb = min_rb; rb <= max_rb; ++rb) - fan_indices[ra][a][rb] = - IndexRange<1>(a, - a+num_detectors_per_ring -1); // I assumed fan size is number of detector per ring - 2 + const int min_rb = ra; + const int max_rb = num_rings-1; // I assumed max ring diff is num_ring_1 + fan_indices[ra].grow(0,half_num_transaxial_crystals_per_block-1); + for (int a = 0; a < half_num_transaxial_crystals_per_block; ++a) + { + // store only 1 half of data as ra,a,rb,b = rb,b,ra,a + fan_indices[ra][a].grow(min_rb, max_rb); + for (int rb = min_rb; rb <= max_rb; ++rb) + fan_indices[ra][a][rb] = + IndexRange<1>(a, + a+num_detectors_per_ring -1); // I assumed fan size is number of detector per ring - 2 + } } + grow(fan_indices); + fill(0); } - grow(fan_indices); - fill(0); -} -GeoData3D& -GeoData3D::operator=(const GeoData3D& other) -{ - base_type::operator=(other); - num_axial_crystals_per_block = other.num_axial_crystals_per_block; - half_num_transaxial_crystals_per_block = other.half_num_transaxial_crystals_per_block; - num_rings = other.num_rings; - num_detectors_per_ring = other.num_detectors_per_ring; - return *this; -} + GeoData3D& + GeoData3D::operator=(const GeoData3D& other) + { + base_type::operator=(other); + num_axial_crystals_per_block = other.num_axial_crystals_per_block; + half_num_transaxial_crystals_per_block = other.half_num_transaxial_crystals_per_block; + num_rings = other.num_rings; + num_detectors_per_ring = other.num_detectors_per_ring; + return *this; + } /*float & GeoData3D::operator()(const int ra, const int a, const int rb, const int b) { @@ -518,186 +518,186 @@ float GeoData3D::operator()(const int ra, const int a, const int rb, const int b : b < half_num_transaxial_crystals_per_block ? (*this)[rb][b%num_detectors_per_ring][ra][a=0); - assert(b>=0); - return - - (*this)[ra][a%num_detectors_per_ring][rb][b=0); + assert(b>=0); + return -float GeoData3D::operator()(const int ra, const int a, const int rb, const int b) const -{ - assert(a>=0); - assert(b>=0); - return + (*this)[ra][a%num_detectors_per_ring][rb][b=0); + assert(b>=0); + return -bool -GeoData3D:: -is_in_data(const int ra, const int a, const int rb, const int b) const -{ - assert(a>=0); - assert(b>=0); - if (rb<(*this)[ra][a].get_min_index() || rb >(*this)[ra][a].get_max_index()) - return false; - if (b>=get_min_b(a)) - return b<=get_max_b(a); - else - return b+num_detectors_per_ring<=get_max_b(a); -} -void GeoData3D::fill(const float d) -{ - base_type::fill(d); -} + (*this)[ra][a%num_detectors_per_ring][rb][b=0); + assert(b>=0); + if (rb<(*this)[ra][a].get_min_index() || rb >(*this)[ra][a].get_max_index()) + return false; + if (b>=get_min_b(a)) + return b<=get_max_b(a); + else + return b+num_detectors_per_ring<=get_max_b(a); + } + void GeoData3D::fill(const float d) + { + base_type::fill(d); + } -int GeoData3D::get_max_ra() const -{ - return base_type::get_max_index(); -} + int GeoData3D::get_min_ra() const + { + return base_type::get_min_index(); + } -int GeoData3D::get_min_a() const -{ - return (*this)[get_min_index()].get_min_index(); -} + int GeoData3D::get_max_ra() const + { + return base_type::get_max_index(); + } -int GeoData3D::get_max_a() const -{ - return (*this)[get_min_index()].get_max_index(); -} + int GeoData3D::get_min_a() const + { + return (*this)[get_min_index()].get_min_index(); + } -int GeoData3D::get_min_rb(const int ra) const -{ - return 0; - // next is no longer true because we store only half the data - //return (*this)[ra][(*this)[ra].get_min_index()].get_min_index(); -} + int GeoData3D::get_max_a() const + { + return (*this)[get_min_index()].get_max_index(); + } -int GeoData3D::get_max_rb(const int ra) const -{ - return (*this)[ra][(*this)[ra].get_min_index()].get_max_index(); -} + int GeoData3D::get_min_rb(const int ra) const + { + return 0; + // next is no longer true because we store only half the data + //return (*this)[ra][(*this)[ra].get_min_index()].get_min_index(); + } -int GeoData3D::get_min_b(const int a) const -{ - return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_min_index(); -} -int GeoData3D::get_max_b(const int a) const -{ - return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_max_index(); -} + int GeoData3D::get_max_rb(const int ra) const + { + return (*this)[ra][(*this)[ra].get_min_index()].get_max_index(); + } + int GeoData3D::get_min_b(const int a) const + { + return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_min_index(); + } -float GeoData3D::sum() const -{ - //return base_type::sum(); - float sum = 0; - for (int ra=get_min_ra(); ra <= get_max_ra(); ++ra) - for (int a = get_min_a(); a <= get_max_a(); ++a) - sum += this->sum(ra,a); - return sum; -} + int GeoData3D::get_max_b(const int a) const + { + return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_max_index(); + } -float GeoData3D::sum(const int ra, const int a) const -{ - //return (*this)[ra][a].sum(); - float sum = 0; - for (int rb=get_min_rb(ra); rb <= get_max_rb(ra); ++rb) - for (int b = get_min_b(a); b <= get_max_b(a); ++b) - sum += (*this)(ra,a,rb,b%num_detectors_per_ring); - return sum; -} -float GeoData3D::find_max() const -{ - return base_type::find_max(); -} + float GeoData3D::sum() const + { + //return base_type::sum(); + float sum = 0; + for (int ra=get_min_ra(); ra <= get_max_ra(); ++ra) + for (int a = get_min_a(); a <= get_max_a(); ++a) + sum += this->sum(ra,a); + return sum; + } -float GeoData3D::find_min() const -{ - return base_type::find_min(); -} + float GeoData3D::sum(const int ra, const int a) const + { + //return (*this)[ra][a].sum(); + float sum = 0; + for (int rb=get_min_rb(ra); rb <= get_max_rb(ra); ++rb) + for (int b = get_min_b(a); b <= get_max_b(a); ++b) + sum += (*this)(ra,a,rb,b%num_detectors_per_ring); + return sum; + } -int GeoData3D::get_num_axial_crystals_per_block() const -{ - return num_axial_crystals_per_block; -} + float GeoData3D::find_max() const + { + return base_type::find_max(); + } -int GeoData3D::get_half_num_transaxial_crystals_per_block() const -{ - return half_num_transaxial_crystals_per_block; -} + float GeoData3D::find_min() const + { + return base_type::find_min(); + } + int GeoData3D::get_num_axial_crystals_per_block() const + { + return num_axial_crystals_per_block; + } -std::ostream& operator<<(std::ostream& s, const GeoData3D& geo_data) -{ - return s << static_cast(geo_data); -} + int GeoData3D::get_half_num_transaxial_crystals_per_block() const + { + return half_num_transaxial_crystals_per_block; + } -std::istream& operator>>(std::istream& s, GeoData3D& geo_data) -{ - s >> static_cast(geo_data); - if (!s) - return s; - geo_data.half_num_transaxial_crystals_per_block = geo_data.get_max_a() - geo_data.get_min_a() + 1; - geo_data.num_axial_crystals_per_block = geo_data.get_max_ra() - geo_data.get_min_ra() + 1; - - - const int max_delta = geo_data[0][0].get_length()-1; - const int num_detectors_per_ring = - geo_data[0][0][0].get_length(); - geo_data.num_detectors_per_ring = num_detectors_per_ring; - geo_data.num_rings = max_delta+1; - - - - for (int ra = 0; ra < geo_data.num_axial_crystals_per_block; ++ra) - + std::ostream& operator<<(std::ostream& s, const GeoData3D& geo_data) + { + return s << static_cast(geo_data); + } + std::istream& operator>>(std::istream& s, GeoData3D& geo_data) { - const int min_rb = ra; - const int max_rb = geo_data.num_rings-1; - for (int a = 0; a < geo_data.half_num_transaxial_crystals_per_block; ++a) - { + s >> static_cast(geo_data); + if (!s) + return s; + geo_data.half_num_transaxial_crystals_per_block = geo_data.get_max_a() - geo_data.get_min_a() + 1; + geo_data.num_axial_crystals_per_block = geo_data.get_max_ra() - geo_data.get_min_ra() + 1; - if (geo_data[ra][a].get_length() != max_rb - max(ra,min_rb) + 1) - { - warning("Reading GeoData3D: inconsistent length %d for rb at ra=%d, a=%d, " - "Expected length %d\n", - geo_data[ra][a].get_length(), ra, a, max_rb - max(ra,min_rb) + 1); - } - geo_data[ra][a].set_offset(max(ra,min_rb)); - for (int rb = geo_data[ra][a].get_min_index(); rb <= geo_data[ra][a].get_max_index(); ++rb) + + const int max_delta = geo_data[0][0].get_length()-1; + const int num_detectors_per_ring = + geo_data[0][0][0].get_length(); + geo_data.num_detectors_per_ring = num_detectors_per_ring; + geo_data.num_rings = max_delta+1; + + + + + for (int ra = 0; ra < geo_data.num_axial_crystals_per_block; ++ra) + + + { + const int min_rb = ra; + const int max_rb = geo_data.num_rings-1; + for (int a = 0; a < geo_data.half_num_transaxial_crystals_per_block; ++a) { - if (geo_data[ra][a][rb].get_length() != num_detectors_per_ring) + + if (geo_data[ra][a].get_length() != max_rb - max(ra,min_rb) + 1) { - warning("Reading GeoData3D: inconsistent length %d for b at ra=%d, a=%d, rb=%d\n" - "Expected length %d\n", - geo_data[ra][a][rb].get_length(), ra, a, rb, num_detectors_per_ring); + warning("Reading GeoData3D: inconsistent length %d for rb at ra=%d, a=%d, " + "Expected length %d\n", + geo_data[ra][a].get_length(), ra, a, max_rb - max(ra,min_rb) + 1); + } + geo_data[ra][a].set_offset(max(ra,min_rb)); + for (int rb = geo_data[ra][a].get_min_index(); rb <= geo_data[ra][a].get_max_index(); ++rb) + { + if (geo_data[ra][a][rb].get_length() != num_detectors_per_ring) + { + warning("Reading GeoData3D: inconsistent length %d for b at ra=%d, a=%d, rb=%d\n" + "Expected length %d\n", + geo_data[ra][a][rb].get_length(), ra, a, rb, num_detectors_per_ring); + } + geo_data[ra][a][rb].set_offset(a); } - geo_data[ra][a][rb].set_offset(a); } } + + return s; } - - return s; -} //////////////////////////////////////////////////// @@ -705,101 +705,101 @@ std::istream& operator>>(std::istream& s, GeoData3D& geo_data) -FanProjData::FanProjData() -{} + FanProjData::FanProjData() + {} -FanProjData::~FanProjData() -{} + FanProjData::~FanProjData() + {} #if 0 -FanProjData::FanProjData(const IndexRange<4>& range) + FanProjData::FanProjData(const IndexRange<4>& range) :base_type(range), num_rings(range.get_length()), num_detectors_per_ring(range[range.get_min_index()].get_length()) { } #endif -FanProjData:: -FanProjData(const int num_rings, const int num_detectors_per_ring, const int max_ring_diff, const int fan_size) -: num_rings(num_rings), num_detectors_per_ring(num_detectors_per_ring), - max_ring_diff(max_ring_diff), half_fan_size(fan_size/2) -{ - assert(num_detectors_per_ring%2 == 0); - assert(max_ring_diff fan_indices; - fan_indices.grow(0,num_rings-1); - for (int ra = 0; ra < num_rings; ++ra) - { - const int min_rb = max(ra-max_ring_diff, 0); - const int max_rb = min(ra+max_ring_diff, num_rings-1); - fan_indices[ra].grow(0,num_detectors_per_ring-1); - for (int a = 0; a < num_detectors_per_ring; ++a) - { - // store only 1 half of data as ra,a,rb,b = rb,b,ra,a - fan_indices[ra][a].grow(max(ra,min_rb), max_rb); - for (int rb = max(ra,min_rb); rb <= max_rb; ++rb) - fan_indices[ra][a][rb] = - IndexRange<1>(a+num_detectors_per_ring/2-half_fan_size, - a+num_detectors_per_ring/2+half_fan_size); - } - } - grow(fan_indices); - fill(0); -} + FanProjData:: + FanProjData(const int num_rings, const int num_detectors_per_ring, const int max_ring_diff, const int fan_size) + : num_rings(num_rings), num_detectors_per_ring(num_detectors_per_ring), + max_ring_diff(max_ring_diff), half_fan_size(fan_size/2) + { + assert(num_detectors_per_ring%2 == 0); + assert(max_ring_diff fan_indices; + fan_indices.grow(0,num_rings-1); + for (int ra = 0; ra < num_rings; ++ra) + { + const int min_rb = max(ra-max_ring_diff, 0); + const int max_rb = min(ra+max_ring_diff, num_rings-1); + fan_indices[ra].grow(0,num_detectors_per_ring-1); + for (int a = 0; a < num_detectors_per_ring; ++a) + { + // store only 1 half of data as ra,a,rb,b = rb,b,ra,a + fan_indices[ra][a].grow(max(ra,min_rb), max_rb); + for (int rb = max(ra,min_rb); rb <= max_rb; ++rb) + fan_indices[ra][a][rb] = + IndexRange<1>(a+num_detectors_per_ring/2-half_fan_size, + a+num_detectors_per_ring/2+half_fan_size); + } + } + grow(fan_indices); + fill(0); + } -float & FanProjData::operator()(const int ra, const int a, const int rb, const int b) -{ - assert(a>=0); - assert(b>=0); - return - ra=0); - assert(b>=0); - return - ra=0); + assert(b>=0); + return + ra=0); - assert(b>=0); - if (rb<(*this)[ra][a].get_min_index() || rb >(*this)[ra][a].get_max_index()) - return false; - if (b>=get_min_b(a)) - return b<=get_max_b(a); - else - return b+num_detectors_per_ring<=get_max_b(a); -} + float FanProjData::operator()(const int ra, const int a, const int rb, const int b) const + { + assert(a>=0); + assert(b>=0); + return + ra=0); + assert(b>=0); + if (rb<(*this)[ra][a].get_min_index() || rb >(*this)[ra][a].get_max_index()) + return false; + if (b>=get_min_b(a)) + return b<=get_max_b(a); + else + return b+num_detectors_per_ring<=get_max_b(a); + } + + void FanProjData::fill(const float d) + { + base_type::fill(d); + } #if 0 -void FanProjData::grow(const IndexRange<4>& range) + void FanProjData::grow(const IndexRange<4>& range) { base_type::grow(range); num_rings =range.get_length(); @@ -807,469 +807,508 @@ void FanProjData::grow(const IndexRange<4>& range) } #endif -int FanProjData::get_min_ra() const -{ - return base_type::get_min_index(); -} + int FanProjData::get_min_ra() const + { + return base_type::get_min_index(); + } -int FanProjData::get_max_ra() const -{ - return base_type::get_max_index(); -} + int FanProjData::get_max_ra() const + { + return base_type::get_max_index(); + } -int FanProjData::get_min_a() const -{ - return (*this)[get_min_index()].get_min_index(); -} + int FanProjData::get_min_a() const + { + return (*this)[get_min_index()].get_min_index(); + } -int FanProjData::get_max_a() const -{ - return (*this)[get_min_index()].get_max_index(); -} + int FanProjData::get_max_a() const + { + return (*this)[get_min_index()].get_max_index(); + } -int FanProjData::get_min_rb(const int ra) const -{ - return max(ra-max_ring_diff, 0); - // next is no longer true because we store only half the data - //return (*this)[ra][(*this)[ra].get_min_index()].get_min_index(); -} + int FanProjData::get_min_rb(const int ra) const + { + return max(ra-max_ring_diff, 0); + // next is no longer true because we store only half the data + //return (*this)[ra][(*this)[ra].get_min_index()].get_min_index(); + } -int FanProjData::get_max_rb(const int ra) const -{ - return (*this)[ra][(*this)[ra].get_min_index()].get_max_index(); -} + int FanProjData::get_max_rb(const int ra) const + { + return (*this)[ra][(*this)[ra].get_min_index()].get_max_index(); + } -int FanProjData::get_min_b(const int a) const -{ - return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_min_index(); -} + int FanProjData::get_min_b(const int a) const + { + return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_min_index(); + } -int FanProjData::get_max_b(const int a) const -{ - return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_max_index(); -} + int FanProjData::get_max_b(const int a) const + { + return (*this)[get_min_index()][a][(*this)[get_min_index()][a].get_min_index()].get_max_index(); + } -float FanProjData::sum() const -{ - //return base_type::sum(); - float sum = 0; - for (int ra=get_min_ra(); ra <= get_max_ra(); ++ra) - for (int a = get_min_a(); a <= get_max_a(); ++a) - sum += this->sum(ra,a); - return sum; -} + float FanProjData::sum() const + { + //return base_type::sum(); + float sum = 0; + for (int ra=get_min_ra(); ra <= get_max_ra(); ++ra) + for (int a = get_min_a(); a <= get_max_a(); ++a) + sum += this->sum(ra,a); + return sum; + } + + float FanProjData::sum(const int ra, const int a) const + { + //return (*this)[ra][a].sum(); + float sum = 0; + for (int rb=get_min_rb(ra); rb <= get_max_rb(ra); ++rb) + for (int b = get_min_b(a); b <= get_max_b(a); ++b) + sum += (*this)(ra,a,rb,b%num_detectors_per_ring); + return sum; + } + + float FanProjData::find_max() const + { + return base_type::find_max(); + } + + float FanProjData::find_min() const + { + return base_type::find_min(); + } + + int FanProjData::get_num_detectors_per_ring() const + { + return num_detectors_per_ring; + } + + + int FanProjData::get_num_rings() const + { + return num_rings; + } + + std::ostream& operator<<(std::ostream& s, const FanProjData& fan_data) + { + return s << static_cast(fan_data); + } + + std::istream& operator>>(std::istream& s, FanProjData& fan_data) + { + s >> static_cast(fan_data); + if (!s) + return s; + fan_data.num_detectors_per_ring = fan_data.get_max_a() - fan_data.get_min_a() + 1; + fan_data.num_rings = fan_data.get_max_ra() - fan_data.get_min_ra() + 1; + + //int max_delta = 0; + //for (int ra = 0; ra < fan_data.num_rings; ++ra) + // max_delta = max(max_delta,fan_data[ra][0].get_length()-1); + const int max_delta = fan_data[0][0].get_length()-1; + const int half_fan_size = + fan_data[0][0][0].get_length()/2; + fan_data.half_fan_size = half_fan_size; + fan_data.max_ring_diff = max_delta; + + for (int ra = 0; ra < fan_data.num_rings; ++ra) + { + const int min_rb = max(ra-max_delta, 0); + const int max_rb = min(ra+max_delta, fan_data.num_rings-1); + for (int a = 0; a < fan_data.num_detectors_per_ring; ++a) + { + if (fan_data[ra][a].get_length() != max_rb - max(ra,min_rb) + 1) + { + warning("Reading FanProjData: inconsistent length %d for rb at ra=%d, a=%d, " + "Expected length %d\n", + fan_data[ra][a].get_length(), ra, a, max_rb - max(ra,min_rb) + 1); + } + fan_data[ra][a].set_offset(max(ra,min_rb)); + for (int rb = fan_data[ra][a].get_min_index(); rb <= fan_data[ra][a].get_max_index(); ++rb) + { + if (fan_data[ra][a][rb].get_length() != 2*half_fan_size+1) + { + warning("Reading FanProjData: inconsistent length %d for b at ra=%d, a=%d, rb=%d\n" + "Expected length %d\n", + fan_data[ra][a][rb].get_length(), ra, a, rb, 2*half_fan_size+1); + } + fan_data[ra][a][rb].set_offset(a+fan_data.num_detectors_per_ring/2-half_fan_size); + } + } + } + + return s; + } + + shared_ptr + get_fan_info(int& num_rings, int& num_detectors_per_ring, + int& max_ring_diff, int& fan_size, + const ProjDataInfo& proj_data_info) + { + const ProjDataInfoCylindricalNoArcCorr * const proj_data_info_ptr = + dynamic_cast(&proj_data_info); + if (proj_data_info_ptr == 0) + { + error("Can only process not arc-corrected data\n"); + } + if (proj_data_info_ptr->get_view_mashing_factor()>1) + { + error("Can only process data without mashing of views\n"); + } + if (proj_data_info_ptr->get_max_ring_difference(0)>0) + { + error("Can only process data without axial compression (i.e. span=1)\n"); + } + num_rings = + proj_data_info.get_scanner_ptr()->get_num_rings(); + num_detectors_per_ring = + proj_data_info.get_scanner_ptr()->get_num_detectors_per_ring(); + const int half_fan_size = + min(proj_data_info.get_max_tangential_pos_num(), + -proj_data_info.get_min_tangential_pos_num()); + fan_size = 2*half_fan_size+1; + max_ring_diff = proj_data_info_ptr->get_max_segment_num(); + + shared_ptr + ret_value(dynamic_cast(proj_data_info_ptr->clone())); + return ret_value; + + } + + void make_fan_data(FanProjData& fan_data, + const ProjData& proj_data) + { + int num_rings; + int num_detectors_per_ring; + int fan_size; + int max_delta; + shared_ptr proj_data_info_ptr = + get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, + *proj_data.get_proj_data_info_sptr()); + + const int half_fan_size = fan_size/2; + fan_data = FanProjData(num_rings, num_detectors_per_ring, max_delta, 2*half_fan_size+1); + + shared_ptr > segment_ptr; + Bin bin; + + for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) + { + segment_ptr.reset(new SegmentBySinogram(proj_data.get_segment_by_sinogram(bin.segment_num()))); + + for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) + for (bin.tangential_pos_num() = -half_fan_size; + bin.tangential_pos_num() <= half_fan_size; + ++bin.tangential_pos_num()) + { + int ra = 0, a = 0; + int rb = 0, b = 0; + + proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); + + fan_data(ra, a, rb, b) = + fan_data(rb, b, ra, a) = + (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()]; + } + } + } + +/// **** This function make fan_data from projecion file while removing the intermodule gaps **** //// +/// *** fan_data doesn't have gaps, proj_data has gaps *** /// + void make_fan_data_remove_gaps(FanProjData& fan_data, + const ProjData& proj_data) + { + int num_rings; + int num_detectors_per_ring; + int fan_size; + int max_delta; + shared_ptr proj_data_info_ptr = + get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, + *proj_data.get_proj_data_info_sptr()); + + const int half_fan_size = fan_size/2; + + + const int virtual_axial_crystals = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_virtual_axial_crystals_per_block(); + + const int virtual_transaxial_crystals = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_virtual_transaxial_crystals_per_block(); + + const int num_transaxial_blocks = + proj_data_info_ptr ->get_scanner_sptr()-> + get_num_transaxial_blocks(); + const int num_axial_blocks = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_axial_blocks(); + const int num_transaxial_crystals_per_block = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_transaxial_crystals_per_block(); + const int num_axial_crystals_per_block = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_axial_crystals_per_block(); + + const int new_num_transaxial_crystals_per_block = num_transaxial_crystals_per_block - virtual_transaxial_crystals; + + const int new_num_axial_crystals_per_block = num_axial_crystals_per_block - virtual_axial_crystals; + + + const int num_transaxial_blocks_in_fansize = fan_size/(num_transaxial_crystals_per_block); + const int new_fan_size = fan_size - num_transaxial_blocks_in_fansize*virtual_transaxial_crystals; + const int new_half_fan_size = new_fan_size/2; + const int num_axial_blocks_in_max_delta = max_delta/(num_axial_crystals_per_block); + const int new_max_delta = max_delta - (num_axial_blocks_in_max_delta)*virtual_axial_crystals-1; + const int new_num_detectors_per_ring = num_detectors_per_ring - num_transaxial_blocks*virtual_transaxial_crystals; + const int new_num_rings = num_rings - (num_axial_blocks-1)*virtual_axial_crystals; +// std::cerr< > segment_ptr; + Bin bin; + for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) { + segment_ptr.reset(new SegmentBySinogram(proj_data.get_segment_by_sinogram(bin.segment_num()))); + + for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring / 2; bin.view_num()++) + for (bin.tangential_pos_num() = -half_fan_size; + bin.tangential_pos_num() <= half_fan_size; + ++bin.tangential_pos_num()) { + int ra = 0, a = 0; + int rb = 0, b = 0; + + proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); + if ((ra == num_rings - 1) || (rb == num_rings - 1) || (a == num_detectors_per_ring - 1) || + (b == num_detectors_per_ring - 1)) { + continue; + } + + int offset = (a%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(a%num_transaxial_crystals_per_block%new_num_transaxial_crystals_per_block):0; + int new_a = a - a/num_transaxial_crystals_per_block*virtual_transaxial_crystals - offset; + offset = (b%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(b%num_transaxial_crystals_per_block%new_num_transaxial_crystals_per_block):0; + int new_b = b - b/num_transaxial_crystals_per_block*virtual_transaxial_crystals - offset; + offset = (ra%num_axial_crystals_per_block)>new_num_axial_crystals_per_block?(ra%num_axial_crystals_per_block%new_num_axial_crystals_per_block):0; + int new_ra = ra - ra/num_axial_crystals_per_block*virtual_axial_crystals - offset; + offset = (rb%num_axial_crystals_per_block)>new_num_axial_crystals_per_block?(rb%num_axial_crystals_per_block%new_num_axial_crystals_per_block):0; + int new_rb = rb - rb/num_axial_crystals_per_block*virtual_axial_crystals - offset; + + + fan_data(new_ra, new_a, new_rb, new_b) = + fan_data(new_rb, new_b, new_ra, new_a) = + (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()]; + + } + } + } + + + void set_fan_data(ProjData& proj_data, + const FanProjData& fan_data) + { + int num_rings; + int num_detectors_per_ring; + int fan_size; + int max_delta; + shared_ptr proj_data_info_ptr = + get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, + *proj_data.get_proj_data_info_sptr()); + + const int half_fan_size = fan_size/2; + assert(num_rings == fan_data.get_num_rings()); + assert(num_detectors_per_ring == fan_data.get_num_detectors_per_ring()); + + Bin bin; + shared_ptr > segment_ptr; + + for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) + { + segment_ptr.reset(new SegmentBySinogram(proj_data.get_empty_segment_by_sinogram(bin.segment_num()))); + + for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) + for (bin.tangential_pos_num() = -half_fan_size; + bin.tangential_pos_num() <= half_fan_size; + ++bin.tangential_pos_num()) + { + int ra = 0, a = 0; + int rb = 0, b = 0; + + proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); + + (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()] = + fan_data(ra, a, rb, b); + } + proj_data.set_segment(*segment_ptr); + } + } + +/// **** This function make proj_data from fan_data while adding the intermodule gaps **** //// +/// *** fan_data doesn't have gaps, proj_data has gaps *** /// + void set_fan_data_add_gaps(ProjData& proj_data, + const FanProjData& fan_data) + { + int num_rings; + int num_detectors_per_ring; + int fan_size; + int max_delta; + shared_ptr proj_data_info_ptr = + get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, + *proj_data.get_proj_data_info_sptr()); -float FanProjData::sum(const int ra, const int a) const -{ - //return (*this)[ra][a].sum(); - float sum = 0; - for (int rb=get_min_rb(ra); rb <= get_max_rb(ra); ++rb) - for (int b = get_min_b(a); b <= get_max_b(a); ++b) - sum += (*this)(ra,a,rb,b%num_detectors_per_ring); - return sum; -} + const int half_fan_size = fan_size/2; + //assert(num_rings == fan_data.get_num_rings()); + //assert(num_detectors_per_ring == fan_data.get_num_detectors_per_ring()); -float FanProjData::find_max() const -{ - return base_type::find_max(); -} -float FanProjData::find_min() const -{ - return base_type::find_min(); -} + // **** Added by Tahereh Nikjenad **** // + const int virtual_axial_crystals = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_virtual_axial_crystals_per_block(); -int FanProjData::get_num_detectors_per_ring() const -{ - return num_detectors_per_ring; -} + const int virtual_transaxial_crystals = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_virtual_transaxial_crystals_per_block(); + const int num_transaxial_crystals_per_block = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_transaxial_crystals_per_block(); + const int num_axial_crystals_per_block = + proj_data_info_ptr->get_scanner_sptr()-> + get_num_axial_crystals_per_block(); -int FanProjData::get_num_rings() const -{ - return num_rings; -} + const int new_num_transaxial_crystals_per_block = num_transaxial_crystals_per_block - virtual_transaxial_crystals; -std::ostream& operator<<(std::ostream& s, const FanProjData& fan_data) -{ - return s << static_cast(fan_data); -} + const int new_num_axial_crystals_per_block = num_axial_crystals_per_block - virtual_axial_crystals; -std::istream& operator>>(std::istream& s, FanProjData& fan_data) -{ - s >> static_cast(fan_data); - if (!s) - return s; - fan_data.num_detectors_per_ring = fan_data.get_max_a() - fan_data.get_min_a() + 1; - fan_data.num_rings = fan_data.get_max_ra() - fan_data.get_min_ra() + 1; - - //int max_delta = 0; - //for (int ra = 0; ra < fan_data.num_rings; ++ra) - // max_delta = max(max_delta,fan_data[ra][0].get_length()-1); - const int max_delta = fan_data[0][0].get_length()-1; - const int half_fan_size = - fan_data[0][0][0].get_length()/2; - fan_data.half_fan_size = half_fan_size; - fan_data.max_ring_diff = max_delta; - - for (int ra = 0; ra < fan_data.num_rings; ++ra) - { - const int min_rb = max(ra-max_delta, 0); - const int max_rb = min(ra+max_delta, fan_data.num_rings-1); - for (int a = 0; a < fan_data.num_detectors_per_ring; ++a) - { - if (fan_data[ra][a].get_length() != max_rb - max(ra,min_rb) + 1) - { - warning("Reading FanProjData: inconsistent length %d for rb at ra=%d, a=%d, " - "Expected length %d\n", - fan_data[ra][a].get_length(), ra, a, max_rb - max(ra,min_rb) + 1); - } - fan_data[ra][a].set_offset(max(ra,min_rb)); - for (int rb = fan_data[ra][a].get_min_index(); rb <= fan_data[ra][a].get_max_index(); ++rb) - { - if (fan_data[ra][a][rb].get_length() != 2*half_fan_size+1) - { - warning("Reading FanProjData: inconsistent length %d for b at ra=%d, a=%d, rb=%d\n" - "Expected length %d\n", - fan_data[ra][a][rb].get_length(), ra, a, rb, 2*half_fan_size+1); - } - fan_data[ra][a][rb].set_offset(a+fan_data.num_detectors_per_ring/2-half_fan_size); - } - } - } + // const int num_axial_detectors = fan_data.get_num_rings(); // Number of ring in fan data (w/o gaps) + // const int num_transaxial_detectors = fan_data.get_num_detectors_per_ring(); // number of detector per ring in fan data w/o gaps - return s; -} + // const int num_axial_crystals_per_block = num_axial_detectors/num_axial_blocks; // number of axial detector per block in fan_data w/o gaps + // const int num_transaxial_crystals_per_block = num_transaxial_detectors/num_transaxial_blocks; // number of transaxial detector per block in fan_data w/o gaps -shared_ptr -get_fan_info(int& num_rings, int& num_detectors_per_ring, - int& max_ring_diff, int& fan_size, - const ProjDataInfo& proj_data_info) -{ - const ProjDataInfoCylindricalNoArcCorr * const proj_data_info_ptr = - dynamic_cast(&proj_data_info); - if (proj_data_info_ptr == 0) - { - error("Can only process not arc-corrected data\n"); - } - if (proj_data_info_ptr->get_view_mashing_factor()>1) - { - error("Can only process data without mashing of views\n"); - } - if (proj_data_info_ptr->get_max_ring_difference(0)>0) - { - error("Can only process data without axial compression (i.e. span=1)\n"); - } - num_rings = - proj_data_info.get_scanner_ptr()->get_num_rings(); - num_detectors_per_ring = - proj_data_info.get_scanner_ptr()->get_num_detectors_per_ring(); - const int half_fan_size = - min(proj_data_info.get_max_tangential_pos_num(), - -proj_data_info.get_min_tangential_pos_num()); - fan_size = 2*half_fan_size+1; - max_ring_diff = proj_data_info_ptr->get_max_segment_num(); - - shared_ptr - ret_value(dynamic_cast(proj_data_info_ptr->clone())); - return ret_value; - -} -void make_fan_data(FanProjData& fan_data, - const ProjData& proj_data) -{ - int num_rings; - int num_detectors_per_ring; - int fan_size; - int max_delta; - shared_ptr proj_data_info_ptr = - get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, - *proj_data.get_proj_data_info_sptr()); - - const int half_fan_size = fan_size/2; - fan_data = FanProjData(num_rings, num_detectors_per_ring, max_delta, 2*half_fan_size+1); - - shared_ptr > segment_ptr; - Bin bin; - - for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) - { - segment_ptr.reset(new SegmentBySinogram(proj_data.get_segment_by_sinogram(bin.segment_num()))); - - for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); - bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); - ++bin.axial_pos_num()) - for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) - for (bin.tangential_pos_num() = -half_fan_size; - bin.tangential_pos_num() <= half_fan_size; - ++bin.tangential_pos_num()) - { - int ra = 0, a = 0; - int rb = 0, b = 0; - - proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); - - fan_data(ra, a, rb, b) = - fan_data(rb, b, ra, a) = - (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()]; - } - } -} + // **** End **** // -/// **** This function make fan_data from projecion file while removing the intermodule gaps **** //// -/// *** fan_data doesn't have gaps, proj_data has gaps *** /// -void make_fan_data_remove_gaps(FanProjData& fan_data, - const ProjData& proj_data) -{ - int num_rings; - int num_detectors_per_ring; - int fan_size; - int max_delta; - shared_ptr proj_data_info_ptr = - get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, - *proj_data.get_proj_data_info_sptr()); - - const int half_fan_size = fan_size/2; + Bin bin; + shared_ptr > segment_ptr; - - // **** Added by me **** // - - - const int num_transaxial_blocks = - proj_data_info_ptr ->get_scanner_sptr()-> - get_num_transaxial_blocks(); - const int num_axial_blocks = - proj_data_info_ptr->get_scanner_sptr()-> - get_num_axial_blocks(); - const int num_transaxial_crystals_per_block = - proj_data_info_ptr->get_scanner_sptr()-> - get_num_transaxial_crystals_per_block(); - const int num_axial_crystals_per_block = - proj_data_info_ptr->get_scanner_sptr()-> - get_num_axial_crystals_per_block(); - - + for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) + { + segment_ptr.reset(new SegmentBySinogram(proj_data.get_empty_segment_by_sinogram(bin.segment_num()))); + + for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) + for (bin.tangential_pos_num() = -half_fan_size; + bin.tangential_pos_num() <= half_fan_size; + ++bin.tangential_pos_num()) + { + int ra = 0, a = 0; + int rb = 0, b = 0; - const int num_transaxial_blocks_in_fansize = fan_size/num_transaxial_crystals_per_block; - const int new_fan_size = fan_size - num_transaxial_blocks_in_fansize; - const int new_half_fan_size = new_fan_size/2; - const int num_axial_blocks_in_max_delta = max_delta/num_axial_crystals_per_block; - const int new_max_delta = max_delta - num_axial_blocks_in_max_delta - 1; - const int new_num_detectors_per_ring = num_detectors_per_ring - num_transaxial_blocks; - const int new_num_rings = num_rings - num_axial_blocks; - - // **** End **** // + proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); - fan_data = FanProjData(new_num_rings, new_num_detectors_per_ring, new_max_delta, 2*new_half_fan_size+1); + (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()] = 0; - - shared_ptr > segment_ptr; - Bin bin; - - for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) - { - segment_ptr.reset(new SegmentBySinogram(proj_data.get_segment_by_sinogram(bin.segment_num()))); - - for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); - bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); - ++bin.axial_pos_num()) - for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) - for (bin.tangential_pos_num() = -half_fan_size; - bin.tangential_pos_num() <= half_fan_size; - ++bin.tangential_pos_num()) - { - int ra = 0, a = 0; - int rb = 0, b = 0; - - proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); - int new_a = a - a/num_transaxial_crystals_per_block; - int new_b = b - b/num_transaxial_crystals_per_block; - int new_ra = ra - ra/ num_axial_crystals_per_block; - int new_rb = rb - rb/num_axial_crystals_per_block; - - if ((ra == num_rings -1) || (rb == num_rings -1) || (a == num_detectors_per_ring-1) || (b == num_detectors_per_ring-1)) continue; - - - fan_data(new_ra, new_a, new_rb, new_b) = - fan_data(new_rb, new_b, new_ra, new_a) = - (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()]; - } - } -} + if ((ra+1)% num_axial_crystals_per_block == 0) continue; + if ((rb+1)% num_axial_crystals_per_block == 0) continue; + if ((a+1)% num_transaxial_crystals_per_block == 0) continue; + if ((b+1)% num_transaxial_crystals_per_block == 0) continue; + if ((ra == num_rings - 1) || (rb == num_rings - 1) ) continue; -void set_fan_data(ProjData& proj_data, - const FanProjData& fan_data) -{ - int num_rings; - int num_detectors_per_ring; - int fan_size; - int max_delta; - shared_ptr proj_data_info_ptr = - get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, - *proj_data.get_proj_data_info_sptr()); - - const int half_fan_size = fan_size/2; - assert(num_rings == fan_data.get_num_rings()); - assert(num_detectors_per_ring == fan_data.get_num_detectors_per_ring()); - - Bin bin; - shared_ptr > segment_ptr; - - for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) - { - segment_ptr.reset(new SegmentBySinogram(proj_data.get_empty_segment_by_sinogram(bin.segment_num()))); - - for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); - bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); - ++bin.axial_pos_num()) - for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) - for (bin.tangential_pos_num() = -half_fan_size; - bin.tangential_pos_num() <= half_fan_size; - ++bin.tangential_pos_num()) - { - int ra = 0, a = 0; - int rb = 0, b = 0; - - proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); - - (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()] = - fan_data(ra, a, rb, b); - } - proj_data.set_segment(*segment_ptr); - } -} -/// **** This function make proj_data from fan_data while adding the intermodule gaps **** //// -/// *** fan_data doesn't have gaps, proj_data has gaps *** /// -void set_fan_data_add_gaps(ProjData& proj_data, - const FanProjData& fan_data) -{ - int num_rings; - int num_detectors_per_ring; - int fan_size; - int max_delta; - shared_ptr proj_data_info_ptr = - get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, - *proj_data.get_proj_data_info_sptr()); - - const int half_fan_size = fan_size/2; - //assert(num_rings == fan_data.get_num_rings()); - //assert(num_detectors_per_ring == fan_data.get_num_detectors_per_ring()); - - - // **** Added by Tahereh Nikjenad **** // - - const int num_transaxial_crystals_per_block = - proj_data_info_ptr->get_scanner_sptr()-> - get_num_transaxial_crystals_per_block(); - const int num_axial_crystals_per_block = - proj_data_info_ptr->get_scanner_sptr()-> - get_num_axial_crystals_per_block(); - - // const int num_axial_detectors = fan_data.get_num_rings(); // Number of ring in fan data (w/o gaps) - // const int num_transaxial_detectors = fan_data.get_num_detectors_per_ring(); // number of detector per ring in fan data w/o gaps - - // const int num_axial_crystals_per_block = num_axial_detectors/num_axial_blocks; // number of axial detector per block in fan_data w/o gaps - // const int num_transaxial_crystals_per_block = num_transaxial_detectors/num_transaxial_blocks; // number of transaxial detector per block in fan_data w/o gaps + int offset = (a%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(a%num_transaxial_crystals_per_block%new_num_transaxial_crystals_per_block):0; + int new_a = a - a/num_transaxial_crystals_per_block*virtual_transaxial_crystals - offset; + offset = (b%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(b%num_transaxial_crystals_per_block)%new_num_transaxial_crystals_per_block:0; + int new_b = b - b/num_transaxial_crystals_per_block*virtual_transaxial_crystals - offset; + offset = (ra%num_axial_crystals_per_block)>new_num_axial_crystals_per_block?(ra%num_axial_crystals_per_block)%new_num_axial_crystals_per_block:0; + int new_ra = ra - ra/num_axial_crystals_per_block*virtual_axial_crystals - offset; + offset = (rb%num_axial_crystals_per_block)>new_num_axial_crystals_per_block?(rb%num_axial_crystals_per_block)%new_num_axial_crystals_per_block:0; + int new_rb = rb - rb/num_axial_crystals_per_block*virtual_axial_crystals - offset; +// int new_a = a - a/num_transaxial_crystals_per_block; +// int new_b = b - b/num_transaxial_crystals_per_block; +// int new_ra = ra - ra/ num_axial_crystals_per_block; +// int new_rb = rb - rb/num_axial_crystals_per_block; - - // **** End **** // - - Bin bin; - shared_ptr > segment_ptr; - - for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) - { - segment_ptr.reset(new SegmentBySinogram(proj_data.get_empty_segment_by_sinogram(bin.segment_num()))); - - for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); - bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); - ++bin.axial_pos_num()) - for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) - for (bin.tangential_pos_num() = -half_fan_size; - bin.tangential_pos_num() <= half_fan_size; - ++bin.tangential_pos_num()) - { - int ra = 0, a = 0; - int rb = 0, b = 0; - - proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); - - (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()] = 1e20; - - - if ((ra+1)% num_axial_crystals_per_block == 0) continue; - if ((rb+1)% num_axial_crystals_per_block == 0) continue; - if ((a+1)% num_transaxial_crystals_per_block == 0) continue; - if ((b+1)% num_transaxial_crystals_per_block == 0) continue; - - int new_a = a - a/num_transaxial_crystals_per_block; - int new_b = b - b/num_transaxial_crystals_per_block; - int new_ra = ra - ra/ num_axial_crystals_per_block; - int new_rb = rb - rb/num_axial_crystals_per_block; - - - (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()] = - fan_data(new_ra, new_a, new_rb, new_b); - } - proj_data.set_segment(*segment_ptr); + + (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()] = + fan_data(new_ra, new_a, new_rb, new_b); + } + proj_data.set_segment(*segment_ptr); + } } -} -void apply_block_norm(FanProjData& fan_data, const BlockData3D& block_data, const bool apply) -{ - const int num_axial_detectors = fan_data.get_num_rings(); - const int num_tangential_detectors = fan_data.get_num_detectors_per_ring(); - const int num_axial_blocks = block_data.get_num_rings(); - const int num_tangential_blocks = block_data.get_num_detectors_per_ring(); - const int num_axial_crystals_per_block = num_axial_detectors/num_axial_blocks; - assert(num_axial_blocks * num_axial_crystals_per_block == num_axial_detectors); - const int num_tangential_crystals_per_block = num_tangential_detectors/num_tangential_blocks; - assert(num_tangential_blocks * num_tangential_crystals_per_block == num_tangential_detectors); - - for (int ra = fan_data.get_min_ra(); ra <= fan_data.get_max_ra(); ++ra) - for (int a = fan_data.get_min_a(); a <= fan_data.get_max_a(); ++a) - // loop rb from ra to avoid double counting - for (int rb = max(ra,fan_data.get_min_rb(ra)); rb <= fan_data.get_max_rb(ra); ++rb) - for (int b = fan_data.get_min_b(a); b <= fan_data.get_max_b(a); ++b) - { - // note: add 2*num_detectors_per_ring to newb to avoid using mod with negative numbers - if (fan_data(ra,a,rb,b) == 0) - continue; - if (apply) - fan_data(ra,a,rb,b) *= - block_data(ra/num_axial_crystals_per_block,a/num_tangential_crystals_per_block, - rb/num_axial_crystals_per_block,b/num_tangential_crystals_per_block); - else - fan_data(ra,a,rb,b) /= - block_data(ra/num_axial_crystals_per_block,a/num_tangential_crystals_per_block, - rb/num_axial_crystals_per_block,b/num_tangential_crystals_per_block); - } -} + void apply_block_norm(FanProjData& fan_data, const BlockData3D& block_data, const bool apply) + { + const int num_axial_detectors = fan_data.get_num_rings(); + const int num_tangential_detectors = fan_data.get_num_detectors_per_ring(); + const int num_axial_blocks = block_data.get_num_rings(); + const int num_tangential_blocks = block_data.get_num_detectors_per_ring(); + const int num_axial_crystals_per_block = num_axial_detectors/num_axial_blocks; + assert(num_axial_blocks * num_axial_crystals_per_block == num_axial_detectors); + const int num_tangential_crystals_per_block = num_tangential_detectors/num_tangential_blocks; + assert(num_tangential_blocks * num_tangential_crystals_per_block == num_tangential_detectors); + + for (int ra = fan_data.get_min_ra(); ra <= fan_data.get_max_ra(); ++ra) + for (int a = fan_data.get_min_a(); a <= fan_data.get_max_a(); ++a) + // loop rb from ra to avoid double counting + for (int rb = max(ra,fan_data.get_min_rb(ra)); rb <= fan_data.get_max_rb(ra); ++rb) + for (int b = fan_data.get_min_b(a); b <= fan_data.get_max_b(a); ++b) + { + // note: add 2*num_detectors_per_ring to newb to avoid using mod with negative numbers + if (fan_data(ra,a,rb,b) == 0) + continue; + if (apply) + fan_data(ra,a,rb,b) *= + block_data(ra/num_axial_crystals_per_block,a/num_tangential_crystals_per_block, + rb/num_axial_crystals_per_block,b/num_tangential_crystals_per_block); + else + fan_data(ra,a,rb,b) /= + block_data(ra/num_axial_crystals_per_block,a/num_tangential_crystals_per_block, + rb/num_axial_crystals_per_block,b/num_tangential_crystals_per_block); + } + } #if 0 -void apply_geo_norm(FanProjData& fan_data, const GeoData& geo_data, const bool apply) + void apply_geo_norm(FanProjData& fan_data, const GeoData& geo_data, const bool apply) { const int num_detectors_per_ring = fan_data.get_num_detectors_per_ring(); const int num_crystals_per_block = geo_data.get_length()*2; for (int a = fan_data.get_min_ra(); a <= fan_data.get_max_ra(); ++a) - for (int b = fan_data.get_min_ra(a); b <= fan_data.get_max_ra(a); ++b) + for (int b = fan_data.get_min_ra(a); b <= fan_data.get_max_ra(a); ++b) { if (fan_data(ra,a,rb,b) == 0) continue; int newa = a % num_crystals_per_block; - int newb = b - (a - newa); + int newb = b - (a - newa); if (newa > num_crystals_per_block - 1 - newa) - { - newa = num_crystals_per_block - 1 - newa; + { + newa = num_crystals_per_block - 1 - newa; newb = - newb + num_crystals_per_block - 1; } // note: add 2*num_detectors_per_ring to newb to avoid using mod with negative numbers @@ -1283,329 +1322,328 @@ void apply_geo_norm(FanProjData& fan_data, const GeoData& geo_data, const bool a } #endif -void apply_geo_norm(FanProjData& fan_data, const GeoData3D& geo_data, const bool apply) -{ - - const int num_axial_detectors = fan_data.get_num_rings(); - const int num_transaxial_detectors = fan_data.get_num_detectors_per_ring(); - const int num_axial_crystals_per_block = geo_data.get_num_axial_crystals_per_block(); - const int num_transaxial_crystals_per_block = geo_data.get_half_num_transaxial_crystals_per_block()*2; - - const int num_transaxial_blocks = num_transaxial_detectors/num_transaxial_crystals_per_block; - const int num_axial_blocks = num_axial_detectors/num_axial_crystals_per_block; - - FanProjData work = fan_data; - work.fill(0); - - - for (int ra = 0; ra < num_axial_crystals_per_block; ++ra) - for (int a = 0; a < num_transaxial_crystals_per_block /2; ++a) - // loop rb from ra to avoid double counting - for (int rb = max(ra,fan_data.get_min_rb(ra)); rb <= fan_data.get_max_rb(ra); ++rb) - for (int b = fan_data.get_min_b(a); b <= fan_data.get_max_b(a); ++b) - { - - - // rotation - - for (int axial_block_num = 0; axial_block_num& data_fan_sums, const FanProjData& fan_data) + { + for (int ra = fan_data.get_min_ra(); ra <= fan_data.get_max_ra(); ++ra) + for (int a = fan_data.get_min_a(); a <= fan_data.get_max_a(); ++a) + data_fan_sums[ra][a] = fan_data.sum(ra,a); + } + + void make_fan_sum_data(Array<2,float>& data_fan_sums, + const ProjData& proj_data) + { + int num_rings; + int num_detectors_per_ring; + int fan_size; + int max_delta; + shared_ptr proj_data_info_ptr = + get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, + *proj_data.get_proj_data_info_sptr()); + const int half_fan_size = fan_size/2; + data_fan_sums.fill(0); + + shared_ptr > segment_ptr; + Bin bin; + + for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) { - if (fan_data(ra,a,rb,b) == 0) - continue; - if (apply) - fan_data(ra,a,rb,b) *= - efficiencies[ra][a]*efficiencies[rb][b%num_detectors_per_ring]; - else - fan_data(ra,a,rb,b) /= - efficiencies[ra][a]*efficiencies[rb][b%num_detectors_per_ring]; - } -} + segment_ptr.reset(new SegmentBySinogram(proj_data.get_segment_by_sinogram(bin.segment_num()))); + + for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) + for (bin.tangential_pos_num() = -half_fan_size; + bin.tangential_pos_num() <= half_fan_size; + ++bin.tangential_pos_num()) + { + int ra = 0, a = 0; + int rb = 0, b = 0; -void make_fan_sum_data(Array<2,float>& data_fan_sums, const FanProjData& fan_data) -{ - for (int ra = fan_data.get_min_ra(); ra <= fan_data.get_max_ra(); ++ra) - for (int a = fan_data.get_min_a(); a <= fan_data.get_max_a(); ++a) - data_fan_sums[ra][a] = fan_data.sum(ra,a); -} + proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); -void make_fan_sum_data(Array<2,float>& data_fan_sums, - const ProjData& proj_data) -{ - int num_rings; - int num_detectors_per_ring; - int fan_size; - int max_delta; - shared_ptr proj_data_info_ptr = - get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, - *proj_data.get_proj_data_info_sptr()); - const int half_fan_size = fan_size/2; - data_fan_sums.fill(0); - - shared_ptr > segment_ptr; - Bin bin; - - for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) - { - segment_ptr.reset(new SegmentBySinogram(proj_data.get_segment_by_sinogram(bin.segment_num()))); - - for (bin.axial_pos_num() = proj_data.get_min_axial_pos_num(bin.segment_num()); - bin.axial_pos_num() <= proj_data.get_max_axial_pos_num(bin.segment_num()); - ++bin.axial_pos_num()) - for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring/2; bin.view_num()++) - for (bin.tangential_pos_num() = -half_fan_size; - bin.tangential_pos_num() <= half_fan_size; - ++bin.tangential_pos_num()) - { - int ra = 0, a = 0; - int rb = 0, b = 0; - - proj_data_info_ptr->get_det_pair_for_bin(a, ra, b, rb, bin); - - const float value = - (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()]; - data_fan_sums[ra][a] += value; - data_fan_sums[rb][b] += value; - } - } -} + const float value = + (*segment_ptr)[bin.axial_pos_num()][bin.view_num()][bin.tangential_pos_num()]; + data_fan_sums[ra][a] += value; + data_fan_sums[rb][b] += value; + } + } + } -void make_fan_sum_data(Array<2,float>& data_fan_sums, - const DetectorEfficiencies& efficiencies, - const int max_ring_diff, const int half_fan_size) -{ - const int num_rings = data_fan_sums.get_length(); - assert(data_fan_sums.get_min_index()==0); - const int num_detectors_per_ring = - data_fan_sums[0].get_length(); + void make_fan_sum_data(Array<2,float>& data_fan_sums, + const DetectorEfficiencies& efficiencies, + const int max_ring_diff, const int half_fan_size) + { + const int num_rings = data_fan_sums.get_length(); + assert(data_fan_sums.get_min_index()==0); + const int num_detectors_per_ring = + data_fan_sums[0].get_length(); - for (int ra = data_fan_sums.get_min_index(); ra <= data_fan_sums.get_max_index(); ++ra) - for (int a = data_fan_sums[ra].get_min_index(); a <= data_fan_sums[ra].get_max_index(); ++a) - { - float fan_sum = 0; - for (int rb = max(ra-max_ring_diff, 0); rb <= min(ra+max_ring_diff, num_rings-1); ++rb) - for (int b = a+num_detectors_per_ring/2-half_fan_size; b <= a+num_detectors_per_ring/2+half_fan_size; ++b) - fan_sum += efficiencies[rb][b%num_detectors_per_ring]; - data_fan_sums[ra][a] = efficiencies[ra][a]*fan_sum; - } -} + for (int ra = data_fan_sums.get_min_index(); ra <= data_fan_sums.get_max_index(); ++ra) + for (int a = data_fan_sums[ra].get_min_index(); a <= data_fan_sums[ra].get_max_index(); ++a) + { + float fan_sum = 0; + for (int rb = max(ra-max_ring_diff, 0); rb <= min(ra+max_ring_diff, num_rings-1); ++rb) + for (int b = a+num_detectors_per_ring/2-half_fan_size; b <= a+num_detectors_per_ring/2+half_fan_size; ++b) + fan_sum += efficiencies[rb][b%num_detectors_per_ring]; + data_fan_sums[ra][a] = efficiencies[ra][a]*fan_sum; + } + } -void make_geo_data(GeoData3D& geo_data, const FanProjData& fan_data) -{ - - const int num_axial_detectors = fan_data.get_num_rings(); - const int num_transaxial_detectors = fan_data.get_num_detectors_per_ring(); - const int num_axial_crystals_per_block = geo_data.get_num_axial_crystals_per_block(); - const int num_transaxial_crystals_per_block = geo_data.get_half_num_transaxial_crystals_per_block()*2; - const int num_transaxial_blocks = num_transaxial_detectors/num_transaxial_crystals_per_block; - const int num_axial_blocks = num_axial_detectors/num_axial_crystals_per_block; - - - // transaxial and axial mirroring - - FanProjData work = fan_data; - work.fill(0); - - for (int ra = fan_data.get_min_ra(); ra <= fan_data.get_max_ra(); ++ra) - for (int a = fan_data.get_min_a(); a <= fan_data.get_max_a(); ++a) - //1// for (int rb = fan_data.get_min_ra(); rb <= fan_data.get_max_ra(); ++rb) - for (int rb = max(ra,fan_data.get_min_rb(ra)); rb <= fan_data.get_max_rb(ra); ++rb) - for (int b = fan_data.get_min_b(a); b <= fan_data.get_max_b(a); ++b) - { - - const int ma = num_transaxial_detectors-1-a; - const int mb = (2*num_transaxial_detectors-1-b)%num_transaxial_detectors; - const int mra = num_axial_detectors-1-ra; - const int mrb = (num_axial_detectors-1-rb); - - - if (ra!=mra && rb !=mrb) - work(ra,a,rb,b)= fan_data(ra,a,rb,b) + fan_data(ra,ma,rb,mb) + fan_data(mra,a,mrb,b) + fan_data(mra,ma,mrb,mb); - else - work(ra,a,rb,b)= fan_data(ra,a,rb,b) + fan_data(ra,ma,rb,mb); - - } - - - geo_data.fill(0); - - - for (int ra = 0; ra < num_axial_crystals_per_block; ++ra) - // for (int a = 0; a <= num_transaxial_detectors/2; ++a) - for (int a = 0; a & data_fan_sums, - const FanProjData& model) -{ - const int num_detectors_per_ring = model.get_num_detectors_per_ring(); - - assert(model.get_min_ra() == data_fan_sums.get_min_index()); - assert(model.get_max_ra() == data_fan_sums.get_max_index()); - assert(model.get_min_a() == data_fan_sums[data_fan_sums.get_min_index()].get_min_index()); - assert(model.get_max_a() == data_fan_sums[data_fan_sums.get_min_index()].get_max_index()); - for (int ra = model.get_min_ra(); ra <= model.get_max_ra(); ++ra) - for (int a = model.get_min_a(); a <= model.get_max_a(); ++a) - { - if (data_fan_sums[ra][a] == 0) - efficiencies[ra][a] = 0; - else - { - float denominator = 0; - for (int rb = model.get_min_rb(ra); rb <= model.get_max_rb(ra); ++rb) - for (int b = model.get_min_b(a); b <= model.get_max_b(a); ++b) - denominator += efficiencies[rb][b%num_detectors_per_ring]*model(ra,a,rb,b); - efficiencies[ra][a] = data_fan_sums[ra][a] / denominator; - } + void iterate_efficiencies(DetectorEfficiencies& efficiencies, + const Array<2,float>& data_fan_sums, + const FanProjData& model) + { + const int num_detectors_per_ring = model.get_num_detectors_per_ring(); + assert(model.get_min_ra() == data_fan_sums.get_min_index()); + assert(model.get_max_ra() == data_fan_sums.get_max_index()); + assert(model.get_min_a() == data_fan_sums[data_fan_sums.get_min_index()].get_min_index()); + assert(model.get_max_a() == data_fan_sums[data_fan_sums.get_min_index()].get_max_index()); + for (int ra = model.get_min_ra(); ra <= model.get_max_ra(); ++ra) + for (int a = model.get_min_a(); a <= model.get_max_a(); ++a) + { + if (data_fan_sums[ra][a] == 0) + efficiencies[ra][a] = 0; + else + { + float denominator = 0; + for (int rb = model.get_min_rb(ra); rb <= model.get_max_rb(ra); ++rb) + for (int b = model.get_min_b(a); b <= model.get_max_b(a); ++b) + denominator += efficiencies[rb][b%num_detectors_per_ring]*model(ra,a,rb,b); + efficiencies[ra][a] = data_fan_sums[ra][a] / denominator; + } + } } -} // version without model -void iterate_efficiencies(DetectorEfficiencies& efficiencies, - const Array<2,float>& data_fan_sums, - const int max_ring_diff, const int half_fan_size) -{ - const int num_rings = data_fan_sums.get_length(); - const int num_detectors_per_ring = data_fan_sums[data_fan_sums.get_min_index()].get_length(); + void iterate_efficiencies(DetectorEfficiencies& efficiencies, + const Array<2,float>& data_fan_sums, + const int max_ring_diff, const int half_fan_size) + { + const int num_rings = data_fan_sums.get_length(); + const int num_detectors_per_ring = data_fan_sums[data_fan_sums.get_min_index()].get_length(); #ifdef WRITE_ALL - static int sub_iter_num = 0; + static int sub_iter_num = 0; #endif - for (int ra = data_fan_sums.get_min_index(); ra <= data_fan_sums.get_max_index(); ++ra) - for (int a = data_fan_sums[ra].get_min_index(); a <= data_fan_sums[ra].get_max_index(); ++a) - { - if (data_fan_sums[ra][a] == 0) - efficiencies[ra][a] = 0; - else - { - float denominator = 0; - for (int rb = max(ra-max_ring_diff, 0); rb <= min(ra+max_ring_diff, num_rings-1); ++rb) - for (int b = a+num_detectors_per_ring/2-half_fan_size; b <= a+num_detectors_per_ring/2+half_fan_size; ++b) - denominator += efficiencies[rb][b%num_detectors_per_ring]; - efficiencies[ra][a] = data_fan_sums[ra][a] / denominator; - } + for (int ra = data_fan_sums.get_min_index(); ra <= data_fan_sums.get_max_index(); ++ra) + for (int a = data_fan_sums[ra].get_min_index(); a <= data_fan_sums[ra].get_max_index(); ++a) + { + if (data_fan_sums[ra][a] == 0) + efficiencies[ra][a] = 0; + else + { + float denominator = 0; + for (int rb = max(ra-max_ring_diff, 0); rb <= min(ra+max_ring_diff, num_rings-1); ++rb) + for (int b = a+num_detectors_per_ring/2-half_fan_size; b <= a+num_detectors_per_ring/2+half_fan_size; ++b) + denominator += efficiencies[rb][b%num_detectors_per_ring]; + efficiencies[ra][a] = data_fan_sums[ra][a] / denominator; + } #ifdef WRITE_ALL - { + { char out_filename[100]; sprintf(out_filename, "MLresult_subiter_eff_1_%d.out", sub_iter_num++); @@ -1623,78 +1661,78 @@ void iterate_efficiencies(DetectorEfficiencies& efficiencies, } } #endif + } } -} -void iterate_geo_norm(GeoData3D& norm_geo_data, - const GeoData3D& measured_geo_data, - const FanProjData& model) -{ - make_geo_data(norm_geo_data, model); - //norm_geo_data = measured_geo_data / norm_geo_data; - - const int num_axial_crystals_per_block = measured_geo_data.get_num_axial_crystals_per_block(); - const int num_transaxial_crystals_per_block = measured_geo_data.get_half_num_transaxial_crystals_per_block()*2; - - const float threshold = measured_geo_data.find_max()/10000.F; - - for (int ra = 0; ra < num_axial_crystals_per_block; ++ra) - for (int a = 0; a =threshold || measured_geo_data(ra,a,rb,b) < 10000*norm_geo_data(ra,a,rb,b))? measured_geo_data(ra,a,rb,b) / norm_geo_data(ra,a,rb,b) - : 0; - } -} + void iterate_geo_norm(GeoData3D& norm_geo_data, + const GeoData3D& measured_geo_data, + const FanProjData& model) + { + make_geo_data(norm_geo_data, model); + //norm_geo_data = measured_geo_data / norm_geo_data; -void iterate_block_norm(BlockData3D& norm_block_data, - const BlockData3D& measured_block_data, - const FanProjData& model) -{ - make_block_data(norm_block_data, model); - //norm_block_data = measured_block_data / norm_block_data; - const float threshold = measured_block_data.find_max()/10000.F; - for (int ra = norm_block_data.get_min_ra(); ra <= norm_block_data.get_max_ra(); ++ra) - for (int a = norm_block_data.get_min_a(); a <= norm_block_data.get_max_a(); ++a) - // loop rb from ra to avoid double counting - for (int rb = max(ra,norm_block_data.get_min_rb(ra)); rb <= norm_block_data.get_max_rb(ra); ++rb) - for (int b = norm_block_data.get_min_b(a); b <= norm_block_data.get_max_b(a); ++b) - { - norm_block_data(ra,a,rb,b) = - (measured_block_data(ra,a,rb,b)>=threshold || - measured_block_data(ra,a,rb,b) < 10000*norm_block_data(ra,a,rb,b)) - ? measured_block_data(ra,a,rb,b) / norm_block_data(ra,a,rb,b) - : 0; - } -} + const int num_axial_crystals_per_block = measured_geo_data.get_num_axial_crystals_per_block(); + const int num_transaxial_crystals_per_block = measured_geo_data.get_half_num_transaxial_crystals_per_block()*2; + const float threshold = measured_geo_data.find_max()/10000.F; -double KL(const FanProjData& d1, const FanProjData& d2, const double threshold) -{ - double sum=0; - for (int ra = d1.get_min_ra(); ra <= d1.get_max_ra(); ++ra) + for (int ra = 0; ra < num_axial_crystals_per_block; ++ra) + for (int a = 0; a =threshold || measured_geo_data(ra,a,rb,b) < 10000*norm_geo_data(ra,a,rb,b))? measured_geo_data(ra,a,rb,b) / norm_geo_data(ra,a,rb,b) + : 0; + } + } + + void iterate_block_norm(BlockData3D& norm_block_data, + const BlockData3D& measured_block_data, + const FanProjData& model) + { + make_block_data(norm_block_data, model); + //norm_block_data = measured_block_data / norm_block_data; + const float threshold = measured_block_data.find_max()/10000.F; + for (int ra = norm_block_data.get_min_ra(); ra <= norm_block_data.get_max_ra(); ++ra) + for (int a = norm_block_data.get_min_a(); a <= norm_block_data.get_max_a(); ++a) + // loop rb from ra to avoid double counting + for (int rb = max(ra,norm_block_data.get_min_rb(ra)); rb <= norm_block_data.get_max_rb(ra); ++rb) + for (int b = norm_block_data.get_min_b(a); b <= norm_block_data.get_max_b(a); ++b) + { + norm_block_data(ra,a,rb,b) = + (measured_block_data(ra,a,rb,b)>=threshold || + measured_block_data(ra,a,rb,b) < 10000*norm_block_data(ra,a,rb,b)) + ? measured_block_data(ra,a,rb,b) / norm_block_data(ra,a,rb,b) + : 0; + } + } + + + double KL(const FanProjData& d1, const FanProjData& d2, const double threshold) { - double asum=0; - for (int a = d1.get_min_a(); a <= d1.get_max_a(); ++a) + double sum=0; + for (int ra = d1.get_min_ra(); ra <= d1.get_max_ra(); ++ra) { - double rbsum=0; - for (int rb = max(ra,d1.get_min_rb(ra)); rb <= d1.get_max_rb(ra); ++rb) + double asum=0; + for (int a = d1.get_min_a(); a <= d1.get_max_a(); ++a) { - double bsum=0; - for (int b = d1.get_min_b(a); b <= d1.get_max_b(a); ++b) - bsum += static_cast(KL(d1(ra,a,rb,b), d2(ra,a,rb,b), threshold)); - rbsum += bsum; + double rbsum=0; + for (int rb = max(ra,d1.get_min_rb(ra)); rb <= d1.get_max_rb(ra); ++rb) + { + double bsum=0; + for (int b = d1.get_min_b(a); b <= d1.get_max_b(a); ++b) + bsum += static_cast(KL(d1(ra,a,rb,b), d2(ra,a,rb,b), threshold)); + rbsum += bsum; + } + asum += rbsum; } - asum += rbsum; + sum += asum; } - sum += asum; + return static_cast(sum); } - return static_cast(sum); -} /*double KL(const GeoData3D& d1, const GeoData3D& d2, const double threshold) { @@ -1727,6 +1765,6 @@ double KL(const FanProjData& d1, const FanProjData& d2, const double threshold) return static_cast(sum); } */ - + END_NAMESPACE_STIR diff --git a/src/utilities/apply_normfactors3D.cxx b/src/utilities/apply_normfactors3D.cxx index ff858845d..71b0ccb99 100644 --- a/src/utilities/apply_normfactors3D.cxx +++ b/src/utilities/apply_normfactors3D.cxx @@ -46,79 +46,86 @@ USING_NAMESPACE_STIR int main(int argc, char **argv) { - if (argc<7 || argc>12) + if (argc<7 || argc>12) { - std::cerr << "Usage: " << argv[0] - << " out_filename in_norm_filename_prefix measured_data apply_or_undo iter_num eff_iter_num [do_eff [ do_geo [ do_block [do_display]]]]\n" - << "apply_or_undo is 1 (multiply) or 0 (divide)\n" - << "do_eff, do_geo, do_block are 1 or 0 and all default to 1\n" - << "do_display is 1 or 0 (defaults to 0)\n"; - return EXIT_FAILURE; + std::cerr << "Usage: " << argv[0] + << " out_filename in_norm_filename_prefix measured_data apply_or_undo iter_num eff_iter_num [do_eff [ do_geo [ do_block [do_display]]]]\n" + << "apply_or_undo is 1 (multiply) or 0 (divide)\n" + << "do_eff, do_geo, do_block are 1 or 0 and all default to 1\n" + << "do_display is 1 or 0 (defaults to 0)\n"; + return EXIT_FAILURE; } - const bool do_display = argc>=11?atoi(argv[10])!=0 : false; - bool do_block = argc>=10?atoi(argv[9])!=0: true; - bool do_geo = argc>=9?atoi(argv[8])!=0: true; - bool do_eff = argc>=8?atoi(argv[7])!=0: true; - - // if (do_geo) - // error("Cannot do geometric factors in 3D yet"); - const int eff_iter_num = atoi(argv[6]); - const int iter_num = atoi(argv[5]); - const bool apply_or_undo = atoi(argv[4])!=0; - shared_ptr measured_data = ProjData::read_from_file(argv[3]); - const std::string in_filename_prefix = argv[2]; - const std::string output_file_name = argv[1]; - const std::string program_name = argv[0]; - shared_ptr out_proj_data_ptr - (new ProjDataInterfile(measured_data->get_exam_info_sptr(), - measured_data->get_proj_data_info_sptr()->create_shared_clone(), - output_file_name)); - - const int num_rings = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_rings(); - const int num_detectors_per_ring = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_detectors_per_ring(); - const int num_transaxial_blocks = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_transaxial_blocks(); - const int num_axial_blocks = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_axial_blocks(); - const int num_transaxial_crystals_per_block = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_transaxial_crystals_per_block(); - const int num_axial_crystals_per_block = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_axial_crystals_per_block(); - - GeoData3D norm_geo_data(num_axial_crystals_per_block, num_transaxial_crystals_per_block/2, num_rings, num_detectors_per_ring); - - BlockData3D norm_block_data(num_axial_blocks, num_transaxial_blocks, - num_axial_blocks-1, num_transaxial_blocks-1); - DetectorEfficiencies efficiencies(IndexRange2D(num_rings, num_detectors_per_ring)); + const bool do_display = argc>=11?atoi(argv[10])!=0 : false; + bool do_block = argc>=10?atoi(argv[9])!=0: true; + bool do_geo = argc>=9?atoi(argv[8])!=0: true; + bool do_eff = argc>=8?atoi(argv[7])!=0: true; + + // if (do_geo) + // error("Cannot do geometric factors in 3D yet"); + const int eff_iter_num = atoi(argv[6]); + const int iter_num = atoi(argv[5]); + const bool apply_or_undo = atoi(argv[4])!=0; + shared_ptr measured_data = ProjData::read_from_file(argv[3]); + const std::string in_filename_prefix = argv[2]; + const std::string output_file_name = argv[1]; + const std::string program_name = argv[0]; + shared_ptr out_proj_data_ptr + (new ProjDataInterfile(measured_data->get_exam_info_sptr(), + measured_data->get_proj_data_info_sptr()->create_shared_clone(), + output_file_name)); + + const int num_transaxial_blocks = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_transaxial_blocks(); + const int num_axial_blocks = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_axial_blocks(); + const int virtual_axial_crystals = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_virtual_axial_crystals_per_block(); + const int virtual_transaxial_crystals = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_virtual_transaxial_crystals_per_block(); + const int num_rings = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_rings() -(num_axial_blocks-1)*virtual_axial_crystals; + const int num_detectors_per_ring = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_detectors_per_ring() -num_transaxial_blocks*virtual_transaxial_crystals; + + const int num_transaxial_crystals_per_block = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_transaxial_crystals_per_block()-virtual_transaxial_crystals; + const int num_axial_crystals_per_block = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_axial_crystals_per_block()-virtual_axial_crystals; + + GeoData3D norm_geo_data(num_axial_crystals_per_block, num_transaxial_crystals_per_block/2, num_rings, num_detectors_per_ring); + + BlockData3D norm_block_data(num_axial_blocks, num_transaxial_blocks, + num_axial_blocks-1, num_transaxial_blocks-1); + DetectorEfficiencies efficiencies(IndexRange2D(num_rings, num_detectors_per_ring)); { - // efficiencies - if (do_eff) - { - char *in_filename = new char[in_filename_prefix.size() + 30]; - sprintf(in_filename, "%s_%s_%d_%d.out", - in_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); - std::ifstream in(in_filename); - in >> efficiencies; - if (!in) - { - warning("Error reading %s, using all 1s instead\n", in_filename); - do_eff = false; - } - delete[] in_filename; - } - - // geo norm + // efficiencies + if (do_eff) + { + char *in_filename = new char[in_filename_prefix.size() + 30]; + sprintf(in_filename, "%s_%s_%d_%d.out", + in_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); + std::ifstream in(in_filename); + in >> efficiencies; + if (!in) + { + warning("Error reading %s, using all 1s instead\n", in_filename); + do_eff = false; + } + delete[] in_filename; + } + + // geo norm if (do_geo) { { @@ -135,45 +142,47 @@ int main(int argc, char **argv) delete[] in_filename; } } - + // block norm - if (do_block) - { - { - char *in_filename = new char[in_filename_prefix.size() + 30]; - sprintf(in_filename, "%s_%s_%d.out", - in_filename_prefix.c_str(), "block", iter_num); - std::ifstream in(in_filename); - in >> norm_block_data; - if (!in) - { - warning("Error reading %s, using all 1s instead\n", in_filename); - do_block = false; - } - delete[] in_filename; - } - } - - { - FanProjData fan_data; - make_fan_data(fan_data, *measured_data); -// make_fan_data_remove_gaps(fan_data, *measured_data); - - if (do_eff) - apply_efficiencies(fan_data, efficiencies, apply_or_undo); - if (do_geo) - apply_geo_norm(fan_data, norm_geo_data, apply_or_undo); - if (do_block) - apply_block_norm(fan_data, norm_block_data, apply_or_undo); - - if (do_display) - display(fan_data, "input*norm"); - set_fan_data(*out_proj_data_ptr, fan_data); -// set_fan_data_add_gaps(*out_proj_data_ptr, fan_data); - - - } + if (do_block) + { + { + char *in_filename = new char[in_filename_prefix.size() + 30]; + sprintf(in_filename, "%s_%s_%d.out", + in_filename_prefix.c_str(), "block", iter_num); + std::ifstream in(in_filename); + in >> norm_block_data; + if (!in) + { + warning("Error reading %s, using all 1s instead\n", in_filename); + do_block = false; + } + delete[] in_filename; + } + } + + { + FanProjData fan_data; + + //make_fan_data(fan_data, *measured_data); + make_fan_data_remove_gaps(fan_data, *measured_data); + + + if (do_eff) + apply_efficiencies(fan_data, efficiencies, apply_or_undo); + if (do_geo) + apply_geo_norm(fan_data, norm_geo_data, apply_or_undo); + if (do_block) + apply_block_norm(fan_data, norm_block_data, apply_or_undo); + + if (do_display) + display(fan_data, "input*norm"); + //set_fan_data(*out_proj_data_ptr, fan_data);apply + set_fan_data_add_gaps(*out_proj_data_ptr, fan_data); + + + } } - return EXIT_SUCCESS; + return EXIT_SUCCESS; } diff --git a/src/utilities/find_ML_normfactors3D.cxx b/src/utilities/find_ML_normfactors3D.cxx index 396fbd8e8..1a060541f 100644 --- a/src/utilities/find_ML_normfactors3D.cxx +++ b/src/utilities/find_ML_normfactors3D.cxx @@ -105,51 +105,58 @@ int main(int argc, char **argv) shared_ptr model_data = ProjData::read_from_file(argv[3]); shared_ptr measured_data = ProjData::read_from_file(argv[2]); const std::string out_filename_prefix = argv[1]; - const int num_rings = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_rings(); - const int num_detectors_per_ring = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_detectors_per_ring(); const int num_transaxial_blocks = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_transaxial_blocks(); - const int num_axial_blocks = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_axial_blocks(); + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_transaxial_blocks(); + const int num_axial_blocks = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_axial_blocks(); + const int virtual_axial_crystals = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_virtual_axial_crystals_per_block(); + const int virtual_transaxial_crystals = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_virtual_transaxial_crystals_per_block(); + const int num_rings = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_rings() -(num_axial_blocks-1)*virtual_axial_crystals; + const int num_detectors_per_ring = + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_detectors_per_ring() -num_transaxial_blocks*virtual_transaxial_crystals; + const int num_transaxial_crystals_per_block = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_transaxial_crystals_per_block(); + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_transaxial_crystals_per_block()-virtual_transaxial_crystals; const int num_axial_crystals_per_block = - measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> - get_num_axial_crystals_per_block(); + measured_data->get_proj_data_info_sptr()->get_scanner_sptr()-> + get_num_axial_crystals_per_block()-virtual_axial_crystals; + - CPUTimer timer; timer.start(); - + FanProjData model_fan_data; FanProjData fan_data; Array<2,float> data_fan_sums(IndexRange2D(num_rings, num_detectors_per_ring)); DetectorEfficiencies efficiencies(IndexRange2D(num_rings, num_detectors_per_ring)); - + GeoData3D measured_geo_data(num_axial_crystals_per_block, num_transaxial_crystals_per_block/2, num_rings, num_detectors_per_ring ); //inputes have to be modified GeoData3D norm_geo_data(num_axial_crystals_per_block, num_transaxial_crystals_per_block/2, num_rings, num_detectors_per_ring ); //inputes have to be modified - + BlockData3D measured_block_data(num_axial_blocks, num_transaxial_blocks, num_axial_blocks-1, num_transaxial_blocks-1); BlockData3D norm_block_data(num_axial_blocks, num_transaxial_blocks, num_axial_blocks-1, num_transaxial_blocks-1); - make_fan_data(model_fan_data, *model_data); + make_fan_data_remove_gaps(model_fan_data, *model_data); { // next could be local if KL is not computed below FanProjData measured_fan_data; float threshold_for_KL; // compute factors dependent on the data { - make_fan_data(measured_fan_data, *measured_data); - + make_fan_data_remove_gaps(measured_fan_data, *measured_data); + /* TEMP FIX */ for (int ra = model_fan_data.get_min_ra(); ra <= model_fan_data.get_max_ra(); ++ra) { @@ -157,13 +164,14 @@ int main(int argc, char **argv) { for (int rb = std::max(ra,model_fan_data.get_min_rb(ra)); rb <= model_fan_data.get_max_rb(ra); ++rb) { - for (int b = model_fan_data.get_min_b(a); b <= model_fan_data.get_max_b(a); ++b) + for (int b = model_fan_data.get_min_b(a); b <= model_fan_data.get_max_b(a); ++b) if (model_fan_data(ra,a,rb,b) == 0) measured_fan_data(ra,a,rb,b) = 0; } } } + threshold_for_KL = measured_fan_data.find_max()/100000.F; //display(measured_fan_data, "measured data");