From c110bd615af1b4c90afbb44793881a893c414ee0 Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Sat, 3 Mar 2018 17:59:18 +0100
Subject: [PATCH] Sinos<->blobs: added support for sub-blob UV shifts
At the moment, slabbing is disabled in the shifted mode for sinos->blobs
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Cxx/include/DctData.h | 7 +-
zUtil_Cxx/include/DctDataAlgorithms.h | 195 ++++++++++++------
zUtil_Cxx/include/DctMatlabData.h | 25 ++-
zUtil_Deformation/Gt6DVolumeToBlobProjector.m | 26 ++-
4 files changed, 189 insertions(+), 64 deletions(-)
diff --git a/zUtil_Cxx/include/DctData.h b/zUtil_Cxx/include/DctData.h
index a571a50e..5cc6adae 100644
--- a/zUtil_Cxx/include/DctData.h
+++ b/zUtil_Cxx/include/DctData.h
@@ -18,7 +18,7 @@
namespace dct {
static const size_t VERSION_MAJOR = 0;
- static const size_t VERSION_MINOR = 0;
+ static const size_t VERSION_MINOR = 1;
// This class is intentionally not memory-safe!! Be careful with it..
template<typename Type, class Alloc = VectorAllocator<Type> >
@@ -158,6 +158,11 @@ namespace dct {
static const char error_wrong_dim_0[] = "Blobs and sinograms should have same length along 1st dimension";
static const char error_wrong_dim_2[] = "Blobs and sinograms should have same length along 3rd dimension";
+ static const char error_incons_use_shifts[] = "Inconsistent use of shifts";
+
+ static const char error_wrong_dim_u_sh[] = "Sinograms should fit into blobs along 1st dimension";
+ static const char error_wrong_dim_v_sh[] = "Sinograms should fit into blobs along 3rd dimension";
+
static const char error_sino_slice_out_of_hlim[] = "Requested a SINOGRAM SLICE outside the range ('sino_offsets' too large)!";
static const char error_blob_slice_out_of_hlim[] = "Requested a BLOB SLICE outside the range ('proj_offsets' too large)!";
diff --git a/zUtil_Cxx/include/DctDataAlgorithms.h b/zUtil_Cxx/include/DctDataAlgorithms.h
index f9206855..f1d504e9 100644
--- a/zUtil_Cxx/include/DctDataAlgorithms.h
+++ b/zUtil_Cxx/include/DctDataAlgorithms.h
@@ -91,6 +91,7 @@ namespace dct {
template<typename Type>
struct DctProjectionTransformCoefficients {
size_t num_projections;
+ bool use_shifts;
std::vector<size_t> dest_slice; // <-- slice in destination sinogram / blob
@@ -99,33 +100,47 @@ namespace dct {
std::vector<Type> coeff; // <-- saxpy coefficient
+ std::vector<size_t> shift_u; // <-- U shift for sinos into blobs
+ std::vector<size_t> shift_v; // <-- V shift for sinos into blobs
+
public:
- DctProjectionTransformCoefficients() : num_projections(0) { };
- DctProjectionTransformCoefficients(const size_t & num_projs)
- : num_projections(num_projs)
+ DctProjectionTransformCoefficients() : num_projections(0), use_shifts(false) { };
+ DctProjectionTransformCoefficients(const size_t & num_projs, const bool & do_use_shifts)
+ : num_projections(num_projs), use_shifts(do_use_shifts)
{
this->resize(num_projs);
}
void resize(const size_t & num_projs) {
num_projections = num_projs;
+
dest_slice.resize(num_projs);
src_vol.resize(num_projs);
src_slice.resize(num_projs);
coeff.resize(num_projs);
+
+ shift_u.resize(num_projs);
+ shift_v.resize(num_projs);
}
void reserve(const size_t & num_projs) {
dest_slice.reserve(num_projs);
src_vol.reserve(num_projs);
src_slice.reserve(num_projs);
coeff.reserve(num_projs);
+
+ shift_u.reserve(num_projs);
+ shift_v.reserve(num_projs);
}
- void emplace_back(const size_t & d_slice, const size_t & s_vol, const size_t & s_slice, const Type & c)
+ void emplace_back(const size_t & d_slice, const size_t & s_vol, const size_t & s_slice, const Type & c, const size_t & d_shift_u = 0, const size_t & d_shift_v = 0)
{
dest_slice.emplace_back(d_slice);
src_vol.emplace_back(s_vol);
src_slice.emplace_back(s_slice);
coeff.emplace_back(c);
+
+ shift_u.emplace_back(d_shift_u);
+ shift_v.emplace_back(d_shift_v);
+
num_projections++;
}
};
@@ -142,12 +157,18 @@ namespace dct {
std::vector<DctProjectionTransformCoefficients<Type> > transforms_to_sinos;
std::vector<DctProjectionTransformCoefficients<Type> > transforms_to_blobs;
- size_t line_length;
- size_t lines_num;
+ size_t blobs_size_u;
+ size_t blobs_size_v;
+ size_t sinos_size_u;
+ size_t sinos_size_v;
+
+ bool use_shifts;
public:
DctProjectionTransform(const size_t & num_threads = std::thread::hardware_concurrency())
- : DctDataProcess<Type, Alloc, vector_size>(num_threads), line_length(0), lines_num(0)
+ : DctDataProcess<Type, Alloc, vector_size>(num_threads)
+ , blobs_size_u(0), blobs_size_v(0), sinos_size_u(0), sinos_size_v(0)
+ , use_shifts(false)
{ }
void set_blobs(const DctProjectionVector<Type, Alloc> & new_blobs);
@@ -170,8 +191,8 @@ namespace dct {
template<typename FunctionType>
void submit_transforms(FunctionType func, const size_t & tot_vols);
- void transform_blobs2sino(const size_t & sino_num, const size_t & min_line, const size_t & max_line);
- void transform_sinos2blob(const size_t & blob_num, const size_t & min_line, const size_t & max_line);
+ void transform_blobs2sino(const size_t & sino_num, const size_t & min_v, const size_t & max_v);
+ void transform_sinos2blob(const size_t & blob_num, const size_t & min_v, const size_t & max_v);
void initialize_transforms_to_blobs();
};
@@ -713,8 +734,8 @@ namespace dct {
this->blobs = new_blobs;
const auto & dims = blobs[0].get_dims();
- this->line_length = dims[0];
- this->lines_num = dims[2];
+ this->blobs_size_u = dims[0];
+ this->blobs_size_v = dims[2];
}
template<typename Type, class Alloc, const size_t vector_size>
@@ -728,8 +749,8 @@ namespace dct {
this->sinos = new_sinos;
const auto & dims = sinos[0].get_dims();
- this->line_length = dims[0];
- this->lines_num = dims[2];
+ this->sinos_size_u = dims[0];
+ this->sinos_size_v = dims[2];
}
template<typename Type, class Alloc, const size_t vector_size>
@@ -738,6 +759,11 @@ namespace dct {
const std::vector<DctProjectionTransformCoefficients<Type> > & new_transforms)
{
this->transforms_to_sinos = new_transforms;
+ this->transforms_to_blobs.clear();
+
+ if (this->transforms_to_sinos.size() > 0) {
+ this->use_shifts = transforms_to_sinos[0].use_shifts;
+ }
}
template<typename Type, class Alloc, const size_t vector_size>
@@ -751,8 +777,8 @@ namespace dct {
this->blobs = std::move(new_blobs);
const auto & dims = blobs[0].get_dims();
- this->line_length = dims[0];
- this->lines_num = dims[2];
+ this->blobs_size_u = dims[0];
+ this->blobs_size_v = dims[2];
}
template<typename Type, class Alloc, const size_t vector_size>
@@ -766,8 +792,8 @@ namespace dct {
this->sinos = std::move(new_sinos);
const auto & dims = sinos[0].get_dims();
- this->line_length = dims[0];
- this->lines_num = dims[2];
+ this->sinos_size_u = dims[0];
+ this->sinos_size_v = dims[2];
}
template<typename Type, class Alloc, const size_t vector_size>
@@ -777,6 +803,10 @@ namespace dct {
{
this->transforms_to_sinos = std::move(new_transforms);
this->transforms_to_blobs.clear();
+
+ if (this->transforms_to_sinos.size() > 0) {
+ this->use_shifts = transforms_to_sinos[0].use_shifts;
+ }
}
template<typename Type, class Alloc, const size_t vector_size>
@@ -798,6 +828,9 @@ namespace dct {
const std::vector<size_t> & blob_offsets = transforms_to_sinos[num_sino].src_vol;
const std::vector<size_t> & proj_offsets = transforms_to_sinos[num_sino].src_slice;
+ const std::vector<size_t> & shifts_u = transforms_to_sinos[num_sino].shift_u;
+ const std::vector<size_t> & shifts_v = transforms_to_sinos[num_sino].shift_v;
+
const size_t * const dims_sino = sinos[num_sino].get_dims();
for (size_t m = 0; m < transforms_to_sinos[num_sino].num_projections; m++) {
@@ -814,13 +847,34 @@ namespace dct {
}
const size_t * const dims_blob = this->blobs[b].get_dims();
- // Sino and Blob have the same size along 1st and 3rd dimensions
-
- if (dims_blob[0] != dims_sino[0]) {
- THROW(WrongArgumentException, DctProjectionTransform, error_wrong_dim_0);
+ if (transforms_to_sinos[num_sino].use_shifts != this->use_shifts) {
+ THROW(WrongArgumentException, DctProjectionTransform, error_incons_use_shifts);
}
- if (dims_blob[2] != dims_sino[2]) {
- THROW(WrongArgumentException, DctProjectionTransform, error_wrong_dim_2);
+
+ if (!transforms_to_sinos[num_sino].use_shifts) {
+ // Sino and Blob have the same size along 1st and 3rd dimensions in initial version
+
+ if (dims_blob[0] != dims_sino[0]) {
+ THROW(WrongArgumentException, DctProjectionTransform, error_wrong_dim_0);
+ }
+ if (dims_blob[2] != dims_sino[2]) {
+ THROW(WrongArgumentException, DctProjectionTransform, error_wrong_dim_2);
+ }
+ } else {
+ // In this case, all the sinos are supposed to fit into the blobs
+ const size_t & su = shifts_u[m];
+ const size_t & sv = shifts_v[m];
+
+ if (dims_blob[0] < (dims_sino[0] + su)) {
+ printf("Num sino: %lu, num blob: %lu, blobs UV size [%lu, %lu], sino UV size [%lu, %lu], sino shift UV [%lu, %lu]\n",
+ num_sino, b, dims_blob[0], dims_blob[2], dims_sino[0], dims_sino[2], su, sv);
+ THROW(WrongArgumentException, DctProjectionTransform, error_wrong_dim_u_sh);
+ }
+ if (dims_blob[2] < (dims_sino[2] + sv)) {
+ printf("Num sino: %lu, num blob: %lu, blobs UV size [%lu, %lu], sino UV size [%lu, %lu], sino shift UV [%lu, %lu]\n",
+ num_sino, b, dims_blob[0], dims_blob[2], dims_sino[0], dims_sino[2], su, sv);
+ THROW(WrongArgumentException, DctProjectionTransform, error_wrong_dim_v_sh);
+ }
}
if (s >= dims_sino[1]) {
@@ -850,6 +904,7 @@ namespace dct {
for (size_t b = 0; b < tot_blobs; b++)
{
transforms_to_blobs[b].reserve(2 * tot_sinos);
+ transforms_to_blobs[b].use_shifts = this->use_shifts;
}
for (size_t num_sino = 0; num_sino < tot_sinos; ++num_sino)
@@ -864,7 +919,13 @@ namespace dct {
const auto & c = metad_s.coeff[m];
// (d_slice, s_vol, s_slice, coefficient)
- transforms_to_blobs[b].emplace_back(p, num_sino, s, c);
+ if (metad_s.use_shifts) {
+ const auto & su = metad_s.shift_u[m];
+ const auto & sv = metad_s.shift_v[m];
+ transforms_to_blobs[b].emplace_back(p, num_sino, s, c, su, sv);
+ } else {
+ transforms_to_blobs[b].emplace_back(p, num_sino, s, c);
+ }
}
}
}
@@ -904,25 +965,32 @@ namespace dct {
for (size_t num_vol = 0; num_vol < complete_vols; num_vol++)
{
- funcs.emplace_back(std::function<void()>(std::bind(func, this, num_vol, size_t(0), lines_num)));
+ funcs.emplace_back(std::function<void()>(std::bind(func, this, num_vol, size_t(0), sinos_size_v)));
}
if (remaining_vols > 0) {
- const size_t tot_theo_slabs = this->num_threads / remaining_vols + ((this->num_threads % remaining_vols) > 0);
- const size_t tot_slabs = std::max(tot_theo_slabs, std::max(lines_num / 2, size_t(1)));
- const size_t lines_per_slab = lines_num / tot_slabs + ((lines_num % tot_slabs) > 0);
+ if (this->use_shifts && func == &DctProjectionTransform<Type, Alloc, vector_size>::transform_sinos2blob) {
+ for (size_t num_vol = complete_vols; num_vol < tot_vols; num_vol++)
+ {
+ funcs.emplace_back(std::function<void()>(std::bind(func, this, num_vol, size_t(0), sinos_size_v)));
+ }
+ } else {
+ const size_t tot_theo_slabs = this->num_threads / remaining_vols + ((this->num_threads % remaining_vols) > 0);
+ const size_t tot_slabs = std::min(tot_theo_slabs, std::max(sinos_size_v / 2, size_t(1)));
+ const size_t lines_per_slab = sinos_size_v / tot_slabs + ((sinos_size_v % tot_slabs) > 0);
- funcs.reserve(complete_vols + remaining_vols * tot_slabs);
+ funcs.reserve(complete_vols + remaining_vols * tot_slabs);
- for (size_t num_blob = complete_vols; num_blob < tot_vols; num_blob++)
- {
- size_t slab_start = 0;
- for (size_t num_slab = 0; num_slab < tot_slabs; ++num_slab) {
- const size_t slab_end = std::min(slab_start + lines_per_slab, lines_num);
+ for (size_t num_vol = complete_vols; num_vol < tot_vols; num_vol++)
+ {
+ size_t slab_start = 0;
+ for (size_t num_slab = 0; num_slab < tot_slabs; ++num_slab) {
+ const size_t slab_end = std::min(slab_start + lines_per_slab, sinos_size_v);
- funcs.emplace_back(std::function<void()>(std::bind(func, this, num_blob, slab_start, slab_end)));
+ funcs.emplace_back(std::function<void()>(std::bind(func, this, num_vol, slab_start, slab_end)));
- slab_start = slab_end;
+ slab_start = slab_end;
+ }
}
}
}
@@ -932,20 +1000,20 @@ namespace dct {
template<typename Type, class Alloc, const size_t vector_size>
void
- DctProjectionTransform<Type, Alloc, vector_size>::transform_blobs2sino(const size_t & sino_num, const size_t & min_line, const size_t & max_line)
+ DctProjectionTransform<Type, Alloc, vector_size>::transform_blobs2sino(const size_t & sino_num, const size_t & min_v, const size_t & max_v)
{
- Type ALIGNED(vector_size) temp_line[line_length];
+ Type ALIGNED(vector_size) temp_line[sinos_size_u];
- const size_t safe_min_line = std::max(min_line, size_t(0));
- const size_t safe_max_line = std::min(max_line, lines_num);
+ const size_t safe_min_v = std::max(min_v, size_t(0));
+ const size_t safe_max_v = std::min(max_v, sinos_size_v);
const DctProjectionTransformCoefficients<Type> & metad = transforms_to_sinos[sino_num];
- for (size_t line_num = safe_min_line; line_num < safe_max_line; line_num++)
+ for (size_t v = safe_min_v; v < safe_max_v; v++)
{
- Type * __restrict const sino_base = sinos[sino_num].get_data() + line_num * sinos[sino_num].get_skip();
+ Type * __restrict const sino_base = sinos[sino_num].get_data() + v * sinos[sino_num].get_skip();
- this->initialize_line(temp_line, line_length);
+ this->initialize_line(temp_line, sinos_size_u);
for (size_t m = 0; m < metad.num_projections; m++)
{
@@ -953,19 +1021,22 @@ namespace dct {
const size_t & b = metad.src_vol[m];
const size_t & s = metad.dest_slice[m];
- OpBinaryPlusScaled<Type> op(metad.coeff[m]);
+ const size_t & su = metad.shift_u[m];
+ const size_t & sv = metad.shift_v[m];
+ const size_t b_shift = this->use_shifts ? (su + sv * blobs[b].get_skip()) : 0;
- const Type * __restrict const blob_data_line = blobs[b].get_data() + line_length * p + line_num * blobs[b].get_skip();
+ OpBinaryPlusScaled<Type> op(metad.coeff[m]);
- this->transform_line_binary(temp_line, blob_data_line, line_length, op);
+ const Type * __restrict const blob_data_line = blobs[b].get_data() + blobs_size_u * p + v * blobs[b].get_skip() + b_shift;
+ this->transform_line_binary(temp_line, blob_data_line, sinos_size_u, op);
const bool is_last = (m+1) == metad.num_projections;
if (is_last || metad.dest_slice[m+1] != s) {
- Type * const __restrict sino_data_line = sino_base + line_length * s;
- this->copy_line(sino_data_line, temp_line, line_length);
+ Type * const __restrict sino_data_line = sino_base + sinos_size_u * s;
+ this->copy_line(sino_data_line, temp_line, sinos_size_u);
if (!is_last) {
- this->initialize_line(temp_line, line_length);
+ this->initialize_line(temp_line, sinos_size_u);
}
}
}
@@ -974,29 +1045,33 @@ namespace dct {
template<typename Type, class Alloc, const size_t vector_size>
void
- DctProjectionTransform<Type, Alloc, vector_size>::transform_sinos2blob(const size_t & blob_num, const size_t & min_line, const size_t & max_line)
+ DctProjectionTransform<Type, Alloc, vector_size>::transform_sinos2blob(const size_t & blob_num, const size_t & min_v, const size_t & max_v)
{
- const size_t safe_min_line = std::max(min_line, size_t(0));
- const size_t safe_max_line = std::min(max_line, lines_num);
+ const size_t safe_min_v = std::max(min_v, size_t(0));
+ const size_t safe_max_v = std::min(max_v, sinos_size_v);
const DctProjectionTransformCoefficients<Type> & metad = transforms_to_blobs[blob_num];
- for (size_t line_num = safe_min_line; line_num < safe_max_line; line_num++)
+ for (size_t v = safe_min_v; v < safe_max_v; v++)
{
- Type * __restrict const base_blob = blobs[blob_num].get_data() + line_num * blobs[blob_num].get_skip();
+ Type * const base_blob = blobs[blob_num].get_data() + v * blobs[blob_num].get_skip();
for (size_t m = 0; m < metad.num_projections; m++)
{
- const size_t & sino_slice = metad.src_slice[m];
- const size_t & sino_vol = metad.src_vol[m];
- const size_t & blob_slice = metad.dest_slice[m];
+ const size_t & s = metad.src_slice[m];
+ const size_t & sino_num = metad.src_vol[m];
+ const size_t & p = metad.dest_slice[m];
+
+ const size_t & su = metad.shift_u[m];
+ const size_t & sv = metad.shift_v[m];
+ const size_t b_shift = this->use_shifts ? (su + sv * blobs[blob_num].get_skip()) : 0;
OpBinaryPlusScaled<Type> op(metad.coeff[m]);
- const Type * __restrict const sino_data_line = sinos[sino_vol].get_data() + line_num * sinos[sino_vol].get_skip() + line_length * sino_slice;
- Type * __restrict const blob_data_line = base_blob + line_length * blob_slice;
+ const Type * __restrict const sino_data_line = sinos[sino_num].get_data() + v * sinos[sino_num].get_skip() + sinos_size_u * s;
+ Type * const blob_data_line = base_blob + blobs_size_u * p + b_shift;
- this->transform_line_binary(blob_data_line, sino_data_line, line_length, op);
+ this->transform_line_binary(blob_data_line, sino_data_line, sinos_size_u, op);
}
}
}
diff --git a/zUtil_Cxx/include/DctMatlabData.h b/zUtil_Cxx/include/DctMatlabData.h
index ae9be3a4..6ec98439 100644
--- a/zUtil_Cxx/include/DctMatlabData.h
+++ b/zUtil_Cxx/include/DctMatlabData.h
@@ -314,6 +314,9 @@ namespace dct {
const mxArray * const proj_offsets_a = mxGetField(coeff, 0, "proj_offsets");
const mxArray * const proj_coeffs_a = mxGetField(coeff, 0, "proj_coeffs");
+ const mxArray * const proj_shifts_uv_a = mxGetField(coeff, 0, "proj_shifts_uv");
+ output.use_shifts = proj_shifts_uv_a != nullptr;
+
if (!sino_offsets_a) {
THROW(WrongArgumentException, DctMatlabData,
"Invalid Coeffs: 'sino_offsets' is missing!");
@@ -344,6 +347,10 @@ namespace dct {
THROW(WrongArgumentException, DctMatlabData,
"Invalid Coeffs: number of 'proj_coeffs' is different from 'blob_offsets'!");
}
+ if (output.use_shifts && (num_elems * 2) != mxGetNumberOfElements(proj_shifts_uv_a)) {
+ THROW(WrongArgumentException, DctMatlabData,
+ "Invalid Coeffs: number of 'proj_shifts_uv' is different from 'blob_offsets'!");
+ }
if (mxDOUBLE_CLASS != mxGetClassID(sino_offsets_a)) {
THROW(WrongArgumentException, DctMatlabData,
@@ -361,6 +368,10 @@ namespace dct {
THROW(WrongArgumentException, DctMatlabData,
"Invalid Coeffs: type of 'proj_coeffs' should be 'double'!");
}
+ if (output.use_shifts && mxDOUBLE_CLASS != mxGetClassID(proj_shifts_uv_a)) {
+ THROW(WrongArgumentException, DctMatlabData,
+ "Invalid Coeffs: type of 'proj_shifts_uv' should be 'double'!");
+ }
output.resize(num_elems);
@@ -376,7 +387,19 @@ namespace dct {
std::transform(projs_pos, projs_pos + num_elems, output.src_slice.begin(), op);
const double * const coeffs = (const double *) mxGetData(proj_coeffs_a);
- std::copy(coeffs, coeffs + num_elems, output.coeff.begin());
+ std::copy_n(coeffs, num_elems, output.coeff.begin());
+
+ if (output.use_shifts) {
+ const double * const shifts_uv = (const double *) mxGetData(proj_shifts_uv_a);
+ for (size_t ind = 0; ind < num_elems; ind++) {
+ // shifts always refer to inside the blobs
+ output.shift_u[ind] = shifts_uv[ind*2];
+ output.shift_v[ind] = shifts_uv[ind*2+1];
+ }
+ } else {
+ std::fill_n(output.shift_u.begin(), num_elems, 0);
+ std::fill_n(output.shift_v.begin(), num_elems, 0);
+ }
return output;
}
diff --git a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
index 9c3266ef..61a9b188 100644
--- a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
+++ b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
@@ -145,13 +145,24 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
blob_offsets = coeffs.blob_offsets(:);
proj_offsets = coeffs.proj_offsets(:);
proj_coeffs = coeffs.proj_coeffs(:);
+
+ if (isfield(coeffs, 'proj_shifts_uv'))
+ shifts_uv = coeffs.proj_shifts_uv;
+ else
+ shifts_uv = zeros(2, numel(proj_offsets));
+ end
+
for m = 1:numel(proj_offsets)
p = proj_offsets(m);
b = blob_offsets(m);
c = proj_coeffs(m);
s = sino_offsets(m);
- blobs{b}(:, p, :) = blobs{b}(:, p, :) + c * sino(:, s, :);
+ blobs_us = (1:size(sino, 1)) + shifts_uv(1, m);
+ blobs_vs = (1:size(sino, 3)) + shifts_uv(2, m);
+
+ blobs{b}(blobs_us, p, blobs_vs) = blobs{b}(blobs_us, p, blobs_vs) + c * sino(:, s, :);
+% blobs{b}(:, p, :) = blobs{b}(:, p, :) + c * sino(:, s, :);
end
catch mexc
fprintf('\n\nError!!\n')
@@ -175,13 +186,24 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
blob_offsets = coeffs.blob_offsets(:);
proj_offsets = coeffs.proj_offsets(:);
proj_coeffs = coeffs.proj_coeffs(:);
+
+ if (isfield(coeffs, 'proj_shifts_uv'))
+ shifts_uv = coeffs.proj_shifts_uv;
+ else
+ shifts_uv = zeros(2, numel(proj_offsets));
+ end
+
for m = 1:numel(proj_offsets)
p = proj_offsets(m);
b = blob_offsets(m);
c = proj_coeffs(m);
s = sino_offsets(m);
- sino(:, s, :) = sino(:, s, :) + c * blobs{b}(:, p, :);
+ blobs_us = (1:sinogram_size(1)) + shifts_uv(1, m);
+ blobs_vs = (1:sinogram_size(3)) + shifts_uv(2, m);
+
+ sino(:, s, :) = sino(:, s, :) + c * blobs{b}(blobs_us, p, blobs_vs);
+% sino(:, s, :) = sino(:, s, :) + c * blobs{b}(:, p, :);
end
catch mexc
fprintf('\n\nError!!\n')
--
GitLab