DigitalBreastTomosynthesis/Reconstruction.m

%% Author: Rodrigo de Barros Vimieiro and Aleks Prša
% Date: April, 2018; 2022
% rodrigo.vimieiro@gmail.com / alex.prsa@gmail.com
% =========================================================================
%{
% 
%     DESCRIPTION:
% 
%     The goal of this software is to present an open source reconstruction 
%     toolbox, which features the four basic types of DBT reconstruction, 
%     with the possibility of different acquisition geometries. This 
%     toolbox is intended for academic usage. Since it is an open source 
%     toolbox, researchers are welcome to contribute to the current version 
%     of the software.
% 
%     Department of Electrical and Computer Engineering, 
%     São Carlos School of Engineering, 
%     University of São Paulo, 
%     São Carlos, Brazil.
%
%     Faculty of Mathematic and Physics
%     University in Ljubljana
%     Ljubljana, Slovenia
% 
%     ---------------------------------------------------------------------
%     Copyright (C) <2018>  <Rodrigo de Barros Vimieiro>
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.    
%  
%}
% =========================================================================
%%                      Reconstruction Code                              %%

close all;clear;clc

%% Global parameters

global showinfo saveinfo animation 

showinfo = uint8(1);        % Show projection animation
saveinfo = uint8(1);        % Save reconstructed volume
animation = uint8(1);       % Graphical animation

%% GUI - Data decision

answer = questdlg(['Load Images or Create a Virtual Phantom?'], ...
    'Data selection', ...
    'Images', 'Virtual', 'Images');

% Handle response
switch answer
    case 'Virtual'
        fprintf('Creating a Virtual Phantom... \n')
        data = 2;
    case 'Images'
        answer = questdlg(['Load Clinical images, Microcalcification images or Noised images?'], ...
            'Image selection', ...
            'Clinical', 'Noised', 'Microcalcifications', 'Microcalcifications');
        switch answer
            case 'Clinical'
                fprintf('Loading Clinical images... \n')
                data = 1;
            case 'Noised'
                fprintf('Loading Noised images... \n')
                data = 1; 
            case 'Microcalcifications'
                fprintf('Loading Microcalcifications images... \n')
                data = 1;
            otherwise
                fprintf('Cancelled by user \n');
    	        return;
        end
end

%% Load components

addpath(genpath('Functions'));
addpath(genpath('Parameters'));

if(~exist('Output','dir'))
    mkdir('Output')
    saveinfo = 1;
end

addpath('Output');

if(data == 1)   %   ** Dicom data **
    
    ParameterSettings_GE

    if(~exist(['Output',filesep,'Clinical'],'dir'))
        mkdir(['Output',filesep,'Clinical'])
    end
           
    % Load projection data    
    path_User = userpath;
    path_ProjData = uigetdir(path_User); 
    
    if(path_ProjData == 0)
        fprintf('Cancelled by user \n');
    	return;
    else
        userpath(path_ProjData)       
        
        fprintf('Reading DICOMs   ')
        [dataProj,infoDicom] = readDicom(path_ProjData,parameter,answer);
        parameter.bitDepth = infoDicom(:,1).BitDepth;
    
        % Pre-process projections
        fprintf('\nPre-processing DICOMs... \n')
        [dataProj,parameter] = dataPreprocess(dataProj,parameter);
        dataProj = flip(dataProj, 2);                               % Source is on the left side

        % Transform intensity image in attenuation coefficients
        dataProj(dataProj == 0) = 1;
        dataProj = -log(dataProj./(2^16-1));
        
        if(numel(find(dataProj==inf)) ~= 0)
            error('Numeric error. Please, double-check if you input data has zeros on it')
        end

    end
    
else    %   ** Virtual data **
    
    ParameterSettings_Phantom;
    
    if(~exist(['Output',filesep,'Virtual'],'dir'))
        mkdir(['Output',filesep,'Virtual'])
    end
    addpath(['Output',filesep,'Virtual']);

    % Create virtual phantom
    data3d = single(phantom(parameter));  
    data3d(data3d<0) = eps;

    data3d = flip(data3d, 1);                               % Phantom in accordance with coordinate system
    data3d = flip(data3d, 2);                               % Source is on the left side

    % Make the projections
    if(animation || saveinfo)

        % Tomosynthesis
        figureScreenSize()
        fprintf('Creating projections   ')
        dataProj = projection(data3d,parameter,[]);
        dataProj = flip(dataProj, 2);                       % Source is on the left side

        % Save phantom projection as DICOM
        filestring = ['Output',filesep,'Virtual',filesep,'Projections'];   
            if(~exist(filestring))
                mkdir(filestring)
            end
        fprintf('\nWriting DICOM projections   ')
        writeDcmProj(dataProj,filestring,parameter);

        if(saveinfo)
            fprintf('\nSaving projections... \n')
            save(['Output',filesep,'Virtual',filesep,'Projections.mat'],'dataProj','-v7.3')
            infoDicom = [];
        else
            load Projections.mat
        end
    else
        load Projections.mat
    end

    % Transform intensity image in attenuation coefficients
    dataProj(dataProj == 0) = 0.00001;
    dataProj = -log(dataProj./(2^16-1));
    dataProj = flip(dataProj, 2);                           % Source is on the left side
       
    if(numel(find(dataProj==inf)) ~= 0)
        error('Numeric error. Please, double-check if you input data has zeros on it')
    end

