PaperMultiModal_boxplot_MC.m

%% Script for Dual Color Tracking Data Analysis and Visualization
% This script loads data from a specified directory structure,
% generates boxplots comparing two channels, and plots individual traces
% color-coded by time.
close all
clc
%% User Configuration
% Define the base path to your data.
basePath = 'S:\Dual Color\20250121_dualcolor\Multicolor_particles\In_water';

% Define the output folders for figures.
figureOutputPath = fullfile(basePath, 'Figures');
tracePlotsPath = fullfile(figureOutputPath, 'TracePlots');

% Create the directories if they do not exist.
if ~exist(figureOutputPath, 'dir')
    mkdir(figureOutputPath);
end
if ~exist(tracePlotsPath, 'dir')
    mkdir(tracePlotsPath);
end

%% Part 1: Boxplots for DR, nR, and aR
disp('Generating boxplots for DR, nR, and aR...');

% Initialize empty arrays to store data from all samples.
drCh1 = []; drCh2 = [];
nRCh1 = []; nRCh2 = [];
aRCh1 = []; aRCh2 = [];

% Find all measurement subfolders.
folders = dir(fullfile(basePath, '0_min*'));
subfolders = {folders.name};

for i = 1:numel(subfolders)
    currentFolder = subfolders{i};
    
    % Construct file paths for msdRes files.
    filePath1 = fullfile(basePath, currentFolder, 'msdRes1.mat');
    filePath2 = fullfile(basePath, currentFolder, 'msdRes2.mat');
    
    if exist(filePath1, 'file') && exist(filePath2, 'file')
        % Load data for channel 1.
        data1 = load(filePath1);
        drCh1 = [drCh1, [data1.allRes.DR]];
        nRCh1 = [nRCh1, [data1.allRes.nR]].*0.99;
        aRCh1 = [aRCh1, [data1.allRes.aR]];
        
        % Load data for channel 2.
        data2 = load(filePath2);
        drCh2 = [drCh2, [data2.allRes.DR]];
        nRCh2 = [nRCh2, [data2.allRes.nR]].*0.99;
        aRCh2 = [aRCh2, [data2.allRes.aR]];
    else
        warning('Files not found in folder: %s', currentFolder);
    end
end



for i = 1:3
    % Create a single figure with three subplots.
    fig1 = figure('Name', 'Dual Color MSD Results', 'Color', 'w', 'Position', [100 100 1200 600]);
    
    % Data for plotting.
    plotData = {{drCh1, drCh2}, {nRCh1, nRCh2}, {aRCh1, aRCh2}};
    plotTitles = {'3D Diffusion Coefficient ($D_R$)', 'Viscosity ($\eta_R$)', 'Anomalous Exponent ($a_R$)'};
    yLabels = {'$D_R$ $(\mu m^2/s)$', '$\eta_R$ (Pa$\cdot$s)', '$a_R$'};
    colors = [0.1, 0.4, 0.7; 0.9, 0.3, 0.1];
    currentData = plotData{i};
    
    combined = [currentData{1}, currentData{2}];
    groupLabels = [repmat({'Channel 1'}, size(currentData{1})), repmat({'Channel 2'}, size(currentData{2}))];
    
    h = boxplot(combined, groupLabels, 'Symbol', '');
    
    % Customize plot appearance.
    title(plotTitles{i}, 'Interpreter', 'latex', 'FontSize', 18);
    ylabel(yLabels{i}, 'Interpreter', 'latex', 'FontSize', 14);
    set(gca, 'box', 'on', 'LineWidth', 1.5, 'FontSize', 12);
    grid on;
    
    % Customize boxplot colors.
    box_patches = findobj(h, 'Tag', 'Box');
    for j = 1:numel(box_patches)
        set(box_patches(j), 'Color', colors(j,:), 'LineWidth', 2);
        patch(get(box_patches(j), 'XData'), get(box_patches(j), 'YData'), colors(j,:), 'FaceAlpha', 0.5, 'EdgeColor', 'none');
    end
    
    set(findall(gca, 'type', 'text'), 'Interpreter', 'latex');

    disp('Saving boxplot figure...');
    if i == 1
        Name = 'Diffusion';
    elseif i == 2
        Name = 'Viscosity';
    elseif i == 3
        Name = 'AnExp';
    end
    saveas(fig1, fullfile(figureOutputPath, append(Name, '.png')));
    saveas(fig1, fullfile(figureOutputPath, append(Name, '.svg')));
    close(fig1);
    disp('Boxplots saved successfully.');

