diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
index 2465838c51e4131ba99030ddb81e6b8bf6320c64..55410bbccad6c17f91bafae69de8c009819d03e9 100644
--- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
+++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
@@ -47,6 +47,39 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
good_orients = algo_params.good_orients;
blobs = algo_params.blobs;
+ % merging all the orientation-space super-sampling projection
+ % coefficients into one structure per orientation. This is needed in
+ % order to fwd/bwd-project all the different sub-orientations
+ % together.
+ if (~isempty(algo_params.geometries_ss))
+ for ii_o = 1:numel(algo_params.offsets_ss)
+ offsets_ss = algo_params.offsets_ss{ii_o};
+ new_offsets_ss = offsets_ss{1};
+ cumulative_offs = numel(new_offsets_ss.proj);
+
+ for ii_ss = 2:numel(offsets_ss)
+ new_offsets_ss.sino_offsets ...
+ = cat(2, new_offsets_ss.sino_offsets, offsets_ss{ii_ss}.sino_offsets + cumulative_offs);
+ new_offsets_ss.blob_offsets ...
+ = cat(2, new_offsets_ss.blob_offsets, offsets_ss{ii_ss}.blob_offsets);
+ new_offsets_ss.proj_offsets ...
+ = cat(2, new_offsets_ss.proj_offsets, offsets_ss{ii_ss}.proj_offsets);
+ new_offsets_ss.proj_coeffs ...
+ = cat(2, new_offsets_ss.proj_coeffs, offsets_ss{ii_ss}.proj_coeffs);
+
+ num_projs = numel(offsets_ss{ii_ss}.proj);
+ for ii_p = 1:num_projs
+ offsets_ss{ii_ss}.proj(ii_p).sino_offset ...
+ = offsets_ss{ii_ss}.proj(ii_p).sino_offset + cumulative_offs;
+ end
+ new_offsets_ss.proj = cat(2, new_offsets_ss.proj, offsets_ss{ii_ss}.proj);
+
+ cumulative_offs = cumulative_offs + num_projs;
+ end
+ algo_params.offsets_ss{ii_o} = new_offsets_ss;
+ end
+ end
+
algo = Gt6DBlobReconstructor(volume_size, blobs, ...
'data_type', self.data_type, ...
'geometries', algo_params.geometries, ...
@@ -139,6 +172,8 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
end
offsets = cat(1, offsets, orient_offsets);
+ % Shifting projection to the concatenated blobs, unless that
+ % blob was shared between parent and twin
if (~isempty(geometries_ss))
orient_offsets_ss = algo_params(ii_g).offsets_ss;
for ii_o = 1:numel(orient_offsets_ss)
diff --git a/zUtil_Deformation/Gt6DVolumeProjector.m b/zUtil_Deformation/Gt6DVolumeProjector.m
index 18fddd9932281235fcf69505f6d16f68ae14aa1c..8ca5b2491a9f872a23179e7b617d5ff5ce6500c4 100644
--- a/zUtil_Deformation/Gt6DVolumeProjector.m
+++ b/zUtil_Deformation/Gt6DVolumeProjector.m
@@ -52,14 +52,6 @@ classdef Gt6DVolumeProjector < handle
if (~self.use_astra_projectors)
self.volume_id = astra_mex_data3d('create', '-vol', self.astra_volume_geometry, 0);
- else
- if (~isempty(self.ss_geometries))
- try
- self.use_astra_oss_projectors = astra_mex_direct_c('supports', 'FP3D_OSS') ...
- && astra_mex_direct_c('supports', 'BP3D_OSS');
- catch
- end
- end
end
try
@@ -143,19 +135,14 @@ classdef Gt6DVolumeProjector < handle
'GPUindex', -1 );
for g_ii = 1:num_ss_geoms
- num_sub_ss_geoms = numel(self.ss_geometries{g_ii});
- self.astra_ss_projection_geometries{g_ii} = cell(num_sub_ss_geoms, 1);
+ ss_geom = cat(1, self.ss_geometries{g_ii}{:});
+
+ self.astra_ss_projection_geometries{g_ii} = astra_create_proj_geom(...
+ 'parallel3d_vec', self.proj_size(2), self.proj_size(1), ss_geom);
+
if (self.use_astra_projectors)
- self.astra_ss_projectors{g_ii} = cell(num_sub_ss_geoms, 1);
- end
- for g_ss_ii = 1:num_sub_ss_geoms
- self.astra_ss_projection_geometries{g_ii}{g_ss_ii} = astra_create_proj_geom(...
- 'parallel3d_vec', self.proj_size(2), self.proj_size(1), self.ss_geometries{g_ii}{g_ss_ii});
-
- if (self.use_astra_projectors)
- self.astra_ss_projectors{g_ii}{g_ss_ii} = astra_create_projector('cuda3d', ...
- self.astra_ss_projection_geometries{g_ii}{g_ss_ii}, self.astra_volume_geometry, opts);
- end
+ self.astra_ss_projectors{g_ii} = astra_create_projector('cuda3d', ...
+ self.astra_ss_projection_geometries{g_ii}, self.astra_volume_geometry, opts);
end
end
end
@@ -171,74 +158,54 @@ classdef Gt6DVolumeProjector < handle
methods (Access = protected)
function ss_sinogram = fwd_project_single_volume_ss(self, volume, ii)
- ss_geoms = self.astra_ss_projection_geometries{ii};
- num_ss_geoms = numel(ss_geoms);
if (self.use_astra_projectors)
if (numel(volume) == 1)
volume = volume * ones(self.volume_geometry, 'single');
end
- ss_projectors = self.astra_ss_projectors{ii};
else
astra_mex_data3d('store', self.volume_id, volume);
end
- if (self.use_astra_oss_projectors)
- ss_sinogram = astra_mex_direct_c('FP3D_OSS', [ss_projectors{:}], volume);
+ num_ss_geoms = numel(self.ss_geometries{ii});
+ volume = volume ./ num_ss_geoms;
+ if (self.use_astra_projectors)
+ ss_sinogram = astra_mex_direct_c('FP3D', self.astra_ss_projectors{ii}, volume);
else
- volume = volume ./ num_ss_geoms;
- ss_sinogram = cell(num_ss_geoms, 1);
- for ii_ss = 1:num_ss_geoms
- if (self.use_astra_projectors)
- ss_sinogram{ii_ss} = astra_mex_direct_c('FP3D', ss_projectors{ii_ss}, volume);
- else
- sino_id = astra_mex_data3d('create', '-proj3d', ss_geoms{ii_ss}, 0);
-
- cfg = astra_struct('FP3D_CUDA');
- cfg.ProjectionDataId = sino_id;
- cfg.VolumeDataId = self.volume_id;
-
- algo_id = astra_mex_algorithm('create', cfg);
- astra_mex_algorithm('iterate', algo_id);
- astra_mex_algorithm('delete', algo_id);
-
- ss_sinogram{ii_ss} = astra_mex_data3d('get_single', sino_id);
- astra_mex_data3d('delete', sino_id);
- end
- end
+ sino_id = astra_mex_data3d('create', '-proj3d', self.astra_ss_projection_geometries{ii}, 0);
+
+ cfg = astra_struct('FP3D_CUDA');
+ cfg.ProjectionDataId = sino_id;
+ cfg.VolumeDataId = self.volume_id;
+
+ algo_id = astra_mex_algorithm('create', cfg);
+ astra_mex_algorithm('iterate', algo_id);
+ astra_mex_algorithm('delete', algo_id);
+
+ ss_sinogram = astra_mex_data3d('get_single', sino_id);
+ astra_mex_data3d('delete', sino_id);
end
end
function volume = bwd_project_single_volume_ss(self, ss_sinogram, ii)
- ss_geoms = self.astra_ss_projection_geometries{ii};
- num_ss_geoms = numel(ss_geoms);
if (self.use_astra_projectors)
- ss_projectors = self.astra_ss_projectors{ii};
- end
-
- if (self.use_astra_oss_projectors)
- volume = astra_mex_direct_c('BP3D_OSS', [ss_projectors{:}], ss_sinogram);
+ volume = astra_mex_direct_c('BP3D', self.astra_ss_projectors{ii}, ss_sinogram);
else
- volume = cell(num_ss_geoms, 1);
- for ii_ss = 1:num_ss_geoms
- if (self.use_astra_projectors)
- volume{ii_ss} = astra_mex_direct_c('BP3D', ss_projectors{ii_ss}, ss_sinogram{ii_ss});
- else
- sino_id = astra_mex_data3d('create', '-proj3d', ss_geoms{ii_ss}, ss_sinogram{ii_ss});
-
- cfg = astra_struct('BP3D_CUDA');
- cfg.ProjectionDataId = sino_id;
- cfg.ReconstructionDataId = self.volume_id;
-
- algo_id = astra_mex_algorithm('create', cfg);
- astra_mex_algorithm('iterate', algo_id);
- astra_mex_algorithm('delete', algo_id);
-
- volume{ii_ss} = astra_mex_data3d('get_single', self.volume_id);
- astra_mex_data3d('delete', sino_id);
- end
- end
- volume = gtMathsSumCellVolumes(volume) ./ num_ss_geoms;
+ sino_id = astra_mex_data3d('create', '-proj3d', self.astra_ss_projection_geometries{ii}, ss_sinogram);
+
+ cfg = astra_struct('BP3D_CUDA');
+ cfg.ProjectionDataId = sino_id;
+ cfg.ReconstructionDataId = self.volume_id;
+
+ algo_id = astra_mex_algorithm('create', cfg);
+ astra_mex_algorithm('iterate', algo_id);
+ astra_mex_algorithm('delete', algo_id);
+
+ volume = astra_mex_data3d('get_single', self.volume_id);
+ astra_mex_data3d('delete', sino_id);
end
+
+ num_ss_geoms = numel(self.ss_geometries{ii});
+ volume = volume ./ num_ss_geoms;
end
function sinogram = fwd_project_single_volume(self, volume, geom_idx)
diff --git a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
index f926af49bf9324a67419d3fe7fe71b2544304a7b..2a16699a01f3b7c59d014068dd35d474b59f95c8 100644
--- a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
+++ b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
@@ -115,14 +115,13 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
c = tic();
ss_offsets_n = self.ss_offsets{n};
- for ii = 1:numel(sinos_ss)
- if (self.blob_sinogram_c_functions)
- blobs = gt6DSingleSinoToBlobs_c(...
- blobs, sinos_ss{ii}, ss_offsets_n{ii}, self.num_threads);
- else
- blobs = self.single_sinogram_to_blobs(...
- blobs, sinos_ss{ii}, ss_offsets_n{ii});
- end
+
+ if (self.blob_sinogram_c_functions)
+ blobs = gt6DSingleSinoToBlobs_c(...
+ blobs, sinos_ss, ss_offsets_n, self.num_threads);
+ else
+ blobs = self.single_sinogram_to_blobs(...
+ blobs, sinos_ss, ss_offsets_n);
end
bs_time = toc(c);
end
@@ -132,24 +131,18 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
% difference that in this case the input is blobs and we want to
% transform them into stack of sinograms
- geoms_ss = self.ss_geometries{n};
ss_offsets_n = self.ss_offsets{n};
- num_sub_ss_geoms = numel(geoms_ss);
- sinos_ss = cell(num_sub_ss_geoms, 1);
-
c = tic();
- for ii = 1:num_sub_ss_geoms
- sinogram_size = [ ...
- self.proj_size(1), size(geoms_ss{ii}, 1), self.proj_size(2)];
-
- if (self.blob_sinogram_c_functions)
- sinos_ss{ii} = gt6DBlobsToSingleSino_c( ...
- sinogram_size, blobs, ss_offsets_n{ii}, self.num_threads);
- else
- sinos_ss{ii} = self.blobs_to_single_sinogram(...
- sinogram_size, blobs, ss_offsets_n{ii});
- end
+ sinogram_size = [ ...
+ self.proj_size(1), numel(ss_offsets_n.proj), self.proj_size(2)];
+
+ if (self.blob_sinogram_c_functions)
+ sinos_ss = gt6DBlobsToSingleSino_c( ...
+ sinogram_size, blobs, ss_offsets_n, self.num_threads);
+ else
+ sinos_ss = self.blobs_to_single_sinogram(...
+ sinogram_size, blobs, ss_offsets_n);
end
bs_time = toc(c);