Select Git revision
SegmentAndTrack.m
Andrew Cohen authored
SegmentAndTrack.m 18.61 KiB
% SegmentAndTrack.m - This is the main program function for the RAB tools
% application.
% /******************************************************************************
%
% This program, part of rab-tools is copyright (C) 2011-2014 Andrew R.
% Cohen and Mark Winter. All rights reserved.
%
% This software may be referenced as:
%
% Clark, B., M. Winter, A.R. Cohen, and B. Link, Generation of Rab-based transgenic
% lines for in vivo studies of endosome biology in zebrafish. Developmental Dynamics,
% 2011. 240(11): p. 2452-65.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions
% are met:
%
% 1. Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
%
% 2. The origin of this software must not be misrepresented; you must
% not claim that you wrote the original software. If you use this
% software in a product, an acknowledgment in the product
% documentation would be appreciated but is not required.
%
% 3. Altered source versions must be plainly marked as such, and must
% not be misrepresented as being the original software.
%
% 4. The name of the author may not be used to endorse or promote
% products derived from this software without specific prior written
% permission.
%
% THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
% OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
% GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
% WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% Andrew R. Cohen acohen@coe.drexel.edu
% RABTools version 1.0 (release) November 2014
%
% ******************************************************************************/
function SegmentAndTrack(varargin)
[ROOT, moviepath, DatasetName, alpha, extMaxH, minPixelSize, maxPixelSize, bParsed] = parseArgs(nargin, varargin);
if ( ~bParsed )
return;
end
moviename = fullfile(moviepath,[DatasetName '_segmented']);
Hulls = Segment(ROOT, moviename, alpha, extMaxH, minPixelSize, maxPixelSize);
TrackedHulls = Track(Hulls);
moviedir = fullfile(moviepath,[DatasetName '_tracked']);
DrawTracks(ROOT, moviedir, 1, TrackedHulls);
statsname = fullfile(moviepath,[DatasetName '_stats']);
stats = WriteStats(statsname, TrackedHulls);
if ( isdeployed )
close all;
exit;
end
end
function [prefix, postfix] = SplitOnLastUScore(instr)
prefix = [];
postfix = [];
uscr = strfind(instr, '_');
if ( isempty(uscr) )
return;
end
if ( uscr(end) == 1 )
postfix = instr(2:end);
return;
end
if ( uscr(end) == length(instr) )
prefix = instr(1:end-1);
return;
end
prefix = instr(1:(uscr(end)-1));
postfix = instr((uscr(end)+1):end);
end
function [numfiles, startidx, endidx, numdigits, fileprefix] = GetFileNames(ROOT)
flist = dir(fullfile(ROOT, '*_*.tif'));
numfiles = 0;
startidx = [];
endidx = [];
numdigits = [];
fileprefix = [];
for i=1:length(flist)
[chkpath, chkname, chkext] = fileparts(flist(i).name);
[chkprefix, chkpostfix] = SplitOnLastUScore(chkname);
if ( isempty(chkprefix) || isempty(chkpostfix) )
continue;
end
chknum = str2double(chkpostfix);
if ( ~isempty(chknum) )
fileprefix = chkprefix;
numdigits = length(chkpostfix);
break;
end
end
if ( i > length(flist) )
return;
end
chkflist = dir(fullfile(ROOT, [fileprefix '_*.tif']));
numfiles = length(chkflist);
startidx = Inf;
endidx = 0;
for i=1:numfiles
[chkpath, chkname, chkext] = fileparts(chkflist(i).name);
[chkprefix, chkpostfix] = SplitOnLastUScore(chkname);
curidx = str2double(chkpostfix);
startidx = min(curidx, startidx); endidx = max(curidx, endidx);
end
if (numfiles ~= (endidx-startidx+1))
disp(sprintf('Error: Noncontiguous input frames\n'));
end
end
function filename = MakeFName(t, ROOT, fileprefix, startidx, numdigits)
frmtstr = ['%' num2str(numdigits, '%02d') 'd'];
filename=fullfile(ROOT, [fileprefix '_' num2str(startidx+t-1,frmtstr) '.tif']);
end
function Hulls = Segment(ROOT, moviename, alpha, extMaxH, minPixelSize, maxPixelSize)
% DatasetName='rab 5c rab 5'
% ROOT='rab 5\';
GENAVI= 0;
if GENAVI
aviobj = avifile(moviename, 'compression', 'none','fps',5);
% else
% mkdir(moviename);
end
Hulls=[];
%frames = dir(fullfile(ROOT, '*_*.tif'));
%frameCount = length(frames);
[frameCount, startidx, endidx, numdigits, fileprefix] = GetFileNames(ROOT);
if ( frameCount < 1 )
disp(sprintf('Error: No video frames found in directory %s\n', ROOT));
if ( isdeployed )
exit;
else
return;
end
end
bWaitbar = 1;
if ( bWaitbar )
hWait = waitbar(0, ['Segmenting Frames(0 of ' num2str(frameCount) ')...']);
end
for t=1:frameCount
if ( bWaitbar )
waitbar(t/frameCount, hWait, ['Segmenting Frames(' num2str(t) ' of ' num2str(frameCount) ')...']);
end
%fname=fullfile(ROOT, ['a00000001_' num2str(t-1,'%02d') '.tif']);
fname = MakeFName(t, ROOT, fileprefix, startidx, numdigits);
im=imread(fname);
im=im2double(im);
im=mat2gray(im);
level=alpha * graythresh(im);
bw=im2bw(im,level);
hold off;imagesc(im);colormap gray;hold on
bwimMax=imextendedmax(im,extMaxH,8);
bwimMax = bwimMax & bw;
% bwimMax=bwmorph(bwimMax,'fill');
%bwimMax=bwmorph(bwimMax,'open');
bwimMax = imopen(bwimMax, strel('disk',1));
[r c]=find(bwimMax & im);
%plot(c,r,'or');
[L num]=bwlabel(bwimMax,4);
for n=1:num
pix=find(L==n); %if length(pix)<minPixelSize,continue,end
%if length(pix)>maxPixelSize,continue,end
nh=[];
bwOrganelle=0*bwimMax;
bwOrganelle(pix)=1;
B=bwboundaries(bwOrganelle);
if ( length(pix)<minPixelSize )
plot(B{1}(:,2),B{1}(:,1),'-b');
continue;
elseif ( length(pix)>maxPixelSize )
plot(B{1}(:,2),B{1}(:,1),'-g');
continue;
else
plot(B{1}(:,2),B{1}(:,1),'-r');
end
[ r c]=find(L==n);
nh.t=t;
nh.xyCenter=[mean(B{1}(:,2)),mean(B{1}(:,1))];
nh.area=bwarea(bwOrganelle);
nh.xyPoints=[c,r];
nh.xyPerim=[B{1}(:,2),B{1}(:,1)];
Hulls=[Hulls;nh];
end
drawnow
if GENAVI
fr=getframe();
aviobj = addframe(aviobj,fr);
% else
% outfname = MakeFName(t, moviename, fileprefix, startidx, numdigits);
% print(gcf,'-dtiff', '-r300', outfname);
end
end
if ( bWaitbar )
close(hWait);
end
if GENAVI
aviobj = close(aviobj);
end
%save ([DatasetName '.mat'], 'Hulls')
end
function TrackedHulls = Track(Hulls)
%load ([DatasetName '.mat'], 'Hulls')
DXY_MAX=25;
%path(path,'assign');
TrackedHulls = Hulls;
tmax=max([TrackedHulls.t]);
for i=1:length(TrackedHulls)
TrackedHulls(i).ID=i;
end
for t=1:tmax-1
th=find([TrackedHulls.t]==t);
th1=find([TrackedHulls.t]==t+1);
if isempty(th) | isempty(th1)
continue
end
d=[];
for i=1:length(th)
for j=1:length(th1)
d(i,j)=Inf;
xy=TrackedHulls(th(i)).xyCenter;
xy1=TrackedHulls(th1(j)).xyCenter;
dxy=norm([xy-xy1]);
if dxy>DXY_MAX
continue
end
area=TrackedHulls(th(i)).area;
area1=TrackedHulls(th1(j)).area;
darea=(max(area,area1)-min(area,area1))/max(area,area1);
lifetime=length(find([TrackedHulls.ID]==TrackedHulls(th(i)).ID));
d(i,j)=25*dxy/DXY_MAX + 10*darea;
d(i,j)=d(i,j)*(t/(lifetime))^2;
end
end
[assign cost]=assignmentoptimal(d);
for i=1:length(assign)
if 0==assign(i)
continue
end
TrackedHulls(th1(assign(i))).ID=TrackedHulls(th(i)).ID;
end
end
%save([DatasetName '_tracked.mat'], 'Hulls')
end
function DrawTracks(ROOT, moviename, bOutTif, Hulls)
%load ([DatasetName '_tracked.mat'])
GENAVI= ~bOutTif;
if GENAVI
aviobj = avifile(moviename,'compression','none');
elseif ( ~exist(moviename, 'dir') )
mkdir(moviename);
end
[frameCount, startidx, endidx, numdigits, fileprefix] = GetFileNames(ROOT);
%renumber
id = unique([Hulls.ID]);
for i = 1 : length(Hulls)
idi = find(id == Hulls(i).ID);
Hulls(i).ID = idi;
end
tmax=max([Hulls.t]);
bWaitbar = 1;
if ( bWaitbar )
hWait = waitbar(0, ['Tracking Frames(0 of ' num2str(tmax) ')...']);
end
for t=1:tmax
if ( bWaitbar )
waitbar(t/tmax, hWait, ['Tracking Frames(' num2str(t) ' of ' num2str(tmax) ')...']);
end
%fname=fullfile(ROOT,['a00000001_' num2str(t-1,'%02d') '.tif']);
fname = MakeFName(t, ROOT, fileprefix, startidx, numdigits);
im=imread(fname); hold off
imagesc(im)
colormap(gray)
axis image
axis off
hold on
th=find([Hulls.t]==t);
text(1,1,['t=' num2str(t)], 'VerticalAlignment','top',...
'fontsize',6,'BackgroundColor',[.7 .9 .7]);
for i=1:length(th)
hi=Hulls(th(i));
oti=find([Hulls.ID]==hi.ID);
if length(oti)<5
plot(hi.xyPerim(:,1),hi.xyPerim(:,2),'-g' )
continue
end
ccc='r';
[xmax imax]=max(hi.xyPerim(:,1));
ymax=hi.xyPerim(imax,2);
plot(hi.xyPerim(:,1),hi.xyPerim(:,2),'-r' )
text(xmax,ymax,[num2str(hi.ID)], 'color',ccc,'fontweight','bold','fontsize',10,...
'VerticalAlignment','top');
end
drawnow
if GENAVI
F = getframe();
aviobj = addframe(aviobj,F);
else
%outfname = fullfile(moviename, ['a00000001_' num2str(t-1,'%02d') '_tracked.tif']);
outfname = MakeFName(t, moviename, fileprefix, startidx, numdigits);
print(gcf,'-dtiff', '-r300', outfname);
end
end
if ( bWaitbar )
close(hWait);
end
if GENAVI
aviobj = close(aviobj);
end
end
function stats = WriteStats(statsname, Hulls)
%load ([DatasetName '_tracked.mat'])
%renumber
id = unique([Hulls.ID]);
for i = 1 : length(Hulls)
idi = find(id == Hulls(i).ID);
Hulls(i).ID = idi;
end
stats=[];
for i = 1 : max([Hulls.ID])
oi=find([Hulls.ID]==i);
if length(oi)<5,continue,end
ss=[];
ss(1)=i;
ss(2)= Hulls(oi(end)).t-Hulls(oi(1)).t+1;
ss(3)= Hulls(oi(1)).t;
ss(4)= Hulls(oi(end)).t; for j=1:length(oi)
idx=(Hulls(oi(j)).t-1)*3+1;
ss(4+idx)= Hulls(oi(j)).area;
ss(5+idx:6+idx)= round(Hulls(oi(j)).xyCenter);
end
ss(60*3+6)=0;
stats=[stats;ss];
end
stats=sortrows(stats,-2);
wrtstat = cell(size(stats) + [3,1]);
wrtstat(4:end,2:end) = num2cell(stats);
wrtstat{2,6} = 'time->';
wrtstat(3,2:5) = {'tracking ID', 'lifetime(frames)', 'first frame', 'last frame'};
for t=1:max([Hulls.t])
wrtstat{2,7+3*(t-1)} = num2str(t);
wrtstat(3,6+3*(t-1):6+3*(t-1)+2) = {'area', 'xcenter', 'ycenter'};
end
bWritten = xlswrite(statsname, wrtstat);
if ( ~bWritten )
disp(sprintf('Warning: Unable to write stats as excel file, writing CSV instead'));
csvwrite(statsname, stats);
end
end
function SegmentFirst(ROOT, moviename, alpha, extMaxH, minPixelSize, maxPixelSize)
[frameCount, startidx, endidx, numdigits, fileprefix] = GetFileNames(ROOT);
if ( frameCount < 1 )
disp(sprintf('Error: No video frames found in directory %s\n', ROOT));
if ( isdeployed )
exit;
else
return;
end
end
fname = MakeFName(1, ROOT, fileprefix, startidx, numdigits);
im=imread(fname);
im=im2double(im);
im=mat2gray(im);
level=alpha * graythresh(im);
bw=im2bw(im,level);
figure;hold off;imagesc(im);colormap gray;hold on
bwimMax=imextendedmax(im,extMaxH,8);
bwimMax = bwimMax & bw;
bwimMax = imopen(bwimMax, strel('disk',1));
[r c]=find(bwimMax & im);
%plot(c,r,'or');
[L num]=bwlabel(bwimMax,4);
for n=1:num
pix=find(L==n);
nh=[];
bwOrganelle=0*bwimMax;
bwOrganelle(pix)=1;
B=bwboundaries(bwOrganelle);
if ( length(pix)<minPixelSize )
plot(B{1}(:,2),B{1}(:,1),'-b', 'LineWidth',2);
elseif ( length(pix)>maxPixelSize )
plot(B{1}(:,2),B{1}(:,1),'-g', 'LineWidth',2); else
plot(B{1}(:,2),B{1}(:,1),'-r', 'LineWidth',2);
end
[ r c]=find(L==n);
end
figure;imagesc(im);colormap gray;
drawnow
end
function [ROOT, moviepath, DatasetName, alpha, extMaxH, minPixelSize, maxPixelSize, bParsed] = parseArgs(callNargin, args)
ROOT = [];
moviepath = [];
DatasetName = [];
alpha = [];
extMaxH = [];
minPixelSize = [];
maxPixelSize = [];
bParsed = 0;
if ( callNargin < 1 )
ROOT = uigetdir();
if ( isempty(ROOT) )
printUsage();
if ( isdeployed )
exit;
else
return;
end
end
else
seps = strfind(args{1}, filesep);
if ( ~isempty(seps) && (seps(end) == length(args{1})) )
ROOT = args{1}(1:end-1);
else
ROOT = args{1};
end
end
if ( ~exist(ROOT,'dir') && exist(ROOT,'file') )
seps = strfind(ROOT, filesep);
if ( isempty(seps) )
disp(sprintf('Error: Cannot find input directory specified: %s\n', ROOT));
return;
end
ROOT = ROOT(1:seps(end)-1);
if ( ~exist(ROOT,'dir') )
disp(sprintf('Error: Cannot find input directory specified: %s\n', ROOT));
return;
end
end
moviepath = [];
% if ( callNargin >= 2 )
% [moviepath,DatasetName,dumpext] = fileparts(args{2});
% else
seps = strfind(ROOT, filesep);
if ( isempty(seps) )
if ( ~isempty(ROOT) )
DatasetName = ROOT;
else
DatasetName = 'default';
end
else
DatasetName = ROOT(seps(end)+1:end);
moviepath = ROOT(1:seps(end)-1);
end
% end
defAlpha = 1.0;
defExtMaxH = 0.35;
defMinPixelSize = 7;
defMaxPixelSize = 100;
if ( callNargin <= 2 )
paramDone = 0;
while ( ~paramDone )
answers = inputdlg({'Alpha:','Extended Max:','Min Pixel Size:','Max Pixel Size:'},'Segmentation Parameters',1,{num2str(defAlpha),num2str(defExtMaxH),num2str(defMinPixelSize),num2str(defMaxPixelSize)});
if ( isempty(answers) )
if ( isdeployed )
exit;
else
return;
end
end
alpha = str2num(answers{1});
extMaxH = str2num(answers{2});
minPixelSize = str2num(answers{3});
maxPixelSize = str2num(answers{4});
moviename = fullfile(moviepath,[DatasetName '_segmented']);
SegmentFirst(ROOT, moviename, alpha, extMaxH, minPixelSize, maxPixelSize);
button = nmquestdlg('Are these parameters correct?','Parameter Check');
close all;
if ( strcmpi(button,'yes') )
bParsed = 1;
return;
elseif ( strcmpi(button,'cancel') )
if ( isdeployed )
exit;
else
return;
end
end
end
end
if ( callNargin >= 2 )
alpha = str2num(args{2});
if ( isempty(alpha) )
printUsage();
if ( isdeployed )
exit;
else
return;
end
end
else
alpha = defAlpha;
end
if ( callNargin >= 3 )
extMaxH = str2num(args{3});
if ( isempty(extMaxH) )
printUsage();
if ( isdeployed )
exit;
else
return;
end
end
else
extMaxH = defExtMaxH;
end
if ( callNargin >= 4 )
minPixelSize = str2num(args{4}); if ( isempty(minPixelSize) )
printUsage();
if ( isdeployed )
exit;
else
return;
end
end
else
minPixelSize = defMinPixelSize;
end
if ( callNargin >= 5 )
maxPixelSize = str2num(args{5});
if ( isempty(maxPixelSize) )
printUsage();
if ( isdeployed )
exit;
else
return;
end
end
else
maxPixelSize = defMaxPixelSize;
end
bParsed = 1;
end
function printUsage()
disp(sprintf('\n'));
disp(sprintf('\n'));
disp(sprintf('Usage: SegmentAndTrack input_dir [alpha_parameter extended_max_parameter minimum_pixel_size maximum_pixel_size]\n'));
disp(sprintf(' If unspecified: alpha_parameter = 1.0\n'));
disp(sprintf(' extended_max_parameter = 0.35\n'));
disp(sprintf(' minimum_pixel_size = 7\n'));
disp(sprintf(' maximum_pixel_size = 100\n'));
disp(sprintf('\n'));
end