end

fprintf('Starting reconstruction');

%% Set specific recon parameters

nIter = [2,3,4,8];           % Iteration to be saved (MLEM or SART)
filterType = 'FBP';          % Filter type: 'BP', 'FBP'
cutoff = 0.7;                % Percentage until cut off frequency (FBP)
filterOrder = 4;             % Butterworth filter order

%% Reconstruction methods

%                       ## Uncomment to use ##
dataRecon3d = FBP(dataProj,filterType,cutoff,filterOrder,parameter);
dataProj = flip(dataProj, 2);
dataRecon3d{1,1} = flip(dataRecon3d{1,1}, 2);
if(saveinfo)
    saveData(dataRecon3d,parameter,answer,[],infoDicom);
end

%                       ## Uncomment to use ##
% dataRecon3d = MLEM(dataProj,nIter,parameter);
% if(saveinfo)
%     saveData(dataRecon3d,parameter,answer,[],infoDicom);
% end

%                       ## Uncomment to use ##
% dataRecon3d = SART(dataProj,nIter,parameter);
% if(saveinfo)
%     saveData(dataRecon3d,parameter,answer,[],infoDicom);
% end

%                       ## Uncomment to use ##
% dataRecon3d = SIRT(dataProj,nIter,parameter);
% if(saveinfo)
%     saveData(dataRecon3d,parameter,answer,[],infoDicom);
% end

fprintf('\n\nFINISHED - Data stored in "Output" folder\n')

%% Show info

% DICOM images
if(animation && data == 1)
%{
    % Projected planes viewer
    figure("Name", 'Projections', NumberTitle='off')
    sliceViewer(dataProj);
    %figure("Name", 'Projections', NumberTitle='off')
    %montage(dataProj, 'Size', [5 5], 'DisplayRange', []);
%
    % Reconstructed slice viewer
    figure("Name", 'Reconstructed images', NumberTitle='off')
    sliceViewer(dataRecon3d{1,1});
    %figure("Name", 'Reconstructed images', NumberTitle='off')
    %montage(dataRecon3d{1,1}, 'Size', [10 round(parameter.nz/10)], 'DisplayRange', []);
%}
end

% Virtual Phantom
if(animation && data == 2)
%
    % Phantom slice viewer
    figure("Name", 'Phantom', NumberTitle='off')
    sliceViewer(data3d);
    %sliceViewer(imrotate(data3d, 90));
    %figure("Name", 'Phantom', NumberTitle='off')
    %orthosliceViewer(imrotate(data3d, 90));
    %figure("Name", 'Phantom', NumberTitle='off')
    %montage(imrotate(data3d, 90), 'Size', [10 parameter.nz/10], 'DisplayRange', []);
%
    % Projected planes viewer
    figure("Name", 'Projections', NumberTitle='off')
    sliceViewer(dataProj);
    %sliceViewer(imrotate(dataProj, 90));
    %figure("Name", 'Projections', NumberTitle='off')
    %montage(imrotate(dataProj, 90), 'Size', [3 3], 'DisplayRange', []);
%
    % Reconstructed slice viewer
    figure("Name", 'Reconstructed phantom', NumberTitle='off')
    sliceViewer(dataRecon3d{1,1});
    %sliceViewer(imrotate(dataRecon3d{1,1}, 90));
    %figure("Name", 'Phantom', NumberTitle='off')
    %orthosliceViewer(imrotate(dataRecon3d{1,1}, 90));
    %figure("Name", 'Reconstructed phantom', NumberTitle='off')
    %montage(imrotate(dataRecon3d{1,1}, 90), 'Size', [10 parameter.nz/10], 'DisplayRange', []);
%{
    % 3D Backprojection Visualization
    figureScreenSize()
    for k=1:parameter.nz
        subplot(1,2,1)
        imshow(data3d(:,:,k))
        title(['Phantom slice ',num2str(k)]);axis on;
        subplot(1,2,2)
        imshow(dataRecon3d{1,end}(:,:,k),[])   
        title(['Reconstructed slice ',num2str(k)]);axis on;
        pause(0.02)
    end
%}
end

beep