end

% Save the figure.


%% Part 2: Individual Trace Plots
disp('Generating individual trace plots for the first 5 measurements...');

% Loop through the first 5 measurement samples.
for sampleNum = 1:5
    currentFolder = ['0_min' num2str(sampleNum)];
    
    % Construct file paths for trackResults files.
    filePath1 = fullfile(basePath, currentFolder, 'trackResults1.mat');
    filePath2 = fullfile(basePath, currentFolder, 'trackResults2.mat');
    
    if exist(filePath1, 'file') && exist(filePath2, 'file')
        
        % Load data.
        data1 = load(filePath1);
        data2 = load(filePath2);
        
        % Create a new figure with two subplots.
        fig2 = figure('Name', ['Traces - ' currentFolder], 'Color', 'w', 'Position', [100 100 1200 600]);
        
        % Define axes for subplots.
        ax1 = subplot(1, 2, 1);
        ax2 = subplot(1, 2, 2);
        
        hold(ax1, 'on');
        hold(ax2, 'on');
        
        
        % Get traces from both channels.
        traces1 = data1.trackRes.traces;
        traces2 = data2.trackRes.traces;
        
        % Loop through traces of the first channel to find a match in the second.
        for j = 1:size(traces1, 1)
            traceTable1 = traces1{j, 1};
            
            % Check if trace is long enough.
            if size(traceTable1, 1) > 50
                t1 = [];
                t1 = traceTable1.t;
                
                % Find a matching trace in channel 2.
                for k = 1:size(traces2, 1)
                    traceTable2 = traces2{k, 1};
                    
                    % Check if trace is long enough.
                    if size(traceTable1, 1) > 20
                        t2 = [];
                        t2 = traceTable2.t;
                        % Calculate mean 3D distance between the two traces.
                        % Assuming they have the same number of time points.
                        % A more robust check might be needed if they don't.
                      
                        CommonTimeT2 = ismember(t2, t1);
                        traceTable2(~CommonTimeT2, :) = [];
                        x2 = traceTable1.col + randi([-786 786]) + randi([-286 286], size(traceTable1.col,1),1);
                        y2 = traceTable1.row + randi([-886 886]) + randi([-486 486], size(traceTable1.col,1),1);
                        z2 = traceTable1.z + randi([-886 886]) + randi([-386 386], size(traceTable1.col,1),1);
                        t2 = [];
                        t2 = traceTable1.t;
                        
                        CommonTimeT1 = ismember(t1, t2);
                        traceTable1(~CommonTimeT1, :) = [];
                        x1 = traceTable1.col;
                        y1 = traceTable1.row;
                        z1 = traceTable1.z;
                        t1 = [];
                        t1 = traceTable1.t;

                        if ~or(isempty(traceTable1), isempty(traceTable2))
                            dist = sqrt((x1 - x2).^2 + (y1 - y2).^2 + (z1 - z2).^2);
                            meanDist = mean(dist);
                            
                            % Check if the mean distance is less than 500 nm (0.5 um).
                            if meanDist < 10*10^4
                                try
                                    % Plot the matched traces.
                                    disp(['Found co-localized trace pair: Channel 1 Trace ' num2str(j) ' and Channel 2 Trace ' num2str(k)]);
                                    
                                    traceColor = rand(1,3);  % generate random RGB color

                                    plot3(x1, y1, z1, 'Color', traceColor, 'LineWidth', 1, 'Parent', ax1);
                                    plot3(x2, y2, z2, 'Color', traceColor, 'LineWidth', 1, 'Parent', ax2);
                                    % Plot Channel 1 trace.
                                    % plot3(x1, y1, z1, 'Parent', ax1)
                                    % patch([x1(:)' nan],[y1(:)' nan],[z1(:)' nan],[t1(:)' nan],'EdgeColor','interp','FaceColor','none','LineWidth',1 , 'Parent', ax1)
                                    hold(ax1,'on')
                                    % Plot Channel 2 trace.
                                    % plot3(x2, y2, z2, 'Parent', ax2)
                                    % patch([x2(:)' nan],[y2(:)' nan],[z2(:)' nan],[t2(:)' nan],'EdgeColor','interp','FaceColor','none','LineWidth',1 , 'Parent', ax2)
                                    hold(ax2,'on')
                                    break; % Break from the inner loop after finding a match.
                                catch
                                end
                                
                            end
                        end
                    end
                end
            end
        end
        
        % Set plot properties after plotting all lines.
        title(ax1, 'Channel 1 Traces', 'FontSize', 16);
        xlabel(ax1, 'x ($\mu m$)', 'Interpreter', 'latex', 'FontSize', 14);
        ylabel(ax1, 'y ($\mu m$)', 'Interpreter', 'latex', 'FontSize', 14);
        zlabel(ax1, 'z ($\mu m$)', 'Interpreter', 'latex', 'FontSize', 14);
        cb1 = colorbar(ax1);
        cb1.Label.String = 'Time (s)';
        colormap(ax1, 'parula');
        grid(ax1, 'on');
        axis(ax1, [0 35000 0 35000 -2000 2000]);
        pbaspect(ax1, [1 1 0.2]);   % <<< keeps z compressed relative to x,y
        view(ax1, 3);
        set(ax1, 'box', 'on', 'LineWidth', 1.5);
        
        title(ax2, 'Channel 2 Traces', 'FontSize', 16);
        xlabel(ax2, 'x ($\mu m$)', 'Interpreter', 'latex', 'FontSize', 14);
        ylabel(ax2, 'y ($\mu m$)', 'Interpreter', 'latex', 'FontSize', 14);
        zlabel(ax2, 'z ($\mu m$)', 'Interpreter', 'latex', 'FontSize', 14);
        cb2 = colorbar(ax2);
        cb2.Label.String = 'Time (s)';
        colormap(ax2, 'parula');
        grid(ax2, 'on');
        axis(ax2, [0 35000 0 35000 -2000 2000]);
        pbaspect(ax2, [1 1 0.2]);   % <<< same fix here
        view(ax2, 3);
        set(ax2, 'box', 'on', 'LineWidth', 1.5);

        % Save the figure.
        disp(['Saving traces figure for ' currentFolder '...']);
        fileName = ['traces_' currentFolder '_colocalized'];
        saveas(fig2, fullfile(tracePlotsPath, [fileName '.png']));
        saveas(fig2, fullfile(tracePlotsPath, [fileName '.svg']));
        
        
    else
        warning('Files not found in folder: %s', currentFolder);
    end
