LoGOF+RelativeAngle_Downstream

All of these functions can be used by a researcher to perform an analysis just using the Filament Displacement tool as well as use this information to further understand the effects of filament orientation on the active forcing environment.

LKxOptFlow_allFrames.m

@@ -0,0 +1,148 @@
+%%  Calculate Optical Flow - all frames at once
+%
+%   Copyright (C) 2020 Leonard Campanello - All Rights Reserved
+%   You may use, distribute, and modify this code under the terms of a 
+%   CC BY 4.0 License. 
+%
+%   If you use, distribute, or modify this code, please cite
+%   Lee*, Campanello*, et al. MBoC 2020
+%   https://doi.org/10.1091/mbc.E19-11-0614
+%
+%   Written by Leonard Campanello
+%   Available on https://github.com/ljcamp1624/OpticalFlowAnalysis
+%
+%   Direct all questions and correspondence to Leonard Campanello and 
+%   Wolfgang Losert at the University of Maryland: 
+%   ljcamp (at) umd (dot) edu
+%   wlosert (at) umd (dot) edu
+%
+%%  Begin function
+%
+%   This function is called by CalculateOpticalFlow.m
+%
+function [vx, vy, rel] = LKxOptFlow_allFrames(images, xySig1, tSig, wSig)
+%%  Convert images
+images = double(images);
+
+%%  Create filters for calculating gradients
+
+%   Process parameters
+xyRange1 = -ceil(3*xySig1):ceil(3*xySig1);
+xySig2 = xySig1/4;
+xyRange2 = -ceil(3*xySig2):ceil(3*xySig2);
+tRange = -ceil(3*tSig):ceil(3*tSig);
+
+%   Build x-gradient filter kernels (along first dimension)
+x = xyRange1;
+y = xyRange2;
+fx1 = exp(-x.*x/2/xySig1/xySig1)/sqrt(2*pi)/xySig1;
+fy1 = exp(-y.*y/2/xySig2/xySig2)/sqrt(2*pi)/xySig2;
+fx1=ones(size(x));
+fy1=ones(size(y));
+gx1 = x/xySig1/xySig1;
+gy1 = 1;
+% xFil1 = (fx1.*gx1)';
+xFil1=gx1';
+% yFil1 = (fy1.*gy1);
+yFil1=gy1;
+%   Build y-gradient filter kernels (along first dimension)
+x = xyRange2;
+y = xyRange1;
+fx2 = exp(-x.*x/2/xySig2/xySig2)/sqrt(2*pi)/xySig2;
+fx2=fx1;
+fy2=fy1;
+fy2 = exp(-y.*y/2/xySig1/xySig1)/sqrt(2*pi)/xySig1;
+gx2 = 1;
+gy2 = y/xySig1/xySig1;
+% xFil2 = (fx2.*gx2)';
+xFil2=gx2';
+yFil2=gy2;
+% yFil2 = (fy2.*gy2);
+
+%   Build t-gradient filter kernels (along first dimension)
+x = xyRange1;
+y = xyRange1;
+t = tRange;
+% fx3 = exp(-x.*x/2/xySig1/xySig1)/sqrt(2*pi)/xySig1;
+% fy3 = exp(-y.*y/2/xySig1/xySig1)/sqrt(2*pi)/xySig1;
+fx3=fx1;
+fy3=fy1;
+% ft3 = exp(-t.*t/2/tSig/tSig)/sqrt(2*pi)/tSig;
+ft3=fx1;
+gx3 = 1;
+gy3 = 1;
+gt3 = t/tSig/tSig;
+% xFil3 = (fx3.*gx3)';
+xFil3=gx3';
+yFil3=gy3;
+% yFil3 = (fy3.*gy3);
+% tFil3 = permute(ft3.*gt3, [3, 1, 2]);
+tFil3 = permute(gt3, [3, 1, 2]);
+
+% %   x-gradient filter (along 1st dimension of the image)
+% [iy, ix] = meshgrid(xyRange2, xyRange1);
+% fx = exp(-ix.*ix/2/xySig1/xySig1)/sqrt(2*pi)/xySig1;
+% fy = exp(-iy.*iy/2/xySig2/xySig2)/sqrt(2*pi)/xySig2;
+% gx = ix.*fx.*fy/xySig1/xySig1;
+% 
+% %   y-gradient filter (along 2nd dimension of the image)
+% [iy, ix] = meshgrid(xyRange1, xyRange2);
+% fx = exp(-ix.*ix/2/xySig2/xySig2)/sqrt(2*pi)/xySig2;
+% fy = exp(-iy.*iy/2/xySig1/xySig1)/sqrt(2*pi)/xySig1;
+% gy = iy.*fx.*fy/xySig1/xySig1;
+% 
+% %   t-gradient filter (along 3rd dimension of the image)
+% [iy, ix, it] = meshgrid(-3:3, -3:3, -3*tSig:3*tSig);
+% fx = exp(-ix.*ix/2)/sqrt(2*pi);
+% fy = exp(-iy.*iy/2)/sqrt(2*pi);
+% ft = exp(-it.*it/2/tSig/tSig)/sqrt(2*pi)/tSig;
+% gt = it.*fx.*fy.*ft/tSig/tSig;
+
+%% Calculate gradients
+% dxI = imfilter(images, gx, 'replicate');
+% dyI = imfilter(images, gy, 'replicate');
+% dtI = imfilter(images, gt, 'replicate');
+dxI = imfilter(imfilter(images, xFil1, 'replicate'), yFil1, 'replicate');
+dyI = imfilter(imfilter(images, xFil2, 'replicate'), yFil2, 'replicate');
+dtI = imfilter(imfilter(imfilter(images, xFil3, 'replicate'), yFil3, 'replicate'), tFil3, 'replicate');
+% tFilTest(:,:,1)=-1;
+% tFilTest(:,:,2)=1;
+% dtI = imfilter(imfilter(imfilter(images, xFil3, 'replicate'), yFil3, 'replicate'), tFilTest, 'replicate');
+
+clear images
+
+%% Calculate structure tensor inputs
+% wdx2 = imgaussfilt(dxI.^2, wSig, 'padding', 'replicate');
+% wdxy = imgaussfilt(dxI.*dyI, wSig, 'padding', 'replicate');
+% wdy2 = imgaussfilt(dyI.^2, wSig, 'padding', 'replicate');
+% wdtx = imgaussfilt(dxI.*dtI, wSig, 'padding', 'replicate');
+% wdty = imgaussfilt(dyI.*dtI, wSig, 'padding', 'replicate');
+wRange = -ceil(3*wSig):ceil(3*wSig);
+x = wRange;
+y = wRange;
+gx = exp(-x.*x/2/wSig/wSig)/sqrt(2*pi)/wSig;
+gy = exp(-y.*y/2/wSig/wSig)/sqrt(2*pi)/wSig;
+xFil4 = (gx)';
+yFil4 = (gy);
+
+wdx2 = imfilter(imfilter(dxI.*dxI, xFil4, 'replicate'), yFil4, 'replicate');
+wdxy = imfilter(imfilter(dxI.*dyI, xFil4, 'replicate'), yFil4, 'replicate');
+wdy2 = imfilter(imfilter(dyI.*dyI, xFil4, 'replicate'), yFil4, 'replicate');
+wdtx = imfilter(imfilter(dxI.*dtI, xFil4, 'replicate'), yFil4, 'replicate');
+wdty = imfilter(imfilter(dyI.*dtI, xFil4, 'replicate'), yFil4, 'replicate');
+clear dxI dyI dtI
+
+%% Calculate Optical Flow
+determinant = (wdx2.*wdy2) - (wdxy.^2);
+vx = ((determinant + eps).^-1).*((wdy2.*-wdtx)+(-wdxy.*-wdty));
+vy = ((determinant + eps).^-1).*((-wdxy.*-wdtx)+(wdx2.*-wdty));
+clear wdtx wdty
+
+trace = wdx2 + wdy2;
+clear wdx2 wdy2 wdxy
+
+e1 = (trace + sqrt(trace.^2 - 4*determinant))/2;
+e2 = (trace - sqrt(trace.^2 - 4*determinant))/2;
+rel = real(min(e1, e2));
+
+end

