-
Nicola Vigano authored
Signed-off-by:Nicola Vigano <nicola.vigano@esrf.fr>
Nicola Vigano authoredSigned-off-by:
do6D.m 15.46 KiB
if (exist('testCreator', 'var'))
delete(testCreator)
clear Gt6DTestCreator
end
if (exist('reconstructor', 'var'))
delete(reconstructor)
clear Gt6DSpotReconstructor
clear Gt6DBlobReconstructor
end
clear Gt6DVolumeProjector
astra_mex_data3d('clear');
astra_mex_algorithm('clear');
% compare_sirt = true;
compare_sirt = false;
if (~exist('mode', 'var'))
% mode = 'synthetic';
% % mode = 'grain576';
%
% mode = 'testPeterSpot2D';
% mode = 'testPeterSpot3D';
mode = 'testPeterSubBlob2D';
% mode = 'testPeterSubBlob3D';
end
% noise_extinction = true;
noise_extinction = false;
switch(mode)
case 'synthetic'
numOfRotations = 32;
testCreator = Gt6DTestCreator(6);
otherOrientations = { ...
% Rotation axis, Initial direction
{[0.5 0 1], [1 0 -0.5]}, ... % - 2 - orientations
{[-0.5 0 1], [1 0 0.5]}, ... % - 3 - orientations
{[0 0.5 1], [1 0 0]}, ... % - 4 - orientations
{[0 -0.5 1], [1 0 0]}, ... % - 5 - orientations
{[-0.25 -0.25 1], [1 1 0.5]}, ... % - 6 - orientations
% {[0.25 0.25 1], [1 1 -0.5]}, ... % - 7-9 - orientations
};
testCreator.createReferenceGeometry(numOfRotations, [40 0 0], ... Position of detector
{[0 1 0], [0 0 1]}, ... Axes of the detector
[0 0 1], ... % Rotation axis
otherOrientations );
% Projects the volumes and creates the theoretical detector images.
[spots, trueDetectorImages] = testCreator.getDetectorImages();
case 'blurred_grain'
grains = cell(1, 1);
grains{1} = load('grain046.mat');
grains{1} = grains{1}.grain046;
sample = GtSample.loadFromFile('sample.mat');
usedGeometry = grains{1}.proj.geom;
rotationAxis = [0 0 1];
angleStep = 0.1; % one degree of integration
testCreator = Gt6DTestCreator(2, perform_permute);
testCreator.createBlurredGeometry(usedGeometry, 5, rotationAxis, angleStep)
% Projects the volumes and creates the theoretical detector images. [spots, trueDetectorImages] = testCreator.getDetectorImages();
case 'testPeterSpot2D'
cd('~/data/testdata/')
% num_orientations = 855;
num_orientations = 341;
% num_orientations = 559;
if (~exist('test_data', 'var'))
fprintf('Loading data..')
% testData = load('testdata_ASTRA_3.mat');
test_data = load('testdata_ASTRA_5'); % 1 deg spread!
fprintf(' Done.\n')
test_data.gr = gtFedTestDiscGrainOri(test_data.fedpars, ...
test_data.bl, test_data.gv.dm, test_data.parameters, ...
'grid', num_orientations);
end
if ( (numel(test_data.gr)-1) ~= num_orientations )
test_data.gr = gtFedTestDiscGrainOri(test_data.fedpars, ...
test_data.bl, test_data.gv.dm, test_data.parameters, ...
'grid', num_orientations);
end
numOfIndexes = 60;
% numOfIndexes = length(testData.fedpars.loadbl);
% indexes = sort(randperm(length(testData.fedpars.loadbl), numOfIndexes));
indexes = 1:numOfIndexes;
indexes = test_data.fedpars.loadbl(indexes);
orientationNum = length(test_data.gr)-1;
testCreator = Gt6DTestCreator(orientationNum);
[spots, trueDetectorImages] ...
= testCreator.createOrientationSpreadGeometry(test_data.gr(2:end), test_data.bl, indexes, '2D');
geometries = testCreator.proj.geometries;
case 'testPeterSpot3D'
cd('~/data/testdata/')
% testData = load('testdata7.mat');
% testData = load('testdata_ASTRA_1.mat');
if (~exist('test_data', 'var'))
fprintf('Loading data..')
test_data = load('testdata_ASTRA_2'); % 1deg spread
test_data.gr = gtFedTestDiscGrainOri(test_data.fedpars, ...
test_data.bl, test_data.gv.dm, test_data.parameters, ...
'grid', 341);
fprintf(' Done.\n')
end
if ( (numel(test_data.gr)-1) ~= 341 )
test_data.gr = gtFedTestDiscGrainOri(test_data.fedpars, ...
test_data.bl, test_data.gv.dm, test_data.parameters, ...
'grid', 341);
end
% numOfIndexes = 60;
% numOfIndexes = length(testData.fedpars.loadbl);
% indexes = sort(randperm(length(testData.fedpars.loadbl), numOfIndexes));
indexes = 1:60;
indexes = test_data.fedpars.loadbl(indexes);
orientationNum = length(test_data.gr)-1;
testCreator = Gt6DTestCreator(orientationNum);
[spots, trueDetectorImages] ...
= testCreator.createOrientationSpreadGeometry(test_data.gr(2:end), test_data.bl, indexes, '3D');
geometries = testCreator.proj.geometries;
case 'testPeterSubBlob2D'
cd('~/data/testdata/')
% num_edge_orients = 4;% num_edge_orients = 5;
num_edge_orients = 6;
% num_edge_orients = 7;
% num_interp = 1;
num_interp = 2;
% num_interp = 4;
% num_interp = 6;
% num_interp = 8;
% num_interp = 12;
% num_interp = 16;
if (~exist('test_data', 'var'))
fprintf('Loading data..')
% testData = load('testdata_ASTRA_6'); % 5 deg spread! scrambled
% testData = load('testdata_ASTRA_9_2D_5deg_nojumps'); % 5 deg spread!
% testData = load('testdata_ASTRA_9_2D_5deg_jumps'); % 5 deg spread! with jumps
test_data = load('up_testdata_ASTRA_5'); % 1 deg spread!
% testData = load('testdata_ASTRA_7_2D'); % 5deg spread with jumps
% testData = load('testdata_ASTRA_8_2D_jumps'); % 5deg spread of only jumps
fprintf(' Done.\n')
sampler = GtOrientationSampling( test_data.bl, test_data.parameters, ...
test_data.gr{1}, 'fedpars', test_data.fedpars, 'verbose', true);
sampler.make_simple_grid('bcc', num_edge_orients, test_data.gv.dm);
% sampler.make_even_simple_grid('cubic', num_edge_orients, test_data.gv.dm);
% sampler.make_oversampling_simple_grid([1 2 3], 2);
end
if ( ~exist('sampler', 'var') || ~sampler.is_valid('bcc', num_edge_orients))
sampler = GtOrientationSampling( test_data.bl, test_data.parameters, ...
test_data.gr{1}, 'fedpars', test_data.fedpars, 'verbose', true);
sampler.make_simple_grid('bcc', num_edge_orients, test_data.gv.dm);
% sampler.make_even_simple_grid('cubic', num_edge_orients, test_data.gv.dm);
% sampler.make_oversampling_simple_grid([1 2 3], 2);
end
if (~isfield(test_data, 'blobondet'))
blobs_on_det = 1:numel(test_data.bl);
else
blobs_on_det = test_data.blobondet;
end
% blobs_on_det = sort(randperm(numel(test_data.bl), 60));
recon_algo_factory = Gt6DReconstructionAlgorithmFactory();
% [reconstructor, blobs] = recon_algo_factory.getBlobReconstructionAlgo( ...
% test_data.bl, blobs_on_det, sampler, '2D');
[reconstructor, blobs] = recon_algo_factory.getSubBlobReconstructionAlgo( ...
test_data.bl, blobs_on_det, sampler, '2D', num_interp);
case 'testPeterSubBlob3D'
cd('~/data/testdata/')
% num_edge_orients = 4;
% num_edge_orients = 5;
num_edge_orients = 6;
% num_edge_orients = 7;
% num_interp = 1;
num_interp = 2;
% num_interp = 16;
if (~exist('test_data', 'var'))
fprintf('Loading data..')
test_data = load('up_testdata_ASTRA_2'); % 1deg spread
% testData = load('testdata_ASTRA_2_3D_5deg_nojumps'); % 5deg spread without jumps
% testData = load('testdata_ASTRA_2_3D_5deg_jumps'); % 5deg spread with jumps
% testData = load('testdata_ASTRA_7'); % 5deg spread with jumps
fprintf(' Done.\n')
sampler = GtOrientationSampling( test_data.bl, test_data.parameters, ...
test_data.gr{1}, 'fedpars', test_data.fedpars, 'verbose', true);
sampler.make_simple_grid('bcc', num_edge_orients, test_data.gv.dm);
% sampler.make_even_simple_grid('cubic', num_edge_orients, test_data.gv.dm);
% sampler.make_oversampling_simple_grid([1 2 3], 2);
end
if ( ~exist('sampler', 'var') || ~sampler.is_valid('bcc', num_edge_orients))
sampler = GtOrientationSampling( test_data.bl, test_data.parameters, ...
test_data.gr{1}, 'fedpars', test_data.fedpars, 'verbose', true);
sampler.make_simple_grid('bcc', num_edge_orients, test_data.gv.dm);
% sampler.make_even_simple_grid('cubic', num_edge_orients, test_data.gv.dm);
% sampler.make_oversampling_simple_grid([1 2 3], 2);
end
if (~isfield(test_data, 'blobondet'))
blobs_on_det = 1:numel(test_data.bl);
else
blobs_on_det = test_data.blobondet;
end
recon_algo_factory = Gt6DReconstructionAlgorithmFactory();
[reconstructor, blobs] = recon_algo_factory.getSubBlobReconstructionAlgo( ...
test_data.bl, blobs_on_det, sampler, '3D', num_interp);
case 'peterDataset'
cd('/mntdirect/_data_id19_degisher/mi1026/id18/AlLi2b_vert_ ')
return
case 'grain576'
cd('/data/id19/graintracking/titanium/ti2007/ti2007_2_AB_/')
return
end
% Initialize reconstruction object
if (~ismember(mode, {'testPeterSubBlob2D', 'testPeterSubBlob3D'}))
zero_vols = testCreator.getInitialVols();
end
if (compare_sirt)
if (ismember(mode, {'testPeterSubBlob2D', 'testPeterSubBlob3D'}))
numBlobs = length(test_data.bl);
sizeBlobs = size(test_data.bl(1).intm);
trueSirtDetect = zeros([sizeBlobs(1:2) numBlobs]);
for ii = 1:numBlobs
trueSirtDetect(:, :, ii) = sum(test_data.bl(ii).intm, 3);
end
sirtDetectImages = permute(trueSirtDetect, [1 3 2]);
else
sirtDetectImages = spots;
end
av_gr = gtFedTestDiscGrainOri(test_data.fedpars, ...
test_data.bl, test_data.gv.dm, test_data.parameters, ...
'mean', num_orientations);
proj = [];
proj.stack = sirtDetectImages;
proj.num_cols = size(sirtDetectImages, 1);
proj.num_rows = size(sirtDetectImages, 3);
proj.geom = double(av_gr{2}.allblobs.proj_geom(1:60, :));
[proj.vol_size_x, proj.vol_size_y, proj.vol_size_z] = size(zero_vols{1});
proj.num_iter = 30;
tic(); reconSIRT_avg = gtAstra3D(0, 0, 0, proj);
fprintf('Reconstructed using SIRT, with %d iterations in %f seconds\n', ...
proj.num_iter, toc());
proj.geom = double(av_gr{1}.allblobs.proj_geom(1:60, :));
tic(); reconSIRT_00 = gtAstra3D(0, 0, 0, proj);
fprintf('Reconstructed using SIRT, with %d iterations in %f seconds\n', ...
proj.num_iter, toc());
end
gtGetMaxDisorientation(test_data.gv.dm, [], 'diameter_average');
gtGetMaxDisorientation(test_data.gv.dm, [], 'diameter_zero');
gtGetMaxDisorientation(test_data.gv.dm, [], 'minimal');
fprintf('Ready to reconstruct\n');
% try
if (ismember(mode, {'testPeterSubBlob2D', 'testPeterSubBlob3D'}))
fprintf('Signal''s Energy %f, Power %f\n', gtGetSignalEnergy(blobs), gtGetSignalPower(blobs))
% semi_angular_deviation = 0
% angular_sampling = 1
% % semi_angular_deviation = max(sampler.gr_stats.proj.max_angle_x, sampler.gr_stats.proj.max_angle_y) / 2;
% semi_angular_deviation_deg = semi_angular_deviation * 180 / pi
% % angular_sampling = 3
% sampler.gr_stats.sample.avg_dist * semi_angular_deviation
% pixel_dist = max(sampler.gr_stats.proj.max_pixel_dist_x, sampler.gr_stats.proj.max_pixel_dist_y);
% fprintf('SubRay distance: %f, total spanned distance %f\n', pixel_dist / angular_sampling, pixel_dist)
%
% reconstructor = Gt6DBlobReconstructor(zero_vols, blobs, ...
% 'geometries', geometries, ...
% 'offsets', offsets, ...
% 'ss_geometries', geometries_ss, ...
% 'ss_offsets', offsets_ss, ...
% 'ss_coeffs', coeffs_ss, ...
% 'angular_sampling', angular_sampling, ...
% 'angular_deviation', semi_angular_deviation * 0.7 );
reconstructor.verbose = true;
else
reconstructor = Gt6DSpotReconstructor(zero_vols, spots, ...
'geometries', geometries, ...
'ss_geometries', geometries_ss, ...
'ss_coeffs', coeffs_ss, ...
'verbose', false);
end
% catch mexc
% gtPrintException(mexc)
% return
% end
if (~noise_extinction)
% Perform reconstruction
t = tic();
% try
if (ismember(mode, {'testPeterSubBlob2D', 'testPeterSubBlob3D'}))
% numIters = 50;
% reconstructor.sirt(numIters, 2e-4);
numIters = 30;
reconstructor.cp(numIters, 1e-3);
% numIters = 500;
% reconstructor.cpTV(numIters, 1e-3);
% numIters = 100;
% reconstructor.cpTVl1(numIters, 1e-5);
else
% numIters = 12;
% reconstructor.cg(numIters, 1e-3);
%
numIters = 20;
reconstructor.cp(numIters, 1e-3);
% numIters = 50;
% reconstructor.sirt(numIters, 1e-3);
end
% catch mexc
% disp(mexc)
% % gtPrintException(mexc)
% return % end
toc(t);
if (strcmpi(mode, 'grain576'))
final_vols = reconstructor.getMergeMultipleJunction(0.7);
end
if (ismember(mode, {'testPeterSpot3D', 'testPeterSpot2D', 'testPeterSubBlob2D', 'testPeterSubBlob3D'}))
result = gtCSAssembleResult(test_data, sampler, reconstructor.getCurrentSolution());
post_result = gtCSPostProcessOrientationSpread(test_data, result, true);
result_viewer = GtOrientationResultView(test_data, result, post_result);
else
result.solution = reconstructor.getCurrentSolution();
end
else
noise_dims = cell(size(blobs));
for ii = 1:numel(blobs)
noise_dims{ii} = size(blobs{ii});
end
signal_energy = gtGetSignalEnergy(blobs);
% n = gtGenerateExtintionNoise(blobs, [0.1 0.15], 80);
n = gtGenerateExtintionNoise(blobs, [0.1 0.15], 0);
% n_white = gtGenerateAdditiveWhiteNoise(noise_dims, signal_energy / 10, 'single');
cryst = test_data.parameters.cryst;
thetatypes = [test_data.bl(blobs_on_det).thetatype];
n_white = gtGenerateNonUniformAdditiveWhiteNoise(noise_dims, thetatypes, cryst, signal_energy / 10, 'single');
fprintf('Signal''s Energy %f, Power %f\n', signal_energy, gtGetSignalPower(blobs))
fprintf('Extinction Noise''s Energy %f, Power %f\n', gtGetSignalEnergy(n), gtGetSignalPower(n))
fprintf('Additive White Noise''s Energy %f, Power %f\n', gtGetSignalEnergy(n_white), gtGetSignalPower(n_white))
noisy_blobs = cell(size(blobs));
for ii = 1:numel(blobs)
noisy_blobs{ii} = blobs{ii} + n{ii} + n_white{ii};
end
reconstructor.reInit(noisy_blobs); % , reconstructor.getCurrentSolution()
numIters = 50;
reconstructor.cp(numIters, 1e-1);
% numIters = 50;
% reconstructor.sirt(numIters, 1e-3);
result_noisy = gtCSAssembleResult(test_data, sampler, reconstructor.getCurrentSolution());
post_result_noisy = gtCSPostProcessOrientationSpread(test_data, result_noisy, true);
result_n_viewer = GtOrientationResultView(test_data, result_noisy, post_result_noisy);
end