% ExVivoHumanLiver_RFA_Decorrelation_ExampleScript.m
% 
% T. Douglas Mast (doug.mast@uc.edu), University of Cincinnati, May 2025
% Example script for reading/processing data from ex-vivo study of
% RFA controlled by echo decorrelation imaging in human liver tissue.
%
% Part of the data set "Real-time control of radiofrequency ablation
%  using three-dimensional ultrasound echo decorrelation imaging in
%  normal and diseased ex vivo human liver" accompanying the article
%  by Ghahramani et al., Physics in Medicine & Biology (2025),
%   https://iopscience.iop.org/article/10.1088/1361-6560/adaacb
%
% Data set is available at the site: https://doi.org/10.7945/2g9q-t294

clear; close all;

% test variables to concatenate over all trials
decorrelation = [];      % cumulative decorrelation within 60^3 voxel
                         %  image volume for each trial, step size 1 mm

ablation = [];           % segmented ablation zones for each trial on
                         %  matching 60^3 grids (ablated voxels = 1,
                         %  unablated = 0)

decorrelation_TCs = [];  % decorrelation at locations of 4 thermocouples
                         %  integrated into the RFA probe, i.e. at
                         %  voxels [22, 54, 45], [34, 49, 45],
                         %  [34, 59, 45], and [30, 54, 39]

temperature_TCs = [];    % temperatures measured by same thermocouples

% loop through all trials' data files to gather/concatenate results
for itrial = 1:109

    % % uncomment to see progress of data loading
    % disp(['Reading data, trial ' num2str(itrial)]);

    load(['output2025/ExVivoHumanLiver_RFA_Decorrelation_RawData_trial' ...
        num2str(itrial) '.mat']);

    % compute log10-scaled cumulative decorrelation per ms for current
    %  trial, as temporal maximum of instantaneous decorrelation maps
    cumulative_decorrelation ...
        = log10(cummax(abs(rawdata.instantaneous_decorrelation),4));

    % last frame of cumulative decorrelation, for use in ROC calculations
    cumdecorr_final = cumulative_decorrelation(:,:,:,end);

    % concatenate cumulative decorrelation and segmented ablation maps
    %  for calculation of an overall ROC curve
    decorrelation = [decorrelation(:); cumdecorr_final(:)];
    ablation = [ablation(:); rawdata.segmented_ablation_zone(:)];

    % concatenate temperatures and cumulative decorrelation measured
    %   at locations of four thermocouples
    if(~isempty(rawdata.thermocouple_temperatures))
        temperature_TCs = [temperature_TCs(:); ...
            rawdata.thermocouple_temperatures(:)];
        decorrelation_TCs = [decorrelation_TCs(:); ...
            cummax(log10(rawdata.decorrelation_at_thermocouples(:)))];
    end 
end

%% example ROC curve across all trials

% compute ROC curve (TPR vs. FPR) and area under curve (AUC)
[FPR,TPR,T,AUC] = perfcurve(ablation,decorrelation,1);

figure(1);
plot(FPR,TPR);
axis square; grid on;
xlabel('False positive rate');
ylabel('True positive rate');
title('ROC curve for local ablation prediction, all trials');
text(0.1,0.05,['AUC = ' num2str(AUC)]);

%% example temperature/decorrelation scatter plot across all trials

figure(2);
plot(decorrelation_TCs,temperature_TCs,'o');
grid on;
xlabel('Cumulative decorrelation, log_{10}(\Delta)/ms');
ylabel('Temperature, ^\circC')
title('Temperature and decorrelation at 4 thermocouples, all trials');

% correlation coefficient and linear regression for temperature vs.
%  decorrelation
[r,p] = corrcoef(decorrelation_TCs,temperature_TCs);
r_decorrelation_temperature = r(2,1);
p_decorrelation_temperature = p(2,1);

tempfit = polyfit(decorrelation_TCs,temperature_TCs,1);
hold on;
plot(get(gca,'XLim'),get(gca,'XLim')*tempfit(1) + tempfit(2),'--r');
hold off;

text(-9.5,135,['r=' num2str(r_decorrelation_temperature)]);
text(-9.5,128,['p=' num2str(p_decorrelation_temperature)]);

%% summary scatter plots of ablation rate and volume vs. ROI decorrelation
load('ExVivoHumanLiver_RFA_Decorrelation_SummaryData.mat');

ablation_volume = vertcat(summary.ablation_volume);
ablation_rate = vertcat(summary.ablation_rate);
ROI_decorrelation = vertcat(summary.ROI_cumulative_decorrelation);

% correlation coefficient and linear regression for ablation volume vs.
%  decorrelation
[r,p] = corrcoef(ROI_decorrelation,ablation_volume);
r_decorrelation_volume = r(2,1);
p_decorrelation_volume = p(2,1);
volumefit = polyfit(ROI_decorrelation,ablation_volume,1);

figure(3);
subplot(211);
plot(ROI_decorrelation,ablation_volume,'o');
grid on;
title('Ablation volume and rate vs. decorrelation, all trials');
xlabel('Cumulative log_{10}-scaled average decorrelation/ms in ROI');
ylabel('Ablation volume (mL)');

hold on;
plot(get(gca,'XLim'),get(gca,'XLim')*volumefit(1) + volumefit(2),'--r');
hold off;
text(-4.25,9,['r=' num2str(r_decorrelation_volume)]);
text(-4.25,8,['p=' num2str(p_decorrelation_volume)]);

% correlation coefficient and linear regression for ablation rate vs.
%  decorrelation
[r,p] = corrcoef(ROI_decorrelation,ablation_rate);
r_decorrelation_rate = r(2,1);
p_decorrelation_rate = p(2,1);
ratefit = polyfit(ROI_decorrelation,ablation_rate,1);

subplot(212);
plot(ROI_decorrelation,ablation_rate,'o');
grid on;
xlabel('Cumulative log_{10}-scaled average decorrelation/ms in ROI');
ylabel('Ablation rate (mL/min)');

hold on;
plot(get(gca,'XLim'),get(gca,'XLim')*ratefit(1) + ratefit(2),'--r');
hold off;
text(-4.25,1.8,['r=' num2str(r_decorrelation_rate)]);
text(-4.25,1.6,['p=' num2str(p_decorrelation_rate)]);