diff --git a/zUtil_Deformation/Gt6DBlobReconstructor.m b/zUtil_Deformation/Gt6DBlobReconstructor.m
index 6b0c4693cfaee9a89055efc41b3f58d294873cf2..b57a06f966d3b4f1f79e6fae894a343a6ea85958 100644
--- a/zUtil_Deformation/Gt6DBlobReconstructor.m
+++ b/zUtil_Deformation/Gt6DBlobReconstructor.m
@@ -115,7 +115,7 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
timing_bp = 0;
for n = 1:numel(self.currentSolution)
cbp = tic();
- v = self.bck_project_single_volumes(p, n);
+ v = self.bwd_project_single_volumes(p, n);
timing_bp = timing_bp + toc(cbp);
[self.currentSolution{n}, nextEnhancedSolution{n}] ...
@@ -396,24 +396,24 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
function initializeWeights(self)
c = tic();
fprintf(' + Projecting ones vols..');
- self.forwardWeights = self.getRowsSum();
+ self.fwd_weights = self.getRowsSum();
fprintf('\b\b (%3.1f s), Fixing it..', toc(c));
c = tic();
- for ii = 1:numel(self.forwardWeights)
- self.forwardWeights{ii} = 1 ./ (self.forwardWeights{ii} + (self.forwardWeights{ii} == 0));
+ for ii = 1:numel(self.fwd_weights)
+ self.fwd_weights{ii} = 1 ./ (self.fwd_weights{ii} + (self.fwd_weights{ii} == 0));
end
fprintf('\b\b (%2.1f s)\n', toc(c));
fprintf(' + Backprojecting ones sinos (fake bproj)..');
c = tic();
for ii = 1:numel(self.geometries)
- self.backwardWeights{ii} = 1 / size(self.geometries{ii}, 1);
+ self.bwd_weights{ii} = 1 / size(self.geometries{ii}, 1);
end
fprintf('\b\b (%2.1f s)\n', toc(c));
end
function [sigma1, sigma1_1, tau] = initializeCPWeights(self)
- sigma1 = self.forwardWeights;
+ sigma1 = self.fwd_weights;
for n = 1:numel(sigma1)
sigma1{n} = sigma1{n};
end
@@ -423,14 +423,14 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
sigma1_1{n} = 1 ./ (1 + sigma1_1{n});
end
- tau = self.backwardWeights;
+ tau = self.bwd_weights;
for n = 1:numel(tau)
tau{n} = - 1 ./ (1 ./ tau{n} + 1);
end
end
function [sigma1, sigma1_1, tau] = initializeCPWeightsTV(self)
- sigma1 = self.forwardWeights;
+ sigma1 = self.fwd_weights;
for n = 1:numel(sigma1)
sigma1{n} = sigma1{n};
end
@@ -440,7 +440,7 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
sigma1_1{n} = 1 ./ (1 + sigma1_1{n});
end
- tau = self.backwardWeights;
+ tau = self.bwd_weights;
for n = 1:numel(tau)
tau{n} = - 1 ./ (1 ./ tau{n} + 1 / numel(tau));
end
diff --git a/zUtil_Deformation/Gt6DSpotReconstructor.m b/zUtil_Deformation/Gt6DSpotReconstructor.m
deleted file mode 100755
index b72e25de66ad9c767233fef78909fe4164438de4..0000000000000000000000000000000000000000
--- a/zUtil_Deformation/Gt6DSpotReconstructor.m
+++ /dev/null
@@ -1,475 +0,0 @@
-classdef Gt6DSpotReconstructor < Gt6DVolumeToSpotProjector
- properties
- spots;
- moderated_spots;
-
- currentSolution = {};
-
- % CG specific variables
- currentResidual;
- currentCondResidual;
-
- currentCorrection;
- currentProjCorrection;
- squareProjCorr;
-
- % FISTA params
- lambda = 1e-3;
- weightFISTA = 1;
-
- % NonNegative param
- non_negative = true;
-
- % Info parameters
- dotResiduals;
- detectorResiduals;
- maxVals;
- minVals;
- thresholds;
- alpha;
- beta;
-
- % Used to display stuff
- verbose = true;
- end
-
- methods (Access = public)
- function self = Gt6DSpotReconstructor(initialVols, spots, varargin)
- fprintf('Initializing SPOT Recontruction:\n - Variables..');
- c = tic();
-
- detectorSize = size(spots);
- self = self@Gt6DVolumeToSpotProjector(initialVols, ...
- [detectorSize(1) detectorSize(3)], varargin{:});
-
- self.statistics.add_task('sirt_apply_sum', 'SIRT Apply time');
- self.statistics.add_task('sirt_fista_sum', 'SIRT Fista time');
- self.statistics.add_task('sirt_bbox_sum', 'SIRT Bbox cut time');
- self.statistics.add_task('sirt_norm_residual', 'SIRT resiual norm time');
-
- self.initProjectionGeometry();
- self.initProjectionGeometrySS();
-
- self.spots = spots;
-
- self.currentSolution = initialVols;
-
- fprintf('\b\b (%2.1f s)\n - Weights:\n', toc(c));
- c = tic();
- self.initializeWeights();
- fprintf(' Tot (%3.1f s)), Done.\n', toc(c));
- end
-
- function cg(self, numIters, lambda)
- if (exist('lambda', 'var'))
- self.lambda = lambda;
- end
- fprintf('Initialization: ')
- self.initializeVariables(numIters);
-
- if (self.verbose)
- fprintf('(lambda: %e) Done.\nReconstruction: ', self.lambda)
- f1 = figure('name', sprintf('CG: Lambda(%e), Iterations(%d)', ...
- self.lambda, numIters));
- ax1 = subplot(2, 2, 1, 'Parent', f1);
- ax2 = subplot(2, 2, 2, 'Parent', f1);
- ax3 = subplot(2, 2, 3, 'Parent', f1);
- drawnow;
- end
-
- t = tic();
- % Initial setup of variables needed for the first iteration
- self.initializeCG();
- % Iterations loop
- for ii = 1:numIters
- fprintf('%05d/%05d', ii, numIters)
-
- try
- self.makeCGIteration(ii, numIters);
- catch mexc
- gtPrintException(mexc, 'Cleaning Gradient Memory');
- self.initializeCG();
- end
-
- if (self.verbose)
- self.visualize(ax1, ax2, ax3, ii);
- drawnow;
- end
-
- fprintf('\b\b\b\b\b\b\b\b\b\b\b')
- end
- fprintf('(%04d) Done in %f seconds.\n', numIters, toc(t) )
- self.printStats();
- end
-
- function sirt(self, numIters, lambda)
- if (exist('lambda', 'var'))
- self.lambda = lambda;
- end
- fprintf('Initialization: ')
- self.initializeVariables(numIters);
-
- if (self.verbose)
- fprintf('(lambda: %e) Done.\nReconstruction: ', self.lambda)
- f1 = figure('name', sprintf('SIRT: Lambda(%e), Iterations(%d)', ...
- self.lambda, numIters));
- ax1 = subplot(2, 2, 1, 'Parent', f1);
- ax2 = subplot(2, 2, 2, 'Parent', f1);
- ax3 = subplot(2, 2, 3, 'Parent', f1);
- drawnow;
- end
-
- t = tic();
- elapsedTime = toc(t);
- % Initial setup of variables needed for the first iteration
- self.currentResidual = self.computeResiduals();
- % Iterations loop
- for ii = 1:numIters
- numchars = fprintf('%03d/%03d (in %f s)', ii, numIters, elapsedTime);
-
- self.makeSIRTIteration(ii);
-
- if (self.verbose)
- self.visualize(ax1, ax2, ax3, ii);
- drawnow;
- end
-
- fprintf('%s', repmat(sprintf('\b'), 1, numchars))
- elapsedTime = toc(t);
- end
- fprintf('(%04d) Done in %f seconds.\n', numIters, toc(t) )
- self.printStats();
- end
-
- function solution = getCurrentSolution(self)
- % solution: x = P^(1/2) u
- solution = internal_cell_div_copy(self.currentSolution, self.backwardWeights);
- end
-
- function finalVols = getMergeMultipleJunction(self, threshold)
- indices = cellfun(@(x){x > threshold}, self.currentSolution);
- finalVols = zeros(size(self.currentSolution{1}));
-
- numVols = length(self.currentSolution);
- for ii = 1:numVols
- first = true;
- for jj = 1:numVols
- if ii ~= jj
- if (first)
- indices{ii} = indices{ii} & (self.currentSolution{ii} >= self.currentSolution{jj});
- else
- indices{ii} = indices{ii} & (self.currentSolution{ii} > self.currentSolution{jj});
- end
- end
- first = false;
- end
- finalVols = finalVols + indices{ii} * ii;
- end
- end
-
- function finalVols = getMergeTwins(self, threshold)
- indices1 = (self.currentSolution{1} > self.currentSolution{2}) ...
- & (self.currentSolution{1} > threshold * max(self.currentSolution{1}(:)));
- indices2 = (self.currentSolution{2} >= self.currentSolution{1}) ...
- & (self.currentSolution{2} > threshold * max(self.currentSolution{2}(:)));
- finalVols = indices1 + 2 * indices2;
- end
-
- function cp(self, numIters, lambda)
- self.initializeVariables(numIters);
- sample_rate = 5;
-
- [sigma1, sigma1_1, tau] = self.initializeCPWeights();
- theta = 1;
-
- p = gtMathsGetSameSizeZeros(self.spots);
- q = gtMathsGetSameSizeZeros(self.currentSolution);
-
- nextEnhancedSolution = self.currentSolution;
-
- fprintf('Reconstruction: ')
-
- c = tic();
- switch (numel(numIters))
- case 1
- iters = 1:numIters;
- case 2
- iters = numIters(1):numIters(2);
- otherwise
- iters = numIters;
- end
- tot_iters = numel(iters);
- for ii = iters
- numchars = fprintf('%03d/%03d (in %f s)', ii-min(iters), tot_iters, toc(c));
-
- l = self.projectVolume(nextEnhancedSolution);
- p = (p + (l - self.spots) .* sigma1) .* sigma1_1;
-
- q = internal_cell_sum_assign(q, nextEnhancedSolution);
- for n = 1:numel(q)
- q{n} = lambda * q{n} ./ max(lambda, abs(q{n}));
- end
-
- nextSolution = self.backprojectVolume(p);
- nextSolution = internal_cell_sum_assign(nextSolution, q);
- nextSolution = internal_cell_prod_assign(nextSolution, tau);
- nextSolution = internal_cell_sum_assign(nextSolution, self.currentSolution);
-
- nextSolution = internal_cell_non_neg_assign(nextSolution);
-
- nextEnhancedSolution = internal_cell_sum_FISTA_scaled_copy(nextSolution, self.currentSolution, theta);
-
- self.currentSolution = nextSolution;
-
- fprintf('%s', repmat(sprintf('\b'), 1, numchars))
-
- if (mod(ii, sample_rate) == 0)
- l = self.projectVolume(self.currentSolution);
- l = l - self.spots;
- self.dotResiduals(ii / sample_rate) = ...
- gtMathsNorm_l2(l) / 2 ...
- + lambda * gtMathsNorm_l1(self.currentSolution) ...
- + gtMathsNorm_l2(p) / 2 ...
- + gtMathsDotProduct(p, self.spots);
- end
- end
- fprintf('(%04d) Done in %f seconds.\n', tot_iters, toc(c) )
-
- self.printStats();
-
- figure;
- subplot(1, 2, 1), plot(self.dotResiduals(1:floor(max(iters) / sample_rate)));
- subplot(1, 2, 2), semilogy(abs(self.dotResiduals));
- end
- end
-
- methods (Access = protected)
- function initializeWeights(self)
- c = tic();
- fprintf(' + Projecting ones vols..');
- self.forwardWeights = self.getRowsSum();
- fprintf('\b\b (%3.1f s), Fixing it..', toc(c));
- c = tic();
- self.forwardWeights = 1 ./ sqrt(self.forwardWeights + (self.forwardWeights == 0));
- fprintf('\b\b (%2.1f s)\n', toc(c));
-
- fprintf(' + Backrojecting ones sinos..');
- self.backwardWeights = self.getColumnsSum();
- fprintf('\b\b (%3.1f s), Fixing it..', toc(c));
- c = tic();
- numGeoms = length(self.backwardWeights);
- for ii = 1:numGeoms
- self.backwardWeights{ii} = 1 ./ sqrt(self.backwardWeights{ii} + (self.backwardWeights{ii} == 0));
- end
- fprintf('\b\b (%2.1f s)\n', toc(c));
-
- self.moderated_spots = self.spots .* self.forwardWeights;
- end
-
- function initializeVariables(self, numIters)
- self.dotResiduals = zeros(numIters, 1);
- self.detectorResiduals = zeros(numIters, 1);
- self.squareProjCorr = zeros(numIters, 1);
- self.maxVals = zeros(numIters, 1);
- self.minVals = zeros(numIters, 1);
- self.alpha = zeros(numIters, 1);
- self.beta = zeros(numIters, 1);
- end
-
- function w = getFISTAWeight(self, t_1)
- w = (self.weightFISTA - 1) / t_1;
- end
-
- function initializeCG(self)
- self.currentResidual = self.computeResiduals(self.currentSolution);
- self.currentCondResidual = self.currentResidual;
-
- self.currentCorrection = internal_cell_prod_scalar_copy(self.currentResidual, -1);
-
- self.currentProjCorrection = self.project(self.currentCorrection);
- end
-
- function [sigma1, sigma1_1, tau] = initializeCPWeights(self)
- sigma1 = self.forwardWeights .^ 2;
- sigma1_1 = 1 ./ (1 + sigma1);
-
- tau = self.backwardWeights;
- for n = 1:numel(tau)
- tau{n} = - 1 ./ (1 ./ (tau{n} .^ 2) + 1);
- end
- end
-
- function [residual, detectorResidual] = computeResiduals(self, newVolumes)
- if (~exist('newVolumes', 'var'))
- newVolumes = self.currentSolution;
- end
- % projectedVols = W * u
- % Where: W = Q^(1/2) A P^(1/2) and u = P^(-1/2) x
- projectedVols = self.project(newVolumes);
-
- % detectorResidual = W u - f
- % Where: f = Q^(1/2) b
- detectorResidual = projectedVols - self.moderated_spots;
-
- % residual = W^T (f - W u)
- residual = self.backproject(detectorResidual);
- end
-
- function makeCGIteration(self, ii, totIters)
- % Square of projected correction -> ||W p_{k}||^2
- backprojCorr = self.backproject(self.currentProjCorrection);
- self.squareProjCorr(ii) = gtMathsDotProduct(self.currentCorrection, backprojCorr);
-% self.squareProjCorr(ii) = gtMathsSquareNorm(self.currentProjCorrection);
-
- % Dot product of current residual and corrected residual
- % -> r_{k}^{T} t_{k}
- self.dotResiduals(ii) = gtMathsSquareNorm(...
- self.currentResidual);
-
- alphaPrefactor = length(self.currentSolution) * numel(self.currentSolution{1}) / totIters;
-% alphaPrefactor = 1;
-
- % Alpha_{k} - Quenched
- self.alpha(ii) = self.dotResiduals(ii) / self.squareProjCorr(ii) * alphaPrefactor;
- self.alpha(ii) = min(self.alpha(ii), 1);
- self.alpha(ii) = max(self.alpha(ii), 0);
-
- % x_{k+1} = r_{k} + alpha_{k} p_{k}
- nextSolution = internal_cell_sum_scaled_copy( ...
- self.currentSolution, self.currentCorrection, self.alpha(ii));
-
- % FISTA - regularization
- nextSolution = self.softThreshold(nextSolution, self.lambda);
-
- % NonNegative constrains
- nextSolution = gtMathsNonNegative(nextSolution);
-
- % FISTA and "non_negative" change the volumes, so we should update
- % the residuals (and corrections) as a consequence.
- nextResidual = self.computeResiduals(nextSolution);
-
- % squareNewRes = ||r_{k+1}||^2
- dotNextResisuals = gtMathsSquareNorm(nextResidual);
-
- % beta_{k+1} = ||r_{k+1}||^2 / ||r_{k}||^2
- self.beta(ii) = dotNextResisuals / self.dotResiduals(ii);
- self.beta(ii) = min(self.beta(ii), 1.0);
- self.beta(ii) = max(self.beta(ii), 0);
-
- nextCorrection = internal_cell_prod_scalar_copy(nextResidual, -1);
- nextCorrection = internal_cell_sum_scaled_assign(nextCorrection, self.currentCorrection, self.beta(ii));
-
- % W p_{k+1}
- self.currentProjCorrection = self.project(nextCorrection);
-
- % Next iteration assignments
- self.currentResidual = nextResidual;
- self.currentCorrection = nextCorrection;
-
- % Save Volumes x_{k}
- self.currentSolution = nextSolution;
-
- self.maxVals(ii) = mean(cellfun(@(x)max(x(:)), nextSolution));
- self.minVals(ii) = mean(cellfun(@(x)min(x(:)), nextSolution));
- end
-
- function visualize(self, ax1, ax2, ax3, ii)
- if ((mod(ii, 10) == 0) || (ii < 25) ...
- || ((ii < 50) && (mod(ii, 2) == 0)) ...
- || ((ii < 80) && (mod(ii, 5) == 0)))
- semilogy(abs(self.dotResiduals(1:ii)), 'Parent', ax1);
- hold(ax1, 'on')
- plot(abs(self.alpha(1:ii)), 'r', 'Parent', ax1);
- plot(self.squareProjCorr(1:ii), 'g', 'Parent', ax1);
- plot(abs(self.beta(1:ii)), 'y', 'Parent', ax1);
- hold(ax1, 'off')
- plot(self.maxVals(1:ii), 'r', 'Parent', ax3);
- hold(ax3, 'on')
- plot(self.minVals(1:ii), 'g', 'Parent', ax3);
- hold(ax3, 'off')
- end
- if (mod(ii, 3) == 0)
- slice_indx = round(size(self.currentSolution{1}, 3) / 2);
- imagesc(self.currentSolution{1}(:, :, slice_indx), 'Parent', ax2)
- end
- end
-
- function makeSIRTIteration(self, ii)
- numCells = numel(self.currentSolution);
-
- % Step length for FISTA
- self.statistics.tic('sirt_norm_residual');
- self.dotResiduals(ii) ...
- = internal_cell_square_norm(self.currentResidual);
- self.statistics.toc('sirt_norm_residual');
-
- self.statistics.tic('sirt_apply_sum');
- if (self.verbose)
- nextCorrection = internal_cell_prod_scalar_copy(self.currentResidual, -1);
- self.currentProjCorrection = self.project(nextCorrection);
- self.squareProjCorr(ii) = gtMathsSquareNorm(self.currentProjCorrection);
- nextSolution = internal_cell_sum_copy(self.currentSolution, nextCorrection);
- else
- nextSolution = internal_cell_sub_copy(self.currentSolution, self.currentResidual);
- end
- self.statistics.toc('sirt_apply_sum');
-
- % fista and nonnegative;
- % For theoretical reasons, step is always 1
- self.statistics.tic('sirt_fista_sum');
- nextSolution = self.applyFISTA(nextSolution, self.currentSolution);
- self.statistics.toc('sirt_fista_sum');
-
- self.statistics.tic('sirt_bbox_sum');
-% if (self.non_negative)
-% nextSolution = gtMathsNonNegative(nextSolution);
-% end
-
- sizes = size(nextSolution{1});
- for n = 1:numCells
- nextSolution{n}([1:3 sizes(1)-2:sizes(1)], :, :) = 0;
- nextSolution{n}(:, [1:3 sizes(2)-2:sizes(2)], :) = 0;
- if (sizes(3) == 8)
- nextSolution{n}(:, :, [1 2 3 6 7 8]) = 0;
- else
- nextSolution{n}(:, :, [1 sizes(3)]) = 0;
- end
- end
- self.statistics.toc('sirt_bbox_sum');
-
- nextResidual = self.computeResiduals(nextSolution);
-
- self.currentSolution = nextSolution;
- self.currentResidual = nextResidual;
-
- self.maxVals(ii) = mean(cellfun(@(x)max(x(:)), nextSolution));
- self.minVals(ii) = mean(cellfun(@(x)min(x(:)), nextSolution));
- end
-
- function nextSolution = softThreshold(self, nextSolution, thresh)
-% waveletLevel = 1;
-% nextSolution = dwt_haar_soft_threshold(nextSolution, waveletLevel, thresh);
- for ii = 1:numel(nextSolution)
- nextSolution{ii} = gtSoftThreshold2(nextSolution{ii}, thresh);
- end
- end
-
- function nextSolution = applyFISTA(self, nextSolution, currentSolution)
- nextSolution = softThreshold(self, nextSolution, self.lambda);
-
- newWeightFISTA = Gt6DSpotReconstructor.getTplus1(self.weightFISTA);
- fistaCurrentWeight = self.getFISTAWeight(newWeightFISTA);
- self.weightFISTA = newWeightFISTA;
-% nextSolution = internal_cell_sum_FISTA_scaled_assign(...
-% nextSolution, currentSolution, fistaCurrentWeight);
- nextSolution = internal_cell_sum_FISTA_scaled_non_neg_assign(...
- nextSolution, currentSolution, fistaCurrentWeight);
- end
- end
-
- methods (Static, Access = public)
- function t_1 = getTplus1(t)
- t_1 = (1 + sqrt(1+ 4 * t * t)) / 2;
- end
- end
-end
diff --git a/zUtil_Deformation/Gt6DVolumeProjector.m b/zUtil_Deformation/Gt6DVolumeProjector.m
index c9cf81a7cf1794c534eedcc1feb675582b3b3288..c26fb640e38a138d64ebd593290775bd969ad65c 100644
--- a/zUtil_Deformation/Gt6DVolumeProjector.m
+++ b/zUtil_Deformation/Gt6DVolumeProjector.m
@@ -9,7 +9,7 @@ classdef Gt6DVolumeProjector < handle
astra_volume_geometries = {};
geometries = {};
- astra_geometries = {};
+ astra_projection_geometries = {};
ss_coeffs = {};
ss_geometries = {};
@@ -29,52 +29,43 @@ classdef Gt6DVolumeProjector < handle
end
methods (Access = public)
- function self = Gt6DVolumeProjector(initialVolumes, detectorSize, varargin)
-% self.checkGPU()
+ function self = Gt6DVolumeProjector(init_volumes, detector_size, varargin)
- if (~exist('initialVolumes', 'var') || isempty(initialVolumes))
+ if (~exist('init_volumes', 'var') || isempty(init_volumes))
error('Gt6DVolumeProjector:wrong_argument', ...
'No initial volumes specified')
end
- numVolumes = length(initialVolumes);
-
- self.statistics.add_task('fproj_raw_sum', 'Raw FP time');
- self.statistics.add_task('bproj_raw_sum', 'Raw BP time');
-
- self.statistics.add_task_partial('fproj_raw_sum', 'copy_input', 'FP input creation/copy time');
- self.statistics.add_task_partial('bproj_raw_sum', 'copy_input', 'BP input creation/copy time');
- self.statistics.add_task_partial('fproj_raw_sum', 'compute', 'FP computation time');
- self.statistics.add_task_partial('bproj_raw_sum', 'compute', 'BP computation time');
- self.statistics.add_task_partial('fproj_raw_sum', 'copy_output', 'FP output get/copy time');
- self.statistics.add_task_partial('bproj_raw_sum', 'copy_output', 'BP output get/copy time');
+ num_volumes = length(init_volumes);
self.statistics.add_task('fproj_raw_ss_single', 'Raw SS FP time (single)');
self.statistics.add_task('bproj_raw_ss_single', 'Raw SS BP time (single)');
- for ii = 2:numVolumes
- self.testVolumeGeometry(initialVolumes{1}, initialVolumes{ii});
+ [vol_size(1), vol_size(2), vol_size(3)] = size(init_volumes{1});
+ for ii = 2:num_volumes
+ if (~all(vol_size == size(init_volumes{ii})))
+ error('PROJECTOR:wrong_argument', ...
+ 'Size mismatch between the old and the new volume');
+ end
end
% volume geometry (x, y, z)
- self.volume_geometries = cell(numVolumes, 1);
- self.astra_volume_geometries = cell(numVolumes, 1);
-
- volSize = [size(initialVolumes{1}, 1), ...
- size(initialVolumes{1}, 2), size(initialVolumes{1}, 3)];
- tempVolGeom = astra_create_vol_geom(volSize(1), volSize(2), volSize(3));
- for n = 1:numVolumes
- self.volume_geometries{n} = volSize;
- self.astra_volume_geometries{n} = tempVolGeom;
+ self.volume_geometries = cell(num_volumes, 1);
+ self.astra_volume_geometries = cell(num_volumes, 1);
+
+ astra_vol_geom = astra_create_vol_geom(vol_size(1), vol_size(2), vol_size(3));
+ for n = 1:num_volumes
+ self.volume_geometries{n} = vol_size;
+ self.astra_volume_geometries{n} = astra_vol_geom;
end
- self.volume_ids = astra_mex_data3d('create', '-vol', self.astra_volume_geometries, initialVolumes);
+ self.volume_ids = astra_mex_data3d('create', '-vol', self.astra_volume_geometries, init_volumes);
- if (~exist('detectorSize', 'var'))
- diagonal = round(norm(size(initialVolumes{1})));
+ if (~exist('detector_size', 'var'))
+ diagonal = round(norm(size(init_volumes{1})));
self.proj_size = [diagonal diagonal];
else
- self.proj_size = detectorSize;
+ self.proj_size = detector_size;
end
self = parse_pv_pairs(self, varargin);
@@ -101,8 +92,9 @@ classdef Gt6DVolumeProjector < handle
% center
% Let's clean up previous geometries
- if (~isempty(self.astra_geometries))
- self.astra_geometries = {};
+ if (~isempty(self.astra_projection_geometries))
+ self.astra_projection_geometries = {};
+
astra_mex_data3d('delete', self.sinogram_ids{:});
self.sinogram_ids = {};
astra_mex_algorithm('delete', self.algo_fproj_ids{:});
@@ -111,25 +103,15 @@ classdef Gt6DVolumeProjector < handle
self.algo_bproj_ids = {};
end
- geometry = self.geometries;
-
- if (iscell(geometry))
- num_geoms = numel(geometry);
- for n = 1:num_geoms
- self.astra_geometries{n} = astra_create_proj_geom(...
- 'parallel3d_vec', self.proj_size(2), self.proj_size(1), geometry{n});
- end
- self.sinogram_ids = astra_mex_data3d('create', '-proj3d', self.astra_geometries);
- else
- self.astra_geometries{end+1} = astra_create_proj_geom(...
- 'parallel3d_vec', self.proj_size(2), self.proj_size(1), geometry);
- self.sinogram_ids = astra_mex_data3d('create', '-proj3d', self.astra_geometries);
+ num_geoms = numel(self.geometries);
+ for n = 1:num_geoms
+ self.astra_projection_geometries{n} = astra_create_proj_geom(...
+ 'parallel3d_vec', self.proj_size(2), self.proj_size(1), self.geometries{n});
end
+ self.sinogram_ids = astra_mex_data3d('create', '-proj3d', self.astra_projection_geometries);
% Let's preallocate the algorithms
cfg = astra_struct('FP3D_CUDA');
-% cfg.option.DetectorSuperSampling = self.angular_sampling;
-% cfg.option.AngularDeviation = self.angular_deviation;
for ii = 1:numel(self.sinogram_ids)
cfg.ProjectionDataId = self.sinogram_ids{ii};
cfg.VolumeDataId = self.volume_ids{ii};
@@ -140,8 +122,6 @@ classdef Gt6DVolumeProjector < handle
end
cfg = astra_struct('BP3D_CUDA');
-% cfg.option.VoxelSuperSampling = self.angular_sampling;
-% cfg.option.AngularDeviation = self.angular_deviation;
for ii = 1:numel(self.sinogram_ids)
cfg.ProjectionDataId = self.sinogram_ids{ii};
cfg.ReconstructionDataId = self.volume_ids{ii};
@@ -171,70 +151,9 @@ classdef Gt6DVolumeProjector < handle
function stats = get_statistics(self)
stats = self.statistics;
end
-
- function forwardWeights = getRowsSum(self)
- numGeoms = numel(self.astra_volume_geometries);
- ones_cell = arrayfun(@(z){1}, 1:numGeoms);
-
- forwardWeights = self.projectVolume(ones_cell);
- end
-
- function backwardWeights = getColumnsSum(self)
- backwardWeights = self.backprojectVolume(1);
- end
end
methods (Access = protected)
- function sinograms = projectVolume(self, volumes)
- numGeoms = numel(self.astra_geometries);
- if (numGeoms ~= numel(volumes))
- error('Gt6DVolumeProjector:wrong_argument', ...
- 'Number of volumes mismatch with allocated buffers')
- end
- self.statistics.tic('fproj_raw_sum');
-
- self.statistics.tic('fproj_raw_sum', 'copy_input');
- astra_mex_data3d('store', self.volume_ids, volumes);
- self.statistics.toc('fproj_raw_sum', 'copy_input');
-
- self.statistics.tic('fproj_raw_sum', 'compute');
- for ii = 1:numGeoms
- astra_mex_algorithm('iterate', self.algo_fproj_ids{ii});
- end
- self.statistics.toc('fproj_raw_sum', 'compute');
-
- self.statistics.tic('fproj_raw_sum', 'copy_output');
- sinograms = astra_mex_data3d('get_single', self.sinogram_ids);
- self.statistics.toc('fproj_raw_sum', 'copy_output');
-
- self.statistics.toc('fproj_raw_sum');
- end
-
- function volumes = backprojectVolume(self, sinograms)
- numGeoms = numel(self.astra_geometries);
- if (numGeoms ~= numel(sinograms))
- error('Gt6DVolumeProjector:backprojVol:wrong_argument', ...
- 'Sinograms should be either a volume, or a cell of volumes')
- end
- self.statistics.tic('bproj_raw_sum');
-
- self.statistics.tic('bproj_raw_sum', 'copy_input');
- astra_mex_data3d('store', self.sinogram_ids, sinograms);
- self.statistics.toc('bproj_raw_sum', 'copy_input');
-
- self.statistics.tic('bproj_raw_sum', 'compute');
- for ii = 1:numGeoms
- astra_mex_algorithm('iterate', self.algo_bproj_ids{ii});
- end
- self.statistics.toc('bproj_raw_sum', 'compute');
-
- self.statistics.tic('bproj_raw_sum', 'copy_output');
- volumes = astra_mex_data3d('get_single', self.volume_ids);
- self.statistics.toc('bproj_raw_sum', 'copy_output');
-
- self.statistics.toc('bproj_raw_sum');
- end
-
function ss_sinogram = projectSingleVolumeSS(self, volumes, ii)
self.statistics.tic('fproj_raw_ss_single');
@@ -293,18 +212,10 @@ classdef Gt6DVolumeProjector < handle
sinogram = astra_mex_data3d('get_single', self.sinogram_ids{geom_idx});
end
- function volume = bck_project_single_volume(self, sinogram, geom_idx)
+ function volume = bwd_project_single_volume(self, sinogram, geom_idx)
astra_mex_data3d('store', self.sinogram_ids{geom_idx}, sinogram);
astra_mex_algorithm('iterate', self.algo_bproj_ids{geom_idx});
volume = astra_mex_data3d('get_single', self.volume_ids{geom_idx});
end
-
- function testVolumeGeometry(~, referenceVol, newVol)
- if (~isequal(size(referenceVol), size(newVol)))
- mexc = MException('PROJECTOR:wrong_argument', ...
- 'Size mismatch between the old and the new volume');
- throw(mexc)
- end
- end
end
end
diff --git a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
index 0696b9c0e0f3c51c0b791c1631512250415bda99..de880eebc2d9de6756e9caf9d180e22b46335a7c 100644
--- a/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
+++ b/zUtil_Deformation/Gt6DVolumeToBlobProjector.m
@@ -5,8 +5,8 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
blobs_depths = [];
% Weights
- forwardWeights = {};
- backwardWeights = {};
+ fwd_weights = {};
+ bwd_weights = {};
blob_sinogram_c_functions = true;
end
@@ -17,25 +17,19 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
self.blobs_depths = blobs_depths;
- self.statistics.add_task('fproj_create_sum', 'FP Create Blob time');
- self.statistics.add_task('bproj_create_sum', 'BP Create Blob time');
-
- self.statistics.add_task('fproj_permute_sum', 'FP Permute Blob time');
- self.statistics.add_task('bproj_permute_sum', 'BP Permute Blob time');
-
self.statistics.add_task('fproj_permute_ss_single', 'SS FP Permute Blob time (single)');
self.statistics.add_task('bproj_permute_ss_single', 'SS BP Permute Blob time (single)');
end
- function forwardWeights = getRowsSum(self)
- forwardWeights = gtMathsGetSameSizeZeros(self.blobs);
+ function fwd_weights = getRowsSum(self)
+ fwd_weights = gtMathsGetSameSizeZeros(self.blobs);
for n = 1:numel(self.geometries);
- forwardWeights = self.fwd_project_single_volumes( ...
- forwardWeights, 1, n);
+ fwd_weights = self.fwd_project_single_volumes( ...
+ fwd_weights, 1, n);
end
end
- function backwardWeights = getColumnsSum(self)
+ function bwd_weights = getColumnsSum(self)
num_blobs = numel(self.blobs_depths);
blobs = cell(1, num_blobs);
for n = 1:num_blobs
@@ -44,50 +38,14 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
end
num_geoms = numel(self.geometries);
- backwardWeights = cell(1, num_geoms);
+ bwd_weights = cell(1, num_geoms);
for n = 1:num_geoms
- backwardWeights{n} = self.bck_project_single_volumes(blobs, n);
+ bwd_weights{n} = self.bwd_project_single_volumes(blobs, n);
end
end
end
methods (Access = protected)
- function blobs = projectVolume(self, volumes)
- % We get volumes in input and we want the blobs as output, but the
- % output from the base function is stacks of sinograms relative to
- % the geometry of the volumes, so we need to actually decompose the
- % output and recompose the slices into blobs
-
- sinograms = self.projectVolume@Gt6DVolumeProjector(volumes);
-
- numBlobs = numel(self.blobs_depths);
- self.statistics.tic('fproj_create_sum');
- blobs = cell(1, numBlobs);
- for m = 1:numBlobs
- blobs{m} = zeros(self.proj_size(1), ...
- self.blobs_depths(m), self.proj_size(2), 'single');
- end
- self.statistics.toc('fproj_create_sum');
-
- self.statistics.tic('fproj_permute_sum');
- for n = 1:numel(sinograms)
- if (self.blob_sinogram_c_functions)
-% blobs = internal_cell_single_sinogram_to_blobs(...
-% blobs, sinograms{n}, self.offsets{n});
- blobs = gt6DSingleSinoToBlobs_c(...
- blobs, sinograms{n}, self.offsets{n});
- else
- blobs = self.single_sinogram_to_blobs( ...
- blobs, sinograms{n}, self.offsets{n});
- end
- end
- self.statistics.toc('fproj_permute_sum');
-
- if (~isempty(self.ss_offsets))
- blobs = self.projectVolumeSS(blobs, volumes);
- end
- end
-
function blobs = projectVolumeSS(self, blobs, volumes)
% We get volumes in input and we want the blobs as output, but the
% output from the base function is stacks of sinograms relative to
@@ -115,42 +73,6 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
blobs = internal_cell_prod_scalar_assign(blobs, renorm_factor);
end
- function volumes = backprojectVolume(self, blobs)
- % Similar problem to the projection function, with the only
- % difference that in this case the input is blobs and we want to
- % transform them into stack of sinograms
-
- numSinograms = numel(self.astra_geometries);
- self.statistics.tic('bproj_create_sum');
- sinograms_sizes = cell(1, numSinograms);
- for n = 1:numSinograms
- sinograms_sizes{n} = [self.proj_size(1), ...
- size(self.geometries{n}, 1), self.proj_size(2)];
- end
- sinograms = cell(1, numSinograms);
- self.statistics.toc('bproj_create_sum');
-
- self.statistics.tic('bproj_permute_sum');
- for n = 1:numSinograms
- if (self.blob_sinogram_c_functions)
-% sinograms{n} = internal_cell_blobs_to_single_sinogram( ...
-% sinograms_sizes{n}, blobs, self.offsets{n});
- sinograms{n} = gt6DBlobsToSingleSino_c( ...
- sinograms_sizes{n}, blobs, self.offsets{n});
- else
- sinograms{n} = self.blobs_to_single_sinogram(...
- sinograms_sizes{n}, blobs, self.offsets{n});
- end
- end
- self.statistics.toc('bproj_permute_sum');
-
- volumes = self.backprojectVolume@Gt6DVolumeProjector(sinograms);
-
- if (~isempty(self.ss_offsets))
- volumes = self.backprojectVolumeSS(volumes, blobs);
- end
- end
-
function volumes = backprojectVolumeSS(self, volumes, blobs)
% Similar problem to the projection function, with the only
% difference that in this case the input is blobs and we want to
@@ -212,7 +134,7 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
end
end
- function volume = bck_project_single_volumes(self, blobs, n)
+ function volume = bwd_project_single_volumes(self, blobs, n)
% Similar problem to the projection function, with the only
% difference that in this case the input is blobs and we want to
% transform them into stack of sinograms
@@ -230,7 +152,7 @@ classdef Gt6DVolumeToBlobProjector < Gt6DVolumeProjector
sinogram_size, blobs, self.offsets{n});
end
- volume = self.bck_project_single_volume(sino, n);
+ volume = self.bwd_project_single_volume(sino, n);
end
end
diff --git a/zUtil_Deformation/Gt6DVolumeToSpotProjector.m b/zUtil_Deformation/Gt6DVolumeToSpotProjector.m
index c2d9ca994d370fc7645aa0908a4ec4b26af74a2a..255bc96bb9a343aeecd3c99065d5681467961662 100644
--- a/zUtil_Deformation/Gt6DVolumeToSpotProjector.m
+++ b/zUtil_Deformation/Gt6DVolumeToSpotProjector.m
@@ -1,68 +1,39 @@
classdef Gt6DVolumeToSpotProjector < Gt6DVolumeProjector
properties
% Weights
- forwardWeights = {};
- backwardWeights = {};
+ fwd_weights = {};
+ bwd_weights = {};
end
methods (Access = public)
function self = Gt6DVolumeToSpotProjector(initialVolumes, detectorSize, varargin)
self = self@Gt6DVolumeProjector(initialVolumes, detectorSize, varargin{:});
-
- self.statistics.add_task('fproj_sum', 'FP Spot Sum time');
- self.statistics.add_task('bproj_dupl', 'BP Spot Duplication time');
-
- self.statistics.add_task('fproj_pre_weight', 'FP Pre-weight time');
- self.statistics.add_task('bproj_pre_weight', 'BP Pre-weight time');
-
- self.statistics.add_task('fproj_post_weight', 'FP Post-weight time');
- self.statistics.add_task('bproj_post_weight', 'BP Post-weight time');
end
- end
- methods (Access = protected)
- function sinograms = projectVolume(self, volumes)
- sinograms = self.projectVolume@Gt6DVolumeProjector(volumes);
-
- self.statistics.tic('fproj_sum');
- sinograms = gtMathsSumCellVolumes(sinograms);
- self.statistics.toc('fproj_sum');
+ function fwd_weights = getRowsSum(self)
+ fwd_weights = gtMathsGetSameSizeZeros(self.sinogram);
+ for n = 1:numel(self.geometries);
+ fwd_weights = self.fwd_project_single_volumes( ...
+ fwd_weights, 1, n);
+ end
end
- function volumes = backprojectVolume(self, sinogram)
- self.statistics.tic('bproj_dupl');
- numGeoms = numel(self.astra_geometries);
- sinograms = cell(numGeoms, 1);
- for n = 1:numGeoms
- sinograms{n} = sinogram;
+ function bwd_weights = getColumnsSum(self)
+ num_geoms = numel(self.geometries);
+ bwd_weights = cell(1, num_geoms);
+ for n = 1:num_geoms
+ bwd_weights{n} = self.bwd_project_single_volumes(1, n);
end
- self.statistics.toc('bproj_dupl');
-
- volumes = self.backprojectVolume@Gt6DVolumeProjector(sinograms);
end
+ end
- function y = project(self, x)
- self.statistics.tic('fproj_pre_weight');
- y = internal_cell_prod_copy(x, self.backwardWeights);
- self.statistics.toc('fproj_pre_weight');
-
- y = self.projectVolume(y);
-
- self.statistics.tic('fproj_post_weight');
- y = y .* self.forwardWeights;
- self.statistics.toc('fproj_post_weight');
+ methods (Access = protected)
+ function sinogram = fwd_project_single_volumes(self, sinogram, volume, n)
+ sinogram = sinogram + self.fwd_project_single_volume(volume, n);
end
- function y = backproject(self, x)
- self.statistics.tic('bproj_pre_weight');
- y = x .* self.forwardWeights;
- self.statistics.toc('bproj_pre_weight');
-
- y = self.backprojectVolume(y);
-
- self.statistics.tic('bproj_post_weight');
- y = internal_cell_prod_assign(y, self.backwardWeights);
- self.statistics.toc('bproj_post_weight');
+ function volume = bwd_project_single_volumes(self, sinogram, n)
+ volume = self.bwd_project_single_volume(sinogram, n);
end
end
end
\ No newline at end of file