LoGOFTool_fixedPad.m

@@ -0,0 +1,75 @@
+% This function allows a researcher to analyze the dynamics of filaments in
+% a dense MT networks
+%
+%%%% USAGE: function [logOF,LogFilt,dsOG] = LoGOFTool_fixedPad(imSeq,paramsLog,paramsOF,padSize)
+%
+%%%% INPUTS
+% for LoG and OF params consults Log_fixedPad_imseq.m and
+% LKxOptFlow_allFrames.m
+%
+%%%% OUTPUTS
+% logOF - the output for quantifying lateral motion
+% LogFilt - the output of LoG only (to be used for further downstream
+% analysis of relative angles)
+% dsOG - the original data adjusted such that overlays of results agree on
+% a pixel-pixel basis
+% DEPENDENCIES
+% Log_fixedPad_imseq.m - generate the anisotropic LoG output
+% LKxOptFlow_allFrames.m - Lucas-Kanade optical flow as developed by Lenny
+%                           Campanello in the Losert Lab (see Lee et. al)
+
+% Original function by Nick Mennona
+% Please direct all correspondence to:
+% nmennona@umd.edu or wlosert@umd.edu
+%%
+function [logOF,LogFilt,dsOG] = LoGOFTool_fixedPad(imSeq,paramsLog,paramsOF,padSize)
+% find motion perpendicular to those pixels identified as filaments in an
+% image
+
+
+%INPUT:
+% imSeq-image sequence for processing/analysis
+% paramsLog:
+%  
+        % logparams.filSig
+        % logparams.numSig
+        % logparams.numAngs
+        % logparams.filterThreshold
+% paramsOF:
+        % xySig
+        % tSig
+        % wSig 
+
+
+    numFrames = size(imSeq,3);
+
+
+    [LogFilt, dsOG] = Log_fixedPad_imseq(imSeq, paramsLog,padSize);
+
+    gaussSig=1.5;
+    filtVid = gaussianTopHat(dsOG,gaussSig);
+
+    xySig = paramsOF.xySig;
+    tSig = paramsOF.tSig;
+    wSig = paramsOF.wSig;
+    [vxMat, vyMat, relMat] = LKxOptFlow_allFrames(filtVid, xySig, tSig, wSig);
+    of_Magnitude = sqrt(vxMat.^2+vyMat.^2);
+    of_Orientation = atan2(vyMat,vxMat);
+
+
+    for tt=1:numFrames
+        tempOr = of_Orientation(:,:,tt);
+        tempOr(tempOr==0)=nan;
+        tempOr(tempOr<0)=pi+tempOr(tempOr<0);
+        tempMag = of_Magnitude(:,:,tt);
+        tempMag(tempMag==0)=nan;
+        tempLOG = LogFilt(:,:,tt);
+
+        % now we have all of the necessary matrices
+        comFrame = abs(tempOr-tempLOG);
+        logOF(:,:,tt)=tempMag.*sin(comFrame);
+        if mod(tt,10)==0
+            disp(['Reading Frame ' num2str(tt)])
+        end
+    end
+end

