% For manuscript "Voxel-wise regression without intensity normalization provides better sensitivity when identifying injury regions after stroke in FDG-PET images"
% Codes generated by Wuxian He
%% Visualize the relationship of voxel intensities between the baseline and MCAO scan.
data1 = data09.BC01;
data2 = data09.BC02;
[img1,~] = imRead3D(fullfile(data1(k).folder,data1(k).name));
[img2,~] = imRead3D(fullfile(data2(k).folder,data2(k).name));

imgLeft1 = img1(41:80,:,:);
imgLeft2 = img2(41:80,:,:);
imgRight1 = img1(1:40,:,:);
imgRight2 = img2(1:40,:,:);
region = pixelLeft{123};
x = imgLeft1(region);
y = imgLeft2(region);
X = [ones(length(x),1) x];
p = X\y;
yCalc = X*p;
R2 = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
figure
scatter(x,y,6,'filled')
grid on
hold on
xlabel('Voxel intensity in the baseline scan','FontSize',15)
ylabel('Voxel intensity in the MCAO scan','FontSize',15)
title('Left striatum','FontSize',15)
t = min(x)-0.1:0.1:max(x)+0.1;
plot(t,t*p(2)+p(1),'LineWidth',1.5)
text(min(x),max(y)-0.1,strcat('R^{2} =',{' '},num2str(R2)),'FontSize',15)
region = pixelRight{123}; % right hemisphere
x = imgRight1(region);
y = imgRight2(region);
X = [ones(length(x),1) x];
p = X\y;
yCalc = X*p;
R2 = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
figure
scatter(x,y,6,'filled')
grid on
hold on
xlabel('Voxel intensity in the baseline scan','FontSize',15)
ylabel('Voxel intensity in the MCAO scan','FontSize',15)
title('Right striatum','FontSize',15)
t = min(x)-0.1:0.1:max(x)+0.1;
plot(t,t*p(2)+p(1),'LineWidth',1.5)
text(min(x),max(y)-0.1,strcat('R^{2} =',{' '},num2str(R2)),'FontSize',15)

%% Stable/unstable regions in each group and data set using box plot.
data1 = data09.normal;
data2 = data09.mcao;
% data1 = random09.A01;
% data2 = random09.A02;
% data1 = data1(1:40);
% data2 = data2(1:40);
% data1 = data09.BC01;
% data2 = data09.BC02;

