%-------------------------------------------------------------------------- % MATLAB code that demonstrates the tracking of a single cell % % © Technical University of Denmark 2022 % % The code was applied to experimental data presented in the paper: % % Acoustic Tethering of Microorganisms, % M. Rode, A. Bioue, F. Miano, H. Bruus, T. Kiørboe and A. Andersen, % Journal of Experimental Biology 225, jeb244089 (2022). % % https://doi.org/10.1242/jeb.244089 %-------------------------------------------------------------------------- % % In the example we track a cell for 10 frames immediately after the % acoustic field is turned on and the cell is brought into the focus plane. % % The tracking is done by a combination of a cross-correlation between % consecutive video frames and an imfindcircle adjustment. % % An image is displayed for each frame when running the script. % % The left panel shows the identified cell position in the previous frame % (star, red) and the current cross-correlation window (square, red). % % The middle panel shows the region (square, purple) within which the % cross-correlation window is allowed, the center of the cross-correlation % window with the maximum cross-correlation (filled circle, orange), and % the circle which best fits the cell (circle, blue). % % The cell position is identified as the center of the circle which best % fits the cell (plus sign, blue) if its distance to the maximum % cross-correlation position is below a threshold, and otherwise the cell % position is identified as the maximum cross-correlation position. % % The right panel shows a table with the identified cell positions. %-------------------------------------------------------------------------- % defining where and in which frame to start the tracking trackSpec.fromFrame = 1; % frame to start tracking in the tiff-file trackSpec.row = 168; % row-position of the cell in the fromFrame trackSpec.column = 866; % column-position of the cell in the fromFrame trackSpec.toFrame = 10; % do the tracking until this frame % controling colors xcorrColor =[0.8500 0.3250 0.0980]; sourceColor =[0.6350 0.0780 0.1840]; targetColor =[0.4940 0.1840 0.5560]; usedcircleColor =[0 0.4470 0.7410]; unusedcircleColor =[0.9290 0.6940 0.1250]; % tiff file with video frames tiffFile = 'particle_tracking_video_frames.tif'; % parameters imfindcircle imfindcircleParam.ObjectPolarity = 'Dark'; imfindcircleParam.Sensitivity = 0.92; imfindcircleParam.Method = 'TwoStage'; % parameters tracking trackingParam.row = trackSpec.row; % actual position trackingParam.column = trackSpec.column; % actual position trackingParam.radius = 20; % initial particle radius to look for trackingParam.min_radius = 8; % don't go below this particle radius trackingParam.bSize = 15; % boundary size, next video frame xcorr-window % table for storing the identified track positions varNames = {'frame','row','column'}; varTypes = {'int16','double','double'}; sz=[trackSpec.toFrame - trackSpec.fromFrame + 1, length(varNames)]; positions = table('Size',sz,'VariableTypes',varTypes,'VariableNames',varNames); % no relative translation in first frame translationRow=0; translationColumn=0; % prepare loop image = imread(tiffFile, trackSpec.fromFrame); % get the first video frame % looking for a circle in the initial frame [circleRow, circleColumn, circleRadius] = ... findCircleCenter(image, trackingParam, imfindcircleParam); if ~isnan(circleColumn) trackingParam.row = circleRow; trackingParam.column = circleColumn; trackingParam.radius = circleRadius; end idx = 1; % index to position table applyFindcircle=true; % apply imfindcircle if useful % saving found position to position-table positions.frame(idx)=trackSpec.fromFrame; positions.row(idx)=trackingParam.row; positions.column(idx)=trackingParam.column; % ready for looping through the frames for frameNbr = trackSpec.fromFrame:trackSpec.toFrame-1 % figure to show the frames involved fig=figure('Position', [100 200 900 300]); sp1=subplot(1,3,1); % source video frame sp2=subplot(1,3,2); % next video frame % show the source video frame in sup-plot 1 imagesc(sp1,image) colormap(sp1,'gray') axis(sp1,'equal') hold(sp1,'on') % asterix - the identified position on the sorce frame scatter(sp1, trackingParam.column, trackingParam.row, 96,... '*', 'MarkerEdgeColor', sourceColor) % limit the image view xlim(sp1, [775 925]) ylim(sp1, [125 275]) % indicate source cross-correlation window rectangle(sp1, 'Position', [floor(trackingParam.column)-trackingParam.radius,... floor(trackingParam.row)-trackingParam.radius,... trackingParam.radius*2+1, trackingParam.radius*2+1],... 'EdgeColor', sourceColor, 'LineWidth', 1) % the video frame shown in subplot 1 title(sp1, sprintf('Frame %i', frameNbr)) % remember last relative translation lastTranslationRow = translationRow; lastTranslationColumn = translationColumn; % do cross-correlation with next frame [image, translationRow, translationColumn] = ... xcorr_frame(tiffFile, image, frameNbr, trackingParam,... lastTranslationRow, lastTranslationColumn); % NB: image has now changed to the next frame in the tiff file % show target cross-correlation window imagesc(sp2, image) colormap(sp2,'gray') axis(sp2,'auto', 'equal', 'tight') hold(sp2,'on') rectangle(sp2, 'Position', [floor(trackingParam.column)-trackingParam.radius-trackingParam.bSize,... floor(trackingParam.row)-trackingParam.radius-trackingParam.bSize,... 2*(trackingParam.radius + trackingParam.bSize) + 1,... 2*(trackingParam.radius + trackingParam.bSize) + 1],... 'EdgeColor', targetColor, 'LineWidth', 1) % adjust the actual position with the cross-correlation trackingParam.row = trackingParam.row + translationRow; trackingParam.column = trackingParam.column + translationColumn; % show the new position found by cross-correlation scatter(sp2, trackingParam.column, trackingParam.row,... 'o','filled',... 'MarkerEdgeColor', xcorrColor,... 'MarkerFaceColor', xcorrColor,... 'LineWidth', 1) % used imfindcircle to keep a more precise center of cell if applyFindcircle % looking for a circle [circleRow, circleColumn, circleRadius] = ... findCircleCenter(image, trackingParam, imfindcircleParam); % check if result from imfindcircle should be used if ~isnan(circleRow) if (sqrt((circleColumn-trackingParam.column)^2+(circleRow-trackingParam.row)^2)... < trackingParam.radius*0.5) % circle center within half a radius - use it trackingParam.row = circleRow; trackingParam.column = circleColumn; trackingParam.radius = ... % use indetified radius in next frame, max(circleRadius,trackingParam.min_radius); % imfindcircle result used circleColor = usedcircleColor; else % imfindcircle result NOT used circleColor = unusedcircleColor; end % show identified circle center scatter(sp2, circleColumn, circleRow, 96,... '+', 'MarkerEdgeColor', circleColor, 'LineWidth', 1) % show identified circle viscircles(sp2, [circleColumn circleRow], circleRadius,... 'LineWidth' , 0.2,... 'Color', circleColor, 'LineWidth', 1); end end % the video frame shown in subplot 2 title(sp2, sprintf('Frame %i', frameNbr + 1)) %limit the image view xlim(sp2, [775 925]) ylim(sp2, [125 275]) % next index in the particle-position table idx = idx + 1; % saving found position to position table positions.frame(idx)=frameNbr+1; positions.row(idx)=trackingParam.row; positions.column(idx)=trackingParam.column; % show the identified positions so far uitable(fig,'Data',positions{:,:},'ColumnName',positions.Properties.VariableNames,... 'RowName',positions.Properties.RowNames,'Units', 'Normalized',... 'Position',[0.70, 0.15, 0.25, 0.70]); drawnow end % cross-correlation between consecutive images function [nextImage, relaticeTranslationRow, relativeTranslationColumn] = ... xcorr_frame(fileName, previousImage, frameNbr,... trackingParam, lastTranslationRow, lastTranslationColumn) % previous video frame cross-correlation window fromRowP = floor(trackingParam.row)... - trackingParam.radius; toRowP = floor(trackingParam.row)... + trackingParam.radius; fromColP = floor(trackingParam.column)... - trackingParam.radius; toColP = floor(trackingParam.column)... + trackingParam.radius; previousWindow... = double(previousImage(fromRowP:toRowP,... fromColP:toColP)); % getting the next frame from the video file nextImage = imread(fileName, frameNbr+1); % next video frame cross-correlation window fromRowS = floor(trackingParam.row+lastTranslationRow)... - trackingParam.radius - trackingParam.bSize; toRowS = floor(trackingParam.row+lastTranslationRow)... + trackingParam.radius + trackingParam.bSize; fromColS = floor(trackingParam.column+lastTranslationColumn)... - trackingParam.radius - trackingParam.bSize; toColS = floor(trackingParam.column+lastTranslationColumn)... + trackingParam.radius + trackingParam.bSize; nextWindow... = double(nextImage(fromRowS:toRowS,... fromColS:toColS)); cc = normxcorr2(previousWindow,nextWindow); % extract the relative translation from the cross-correlation [~, imax2] = max(abs(cc(:))); [rPeak2, cPeak2] = ind2sub(size(cc),imax2(1)); relaticeTranslationRow = rPeak2-size(previousWindow,1)... - trackingParam.bSize+lastTranslationRow; relativeTranslationColumn = cPeak2-size(previousWindow,2)... - trackingParam.bSize+lastTranslationColumn; end % use imfindcircle to find possible center of the particle function [circleRow, circleColumn, circleRadius] = ... findCircleCenter(image, trackingParam, imfindcircleParam) % defining window for imfindcricle fromRow = floor(trackingParam.row)... - trackingParam.radius - trackingParam.bSize; toRow = floor(trackingParam.row)... + trackingParam.radius + trackingParam.bSize; fromCol = floor(trackingParam.column)... - trackingParam.radius - trackingParam.bSize; toCol = floor(trackingParam.column)... + trackingParam.radius + trackingParam.bSize; circleWindow = double(image(fromRow:toRow,... fromCol:toCol)); % look for circle [centers, radii, metric] = imfindcircles(circleWindow,... [floor(trackingParam.radius*0.75) ceil(trackingParam.radius*1.25)],... 'ObjectPolarity',imfindcircleParam.ObjectPolarity,... 'Sensitivity',imfindcircleParam.Sensitivity,... 'Method' ,imfindcircleParam.Method); cSize=size(centers); if cSize(1)>0 % returning the best circle index = find(metric==max(metric)); circleRow = fromRow-1+centers(index(1),2); circleColumn = fromCol-1+centers(index(1),1); circleRadius = radii(index(1)); else circleRow = NaN; circleColumn = NaN; circleRadius = NaN; end end