createManyBulkMasks.m

@@ -0,0 +1,34 @@
+function regions = createManyBulkMasks(bulkMask,originalBoundary,numBulkMasks)
+% goal of this function is to take a bulkMask that has already been
+% generated from findCellandBulkBoundary.m and create smaller versions of
+% the same bulkMask by erosion to deduce regional specific differences in
+% dynamics 
+
+% as a geometric representation, the original bulkMask and all others until
+% the numBulkMasks (final) mask are all ring like
+
+% INPUT:
+% bulkMask: original bulkMask to be eroded
+% numBulkMasks: number of regions to be generated
+
+% OUTPUT:
+% regions: struct which contains all of the masks of interest
+randomSEval=20;
+regions(1).bulkMask = originalBoundary;
+for i=1:numBulkMasks-1
+    if i==1
+        tempNanMask = im2double(bulkMask);
+        tempNanMask(tempNanMask==0)=NaN;
+        [newbulkMask,bndryMask]=getBndryBulk(bulkMask,randomSEval,tempNanMask);
+        regions(i+1).bulkMask = bndryMask;
+        tempBulk = newbulkMask;
+    else
+        tempNanMask = im2double(tempBulk);
+        tempNanMask(tempNanMask==0)=NaN;
+        [newbulkMask,bndryMask]=getBndryBulk(tempBulk,randomSEval,tempNanMask);
+        regions(i+1).bulkMask = bndryMask;
+        tempBulk = newbulkMask;
+    end
+end
+regions(numBulkMasks+1).bulkMask=newbulkMask;
+end

findCellandBulkBoundary.m

