Code covered by the BSD License

# Simple cell tracking/contour for extensional flow

### Kevin Roth (view profile)

This script tracks two cells in an extensional flow field while measuring the contours.

ExtensionalFlowCellContour.m
%Image Processing: Full .avi file analysis
%Part 2 of 2 of full .avi analysis
%Kevin Roth
%2012-02-01

%This script is my first foray into MATLAB and image processing in
%general.  The particle tracking is not too complex, and the script relies
%on only having two cells in the frame at a time.

%If you have any issues, feel free to email me and I'd be glad to help.

%Instructions:

%First: run 'ReadAndCropVideo.m' to decide how to crop video for proper
%analysis (need to crop the video such that only two cells are in the frame
%at all times)

%Before running 'IPFullAviAnalysisUpdate2.m,' RE-SAVE IT IN THE NEW FOLDER.

%And then...Change the values for cropping to the correct values from
%'ReadAndCropVideo.m' and run program (line **80**).

%Specify if the incoming cell is the top cell or the bottom cell by
%defining z, line **69**.

%And then...Watch it go.

%Change line **1978** to delta_yp < ___. (use judgement, pick a value that
%represents when the cells come into contact and begin deforming)

%The equivalent diameter is used for the deformation
%values, as well as for calculating the relative velocity of cells incoming
%for the drag force

%And then...see the bottom of this file, open cftool and analyze data.

%And then...change the values for k and \eta at the bottom and re-run the
%last bit of code to get a usable graph with data and function.  Then SAVE.

%It should be noted that the resolution of the objective used for imaging
%was 10 pixels/micron.  If the resolution needs to be changed, find the
%spots in the code where the units are converted from pixels to microns.

%Stokes' drag is used to calculate the drag force on the incoming cell.

%The Kelvin-Voigt model is the viscoelastic model described below in the
%curve fitting section at the end of the script.

close all
clearvars
clc

nFrames = cellmov.NumberOfFrames;
vidHeight = cellmov.Height;
vidWidth = cellmov.Width;
fRate = cellmov.FrameRate;
cnt = 1;
cnt2 = 1;
cnt3 = 1;
%Name variables

dragCell = 6;
%dragCell = 4  bottom cell is the incoming cell
%dragCell = 6  top cell is the incoming cell

for j = 1:nFrames

%First, perform transformations on images

cropimage = imcrop(imageData, [_ _ _ _]);
%Crop image
grayData = rgb2gray(cropimage);
%Converts RGB image to grayscale
level = graythresh(grayData);
%Automatically computes an appropriate threshold to convert
%grayscale image to binary
bw = im2bw(grayData,level);
%Converts image to binary image based on calculated threshold
%'level'
bw = bwareaopen(bw,50);
%Removes background noise
I3 = edge(bw);
%Find edges of a binary image (1's and 0's)
I4 = imfill(I3,'holes');
%Fill image regions and holes
I5 = bwareaopen(I4,2000);
%Morphologically open binary image (remove small objects)
BWoutline = bwperim(I5);
%Outline on original image with the perimeter generated from the
%image processing
Segout = cropimage;
Segout(BWoutline) = 255;

figure(1), imshow(Segout);
%Show the outline of the image processing function on the original
%images while the program is running

cc = bwconncomp(I5,8);
%Find connected components in a binary image
cc.NumObjects;
%Stores the number of connected components in that frame in the
%structure cc

%IF there are 2 cells per image, then perform and store numerical
%analysis on frame

if cc.NumObjects > 1
t2(cnt) = j;
cnt = cnt+1;
end

%This short if statement is for defining a vector containing only
%the frames where both cells get detected.

if cc.NumObjects > 1
cell2 = false(size(I5));
cell2(cc.PixelIdxList{2}) = true;
%figure(2), imshow(cell2);
%Show the second cell when present

labeled = labelmatrix(cc);
%Creates a label matrix from bwconncomp structure
graindata = regionprops(cc,'Area','EquivDiameter',...
'MajorAxisLength','MinorAxisLength','Centroid','Orientation');
%Measures properties of image regions AND SAVES THEM!!!

if graindata(1).Centroid(2) < graindata(2).Centroid(2)
topCell(1,j) = graindata(1,1);
bottomCell(1,j) = graindata(2,1);
end

%This if statement says that if the y position of the
%centroid of the area that MATLAB labels as #1 is less
%(higher on the screen) than the y position of area labeled
%as #1, then the data will save the #1 area as
%the top cell and the #2 area the bottom cell.

if graindata(1).Centroid(2) > graindata(2).Centroid(2)
topCell(1,j) = graindata(2,1);
bottomCell(1,j) = graindata(1,1);
end

