s%"Flow analyzer." Given a movie showing a signal that changes over time, %this will perform an analysis that gives you the position and vector %scomponents of the optical flow throughout the movie. Plotting the flow %heatmap is optional but useful for visualization. For large movies, the %grid size (gSz) may be increased so that, instead of sampling every pixel, %a gSz-by-gSz window is averaged and assigned to the center point of the %window. Finally, if a labeled image (in which the region belonging to cell %1 has a pixel value of 1, the region of cell 2 has pixel value 2, etc.) is %present, this will record a separate structure that contains the flow %results within each individual cell. It will also filter and plot the mean %optical flow within each cell over time and identify the point of %activation by calculating a differential based on that signal and finding %the maximum. Written by Erik N. Schaumann at the University of Chicago in %the groups of Margaret L. Gardel and Bozhi Tian. %Citation: Rotenberg, MY; Yamamoto, N; Schaumann, EN; Matino, L; Santoro, F; Tian, B. Living myofibroblast-silicon composites for probing in vivo electrical coupling in cardiac systems. Proc. Natl. Acad. Sci. 2019. DOI: 10.1073/pnas.1913651116 %Hardcoded variables fName = 'dFoF'; %Name of the .avi file you want to analyze fLabs = 'Labs.tif'; %Name of the image with uniquely and sequentially labeled cells tInt = 1; %Time interval between frames (in s) tStim = 200; %The time point at which stimulation occurred - throw out time points beneath this gSz = 1; %Gridsize - gSz^2 is the size of the region over which averaging will be performed. Example: gSz = 2 samples every second pixel, and averages the 2-by-2 region thus formed and records that data for the center point. plotHeat = 1; %Do you want to plot the heatmap of optical flow? Set to 0 if you want to save a trivial amount of runtime. diff_window = 1; %This should be ~1/2 of the time between the start and end of the initial Ca wave. save('FlanalyzerSettings','tStim','tInt','gSz','fName','fLabs','plotHeat','diff_window') %Load in video and image info v = VideoReader(strcat(fName,'.avi')); %Set up the template x1 = 1+0.5*(gSz-1):gSz:v.Width-mod(v.Width,gSz); y1 = 1+0.5*(gSz-1):gSz:v.Height-mod(v.Height,gSz); %Centerpoints for averaging [xi,yi] = meshgrid(y1,x1); pos_0 = [xi(:),yi(:)]; vec_0 = zeros(length(pos_0),2); sPtsx = 1:gSz:v.Width-mod(v.Width,gSz); sPtsy = 1:gSz:v.Height-mod(v.Height,gSz); %Set up the right points to sample [sxi,syi] = meshgrid(sPtsx,sPtsy); sPos = [sxi(:),syi(:)]; %Looks for the imgage specified in the settings - if not present, only does %the full field of view try labs = imread(fLabs); labTag = 1; numCells = max(labs(:)); catch warning('No labeled image - analysis on full field of view only') labTag = 0; end %Set up structures opticFlow = opticalFlowLK(); grandFlow = struct('pos',[],'vec',[],'max',0); cellData = struct('Time',struct('pos',[],'vec',[],'max',[],'mean',[]),'Centroid',[]); %Initialize our counter ctr = 1; %tic %The base loop that estimates the optical flow and extracts the relevant %information from individual cells, if a labeled image is available while (v.CurrentTime + 1/v.FrameRate) <= v.Duration && v.CurrentTime < 1/v.FrameRate*ctr %Read in the frames frameRGB = readFrame(v); frameGray = rgb2gray(frameRGB); %Compute the estimated optical flow and add to the compilation flowRaw = estimateFlow(opticFlow,frameGray); if plotHeat == 1 %Optional plotting of the optical flow heatmap (multiplied by 1000 and converted into 16 bits for visibility)! Im = uint16(flowRaw.Magnitude*1000); iName = strcat(fName,'_Heatmap.tif'); imwrite(Im,iName,'WriteMode','append') end %Read in the maximum flow in the entire frame grandFlow(ctr).max = max(flowRaw.Magnitude(:)); vec = vec_0; pos = pos_0; for i = 1:length(pos) vxtb = flowRaw.Vx(sPos(i,2):sPos(i,2)+gSz-1,sPos(i,1):sPos(i,1)+gSz-1);%Window of Vx's to be averaged vytb = flowRaw.Vy(sPos(i,2):sPos(i,2)+gSz-1,sPos(i,1):sPos(i,1)+gSz-1);%" but for Vy %Assign the value of the average throughout the window to the %center of the window vxtb = reshape(vxtb,[],1); vec(i,1) = mean(vxtb); vytb = reshape(vytb,[],1); vec(i,2) = mean(vytb); end %Save space by removing vectors with 0 for both components delMe = find(all(vec == 0,2)); pos(delMe,:) = []; vec(delMe,:) = []; grandFlow(ctr).pos = pos; grandFlow(ctr).vec = vec; if labTag == 1 %Get the data for individual cells and put it where it needs to be for i = 1:numCells %IndyFlow is the structure containing relevant flow information %for each individual cell. This gets concatenated to the final %cellData structure. indyFlow = struct('pos',[],'vec',[],'max',[],'mean',[]); %Calculating things for the individual cells cellIm = single(roicolor(labs,i)); %Only look at the region with the correct index indyFlow.pos = pos_0; indyX = cellIm.*flowRaw.Vx;%Set all vectors outside the cell to 0 indyX = reshape(indyX,[],1); indyY = cellIm.*flowRaw.Vy; indyY = reshape(indyY,[],1); indyFlow.vec = horzcat(indyX,indyY); delMe = find(all(indyFlow.vec == 0,2));%Remove 0 vectors indyFlow.pos(delMe,:) = []; indyFlow.vec(delMe,:) = []; indyMag = cellIm.*flowRaw.Magnitude; indyFlow.max = max(indyMag(:));%Add the maximum and mean magnitudes to the structure indyFlow.mean = mean(indyMag(:)); %Data for future cell-cell transmission cellData(i).Time(ctr) = indyFlow; if any(cellIm(:)) ~= 0 regStats = regionprops(cellIm); cellData(i).Centroid = regStats.Centroid; end end end ctr = ctr + 1; %toc end %Save the results sName = strcat(fName,'_Flowx',num2str(gSz),date); save(sName,'grandFlow','cellData') tRange = tStim:size(grandFlow,2); %Compute and plot cell-to-cell transmission, if a labeled image is available. nms = {}; ncnt = 1; for i = 1:length(cellData) if any([cellData(i).Time(tRange).mean]) %Pass the mean flow through a 1D Gaussian filter, then make the %"differential" signal given by d(i) = f(i+w) - f(i-w) where d is %the differential, f is the Gaussian filtered flow signal, and w is %the diff_window variable. Meanz = imgaussfilt(double([cellData(i).Time.mean]),diff_window); filt = Meanz(tRange); nFilt = filt/max(filt); temp_filt = [zeros(1,diff_window,1),nFilt,zeros(1,diff_window)]; dFilt = circshift(temp_filt,-diff_window) - circshift(temp_filt,diff_window); dFilt(end:-1:end-diff_window+1) = []; dFilt(1:diff_window) = []; try %Find Tmax, a.k.a the highest peak. [Pks,Ind] = findpeaks(dFilt); Tmax = find(dFilt == max(Pks)); %Plot the Gaussian filtered flow signal subplot(2,1,1) gf1(ncnt) = plot(tRange*tInt,nFilt); ga = gca; gf1(ncnt).LineWidth = 1.2; axis([0 Inf 0 1]) ga.FontSize = 14; hold on %Plot the "differential" signal subplot(2,1,2) gf2(ncnt) = plot(tRange*tInt,dFilt); ga = gca; gf2(ncnt).LineWidth = 1.2; axis([0 Inf -1 1]) ga.FontSize = 14; %(Optional) Text to indicate the location of Tmax txt = strcat('t = ',num2str((Tmax-1)*tInt),' s'); text((Tmax-1)*tInt,double(max(dFilt)),txt) %Making sure that the legend is accurate nme = strcat('Cell ',num2str(i)); nms{ncnt} = nme; ncnt = ncnt + 1; catch %If Labs.tif specifies regions that don't display any optical %flow, there won't be any peaks, so the Tmax assignment throws %an error. The second iteration should omit those cells, but in %the meantime the rest of the analysis will keep running. warning('No peaks in cell %d!!!',i) end end hold on end %Axis labels and legend xlabel('Time (s)') subplot(2,1,1) legend(nms,'Location','northwest') legend('boxoff') ylabel('$<\nu>$','Interpreter','Latex') subplot(2,1,2) ylabel('$\Delta <\nu>$','Interpreter','Latex') indName = strcat(fName,'_Cell_Flow');