
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