%This if statement says that if MATLAB flips the labeling
%of areas #1 and #2, to save the area with the lower value
%of y-position (in this case, area #2) as topCell.  This
%should eliminate any switching of the labeling for the
%tracking...

topCellMinD(j) = topCell(1,j).MinorAxisLength/10;   %um
topCellMajD(j) = topCell(1,j).MajorAxisLength/10;   %um
topCellArea(j) = topCell(1,j).Area/100;             %um^2
topCellEquivD(j) = topCell(1,j).EquivDiameter/10;   %um
topCelly(j) = topCell(1,j).Centroid(2)/10;          %um
topCellx(j) = topCell(1,j).Centroid(1)/10;          %um
topCellDeformation(j) = (topCellMajD-topCellMinD)/...
(topCellMajD+topCellMinD);                      %unitless
topCellOrient(j) = topCell(1,j).Orientation;        %degrees

bottomCellMinD(j) = bottomCell(1,j).MinorAxisLength/10;
bottomCellMajD(j) = bottomCell(1,j).MajorAxisLength/10;
bottomCellArea(j) = bottomCell(1,j).Area/100;
bottomCellEquivD(j) = bottomCell(1,j).EquivDiameter/10;
bottomCelly(j) = bottomCell(1,j).Centroid(2)/10;
bottomCellx(j) = bottomCell(1,j).Centroid(1)/10;
bottomCellDeformation(j) = (bottomCellMajD-bottomCellMinD)/...
(bottomCellMajD+bottomCellMinD);
bottomCellOrient(j) = bottomCell(1,j).Orientation;

time(j) = j/fRate; %s

%This sequence of commands tells MATLAB to save each part
%of the structure as it's own array-so I only have to deal
%with 1x(number of frames analyzed) arrays instead of
%structures.

%With the current optical setup (40x objectives), there are
%10 pixels/micrometer.  In order to put values in microns,
%each length must be divided by 10 (pixels).  For area, one
%must divide by 100: (100 pixels^2/micron^2).

%The time for each frame is the number of frame divided by
%the frame rate (fRate).

if (bottomCelly(j)-topCelly(j)) < 7.5
t3(cnt2) = j;
cnt2 = cnt2+1;
end

%This if statement tells MATLAB that I want to know which
%frames the y position difference between the cells is less
%than 10 um
%t3 is used for deciding when the collision begins (for the
%curve fitting)

if (bottomCelly(j)-topCelly(j)) < 14;
%Use 14um here, standardized across videos
t4(cnt3) = j;
cnt3 = cnt3+1;
end

%This if statement tells MATLAB to define a new time scale when
%the difference between the two cells is between between
%(Equivalent Diameter of stagnant cell)+(Equivalent Diameter of
%incoming cell).  This is for the drag force calculation.
%Also, going to get used for the calculation of the average
%diameter of cells, for the first ten frames deal (replacing
%t2).

else
%figure(2), imshow(cropimage)

%Use this to see when the second cell is not present so you can
%tell the difference between frames (only for watching while
%program runs).

end

end

%It should be noted that MATLAB will save data in the column number with
%the same frame number it is taken from.  If two cells are not detected in
%a frame, then MATLAB does not save this frame for data, and will fill the
%cell with a zero. If there are no more 2 cell frames detected, the program
%will not fill those cells in with zeros.  For example, using m19 from
%2011-01-17, there are 452 frames, and MATLAB only records data from 426 of
%them because after frame 426 a second cell is not detected, so no more
%data is taken.  But looking through the data, there are zeros where there
%wasn't a cell detected, so MATLAB is NOT compressing the data and ignoring
%frames where two cells are not present.

%Plotting Data

%figure(3), plot(time,topCellMinD,'.',time,topCellMajD,'x');
%xlabel('Time (s)'),ylabel('Cell Diameter (\mum)'),title('Top Cell');
%legend('Minor Axis','Major Axis'),axis([0 nFrames/fRate 5.5 8.0]);

%figure(4), plot(time,bottomCellMinD,'.',time,bottomCellMajD,'x');
%xlabel('Time (s)'),ylabel('Cell Diameter (\mum)'),title('Bottom Cell');
%legend('Minor Axis','Major Axis'),axis([0 nFrames/fRate 5.5 8.0]);

figure(5), plot(time,topCelly,'.',time,bottomCelly,'x');
xlabel('Time (s)'),ylabel('y-Coordinate (/mum)'),title('Cell Position');
legend('Top Cell','Bottom Cell');

figure(6), plot(time,topCellEquivD,'.',time,bottomCellEquivD,'x');
xlabel('Time (s)'),ylabel('Equiv D (/mum)'),title('Cell Equiv D');
legend('Top Cell','Bottom Cell');

%Data Analysis

%By adding the if statement in the for loop above defining t2 to be the
%frames when there are two cells only, we are able to eliminate the
%need to look for the first ten good frames in the video to take the
%mean.  It's now automatic.

topCellEquivD_mean = mean(topCellEquivD(t4(1:10)));
bottomCellEquivD_mean = mean(bottomCellEquivD(t4(1:10)));
%units: um
%changed this mean calculation time scale, but the rest of the data
%should be t2 dependent, because that is the time scale where both
%cells are present in the frame.  This will just help eliminate the
%funky measurements from cell rotation.  Now, it will be the ten frames
%as the drag is being calculated, so shouldn't be any issue with
%deforming in that range as well.

%Calculate the mean value of the first ten measurements of the major
%and minor axis diameters of both cells

%Percent Deformation

topCellEquivDef = (topCellEquivD(t2)-topCellEquivD_mean)/...
topCellEquivD_mean*100;
topCellEquivDef = (topCellEquivD(t2)-topCellEquivD_mean)/...
topCellEquivD_mean*100;

time2 = t2/fRate;
t3_shift = t3-t3(1);
time3=  t3_shift/fRate;
%t2 is the frame number, this will give time2 units of time (s)
%t2 defines the frame numbers where both cells get detected
%t3 defines the frames where the cell's y position difference is less
%than 10 um

%Percent Deformation = [a(0)-a(t)]/a(0)*100%

%figure(6), plot(time2,topCellMinDef,'.',time2,topCellMajDef,'x');
%xlabel('Time (s)'),ylabel('Percent Deformation');
%title('Top Cell Deformation'),legend('Minor Axis','Major Axis');
%axis([0 nFrames/fRate -15 15]);

%figure(7), plot(time2,bottomCellMinDef,'.',time2,bottomCellMajDef,'x');
%xlabel('Time (s)'),ylabel('Percent Deformation');
%title('Bottom Cell Deformation'),legend('Minor Axis','Major Axis');
%axis([0 nFrames/fRate -15 15]);

%Alternative Cell Deformation Measurement

figure(8), plot(time,topCellDeformation,'.',time,...
bottomCellDeformation,'x');
xlabel('Time (s)'),ylabel('[long-short]/[long+short]');
title('Cell Deformation/Strain?'),legend('Top Cell','Bottom Cell');
%axis([0 nFrames/fRate 0 0.1]);

%Orientation

figure(9), plot(time,topCellOrient,'.',time,bottomCellOrient,'x');
xlabel('Time (s)'),ylabel('Angle (Degrees)');
title('Major Axis Orientation vs. Horizontal Line');
legend('Top Cell','Bottom Cell');
%axis([0 nFrames/fRate -90 90]);

%Delta x; like spring

topCellDeltaX = (topCellEquivD_mean-topCellEquivD(t2))*10^-6;
bottomCellDeltaX = (bottomCellEquivD_mean-bottomCellEquivD(t2))*10^-6;
topCellD_mean = topCellEquivD_mean;
bottomCellD_mean = bottomCellEquivD_mean;
topCellD = topCellEquivD;
bottomCellD = bottomCellEquivD;
%units of the DeltaX: meters
%units of the mean and D: um
%Only need to change the variable from Min diameter to Maj diameter in
%these variables, they get referenced throughout the rest of the
%program.

figure(10), plot(time2,topCellDeltaX,'.',time2,bottomCellDeltaX,'x');
xlabel('Time(s)'),ylabel('delta X (m)');
title('Spring-like behavior of a super sweet RBC');
legend('Top Cell','Bottom Cell');

%Percent Deformation for curve fitting

%AbsbottomCellMinDef = abs(bottomCellMinDef);
%AbsbottomCellMajDef = abs(bottomCellMajDef);
%AbstopCellMinDef = abs(topCellMinDef);
%AbstopCellMajDef = abs(topCellMajDef);

%Drag Force Calculation

if dragCell < 5
bottomCellVelocity = abs(bottomCelly(t4(10))-bottomCelly(t4(1)))/...
((t4(10)-t4(1))/fRate)*10^-6;
%Convert cell velocity to units of m/s from um/s
%This is only the y velocity of the cell, it's the important value.
%And there is very little x movement.  Take absolute value because the
%cell is travelling in the negative y direction

viscosity = 0.001;
%units: kg/m*s or Pa*s, it's 0.01 P = 1 cP = 0.001 Pa*s

rcell = (bottomCellEquivD_mean/2)*10^-6;
%units of m from um
%Average cell diameter over first ten frames
%Actually a RADIUS, divide the equivalent diameter by 2

dragForce = 6*pi*viscosity*rcell*bottomCellVelocity;
%units of N
end

if dragCell > 5
topCellVelocity = abs(topCelly(t4(10))-topCelly(t4(1)))/...
((t4(10)-t4(1))/fRate)*10^-6;
%Convert cell velocity to units of m/s from um/s
%This is only the y velocity of the cell, it's the important value.
%And there is very little x movement.  Take absolute value because the
%cell is travelling in the negative y direction

viscosity = 0.001;
%units: kg/m*s or Pa*s, it's 0.01 P = 1 cP = 0.001 Pa*s

rcell = (topCellEquivD_mean/2)*10^-6;
%units of m from um
%Average cell diameter over first ten frames
%Actually a RADIUS, divide the equivalent diameter by 2

dragForce = 6*pi*viscosity*rcell*topCellVelocity;
%units of N
end

%Curve Fitting

topCellDeltaX_Fd = topCellDeltaX/dragForce;
bottomCellDeltaX_Fd = bottomCellDeltaX/dragForce;
%Normalize data to eliminate one variable

yPositionDifference = -topCelly(t2)+bottomCelly(t2);
yPosDif = -topCelly(t3)+bottomCelly(t3);

%figure(11), plot(time2,yPositionDifference,'.',time3,yPosDif,'x');
%xlabel('Time (s)'),ylabel('topCelly-bottomCelly');

topDeltaX = (topCellD_mean-topCellD(t3))*10^-6;
bottomDeltaX = (bottomCellD_mean-bottomCellD(t3))*10^-6;
%Units of m
topDeltaX_Fd = topDeltaX/dragForce;
bottomDeltaX_Fd = bottomDeltaX/dragForce;

figure(12), plot(time3,topDeltaX,'.',time3,bottomDeltaX,'x');
xlabel('Time(s)'),ylabel('Delta X (\mum)')
title('Delta X for \DeltaYpos<9')

%Instructions for cftool

%Start by opening the curve fitting tool, type cftool in command
%window, press enter.

%Create data sets using time3 for the x and (bottom/top)DeltaX_Fd
%for the y.

%Choose "Fitting," and create a new fit.  Under "Type of fit" choose
%"Custom Equations."  Use 1/k*(1-exp(-k/n*x)).  Need good guesses for the
%constants.  Use: k = 1e-6 and n = 1e-6.  Apply fit.

%NEED TO CHANGE VALUES DEPENDING ON VALUES GIVEN BY CFTOOL

%TOP CELL
k_top = 0;%Units: N/m
k_top_upper = 0; %This is for the 95% confidence bounds
k_top_lower = 0;
n_top = 0; %Units: N*s/m
n_top_upper = 0;
n_top_lower = 0;
R2_top = 0;
SSE_top = 0;
RMSE_top = 0;
expFit_top = 1/k_top*(1-exp(-k_top/n_top*time3));
expFit_top_upper = 1/k_top_upper*(1-exp(-k_top_upper/n_top_upper*time3));
expFit_top_lower = 1/k_top_lower*(1-exp(-k_top_lower/n_top_lower*time3));

figure(13), plot(time3,topDeltaX_Fd,'.',time3,expFit_top,'-',time3,...
expFit_top_upper,'--',time3,expFit_top_lower,'-.');
xlabel('Time (s)'),ylabel('\Delta(Equivalent Diameter)/Drag Force (m/N)');
title('Top Cell Equivalent Diameter in Cell Collision');
text(0.75,4e4,'\it1/k[1-\ite^{-k/\etat}]','FontSize',18);
legend('Experimental Data','Exponential Fit','95% Upper','95% Lower');

%BOTTOM CELL
k_bottom = 0;
k_bottom_upper = 0;
k_bottom_lower = 0;
n_bottom = 0;
n_bottom_upper = 0;
n_bottom_lower = 0;
R2_bottom = 0;
SSE_bottom = 0;
RMSE_bottom = 0;
expFit_bottom = 1/k_bottom*(1-exp(-k_bottom/n_bottom*time3));
expFit_bottom_upper = 1/k_bottom_upper*(1-exp(-k_bottom_upper/...
n_bottom_upper*time3));
expFit_bottom_lower = 1/k_bottom_lower*(1-exp(-k_bottom_lower/...
n_bottom_lower*time3));

figure(14), plot(time3,bottomDeltaX_Fd,'.',time3,expFit_bottom,'-',...
time3,expFit_bottom_upper,'--',time3,expFit_bottom_lower,'-.');
xlabel('Time (s)'),ylabel('\Delta(Equivalent Diameter)/Drag Force (m/N)');
title('Bottom Cell Equivalent Diameter in Cell Collision');
text(0.75,4e4,'\it1/k[1-\ite^{-k/\etat}]','FontSize',18);
legend('Expermental Data','Exponential Fit','95% Upper','95% Lower');