clear all; clc; 
%% Author Statement
%Author: Scott Black, University of Strathclyde, Department of Biomedical Engineering
%Date: 26/12/2021

%Purpose: The purpose of this code is to combine several DICOM stacks on a slice-by-slice 
%basis to create a new dataset. Specifically, this script was used to combine 3 DICOM 
%stacks of human aortae which were created from 4D Flow-MRI data at discrete time steps 
%(systolic acceleration, peak systole, systolic deceleration). The resultant dataset is 
%therefore a temporal composite which combines each of these time-steps. 
%This code can be altered to combine any number of time steps.

%% Read in DICOM files
%Load in the working directory for the location of each DICOM stack of images
%Systolic acceleration
projectdir_SysAcc = 'Insert full path where DICOM stack is located'; 
%Peak systole
projectdir_Systole = 'Insert full path where DICOM stack is located';    
%Systolic Deceleration
projectdir_SysDec = 'Insert full path where DICOM stack is located';    

%Create a variable for all of the data from each DICOM stack
dicomlist_Systole = dir(fullfile(projectdir_Systole, '*.dcm') );
dicomlist_SysAcc = dir(fullfile(projectdir_SysAcc, '*.dcm') );
dicomlist_SysDec = dir(fullfile(projectdir_SysDec, '*.dcm') );

%Determine the number of images in each stack
nfile_Systole = length(dicomlist_Systole);
nfile_SysAcc = length(dicomlist_SysAcc);
nfile_SysDec = length(dicomlist_SysDec);

%Load in the images themselves
%All time steps should have the same number of images, so it does not
%matter which nfile_timestep is used in the for loop
for K = 1 : nfile_Systole 
     dimages_Systole{K} = dicomread(fullfile(projectdir_Systole, dicomlist_Systole(K).name));
     dimages_SysAcc{K} = dicomread(fullfile(projectdir_SysAcc, dicomlist_SysAcc(K).name));
     dimages_SysDec{K} = dicomread(fullfile(projectdir_SysDec, dicomlist_SysDec(K).name));
end

%Add the path ofor each DICOM stack into the workspace again
addpath('Insert full path where DICOM stack is located');    
addpath('Insert full path where DICOM stack is located');    
addpath('Insert full path where DICOM stack is located');   
%%
%To read the resultant composite image in segmentation software, each image 
%in the stack must have the same matadata. This will not change anything of
%the image itself, it just meakes it readable (i.e. in SimVascular and ITKSnap)
%I have chosen to copy the metadata from systolic acceleration into the
%final composite. However, it is not important what time step is used.
 uid1=dicomuid;
 uid2=dicomuid;

 %Alpha is the blending factor when merging images
 alpha=0.5;
for K=1:nfile_Systole
%Create new metadata
    var = dicomlist_Systole(K).name;
    metadata = dicominfo(var);
     info.InstanceNumber=K;
%Check to ensure the metadata of a few slices are the same
     info.SeriesInstanceUID=uid1;
     info.StudyInstanceUID=uid2;
%Create the file name of the composite DICOM stack
    fileName = sprintf('TemporalComposite_%d.dcm',K);
%With Imadjust(image,[]), the lower and upper limits are determined by the
%user. Althernatively, they can be removed, again at the discretion of the
%user
%Merge Systolic acceleration and peak systole together using alpha blending
%Similar to how you would use imfuse(). However, with the equations below,
%you can control the blending factor(alpha, defined above as 0.5)
    x = alpha * imadjust(dimages_Systole{K},[0.001 0.008]) + (1 - alpha) * imadjust(dimages_SysAcc{K},[0.001 0.008]);
    y=  alpha * imadjust(dimages_SysDec{K},[0.001 0.008]) + (1 - alpha) * x;
%To even out the contrast and create uniform signal intensity, the fused
%images should be averaged
    averageImage = uint16(double(x) + double(y))/2;
%Write the temporal composite image to the working directory
    dicomwrite(averageImage,fileName,metadata, 'CreateMode', 'copy');
    EnhancedImages{K}=dicomread(fileName);
end

%%
%To visualise the temporal composite images, for example at DICOM slices
%845, 890, and 950
figure
subplot(1,3,1)
imshow(EnhancedImages{845},[])
title('Slice 845')
subplot(1,3,2)
imshow(EnhancedImages{890},[])
title('Slice 890')
subplot(1,3,3)
imshow(EnhancedImages{950},[])
title('Slice 950')
set(gcf,'color','white')