%% Author: Rodrigo de Barros Vimieiro / Aleks Prša
% Date: April, 2018 / May 2022
% rodrigo.vimieiro@gmail.com / alex.prsa@gmail.com
% =========================================================================
%{
% -------------------------------------------------------------------------
%    writeDicom(dataRecon3d,infoDicomRecon,filestring,reconmeth,parameter)
% -------------------------------------------------------------------------
%     DESCRIPTION:
%     This function write a set of Dicom images with respective headers.
%  
%     INPUT:
%     - dataRecon3d = Reconstructed data
%     - infoDicomProjs = Dicom header from projections
%     - filestring = Path to save
%     - reconmeth = Reconstructed method
%     - parameter = Parameter of all geometry
% 
%     ---------------------------------------------------------------------
%     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/>.
%}
% =========================================================================
%% Write Dicom images
function writeDicom(dataRecon3d,infoDicomRecon,filestring,reconmeth,parameter)

% Get reconstructed dicom headers
 %uiwait(msgbox('Select reconstructed image Dicom file to get headers.','Dicom','Warn'));
 %[filename,path] = uigetfile('*.dcm');
 %infoDicomRecon = dicominfo([path,filename]);

% Weather is Hologic or GE
dc=0;
if(strcmp(parameter.type,'vct')||strcmp(parameter.type,'hologic'))
    dc = -1;
end

% Adjust data to save, based on the recon method
if(strcmp(reconmeth,'BP')||strcmp(reconmeth,'MLEM')||strcmp(reconmeth,'SART')||strcmp(reconmeth,'SIRT'))
    for k=1:size(dataRecon3d,2)
        dataRecon3d{1,k} = (2^16-1).*mat2gray(dataRecon3d{1,k});
        dataRecon3d{1,k} = uint16(dataRecon3d{1,k});
    end
else
    if(strcmp(reconmeth,'FBP'))
        fprintf('\nAdjusting grayscale - windowing   ')
        for k=1:size(dataRecon3d,2)
            dataRecon3d{1,k} = dataRecon3d{1,k} + abs(min(dataRecon3d{1,k}(:)));
            dataRecon3d{1,k} = (2^16-1).*mat2gray(dataRecon3d{1,k});
            dataRecon3d{1,k} = uint16(dataRecon3d{1,k});

            % Adjusting grayscale window - windowing            
            data = dataRecon3d{1,k};
            data = reshape(data,1,numel(data));
            data = data ./ 256;
            data = fix(data);
            nmin = 0;
            nmax = 0;
            binmin = 0;
            binmax = 0;
            
            % Determining lower grayscale limit
            fprintf('\n   Lower limit')
            for bin=1:256
                for i=1:numel(data)
                    if(data(i) == (bin-1))
                        nmin = nmin + 1;
                    end
                    if(nmin > 14739305)
                        binmin = bin;
                        break
                    end
                end
                if(nmin > 14739305)
                    binmin = bin;
                    break
                end
            end

            % Determining upper grayscale limit
            fprintf('\n   Upper limit')
            for bin=1:256
                for i=1:numel(data)
                    if(data(i) == (257-bin))
                        nmax = nmax + 1;
                    end
                    if(nmax > 37331)
                        binmax = bin;
                        break
                    end
                end
                if(nmax > 37331)
                    binmax = bin;
                    break
                end
            end

            wmin = binmin*256;                                  % Lover grayscale limit
            wmax = (256-binmax)*256;                            % Upper grayscale limit

            dataRecon3d{1,k}(dataRecon3d{1,k} < wmin) = wmin;
            dataRecon3d{1,k}(dataRecon3d{1,k} > wmax) = wmax;

            %dataRecon3d{1,k}(dataRecon3d{1,k} < 32325) = 32325;
            %dataRecon3d{1,k}(dataRecon3d{1,k} > 58600) = 58600;          

            dataRecon3d{1,k} = (2^16-1).*mat2gray(dataRecon3d{1,k});
            dataRecon3d{1,k} = uint16(dataRecon3d{1,k});
        end
    else
        error('Invalid recon method');
    end
end

if(~exist(filestring))
    mkdir(filestring)
end

% Write Dicom files in respective folder
fprintf('\nWriting DICOMs   ')
fprintf('%3d/%3d', 1, parameter.nz)
for i=1:size(dataRecon3d{1,1},3)
    fprintf('\b\b\b\b\b\b\b%3d/%3d', i, parameter.nz)
    if(i<10)
        fileimg = [filestring,filesep,'0',num2str(i+dc),'.dcm'];
    else
        fileimg = [filestring,filesep,num2str(i+dc),'.dcm'];
    end
    dicomwrite(dataRecon3d{1,1}(:,:,i), fileimg, infoDicomRecon(:,1), 'CreateMode', 'copy');
end

end