l2 = zeros(length(data1),length(pixelLeft));
for i = 1:length(data1)
    [img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
    [img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
    imgLeft1 = img1(41:80,:,:);
    imgLeft2 = img2(41:80,:,:);
    for j = 1:length(pixelLeft)
        x = imgLeft1(pixelLeft{j});
        y = imgLeft2(pixelLeft{j});
        X = [ones(length(x),1) x];
        p = X\y;
        yCalc = X*p;
        l2(i,j) = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
    end
end
r2 = zeros(length(data1),length(pixelRight)); % right brain
for i = 1:length(data1)
    [img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
    [img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
    imgRight1 = img1(1:40,:,:);
    imgRight2 = img2(1:40,:,:);
    for j = 1:length(pixelRight)
        x = imgRight1(pixelRight{j});
        y = imgRight2(pixelRight{j});
        X = [ones(length(x),1) x];
        p = X\y;
        yCalc = X*p;
        r2(i,j) = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
    end
end

% stable regions
indLeft = [];
for i = 1:length(pixelRight)
    if median(l2(:,i)) >= 0.9
        if prctile(l2(:,i),25) >= 0.8
            if length(pixelLeft{i}) >= 15
                indLeft = [indLeft i];
            end
        end
    end
end

indRight = [];
for i = 1:length(pixelRight)
    if median(r2(:,i)) >= 0.9
        if prctile(r2(:,i),25) >= 0.8
            if length(pixelRight{i}) >= 15
                indRight = [indRight i];
            end
        end
    end
end

stable09.Aleft = indLeft;
stable09.Aright = indRight;

XX = {stable09.Aleft,stable09.Bleft,stable09.Cleft};
VN=venn(XX{1:3}); % Venn plot to show the consistency
VN=VN.labels('1st subgroup','2nd subgroup','3rd subgroup');
VN=VN.draw();
colorList=[5 185 226;
          239 122 109;
          184 131 212]./255;
for i=1:3
    VN.setPatchN(i,'FaceColor',colorList(i,:),'EdgeColor',colorList(i,:))
end
set(gcf,'Position',[535,245,746,590])
set(gca,'Position',[0.077,0.004,0.836,0.991])
VN.setFont('Color',[0,0,0],'FontSize',20)
VN.setLabel('Color',[0,0,0],'FontSize',20,'FontName','Arial')
% 0.47,0.00169,0.068,0.089

% Unstable regions
data1 = data09.normal;
data2 = data09.mcao;
% data1 = random09.A01;
% data2 = random09.A02;
% data1 = data1(1:40);
% data2 = data2(1:40);
% data1 = random15.A01;
% data2 = random15.A02;
% data1 = data1(1:32);
% data2 = data2(1:32);
% data1 = [random09.A01;random09.B01;random09.C01];
% data2 = [random09.A02;random09.B02;random09.C02];
% data1 = [random15.A01;random15.B01;random15.C01];
% data2 = [random15.A02;random15.B02;random15.C02];

l2 = zeros(length(data1),length(pixelLeft));
for i = 1:length(data1)
    [img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
    [img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
    imgLeft1 = img1(41:80,:,:);
    imgLeft2 = img2(41:80,:,:);
    for j = 1:length(pixelLeft)
        x = imgLeft1(pixelLeft{j});
        y = imgLeft2(pixelLeft{j});
        X = [ones(length(x),1) x];
        p = X\y;
        yCalc = X*p;
        l2(i,j) = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
    end
end
r2 = zeros(length(data1),length(pixelRight)); % right brain
for i = 1:length(data1)
    [img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
    [img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
    imgRight1 = img1(1:40,:,:);
    imgRight2 = img2(1:40,:,:);
    for j = 1:length(pixelRight)
        x = imgRight1(pixelRight{j});
        y = imgRight2(pixelRight{j});
        X = [ones(length(x),1) x];
        p = X\y;
        yCalc = X*p;
        r2(i,j) = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
    end
end

a = zeros(222,1);
for i = 1:222
a(i) = length(pixelRight{i});
end
b = a>=15; % region with at least 15 voxels

ind = median(l2) - median(r2);
ind = ind.*b';
ind = find(ind >= 0.2);
unstable09.A = ind;

 XX = {unstable09.A,unstable09.B,unstable09.C};
VN=venn(XX{1:3}); % Venn plot to show the consistency
VN=VN.labels('1st subgroup','2nd subgroup','3rd subgroup');
VN=VN.draw();
colorList=[5 185 226;
          239 122 109;
          184 131 212]./255;
for i=1:3
    VN.setPatchN(i,'FaceColor',colorList(i,:),'EdgeColor',colorList(i,:))
end
set(gcf,'Position',[535,245,746,590])
set(gca,'Position',[0.077,0.004,0.836,0.991])
VN.setFont('Color',[0,0,0],'FontSize',20)
VN.setLabel('Color',[0,0,0],'FontSize',20,'FontName','Arial')

%% Bar plot showing the R^2 difference between the right and left brain.
data1 = data09.normal;
data2 = data09.mcao;
% data1 = random09.A01;
% data2 = random09.A02;
% data1 = data1(1:40);
% data2 = data2(1:40);
% data1 = [data1;random09.B01];
% data2 = [data2;random09.B02];
% data1 = data1(1:80);
% data2 = data2(1:80);
% data1 = [data1;random09.C01];
% data2 = [data2;random09.C02];
% data1 = data1(1:120);
% data2 = data2(1:120);

l2 = zeros(length(data1),length(pixelLeft));
for i = 1:length(data1)
    [img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
    [img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
    imgLeft1 = img1(41:80,:,:);
    imgLeft2 = img2(41:80,:,:);
    for j = 1:length(pixelLeft)
        x = imgLeft1(pixelLeft{j});
        y = imgLeft2(pixelLeft{j});
        X = [ones(length(x),1) x];
        p = X\y;
        yCalc = X*p;
        l2(i,j) = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
    end
end
r2 = zeros(length(data1),length(pixelRight)); % right brain
for i = 1:length(data1)
    [img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
    [img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
    imgRight1 = img1(1:40,:,:);
    imgRight2 = img2(1:40,:,:);
    for j = 1:length(pixelRight)
        x = imgRight1(pixelRight{j});
        y = imgRight2(pixelRight{j});
        X = [ones(length(x),1) x];
        p = X\y;
        yCalc = X*p;
        r2(i,j) = 1 - sum((y - yCalc).^2)/sum((y - mean(y)).^2);
    end
end

a = zeros(222,1);
for i = 1:222
a(i) = length(pixelRight{i});
end
b = a>15; % region with at least 15 voxels
ind = median(l2) - median(r2);
ind = ind.*b';
ind = ind(ind~=0);
ind_09 = ind(~isnan(ind)); % remove NaN

figure
bar(ind_09)
grid on
hold on
bar(ind_15,'FaceColor',[250 127 111]/255)
set(gcf,'Position',[156,136,1027,719])
legend('1st data set','2nd data set','FontSize',20)
xlabel('Region indices','FontSize',22)
ylabel('R^{2} differences','FontSize',22)
%% Reference region with dice similarity coefficient
brainmask=spm_vol('C:\Users\henry\Documents\MATLAB\Code\19Shenzhen\Template\brainmask.img');
brainmask = spm_read_vols(brainmask);
brainloc=brainmask>=0.8;
clear brainmask

data1 = random09.A01;
data2 = random09.A02;
[y1, y2, ~, ~] = normalization_clus(data1, data2, brainloc, cerebellum);
volume.hypo09A = y1;
data1 = random09.B01;
data2 = random09.B02;
[y1, y2, ~, ~] = normalization_clus(data1, data2, brainloc, cerebellum);
volume.hypo09B = y1;
data1 = random09.C01;
data2 = random09.C02;
[y1, y2, ~, ~] = normalization_clus(data1, data2, brainloc, cerebellum);
volume.hypo09C = y1;

data1 = random15.A01;
data2 = random15.A02;
[y1, y2, ~, ~] = normalization_clus(data1, data2, brainloc, cerebellum);
volume.hypo15A = y1;
data1 = random15.B01;
data2 = random15.B02;
[y1, y2, ~, ~] = normalization_clus(data1, data2, brainloc, cerebellum);
volume.hypo15B = y1;
data1 = random15.C01;
data2 = random15.C02;
[y1, y2, ~, ~] = normalization_clus(data1, data2, brainloc, cerebellum);
volume.hypo15C = y1;

DSC.cere09AB = simi_index(volume.hypo09A,volume.hypo09B);
DSC.cere09AC = simi_index(volume.hypo09A,volume.hypo09C);
DSC.cere09BC = simi_index(volume.hypo09B,volume.hypo09C);

DSCinter.cereA1 = simi_index(volume.hypo09A,volume.hypo15A);
DSCinter.cereA2 = simi_index(volume.hypo09A,volume.hypo15B);
DSCinter.cereA3 = simi_index(volume.hypo09A,volume.hypo15C);
DSCinter.cereB1 = simi_index(volume.hypo09B,volume.hypo15A);
DSCinter.cereB2 = simi_index(volume.hypo09B,volume.hypo15B);
DSCinter.cereB3 = simi_index(volume.hypo09B,volume.hypo15C);
DSCinter.cereC1 = simi_index(volume.hypo09C,volume.hypo15A);
DSCinter.cereC2 = simi_index(volume.hypo09C,volume.hypo15B);
DSCinter.cereC3 = simi_index(volume.hypo09C,volume.hypo15C);

%% Hypo-uptake volume using existing methods on each individual rat.
V = spm_vol('C:\Users\henry\Documents\MATLAB\Code\Template\template_T2.img');
V.fname = 'region.img';
data1 = [random09.A01;random09.B01;random09.C01];
data2 = [random09.A02;random09.B02;random09.C02];
% Grand mean normalization
[~, ~, ~, t] = normalization_GM(data1, data2, brainloc);
[img1,~] = imRead3D(fullfile(data1(1).folder,data1(1).name));
[img2,~] = imRead3D(fullfile(data2(1).folder,data2(1).name));
img1 = img1/mean(img1(brainloc));
img2 = img2/mean(img2(brainloc));
a = img1-img2; 
b = a>=0.15; % t
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; 
spm_write_vol(V,img);
% -------- Automatic screenshot in python -------- %
% Reference region normalization
[~,~,~, t] = normalization_clus(data1, data2, brainloc, cerebellum);
[img1,~] = imRead3D(fullfile(data1(1).folder,data1(1).name));
[img2,~] = imRead3D(fullfile(data2(1).folder,data2(1).name));
img1 = img1/mean(img1(cerebellum));
img2 = img2/mean(img2(cerebellum));
a = img1-img2; 
b = a>=0.18; % t
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; 
spm_write_vol(V,img);
% -------- Automatic screenshot in python -------- %
% Histogram-based normalization
[~,~,~, t] = normalization_HB(data1, data2, template, brainloc);
[img1,~] = imRead3D(fullfile(data1(1).folder,data1(1).name));
[img2,~] = imRead3D(fullfile(data2(1).folder,data2(1).name));
img1 = img1/hist_based_norm(img1,template,brainloc);
img2 = img2/hist_based_norm(img2,template,brainloc);
a = img1-img2; 
b = a>=0.7; % t
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; 
spm_write_vol(V,img);
% -------- Automatic screenshot in python -------- %
% Yakushev's method
[~, ~, ~, t] = normalization_GM(data1, data2, brainloc);
[img1,~] = imRead3D(fullfile(data1(1).folder,data1(1).name));
[img2,~] = imRead3D(fullfile(data2(1).folder,data2(1).name));
img1 = img1/mean(img1(brainloc));
img2 = img2/mean(img2(brainloc));
a = img1-img2;
b = find(a<=-t); % Hyper-uptake clusters
[img1,~] = imRead3D(fullfile(data1(1).folder,data1(1).name));
[img2,~] = imRead3D(fullfile(data2(1).folder,data2(1).name));
img1 = img1/mean(img1(b));
img2 = img2/mean(img2(b));
a = img1-img2; 
b = a>=0.35; % t
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; 
spm_write_vol(V,img);
% -------- Automatic screenshot in python -------- %
% Unstable region
unstable_single(img1,img2,pixelLeft,pixelRight,0.2,V);

%% ROC curve
i = 1;
% GMN
[~,~,~, t] = normalization_GM(data1, data2, brainloc);
t = [-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2 2.5 3 3.5]*t;
roc = [1 1];
k = 1;
[img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
[img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
img1 = img1/mean(img1(brainloc));
img2 = img2/mean(img2(brainloc)); % intensity normalization grand mean
% ------ threshold iteration ------ %
a = img1-img2;
b = a>=t(k);
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; % remove extracranial part
spm_write_vol(V,img);
k = k+1;
% Show the region using Mango.exe and save the screenshot
name = data(i).name; % the name of the TTC segmentation results.
a = imread(name);
a = rgb2gray(a);
b = imread(strcat(name(1:5),'pet.png'));
b = rgb2gray(b);
b = imresize(b,[size(a,1) size(a,2)]);
x = b==76;
y = bwconncomp(x);
y = y.PixelIdxList;
c = b;
c(:,1:174)=0;
z = find(c~=0);
idx = cellfun(@numel,y);
idx = find(idx>=50);
ind_b = [];
for j = 1:length(idx)
ind_b = [ind_b;y{idx(j)}];
end
eval(strcat('ind_a = ind.',name(1:4),';'));
FP = setdiff(ind_b,ind_a);
TP = intersect(ind_b,ind_a);
FN = setdiff(ind_a,ind_b);
TN = setdiff(z,union(ind_b,ind_a));
roc = [roc;length(FP)/length(union(FP,TN)) length(TP)/length(union(TP,FN))];
% ------ after all thresholds all evaluated ------ %
roc = sort(roc,'descend');
eval(strcat('ROC.',name(1:4),'GMN=roc;'))

% RRN
[~,~,~, t] = normalization_clus(data1, data2, brainloc, cerebellum);
t = [-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2 2.5 3 3.5]*t;
roc = [1 1];
k = 1;
[img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
[img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
img1 = img1/mean(img1(cerebellum));
img2 = img2/mean(img2(cerebellum)); % intensity normalization using the cerebellum (indices for 80*108*63)
% ------ threshold iteration ------ %
a = img1-img2;
b = a>=t(k);
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; % remove extracranial part
spm_write_vol(V,img);
k = k+1;
% Show the region using Mango.exe and save the screenshot
name = data(i).name; % the name of the TTC segmentation results.
a = imread(name);
a = rgb2gray(a);
b = imread(strcat(name(1:5),'pet.png'));
b = rgb2gray(b);
b = imresize(b,[size(a,1) size(a,2)]);
x = b==76;
y = bwconncomp(x);
y = y.PixelIdxList;
c = b;
c(:,1:174)=0;
z = find(c~=0);
idx = cellfun(@numel,y);
idx = find(idx>=50);
ind_b = [];
for j = 1:length(idx)
ind_b = [ind_b;y{idx(j)}];
end
eval(strcat('ind_a = ind.',name(1:4),';'));
FP = setdiff(ind_b,ind_a);
TP = intersect(ind_b,ind_a);
FN = setdiff(ind_a,ind_b);
TN = setdiff(z,union(ind_b,ind_a));
roc = [roc;length(FP)/length(union(FP,TN)) length(TP)/length(union(TP,FN))];
% ------ after all thresholds all evaluated ------ %
roc = sort(roc,'descend');
eval(strcat('ROC.',name(1:4),'RRN=roc;'))

% HN
[~,~,~, t] = normalization_HB(data1, data2, template, brainloc);
t = [-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2 2.5 3 3.5]*t;
roc = [1 1];
k = 1;
[img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
[img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
img1 = img1/hist_based_norm(img1,template,brainloc);
img2 = img2/hist_based_norm(img2,template,brainloc);
% ------ threshold iteration ------ %
a = img1-img2;
b = a>=t(k);
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; % remove extracranial part
spm_write_vol(V,img);
k = k+1;
% Show the region using Mango.exe and save the screenshot
name = data(i).name; % the name of the TTC segmentation results.
a = imread(name);
a = rgb2gray(a);
b = imread(strcat(name(1:5),'pet.png'));
b = rgb2gray(b);
b = imresize(b,[size(a,1) size(a,2)]);
x = b==76;
y = bwconncomp(x);
y = y.PixelIdxList;
c = b;
c(:,1:174)=0;
z = find(c~=0);
idx = cellfun(@numel,y);
idx = find(idx>=50);
ind_b = [];
for j = 1:length(idx)
ind_b = [ind_b;y{idx(j)}];
end
eval(strcat('ind_a = ind.',name(1:4),';'));
FP = setdiff(ind_b,ind_a);
TP = intersect(ind_b,ind_a);
FN = setdiff(ind_a,ind_b);
TN = setdiff(z,union(ind_b,ind_a));
roc = [roc;length(FP)/length(union(FP,TN)) length(TP)/length(union(TP,FN))];
% ------ after all thresholds all evaluated ------ %
roc = sort(roc,'descend');
eval(strcat('ROC.',name(1:4),'HN=roc;'))

% CN
[~, ~, ~, t] = normalization_GM(data1, data2, brainloc);
c = t;
t = [-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2 2.5 3 3.5]*(t+0.16);
roc = [1 1];
k = 1;
[img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
[img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
img1 = img1/mean(img1(brainloc));
img2 = img2/mean(img2(brainloc));
a = img1-img2;
b = find(a<=-c); % Hyper-uptake clusters
[img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
[img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
img1 = img1/mean(img1(b));
img2 = img2/mean(img2(b));
% ------ threshold iteration ------ %
a = img1-img2;
b = a>=t(k);
b = double(b);
img = [b(1:40,:,:);zeros(40,108,63)]; % only show the right brain hypo-metabolism
img = img.*brainloc; % remove extracranial part
spm_write_vol(V,img);
k = k+1;
% Show the region using Mango.exe and save the screenshot
name = data(i).name; % the name of the TTC segmentation results.
a = imread(name);
a = rgb2gray(a);
b = imread(strcat(name(1:5),'pet.png'));
b = rgb2gray(b);
b = imresize(b,[size(a,1) size(a,2)]);
x = b==76;
y = bwconncomp(x);
y = y.PixelIdxList;
c = b;
c(:,1:174)=0;
z = find(c~=0);
idx = cellfun(@numel,y);
idx = find(idx>=50);
ind_b = [];
for j = 1:length(idx)
ind_b = [ind_b;y{idx(j)}];
end
eval(strcat('ind_a = ind.',name(1:4),';'));
FP = setdiff(ind_b,ind_a);
TP = intersect(ind_b,ind_a);
FN = setdiff(ind_a,ind_b);
TN = setdiff(z,union(ind_b,ind_a));
roc = [roc;length(FP)/length(union(FP,TN)) length(TP)/length(union(TP,FN))];
% ------ after all thresholds all evaluated ------ %
roc = sort(roc,'descend');
eval(strcat('ROC.',name(1:4),'CN=roc;'))

% Unstable region
[img1,~] = imRead3D(fullfile(data1(i).folder,data1(i).name));
[img2,~] = imRead3D(fullfile(data2(i).folder,data2(i).name));
roc = [1 1];
t = [-1 -0.8 -0.6 -0.4 -0.2 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.4 0.5 0.6 0.8 1.0];
k = 1;
% ------ threshold iteration ------ %
unstable_single(img1,img2,pixelLeft,pixelRight,t(k),V);
k = k+1;
% Show the region using Mango.exe and save the screenshot
name = data(i).name; % the name of the TTC segmentation results.
a = imread(name);
a = rgb2gray(a);
b = imread(strcat(name(1:5),'pet.png'));
b = rgb2gray(b);
b = imresize(b,[size(a,1) size(a,2)]);
x = b==76;
y = bwconncomp(x);
y = y.PixelIdxList;
c = b;
c(:,1:174)=0;
z = find(c~=0);
idx = cellfun(@numel,y);
idx = find(idx>=50);
ind_b = [];
for j = 1:length(idx)
ind_b = [ind_b;y{idx(j)}];
end
eval(strcat('ind_a = ind.',name(1:4),';'));
FP = setdiff(ind_b,ind_a);
TP = intersect(ind_b,ind_a);
FN = setdiff(ind_a,ind_b);
TN = setdiff(z,union(ind_b,ind_a));
roc = [roc;length(FP)/length(union(FP,TN)) length(TP)/length(union(TP,FN))];
% ------ after all thresholds all evaluated ------ %
roc = sort(roc,'descend');
eval(strcat('ROC.',name(1:4),'UR=roc;'))
i = i+1;


i = 1;
b = zeros(5,1); % AUC of an ROC curve
name = data(i).name;
figure
hold on
grid on
eval(strcat('a = ROC.',name(1:4),'GMN;'))
b(1) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
eval(strcat('a = ROC.',name(1:4),'RRN;'))
b(2) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
eval(strcat('a = ROC.',name(1:4),'HN;'))
b(3) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
eval(strcat('a = ROC.',name(1:4),'CN;'))
b(4) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
eval(strcat('a = ROC.',name(1:4),'UR;'))
b(5) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
b = round(b,2);
legend(strcat('GMN: ',num2str(-b(1))),strcat('RRN: ',num2str(-b(2))),strcat('HN: ',num2str(-b(3))),strcat('CN: ',num2str(-b(4))),strcat('UR: ',num2str(-b(5))),'Location','Best','FontSize',12,'FontWeight','bold')
xlabel('False positive rate','FontSize',15,'FontWeight','bold')
ylabel('True positive rate','FontSize',15,'FontWeight','bold')
if i == 1
    title(['The ',num2str(i),'st rat'],'FontSize',15,'FontWeight','bold')
elseif i == 2
    title(['The ',num2str(i),'nd rat'],'FontSize',15,'FontWeight','bold')
elseif i == 3
    title(['The ',num2str(i),'rd rat'],'FontSize',15,'FontWeight','bold')
else
    title(['The ',num2str(i),'th rat'],'FontSize',15,'FontWeight','bold')
end
i = i+1;

% Get the average ROC curve using thresholding averaging.
figure
hold on
grid on
b = zeros(5,1);
a = 0;
for i = 1:11
    name = data(i).name;
    eval(strcat('a = a + ROC.',name(1:4),'GMN;'))
end
a = a/11; % Average ROC for GMN
b(1) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
a = 0;
for i = 1:11
    name = data(i).name;
    eval(strcat('a = a + ROC.',name(1:4),'RRN;'))
end
a = a/11; % Average ROC for RRN
b(2) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
a = 0;
for i = 1:11
    name = data(i).name;
    eval(strcat('a = a + ROC.',name(1:4),'HN;'))
end
a = a/11; % Average ROC for HN
b(3) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
a = 0;
for i = 1:11
    name = data(i).name;
    eval(strcat('a = a + ROC.',name(1:4),'CN;'))
end
a = a/11; % Average ROC for CN
b(4) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
a = 0;
for i = 1:11
    name = data(i).name;
    eval(strcat('a = a + ROC.',name(1:4),'UR;'))
end
a = a/11; % Average ROC for UR
b(5) = trapz(a(:,1),a(:,2));
plot(a(:,1),a(:,2),'LineWidth',2)
b = round(b,2);
legend(strcat('GMN: ',num2str(-b(1))),strcat('RRN: ',num2str(-b(2))),strcat('HN: ',num2str(-b(3))),strcat('CN: ',num2str(-b(4))),strcat('UR: ',num2str(-b(5))),'Location','Best','FontSize',12,'FontWeight','bold')
xlabel('False positive rate','FontSize',15,'FontWeight','bold')
ylabel('True positive rate','FontSize',15,'FontWeight','bold')
title('Averaged ROC','FontSize',15,'FontWeight','bold')

%% Single-subject compares with two-sample t-test
% CN for the 11 rats with TTC
data1 = dir(fullfile('C:\Users\henry\OneDrive\Desktop\Second paper\TTC\PET','swr*1.img'));
data2 = dir(fullfile('C:\Users\henry\OneDrive\Desktop\Second paper\TTC\PET','swr*2.img'));
[~, a, ~, ~] = normalization_GM(data1, data2, brainloc);
b = find(a==1);
[c, ~, ~, ~] = normalization_clus(data1, data2, brainloc, b);
c = -c;
c(41:80,:,:) = 0;
V.fname = 'hypo-uptake.img';
spm_write_vol(V,c);

t = zeros(11,4);
for i = 1:11
    name = data(i).name; % the name of the TTC segmentation results.
    a = imread(name);
    b = imread(strcat(name(1:5),'pet.png'));
    b = rgb2gray(b);
    b = imresize(b,[size(a,1) size(a,2)]);
    x = b==76;
    y = bwconncomp(x);
    y = y.PixelIdxList;
    c = b;
    c(:,1:174)=0;
    z = find(c~=0);
    idx = cellfun(@numel,y);
    idx = find(idx>=50);
    ind_b = [];
    for j = 1:length(idx)
        ind_b = [ind_b;y{idx(j)}];
    end
    eval(strcat('ind_a = ind.',name(1:4),';'));
    FP = setdiff(ind_b,ind_a);
    TP = intersect(ind_b,ind_a);
    FN = setdiff(ind_a,ind_b);
    TN = setdiff(z,union(ind_b,ind_a));
    t(i,1) = length(FP)/length(union(FP,TN));
    t(i,2) = length(TP)/length(union(TP,FN));
    b = imread(strcat(name(1:5),'cn.png'));
    b = rgb2gray(b);
    b = imresize(b,[size(a,1) size(a,2)]);
    x = b==76;
    y = bwconncomp(x);
    y = y.PixelIdxList;
    c = b;
    c(:,1:174)=0;
    z = find(c~=0);
    idx = cellfun(@numel,y);
    idx = find(idx>=50);
    ind_b = [];
    for j = 1:length(idx)
        ind_b = [ind_b;y{idx(j)}];
    end
    FP = setdiff(ind_b,ind_a);
    TP = intersect(ind_b,ind_a);
    FN = setdiff(ind_a,ind_b);
    TN = setdiff(z,union(ind_b,ind_a));
    t(i,3) = length(FP)/length(union(FP,TN));
    t(i,4) = length(TP)/length(union(TP,FN));
end




