Code covered by the BSD License

# Blood cells tracking and measurement by using spatiotemporal images analysis

### Chen yuan (view profile)

01 Feb 2011 (Updated )

Three demos for blood cells tracking.

demo3.m
```%-------------------------------------------------------------------------
% This is a simulation demo to show multiple trackpaths for tracking the
% blood cell motion. In this demo, the diameter of the blood cell is 16
% pixels, so we use 9 trackpaths(in row 23) and the distance between each
% trackpath is 8 pixels(in row 25).
% Written by Yuan Chen,Nanjing university of aeronautic and astronautic(2010)
%-------------------------------------------------------------------------
clear all;
%% read the avi file, and calculate the mean image;
video = {video.cdata};                  % exracte the avi data(matrix);
z = 0;
for i=1:length(video);                  % change 'color' video to gray;
video{i}=rgb2gray(video{i});
z = z + double(video{i});
end
z = z/length(video);                    % the mean image;
[mz,nz] = size(z);
[mu1,mu2,v1x,v1y,v2x,v2y] = eigfunction(z,7,5);  % calculate the eigen values and eigen vectors of the mean image;

%% the trackpaths generation;
N = 9;                                  % THE number of trackpath can be set within [5,65];NOTE that,larger N will consume
% much more processing times,we suggrest that N = 9, and no larger than 16;
t = 64/(N-1);                           % the distance between each trakcpath;
T = 64:t:128;                           % T are the radius of the trackpaths(the center(148,150)in the mean image z);
for k = 1:length(T)
tmp = zeros(size(z));
for i = 1:138
for j = 1:300
if sqrt((148-i)^2 + (150-j)^2 ) < T(k)+1 && sqrt((148-i)^2 + (150-j)^2 ) >= T(k);
tmp(i,j) = 255;
end
end
end
tmp = im2bw(tmp);
tmp = bwmorph(tmp,'thin',inf);
trackpath{k} = tmp;
end
display([num2str(N),' trackpaths are selected to generate ',num2str(N),' ST images...Processing them may take several minutes'])
%% ST image generation, processing and traces extraction;
for i=1:size(trackpath,2)                                   % the number of the generated trackpath;
map = [];
[xx,yy] = startpoint(trackpath{i});                     % finding the starting point of a trackpath;
ind_trackpath = ord_line_indx(trackpath{i},xx,yy);      % order the trackpath points;
for j = 1:size(video,2)                                 % ST image generation;
line = video{j}(ind_trackpath);
map = [map;line];
end
map = map';
Map{i} = map;                                           % store the ST image;
display(['processing ST image',int2str(i)])
[F,Theta,angle] = J4(map,2.5);                          % raw ST image enhanced by Jacob filter;
R = nsf(F);                                             % applying noise supression function;
K = angfilter(R,Theta,angle,20);                        % applying orientation filter function;
th = graythresh(K);                                     % threshold for cell tracked trackpath;
if th <= 0.01;                                          % threshold for none cell tracked trackpath;
th = 0.01;
end
bw = im2bw(K,th);
thin = bwmorph(bw,'thin',inf);                          % thinning extracted traces;
trace{i} = bwareaopen(thin,10);                         % remove noise trace if its length smaller than 10 pixels;
end

%% calculate the trackpaths coordinate and grayscale;
h = fspecial('gaussian',7);
for i = 1:length(trackpath)                                 % denoise;
Map_denoise{i} = conv2(Map{i},h,'same');
Map_denoise{i} = conv2(Map_denoise{i},h,'same');
end
% -------find out the coordinate(X_trackpath,Y_trackpath) and grayscale of the trackpaths;
for i = 1:length(trackpath)
[xx,yy] = startpoint(trackpath{i});                     % find the start coordinate of a trackpath;
ind_trackpath = ord_line_indx(trackpath{i},xx,yy);      % search the trackpath's point from the start point orderly;
for j = 1:length(ind_trackpath)
[X_trackpath(j),Y_trackpath(j)] = ind2sub([mz,nz],ind_trackpath(j)); % change the index to (X,Y) coordinate;
end
Trackpath_value{i}(:,1) = X_trackpath';
Trackpath_value{i}(:,2) = Y_trackpath';
Trackpath_value{i}(:,3) = 0;
[Y_trace,X_trace] = find(trace{i}==1);                  % find extracted trace coordinate;
for k = 1:length(Y_trace)                               % store the grayscale of the trace;
Trackpath_value{i}(Y_trace(k),3) = Map_denoise{i}(Y_trace(k),X_trace(k));
end
clear X_trackpath;clear Y_trackpath;
end
% -------backward mapping the extracted traces to the trackpaths;
for i = 1:length(Trackpath_value)
for j = 1:size(Trackpath_value{i},1)
z(Trackpath_value{i}(j,1),Trackpath_value{i}(j,2)) = Trackpath_value{i}(j,3);
end
end
% figure,imshow(z,[]),title('The tracking state of trackpaths')

%% trackpaths alignment
k = round(length(T)/2);                                     % the centerline trackpath;
% -------turn all trackpaths and all ST images to the same length(or size) as the centreline trackpath and centreline ST image
% according to the eigen victor v1x,v1y;
for i = 1:length(T)
p = findpoints(trackpath{k},trackpath{i},v1x,v1y);      % the corresponding relationship between the k and the i trackpaths.
P{i} = round(smooth(p));
if i~=k
for j = 1:length(P{i})
trace_aligned{i}(j,:) = trace{i}(P{i}(j),:);
Points_aligned{i}(j,:) = Trackpath_value{i}(P{i}(j),:);
end
else
trace_aligned{i} = trace{i};
Points_aligned{i} = Trackpath_value{i};
end
end

% -------base trace generation;
B_trace = 0;
for i = 1:length(trace_aligned)
B_trace = B_trace + trace_aligned{i};
end
% figure,imshow(B_trace),title('the result of traces mapping')
B_trace = im2bw(B_trace);
B_trace = bwmorph(B_trace,'dilate');
B_trace = bwmorph(B_trace,'thin',inf);
% figure,imshow(B_trace),title('Base-trace')

%%  fusion of aligned points in trackpaths to calculate the trajectory;
l = 1;
for j = 1:size(Points_aligned{1},1)             % the number of the aligned points in the trackpath;
Fusion_points = [];
for i = 1:size(Points_aligned,2)            % the number of the trackpaths;
if Points_aligned{i}(j,3) > 0           % >0 means the cell is tracked;
X = Points_aligned{i}(j,1);         % X coordinate;
Y = Points_aligned{i}(j,2);         % Y coordinate;
V = Points_aligned{i}(j,3);         % grayvalue of (X,Y);
Fusion_points = [Fusion_points,[X;Y;V]];
end
end
if size(Fusion_points,2) > 2                % the cell is tracked by three or more than three trackpths;
X_fused(l,1) = mean(Fusion_points(1,:));% the fused position is the center among the aligned points;
Y_fused(l,1) = mean(Fusion_points(2,:));
l = l+1;
elseif size(Fusion_points,2) == 2           % the cell is tracked by two trackpaths(Eq.(15)in the paper);
x1 = Fusion_points(1,1); x2 = Fusion_points(1,2);
y1 = Fusion_points(2,1); y2 = Fusion_points(2,2);
v1 = Fusion_points(3,1); v2 = Fusion_points(3,2);
[X_fused(l,1),Y_fused(l,1)] = fusion2p(x1,y1,x2,y2,v1,v2,0.4); % the fused position calculation;
l = l+1;
elseif size(Fusion_points,2) ==1            % the cell is tracked by only one trackpath;
X_fused(l,1) = Fusion_points(1,1);      % no fusion is needed;
Y_fused(l,1) = Fusion_points(2,1);
l = l+1;
end
end

X_trajectory = round(X_fused);               % calculated trajectory;
Y_trajectory = round(Y_fused);
X_smooth = round(smooth(X_fused,10));        % smooth process to prevent jagged trajectory;
Y_smooth = round(smooth(Y_fused,10));
for i = 1:length(X_fused)
z(X_trajectory(i),Y_trajectory(i)) = 255;
end
figure,imshow(z,[]);title(['Blood cell tracking by ',num2str(N), ' trackpaths(black line)'])
hold on;plot(Y_smooth,X_smooth,'r')

% -------the real trajectory of the cell;
x = 76:228;
y = 124*sin(pi/162*x+80.32);
y = round(-y+151);
hold on,plot(x,y);
h = legend('Tracked trajectory','Real trajectory',2);

%% calculating the error and variance;
for i = 1:length(X_trajectory)                  % by finding the closest points in the real trajectory to calculate the error;
for j = 1:length(x)
dis(j) = sqrt((X_trajectory(i) - y(j))^2 + (Y_trajectory(i) - x(j))^2);
Error(i) = min(dis);
end
end
ES = [mean(Error),std2(Error)];                 % mean error and standard deviation;
display(['The Error and Variance of the tracked tracjectory by using ',num2str(N),' trackpaths are'])
display([num2str(ES(1)) ' pixels and ' num2str(ES(2)),' pixels'])```