@@ -0,0 +1,132 @@
+% This function is not necessary for the use of LoGOF generally, but can be used 
+% to help researchers segment specific cells by hand in a FOV
+
+function [bulkReg, boundaryReg,gaussianMask] = findCellandBulkBoundary(rawImSeq)
+% INPUT: rawImSeq-take in the raw image sequence...output from LOGOFToolv3
+% this function takes the first frame from that raw image and then tries to
+% get a cellMask for analysis
+% OUTPUT:
+% bulkReg--bulk region that is (hardcoded at this point) 85% of the cell
+% mask area
+% boundaryReg--the rest to be composed of as the boundary
+
+% PLEASE FOLLOW INSTRUCTIONS ON HOW TO DRAW THE SEGMENTING POLYGON
+    close all
+    firstFrame = rawImSeq(:,:,1);
+    figure;
+    imagesc(firstFrame);truesize
+    disp('Look at the image!')
+    if input('Do you need to draw a polygon? 1 for Yes, 0 for No ')
+        if input('Does it take up the full Screen? 1 for Yes 0 for No ')
+            lineROI = drawline;
+            numX = floor(lineROI.Position(1));% draw line left to right
+            numXf = floor(lineROI.Position(2));
+            tempGetRid = NaN(size(firstFrame));
+            tempGetRid(1:end,numX:numXf)=1;
+            segmentedCellRegion = tempGetRid.*firstFrame;
+            imagesc(segmentedCellRegion)
+%             disp('Does this look good?')
+%             pause(10)
+            %need to start manually changing parameters HERE
+            [A,~,~]=findCellBoundary(segmentedCellRegion,50,20000);
+            figure(2); imagesc(A);truesize
+            X = input('Does this look good?');
+            while X~=1
+                newThresh = input('Enter new threshold (default was 50) ');
+                newArea = input('Enter new obj area (default was 20000) ');
+                [A,~,~]=findCellBoundary(segmentedCellRegion,newThresh,newArea);
+                imagesc(A);truesize
+                X = input('Does this look good?');
+            end
+            imagesc(A)
+%             pause(10)
+        else
+            cellROI = drawpolygon;
+            tempGetRid = NaN(size(firstFrame));
+            % these values are based off of choosing top left, top right, bottom right,
+            % bottom left in that sequence for the ROI!
+            numX=floor(cellROI.Position(5)); numXf = floor(cellROI.Position(8));
+            numY = floor(cellROI.Position(1));numYf = floor(cellROI.Position(2));
+            close all
+            tempGetRid(numX:numXf,numY:numYf)=1;
+            segmentedCellRegion = tempGetRid.*firstFrame;
+            [A,~,~]=findCellBoundary(segmentedCellRegion,50,20000);
+            X = input('Does this look good?');
+            while X~=1
+                newThresh = input('Enter new threshold (default was 50) ');
+                newArea = input('Enter new obj area (default was 20000) ');
+                [A,~,~]=findCellBoundary(segmentedCellRegion,newThresh,newArea);
+                X = input('Does this look good?');
+            end
+            imagesc(A)
+        end
+    else
+        segmentedCellRegion = firstFrame;
+        [bigTemp,~,~]=findCellBoundary(segmentedCellRegion,50,20000);
+        bigTemp(1,:)=1; bigTemp=imfill(bigTemp,'holes');bigTemp(1,:)=0;
+        bigTemp(end,:)=1;bigTemp=imfill(bigTemp,'holes');bigTemp(end,:)=0;
+        A=bigTemp;
+    end
+
+    gaussianMask = imgaussfilt(im2double(A),10)>.15;
+    disp('CellMask generated and smoothed! Check it out!')
+    imagesc(gaussianMask)
+    maskGood=0;
+    while ~maskGood
+        maskGood=input('Is this okay? 1 for Yes, 0 for No ');
+        if maskGood==0
+            filterThresh = input('Enter new threshhold ');
+            sizeThresh = input('Enter new Size Threshold (orig = 20000) '); 
+            [A,~,~]=findCellBoundary(segmentedCellRegion,filterThresh,sizeThresh);
+            gaussianMask = imgaussfilt(im2double(A),10)>.15;
+        end
+        imagesc(gaussianMask)
+        figure
+        imagesc(firstFrame)
+    end
+
+    percentage = .25;
+    [bulkReg,boundaryReg] = findPercentageBoundary(gaussianMask,percentage);
+end
+
+
+% 
+% test = params(2).logOF;
+% figure;
+% imagesc(test(:,:,3))
+% cellROI = drawpolygon;
+% %%
+% tempGetRid = NaN(size(test(:,:,1)));
+% % these values are based off of choosing top left, top right, bottom right,
+% % bottom left in that sequence for the ROI!
+% numX=cellROI.Position(5); numXf = cellROI.Position(8);
+% numY = cellROI.Position(1);numYf = cellROI.Position(2);
+% tempGetRid(cellROI.Position(5):cellROI.Position(8),cellROI.Position(1):cellROI.Position(2))=1;
+% imagesc(tempGetRid.*testGetRid)
+% [A,B,C]=findCellBoundary(tempGetRid.*testGetRid);imagesc(A)
+% %%
+% % for i=1:size(dsOG,3)
+% %     [bigTemp,~,~] = findCellBoundary(dsOG(:,:,i));
+% %     % 4/26/22--interesting workaround I found
+% % bigTemp(1,:)=1; bigTemp=imfill(bigTemp,'holes');bigTemp(1,:)=0;
+% % bigTemp(end,:)=1;bigTemp=imfill(bigTemp,'holes');bigTemp(end,:)=0;
+% %     labelledimages{i,1} = bigTemp;
+% % end
+% % disp('Found Boundary and Bulk regions')
+% %% GAUSSIAN SMOOTH THE BINARIZED IMAGE TO GET A NICE SMOOTH MASK
+% bigTemp=imgaussfilt(im2double(A),10)>.15;
+% imagesc(bigTemp)
+% %
+% percentage = .25;
+% [bulkReg,boundaryReg] = findPercentageBoundary(bigTemp,percentage);
+% imshowpair(bulkReg,boundaryReg)
+
+
+
+
+
+
+
+
+
+