Potential Bug? Freezing while making plots iteratively

8 views (last 30 days)
Hi community,
I am at a loss for why my matlab code is freezing up, I am not using that much of matlab's capabilities...I am attempting to read through a directory, and create a tiledlayout 2x2 plot for each folder in the directory. In this, I create an exponential trendline to fit some of my data, and once I try plotting all of the stuff in the tiled layout, it can cause the system to freeze up right before saving the picture, preventing the for loop from reaching the next file. It shouldn't be a code failure, as this code has worked in pieces before, but causing the entire app become entirely unresponsive for hours requiring a force quit on a 2022 computer makes me relatively surprised.
For the purposes of ease of viewing, I will attach the code in a .m, and paste the file's full code here, and I will also attach a link to a compressed folder of a few files that will allow one to run this code.
To be clear the freeze and crash happens sporadically, sometimes on the first file, sometimes on the 5 file, which makes me believe there is a bug of some sorts! I see it more frequently occur when I have a debug run-stop at the last end line after the saveas function, so that I can view the figure for any issues before moving to the next folder. I appreciate the help!
clear all
close all
clc
%% Change filepath and timestep frequency
filepath = '/Users/nickscott/Library/CloudStorage/Box-Box/Parekh Lab Public/Nick/CellProfiler/ProperFPS_NikonExp';
cd '/Users/nickscott/Library/CloudStorage/Box-Box/Parekh Lab Public/Nick/CellProfiler/ProperFPS_NikonExp/'
topLevelFolder = pwd;
files = dir(topLevelFolder);
dirFlags = [files.isdir];
% Extract only those that are directories.
subFolders = files(dirFlags); % A structure with extra info.
% Get only the folder names into a cell array.
folderList = {subFolders(3:end).name}; % Start at 3 to skip . and ..
for i= 1: length(folderList)
thisFolder = folderList{i};
fprintf('Processing folder %s\n', thisFolder);
% Get txt files.
filePattern = sprintf('%s/*.txt', thisFolder);
baseFileNames = dir(filePattern);
bFN = struct2cell(baseFileNames)
dt{i}= bFN{1}
end
%% Constants
FigList = {'Median Trendline of Perinuclear Intensity', 'Median Trendline of Nuclei Intensity'...
'Best Fit Using Exponential Function', 'Median Trendline of Nuclei/Perinuclear Intensity'};
PeriRing = '/ShNucDiffPipe_Perinuclear_Ring.csv'; %Perinuclear Ring file name
Nuc = '/ShNucDiffPipe_ShrunkenNuclei.csv'; %Stained Nuclei file name
%% Giant 'for loop' for every folder in the batch
for p= 1:length(dt)
%% Naming Extraction
newstr= extractBetween([filepath '/' folderList{p} '/'], "NikonExp/","/"); %Pull Metadata folder name
FigTitle = regexprep(newstr, '_', ' '); %Change _ to spaces
%% Read Table, Extract Variables of Interest
%Automatically extract the name
TPR = readtable([filepath '/' folderList{p} PeriRing], 'Delimiter',','); %perinuclear table
TN = readtable([filepath '/' folderList{p} Nuc], 'Delimiter',','); %nuclei table
TT = readtable([filepath '/' folderList{p} '/' char(dt{1,p})], 'Delimiter', 'tab');
% Modify Time frame to start out at 0 instead of being image set index from
% ... Cellprofiler, may need to change if pipeline can be adjusted
TPR{:,1}=TPR{:,1}-TPR{1,1};
TN{:,1}=TN{:,1}-TN{1,1};
%Extract Unique nuclei objects (must be in column 2)
[~,~,j]=unique(TPR(:,2));
ArrayPR = accumarray(j,(1:length(j)).', [], @(k) {TPR(k,:)});
[~,~,l]=unique(TN(:,2));
ArrayN = accumarray(l,(1:length(l)).', [], @(k) {TN(k,:)});
table_heightsPR = cellfun(@height,ArrayPR);
table_heightsN = cellfun(@height,ArrayN);
height_thresholdPR = 0.75*max(table_heightsPR);
height_thresholdN = 0.75*max(table_heightsN);
LineNumb = 0;
if contains(TT.Properties.VariableNames{2}, 'Time_s_ms_')
for i= 1:length(table_heightsN)
if height(ArrayN{i}) >= height_thresholdN
ArrayN{i,1}.Timestamp = zeros(table_heightsN(i),1);
[row,col] = size(ArrayN{i,1});
LineNumb = LineNumb + 1;
for j=1:row
FrameNum = ArrayN{i,1}.ImageNumber(j);
FrameNum = FrameNum+1;
ind1 = find(TT.Index == FrameNum);
ArrayN{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_ms_(ind1)),'%f.%f').*[1;0.0001]);
end
else
end
end
for i= 1:length(table_heightsPR)
if height(ArrayPR{i}) >= height_thresholdPR
ArrayPR{i,1}.Timestamp = zeros(table_heightsPR(i),1);
[row,col] = size(ArrayPR{i,1});
for j=1:row
FrameNum = ArrayPR{i,1}.ImageNumber(j);
FrameNum = FrameNum+1;
ind2 = find(TT.Index == FrameNum);
ArrayPR{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_ms_(ind2)),'%f.%f').*[1;0.0001]); %sum(sscanf(char(TT.Time_m_s_ms_(ind1)),'%f:%f').*[60;1]);
end
else
end
end
elseif contains(TT.Properties.VariableNames{2}, 'Time_m_s_ms_')
for i= 1:length(table_heightsN)
if height(ArrayN{i}) >= height_thresholdN
ArrayN{i,1}.Timestamp = zeros(table_heightsN(i),1);
[row,col] = size(ArrayN{i,1});
LineNumb = LineNumb + 1;
for j=1:row
FrameNum = ArrayN{i,1}.ImageNumber(j);
FrameNum = FrameNum+1;
ind1 = find(TT.Index == FrameNum);
ArrayN{i,1}.Timestamp(j) = sum(sscanf(char(TT.Time_m_s_ms_(ind1)),'%f:%f').*[60;1]);
end
else
end
end
for i= 1:length(table_heightsPR)
if height(ArrayPR{i}) >= height_thresholdPR
ArrayPR{i,1}.Timestamp = zeros(table_heightsPR(i),1);
[row,col] = size(ArrayPR{i,1});
for j=1:row
FrameNum = ArrayPR{i,1}.ImageNumber(j);
FrameNum = FrameNum+1;
ind2 = find(TT.Index == FrameNum);
ArrayPR{i,1}.Timestamp(j) = sum(sscanf(char(TT.Time_m_s_ms_(ind2)),'%f:%f').*[60;1]);
end
else
end
end
elseif contains(TT.Properties.VariableNames{2}, 'Time_s_')
for i= 1:length(table_heightsN)
if height(ArrayN{i}) >= height_thresholdN
ArrayN{i,1}.Timestamp = zeros(table_heightsN(i),1);
[row,col] = size(ArrayN{i,1});
LineNumb = LineNumb + 1;
for j=1:row
FrameNum = ArrayN{i,1}.ImageNumber(j);
FrameNum = FrameNum+1;
ind1 = find(TT.Index == FrameNum);
ArrayN{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_(ind1)),'%f').*[1]);
end
else
end
end
for i= 1:length(table_heightsPR)
if height(ArrayPR{i}) >= height_thresholdPR
ArrayPR{i,1}.Timestamp = zeros(table_heightsPR(i),1);
[row,col] = size(ArrayPR{i,1});
for j=1:row
FrameNum = ArrayPR{i,1}.ImageNumber(j);
FrameNum = FrameNum+1;
ind2 = find(TT.Index == FrameNum);
ArrayPR{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_(ind2)),'%f').*[1]);
end
else
end
end
end
%% Plot to the tiled subplots *if* the nuclei doesn't float in and out of frame significantly
Tile= tiledlayout(2,2);
set(gcf, 'WindowState', 'maximized');
nexttile
hold on
for i= 1:length(ArrayPR) %i=Nucleus number
if height(ArrayPR{i}) >= height_thresholdPR
plot(ArrayPR{i,1}.Timestamp,ArrayPR{i,1}.Intensity_MedianIntensity_Alexa488,'color',[0,0,0]+0.75)
else
end
end
% median intensity for each image number calculation:
% first, concatenate all the tables into one big table:
% t = vertcat(ArrayPR{:});
% or concatenate just the tables whose height is >= height_threshold:
t = vertcat(ArrayPR{table_heightsPR >= height_thresholdPR});
% get all the image numbers:
all_image_time = t.Timestamp;
% get all the intensities:
all_intensity = t.Intensity_MedianIntensity_Alexa488;
mask = isnan(all_intensity);
all_intensity(mask) = 0;
% get the unique image numbers:
u_image_time = unique(all_image_time);
% calculate the median intensity for each image number:
% number of unique image numbers:
n_images = numel(u_image_time);
% pre-allocating the median intensity array:
median_intensity_pr = zeros(1,n_images);
% for each unique image number:
for f = 1:n_images
% get a logical index saying whether each element of all_image_number
% is equal to this image number (idx is true where
% all_image_number == u_image_number(i) and false elsewhere):
idx_pr = all_image_time == u_image_time(f);
% use that logical index to calculate the median intensity for this
% image number:
median_intensity_pr(f) = median(all_intensity(idx_pr));
end
% plot
plot(u_image_time,median_intensity_pr,'k','LineWidth',2); %Black median line
title(FigList{1})
hold off
nexttile
hold on
for i= 1:length(ArrayN) %i=Nucleus number
if height(ArrayN{i}) >= height_thresholdN
plot(ArrayN{i,1}.Timestamp,ArrayN{i,1}.Intensity_MedianIntensity_Alexa488,'color',[0,0,0]+0.75)
else
end
end
% median intensity for each image number calculation:
% first, concatenate all the tables into one big table:
% t = vertcat(ArrayPR{:});
% or concatenate just the tables whose height is >= height_threshold:
N = vertcat(ArrayN{table_heightsN >= height_thresholdN});
% get all the image numbers:
all_image_time_N = N.Timestamp;
% get all the intensities:
all_intensity_N = N.Intensity_MedianIntensity_Alexa488;
% get the unique image numbers:
u_image_time_N = unique(all_image_time_N);
% calculate the median intensity for each image number:
% number of unique image numbers:
n_images_N = numel(u_image_time_N);
% pre-allocating the median intensity array:
median_intensity_N = zeros(1,n_images);
% for each unique image number:
for f = 1:n_images
% get a logical index saying whether each element of all_image_number
% is equal to this image number (idx is true where
% all_image_number == u_image_number(i) and false elsewhere):
idx = all_image_time_N == u_image_time_N(f);
% use that logical index to calculate the median intensity for this
% image number:
median_intensity_N(f) = median(all_intensity_N(idx));
end
% plot
plot(u_image_time_N,median_intensity_N,'k','LineWidth',2); %Black median line
title(FigList{2})
hold off
%plot(t(indexOfInterest),signal(indexOfInterest)) % plot on new axes
%if height(ArrayN{i}) >= height_thresholdN
indexOfInterest = (ArrayN{1,1}.Timestamp <= 12.5) & (ArrayN{1,1}.Timestamp >= 2.4);
IoI1 = find(indexOfInterest,1,'first');
nexttile
hold on
xdata = u_image_time(IoI1:127);
ydata = transpose(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127));
fun = @(cs,xdata) cs(1)-cs(2)*exp(-cs(3)*xdata);
x0 = [1.5, 2.5, 0.3];
lb = [-Inf, 0,0];
ub = [Inf, Inf, Inf];
x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)
times = linspace(xdata(1),xdata(end));
plot(xdata,ydata,'ko',times,fun(x,times),'b-')
txt = {'a-b*exp(-c*t)', ['a= ' num2str(x(1))], ['b= ' num2str(x(2))], ['c= ' num2str(x(3))]};
text(max(u_image_time(IoI1:127))*2/3, (max(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127))-min(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127)))/5+min(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127)), txt)
title(FigList{3})
hold off
nexttile
hold on
for i= 1:length(ArrayN) %i=Nucleus number
if height(ArrayN{i}) >= 0.75*max(table_heightsN)
plot(ArrayN{i,1}.Timestamp,ArrayN{i,1}.Intensity_MedianIntensity_Alexa488./ArrayPR{i,1}.Intensity_MedianIntensity_Alexa488,'color',[0,0,0]+0.75)
else
end
end
plot(u_image_time,(median_intensity_N./median_intensity_pr),'k','LineWidth',2); %Black median line
title(FigList{4})
hold off
a = string(LineNumb);
tite = strcat(FigTitle, ' N=', a);
title(Tile,tite, 'fontsize', 16)
xlabel(Tile,'Time Scale (sec)', 'fontsize', 16)
ylabel(Tile,'Median Fluorescence Intensity', 'fontsize', 16)
Tile.TileSpacing = 'compact';
Tile.Padding = 'normal';
set(findall(Tile, '-property', 'fontsize'), 'fontsize', 16)
picName = char(strcat(filepath, '/', folderList{p}, '/NormalizedTrendFit', newstr, '_', datestr(now, 'yyyymmddHHMM'), '.png'))
saveas(gcf, picName);
end
  5 Comments
dpb
dpb on 19 Sep 2022
Edited: dpb on 19 Sep 2022
See what
rendererinfo
returns to see what it says.
Then
opengl software|hardware
swaps back/forth. See the doc for all the gory details, of course...

Sign in to comment.

Answers (1)

Matt J
Matt J on 19 Sep 2022
Edited: Matt J on 19 Sep 2022
Does the freeze always happen right before the saveas? If so, it might be worth issuing a drawnow() command right before the save step.
picName = char(strcat(filepath, '/', folderList{p}, '/NormalizedTrendFit', newstr, '_', datestr(now, 'yyyymmddHHMM'), '.png'))
drawnow
saveas(gcf, picName);
My suspicion is that the figure cannot be saved because the drawing of the figure is not fully complete, resulting in some sort of tug of war between saveas and the renderer.
  1 Comment
Nicholas Scott
Nicholas Scott on 19 Sep 2022
it always occurs right before the saveas and the figure is not at full size. I'll add that and see if it alleviates the issue. Thanks Matt J!

Sign in to comment.

Categories

Find more on Graphics Object Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!