Code covered by the BSD License  

Highlights from
Simple cell tracking/contour for extensional flow

image thumbnail

Simple cell tracking/contour for extensional flow



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

%Image Processing: Full .avi file analysis
%Part 2 of 2 of full .avi analysis
%Kevin Roth

%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.


%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

cellmov = mmreader('*.avi');      
    %Read in video
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
    imageData = read(cellmov,j);
        %Read in image
    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
    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
        %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;
        %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',...
            %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);            
                %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);
                %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
            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)/...
            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
                %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;
            %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;
            %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
        %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).

%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)/...
    topCellEquivDef = (topCellEquivD(t2)-topCellEquivD_mean)/...

    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,...
xlabel('Time (s)'),ylabel('[long-short]/[long+short]');
title('Cell Deformation/Strain?'),legend('Top Cell','Bottom Cell');
%axis([0 nFrames/fRate 0 0.1]);

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
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)))/...
        %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

if dragCell > 5
    topCellVelocity = abs(topCelly(t4(10))-topCelly(t4(1)))/...
        %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

%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.

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;
AdjustedR_Square_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,...
xlabel('Time (s)'),ylabel('\Delta(Equivalent Diameter)/Drag Force (m/N)');
title('Top Cell Equivalent Diameter in Cell Collision');
legend('Experimental Data','Exponential Fit','95% Upper','95% Lower');

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;
AdjustedR_Square_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/...
expFit_bottom_lower = 1/k_bottom_lower*(1-exp(-k_bottom_lower/...

figure(14), plot(time3,bottomDeltaX_Fd,'.',time3,expFit_bottom,'-',...
xlabel('Time (s)'),ylabel('\Delta(Equivalent Diameter)/Drag Force (m/N)');
title('Bottom Cell Equivalent Diameter in Cell Collision');
legend('Expermental Data','Exponential Fit','95% Upper','95% Lower');

Contact us