%% Load file
% Loading PD-L1.tif
[PDL1_file, PDL1_path] = uigetfile('*.tif'); % PD-L1_1.tif
tiffread(PDL1_file)
PDL1 = ans.data;
clear ans

% Loading masking file
[neuron_file, neuron_path] =  uigetfile('*.tif'); % Neuron_2.tif
tiffread(neuron_file)
neuron = ans.data;
clear ans

[astrocyte_file, astrocyte_path] =  uigetfile('*.tif'); % Astrocyte_2.tif
tiffread(astrocyte_file)
astrocyte = ans.data;
clear ans


figure(1), 
subplot(1,3,1), imshow(PDL1, [] ), title('PD-L1'),
subplot(1,3,2), imshow(neuron), title('Neuron masking'),
subplot(1,3,3), imshow(astrocyte), title('Astrocyte masking');

neuron(neuron<255) = 0;
neuron(neuron==255) = 1;

astrocyte(astrocyte<255) = 0;
astrocyte(astrocyte==255) = 1;

%% masking
masked_neuron = PDL1.*neuron;
masked_astrocyte = PDL1.*astrocyte;
figure(2),
subplot(1,2,1), imshow(masked_neuron, []), title('PD-L1 expressing in neuron')
subplot(1,2,2), imshow(masked_astrocyte, []), title('PD-L1 expressing in astrocyte')

neuron_number = length( find(neuron==1) );
astrocyte_number = length( find(astrocyte==1) );

average_neuron = sum( sum(masked_neuron) ) / neuron_number;
average_astrocyte = sum( sum(masked_astrocyte) ) / astrocyte_number;

x = categorical({'neuron', 'astrocyte'});
x = reordercats(x,{'neuron', 'astrocyte'});
y = [average_neuron average_astrocyte];
figure(3), bar(x, y)
saveas(gcf,'PDL1.jpg')


%% background correction
set = neuron+astrocyte;
background_index = find(set==0);
PDL1_array = reshape(PDL1, [], 1);
background = PDL1_array(background_index);
background_signal = mean(background);

correction_neuron = average_neuron/background_signal;
correction_astrocyte = average_astrocyte/background_signal;
X = categorical({'neuron', 'astrocyte'});
X = reordercats(x,{'neuron', 'astrocyte'});
Y = [correction_neuron-1 correction_astrocyte-1];
figure(4), bar(X, Y), ylim([0 1])
saveas(gcf,'PDL1_background_normalize.jpg')