% Plot one 3D PET concentration frame.
clear; close all; clc;

selected_frame = 7; % choose 1--19

load('PET_3D_concentration.mat', 'C3D');
V = squeeze(C3D(selected_frame,:,:,:)); % [x,z,y]
V = flip(permute(V,[3 2 1]),1);         % manuscript orientation

voxel_xy = 2.329; % mm
voxel_z  = 2.388; % mm
x = ((1:size(V,1))-0.5)*voxel_xy;
z = ((1:size(V,2))-0.5)*voxel_z;
y = ((1:size(V,3))-0.5)*voxel_xy;

blue  = [0.0784 0.1686 0.5490];
white = [0.8000 0.8000 0.8000];
red   = [0.5020 0.0000 0.0000];
t = (0:21)'/21;
blue_to_white = round((1-t).*blue + t.*white,4);
white_to_red = round((1-t).*white + t.*red,4);
cmap = [repmat(blue,8,1); blue_to_white(2:end,:); ...
        repmat(white,6,1); white_to_red(2:end,:); ...
        repmat(red,6,1)];

figure('Color','w');
h = PATCH_3Darray(V,x,z,y,cmap,[0 0.15],'col');
if iscell(h)
    for k = 1:numel(h)
        set(h{k},'EdgeColor','none');
    end
else
    set(h,'EdgeColor','none');
end

colormap(cmap);
caxis([0 0.15]);
axis equal tight off;
view(305,30);
cb = colorbar;
cb.Label.String = 'C/C_0 [-]';
title(sprintf('PET frame %d', selected_frame));