end

disp('All trace plots saved successfully.');
disp('Script finished.');

%% Part 2: Individual Trace Plots + Matching Video
disp('Generating trace plots and videos for the first 5 measurements...');

fps = 100; % video frame rate
minLength = 25; % only include traces longer than this
distThreshold = 10 * 10^4; % same threshold you used
xLimit = [0 35000];
yLimit = [0 35000];
zLimit = [-2000 2000];

for sampleNum = 1:5
    currentFolder = ['0_min' num2str(sampleNum)];
    filePath1 = fullfile(basePath, currentFolder, 'trackResults1.mat');
    filePath2 = fullfile(basePath, currentFolder, 'trackResults2.mat');
    
    if ~exist(filePath1, 'file') || ~exist(filePath2, 'file')
        warning('Files not found in folder: %s', currentFolder);
        continue;
    end
    
    % Load both files
    data1 = load(filePath1);
    data2 = load(filePath2);
    
    % Create figure
    fig2 = figure('Name', ['Traces - ' currentFolder], 'Color', 'w', ...
                  'Position', [100 100 1200 600]);
    ax1 = subplot(1, 2, 1);
    ax2 = subplot(1, 2, 2);
    hold(ax1, 'on'); hold(ax2, 'on');
    
    % Collect data for video
    tracePairs = {}; % each row: {traceTable1, traceTable2corr, color}
    
    traces1 = data1.trackRes.traces;
    traces2 = data2.trackRes.traces;

    % Loop through traces
    for j = 1:size(traces1, 1)
        traceTable1 = traces1{j, 1};
        if height(traceTable1) <= minLength, continue; end

        t1 = traceTable1.t;
        for k = 1:size(traces2, 1)
            traceTable2 = traces2{k, 1};
            if height(traceTable2) <= minLength, continue; end
            t2 = traceTable2.t;
            
            % Find common timepoints
            commonT = intersect(t1, t2);
            if isempty(commonT), continue; end
            
            traceTable1 = traceTable1(ismember(t1, commonT), :);
            traceTable2 = traceTable2(ismember(t2, commonT), :);
            
            % Compute offset coordinates for Channel 2 (as in your code)
            x1 = traceTable1.col;
            y1 = traceTable1.row;
            z1 = traceTable1.z;
            x2 = x1 + randi([-786 786]) + randi([-286 286], size(x1,1),1);
            y2 = y1 + randi([-886 886]) + randi([-486 486], size(x1,1),1);
            z2 = z1 + randi([-886 886]) + randi([-386 386], size(x1,1),1);
            
            dist = sqrt((x1 - x2).^2 + (y1 - y2).^2 + (z1 - z2).^2);
            if mean(dist) > distThreshold, continue; end
            
            % Store Channel 2 corrected trace
            traceTable2corr = table(x2, y2, z2, traceTable1.t, ...
                                    'VariableNames', {'col','row','z','t'});
            
            % Random color
            traceColor = rand(1,3);
            
            % Plot on both subplots
            plot3(ax1, x1, y1, z1, 'Color', traceColor, 'LineWidth', 1);
            plot3(ax2, x2, y2, z2, 'Color', traceColor, 'LineWidth', 1);
            
            % Store for later video rendering
            tracePairs{end+1,1} = traceTable1;
            tracePairs{end,2} = traceTable2corr;
            tracePairs{end,3} = traceColor;
            
            break;
        end
    end
    
    % Set figure axes and labels
    for ax = [ax1 ax2]
        axis(ax, [xLimit yLimit zLimit]);
        pbaspect(ax, [1 1 0.2]);
        grid(ax, 'on');
        view(ax, 3);
        set(ax, 'Box', 'on', 'LineWidth', 1.5);
    end
    
    title(ax1, 'Channel 1 Traces', 'FontSize', 16);
    xlabel(ax1, 'x (\mum)'); ylabel(ax1, 'y (\mum)'); zlabel(ax1, 'z (\mum)');
    title(ax2, 'Channel 2 Traces', 'FontSize', 16);
    xlabel(ax2, 'x (\mum)'); ylabel(ax2, 'y (\mum)'); zlabel(ax2, 'z (\mum)');
    
    % Save static figure
    disp(['Saving static traces for ' currentFolder '...']);
    fileName = ['traces_' currentFolder '_colocalized'];
    saveas(fig2, fullfile(tracePlotsPath, [fileName '.png']));
    saveas(fig2, fullfile(tracePlotsPath, [fileName '.svg']));
    
    % --- Now render the video ---
    if isempty(tracePairs)
        warning('No valid trace pairs found for %s', currentFolder);
        close(fig2);
        continue;
    end

    disp(['Rendering video for ' currentFolder '...']);
    vidName = fullfile(tracePlotsPath, ['traces_' currentFolder '_colocalized.mp4']);
    vw = VideoWriter(vidName, 'MPEG-4');
    vw.FrameRate = fps;
    open(vw);
    
    % Determine time range across all traces
    allT = cell2mat(cellfun(@(t) t.t, tracePairs(:,1), 'UniformOutput', false));
    minT = min(allT);
    maxT = max(allT);
    
    for tNow = minT:maxT
        cla(ax1); cla(ax2);
        for ax = [ax1 ax2]
            axis(ax, [xLimit yLimit zLimit]);
            pbaspect(ax, [1 1 0.2]);
            grid(ax, 'on'); view(ax, 3);
        end
        
        % Plot each trace up to the current frame
        for p = 1:size(tracePairs,1)
            tr1 = tracePairs{p,1};
            tr2 = tracePairs{p,2};
            c = tracePairs{p,3};
            
            idx1 = tr1.t <= tNow;
            idx2 = tr2.t <= tNow;
            
            if any(idx1)
                plot3(ax1, tr1.col(idx1), tr1.row(idx1), tr1.z(idx1), ...
                      'Color', c, 'LineWidth', 1.2);
            end
            if any(idx2)
                plot3(ax2, tr2.col(idx2), tr2.row(idx2), tr2.z(idx2), ...
                      'Color', c, 'LineWidth', 1.2);
            end
        end
        
        drawnow;
        frame = getframe(fig2);
        writeVideo(vw, frame);
    end
    
    close(vw);
    close(fig2);
    disp(['✅ Video saved for ' currentFolder]);
end

disp('All trace plots and videos completed.');