Create a gist now

Instantly share code, notes, and snippets.

Embed
What would you like to do?
MATLAB function to autodetect cells and perform image processing to quantify changes in fluorescence intensity
function [varargout] = CellFinder(varargin)
%% CellFinder - FUNCTION TO AUTODETECT CELLS
clc; close all; clear;
% Change the current folder to the folder of this m-file.
cd(fileparts(which(mfilename('fullpath'))));
%% -- DEAL OR SET INPUT ARGS
% if exist('varargin','var') && nargin > 0
% [v] = deal(varargin{:});
% end
cellSizeMin = 60;
%% -- GET ALL THE THE MEDIA FILES LOCATED IN SPECIFIED FOLDER
regexpStr = '((\S)+(\.tiff|\.tif|\.jpg|\.png+))';
filedir = uigetdir();
allfileinfo = dir(filedir);
allfilenames = {allfileinfo.name};
filepaths = fullfile(filedir,allfilenames(~cellfun('isempty',regexp(allfilenames,regexpStr))));
ImageFiles=filepaths;
%% -- RESIZE IMAGES AND STORE IMAGE PIXEL MATRIX VALUES INTO A CELL ARRAY
% THE FINAL CELL ARRAY OF INTEREST WILL BE CALLED 'iDUBs'
% To access the pixel value matrix from a single image from this cell array,
% for example the first image of the z-stack, simply use:
numF = numel(ImageFiles);
fprintf('\n There are %s total frames \n',num2str(numF))
for nT = 1:numF
[f,map] = imread(ImageFiles{nT}); % get image data from file
%infoto = imfinfo(ImageFiles{nT}); % image meta info
I = f; % store copy of image data
% colormap will be a 512x512 matrix of class uint8 (values range from 0-255)
iIND = gray2ind(I);
iIND = imresize(iIND, [512 NaN]);
iINDs{nT} = iIND;
% colormap will be a 512x512 matrix of class double (values range from 0-1)
iDUB = im2double(I);
iDUB = imresize(iDUB, [512 NaN]);
iDUB(iDUB > 1) = 1;
iDUBs{nT} = iDUB;
imagesc(iDUB)
pause(.05)
end
%% PROMPT TO GET STIMULATION FRAME
promptTXT = {'Enter Stimulation Frame(s): '};
dlg_title = 'Input'; num_lines = 1;
presetval = {'29:31'};
dlgOut = inputdlg(promptTXT,dlg_title,num_lines,presetval);
stimFrames = str2num(dlgOut{:});
% Remove stim frames
iDUBs(stimFrames) = [];
numF = size(iDUBs,2);
fprintf('\n There are now %s total frames \n',num2str(numF))
%% -- SUM OVER Z-STACK OF IMAGES
IMGsum = iDUBs{1};
IMG = iDUBs{1};
for nT = 2:numF
IMG(:,:,nT) = iDUBs{nT};
end
IMGmu = mean(IMG,3);
%% PLOT MEAN IMAGE VS FIRST IMAGE OF STACK
close all;
fh1 = figure('Units','normalized','OuterPosition',[.1 .1 .8 .6],'Color','w','MenuBar','none');
hax1 = axes('Position',[.05 .05 .4 .9]);
hax2 = axes('Position',[.52 .05 .4 .9]);
axes(hax1);
imagesc(IMG(:,:,1))
axes(hax2);
imagesc(IMGmu)
pause(2)
%% USE MOUSE TO DRAW BOX AROUND BACKGROUND AREA
close all;
fh4=figure('Units','normalized','OuterPosition',[.1 .1 .6 .8],'Color','w','MenuBar','none');
hax4 = axes('Position',[.05 .05 .9 .9],'Color','none');
imagesc(IMGmu);
title('USE MOUSE TO DRAW BOX AROUND BACKGROUND AREA THEN CLOSE FIGURE WINDOW')
disp('DRAW BOX AROUND A BACKGROUND AREA THEN CLOSE FIGURE WINDOW')
h1 = imrect;
pos1 = round(getPosition(h1)); % [xmin ymin width height]
%% GET FRAME COORDINATES AND CREATE XY MASK
MASKTBLR = [pos1(2) (pos1(2)+pos1(4)) pos1(1) (pos1(1)+pos1(3))];
% Background
mask{1} = zeros(size(iDUB));
mask{1}(MASKTBLR(1):MASKTBLR(2), MASKTBLR(3):MASKTBLR(4)) = 1;
mask1 = mask{1};
%% CHECK THAT MASK(S) ARE CORRECT
close all;
fh5=figure('Units','normalized','OuterPosition',[.1 .1 .6 .8],'Color','w','MenuBar','none');
hax5 = axes('Position',[.05 .05 .9 .9],'Color','none');
imagesc(IMGmu.*mask{1});
pause(1)
%% -- GET MEAN OF BACKGROUND PIXELS & SUBTRACT FROM IMAGE
IMGS = IMG;
f1BACKGROUND = IMGmu .* mask1;
meanBG = mean(f1BACKGROUND(f1BACKGROUND > 0));
meanALL = mean(IMGmu(:));
IMU = IMGmu - meanBG;
IMU(IMU <= 0) = 0;
IMG = IMGS - meanBG;
IMG(IMG <= 0) = 0;
%% -- MESH SURFACE PLOT
close all;
fh5=figure('Units','normalized','OuterPosition',[.1 .1 .6 .8],'Color','w','MenuBar','none');
hax5 = axes('Position',[.05 .05 .9 .9],'Color','none');
axes(hax5);
mesh(IMU)
for vv = 1:90
view([vv 90-vv])
pause(.02)
end
pause(2)
%% check vague SNR for masks.
hist(IMU(:),80);
pause(.5)
%% SET THRESHOLD BASED ON SNR HISTOGRAM
promptTXT = {'Enter Threshold Mask Values:'};
dlg_title = 'Input'; num_lines = 1;
presetval = {num2str(0.015)};
dlgOut = inputdlg(promptTXT,dlg_title,num_lines,presetval);
threshmask = str2num(dlgOut{:});
%% -- REMOVE PIXELS BELOW THRESHOLD
IMU(IMU <= threshmask) = 0;
IMG(IMG <= threshmask) = 0;
%% -- VIEW MEAN IMAGE VS FIRST IMAGE OF STACK AFTER BG SUBTRACTION
close all;
GUIfh=figure('Units','normalized','OuterPosition',[.1 .1 .8 .6],'Color','w','MenuBar','none');
hax1 = axes('Position',[.05 .05 .4 .9]);
hax2 = axes('Position',[.52 .05 .4 .9]);
axes(hax1);
imagesc(IMG(:,:,1))
axes(hax2);
imagesc(IMU)
pause(2)
%% --- PERFORM CONVOLUTION WITH GAUSSIAN FILTER ---
doConvnFilter = 1;
if doConvnFilter
% MaskMaker([PEAK HEIGHT] [MASK SIZE] [SLOPE SD] [RESOLUTION])
Mask = MasksMaker(2.5, 5, .17, .1);
fIMU = convn(IMU,Mask,'same');
fIMG = IMG(:,:,1);
for nn = 1:size(IMG,3)
fIMG(:,:,nn) = convn(IMG(:,:,nn),Mask,'same');
end
else
fIMU = IMU;
fIMG = IMG;
end
%% PLOT MEAN IMAGE VS FIRST IMAGE OF STACK
close all;
fh1 = figure('Units','normalized','OuterPosition',[.1 .1 .8 .6],'Color','w','MenuBar','none');
hax1 = axes('Position',[.05 .05 .4 .9]);
hax2 = axes('Position',[.52 .05 .4 .9]);
axes(hax1);
imagesc(fIMG(:,:,1))
axes(hax2);
imagesc(fIMU)
pause(2)
%% --- PERFORM CELL AUTOTRACE ---
% Pad the image with zeros so no cell "edge" touches the image outer boarder
fIMU(1:5,:) = 0;
fIMU(end-5:end,:) = 0;
fIMU(:,1:5) = 0;
fIMU(:,end-5:end) = 0;
dim = size(fIMU);
col = round(dim(2)/2)-90;
row = min(find(fIMU(:,col)));
fIMG_filled = imfill(fIMU,'holes');
fbIMG = bwboundaries(fIMG_filled);
for k=1:numel(fbIMG)
fbkIMG(k) = numel(fbIMG{k});
end
% cellSizeMin = 60;
tIMG = fbkIMG < cellSizeMin;
fbIMG(tIMG) = [];
%% --- PLOT AUTOTRACED EDGES ---
close all
fh1=figure('Units','normalized','OuterPosition',[.1 .1 .8 .6],'Color','w','MenuBar','none');
hax1 = axes('Position',[.05 .05 .45 .9],'Color','none');
hax2 = axes('Position',[.52 .05 .45 .9],'Color','none');
axes(hax1)
imagesc(fIMU)
colormap(hax1,hot)
title('Convolution Processed Image Mean')
hold on
for k=1:numel(fbIMG)
plot(hax1,fbIMG{k}(:,2),fbIMG{k}(:,1),'g','LineWidth',2);
end
hold on
pause(.5)
%axes(GUIfh.Children(1).Children(2))
axes(hax2)
imagesc(fIMG(:,:,1))
colormap(hax2,hot)
title('Convolution Processed First Image of Stack')
hold on
for k=1:numel(fbIMG)
plot(hax2,fbIMG{k}(:,2),fbIMG{k}(:,1),'w','LineWidth',2);
end
hold on
pause(2)
%% CREATE BINARY IMAGE AND GET LABELS FOR EACH CELL
% Create black and white (BW) binary image
fIMG_BW = im2bw(fIMG_filled, .5);
% Returns labels for connected components of a binary image
[L, NUM] = bwlabeln(fIMG_BW);
fprintf('\n There are %s discrete cells \n',num2str(NUM))
axes(hax1)
imagesc(fIMG_BW);
title('BINARY IMAGE')
axes(hax2)
imagesc(L);
title('LABELS (EACH CELL HAS A UNIQUE COLOR USED AS LABEL)')
pause(2)
%% -- CALCULATE STATISTICS FROM REGIONPROPS AT LOCATIONS (L) ON IMAGE STACK
STATS = {regionprops(L,fIMG(:,:,1), 'MeanIntensity','Perimeter','PixelValues')};
for nn = 1:size(IMG,3)
STATS{nn} = regionprops(L,fIMG(:,:,nn), 'MeanIntensity','Perimeter','PixelValues');
% STATS = regionprops(L,fIMG(:,:,1), 'all')
end
% STATS{1}(1).MeanIntensity
% STATS{1}(1).Perimeter
% STATS{1}(1).PixelValues
PerimSizes = [STATS{1}(:).Perimeter];
TooSmall = PerimSizes < cellSizeMin/4;
for nn = 1:size(STATS,2)
STATS{nn}(TooSmall) = [];
end
numIMGS = size(STATS,2);
numCells = sum(~TooSmall);
fprintf('\n There are %s cells that passed size threshold \n', num2str(numCells) )
%% GET DESIRED STATS PER STACK IMAGE
% Each column is a cell
% Each row is a frame
for mm = 1:numIMGS
for nn = 1:numCells
PERIMETER(mm,nn) = STATS{mm}(nn).Perimeter;
end
end
for mm = 1:numIMGS
for nn = 1:numCells
INTENSITY(mm,nn) = STATS{mm}(nn).MeanIntensity;
end
end
% for mm = 1:numIMGS
% for nn = 1:numCells
% PIXELS(mm,nn) = STATS{mm}(nn).PixelList;
% end
% end
%% CALCULATE DELTA F OVER F
BL = repmat( mean(INTENSITY(1:10,:)) , size(INTENSITY,1) , 1);
IMGdF = (INTENSITY - BL) ./ BL;
%% PLOT DELTA F OVER F VALUES
fh1=figure('Units','normalized','OuterPosition',[.1 .1 .8 .7],'Color','w','MenuBar','none');
hax1 = axes('Position',[.05 .05 .9 .9],'Color','none');
axes(hax1)
ph1 = plot(IMGdF,'LineWidth',3);
%%
varargout = {IMGdF};
end
%====================================================================
function Mask = MasksMaker(varargin)
%%
% Mask = MaskMaker(2.5, 11, .18, .1)
% MaskMaker([HIGHT OF PEAK] [SIZE OF MASK] [STDEV OF SLOPE] [RESOLUTION])
GNpk = 2.5; % HIGHT OF PEAK
GNnum = 11; % SIZE OF MASK
GNsd = 0.18; % STDEV OF SLOPE
GNres = 0.1; % RESOLUTION
%% -- DEAL ARGS
if nargin < 1
GNpk = 2.5; % HIGHT OF PEAK
GNnum = 11; % SIZE OF MASK
GNsd = 0.18; % STDEV OF SLOPE
GNres = 0.1; % RESOLUTION
elseif nargin == 1
v1 = varargin{1};
GNpk = v1; % HIGHT OF PEAK
GNnum = 11; % SIZE OF MASK
GNsd = 0.18; % STDEV OF SLOPE
GNres = 0.1; % RESOLUTION
elseif nargin == 2
[v1, v2] = deal(varargin{:});
GNpk = v1; % HIGHT OF PEAK
GNnum = v2; % SIZE OF MASK
GNsd = 0.18; % STDEV OF SLOPE
GNres = 0.1; % RESOLUTION
elseif nargin == 3
[v1, v2, v3] = deal(varargin{:});
GNpk = v1; % HIGHT OF PEAK
GNnum = v2; % SIZE OF MASK
GNsd = v3; % STDEV OF SLOPE
GNres = 0.1; % RESOLUTION
elseif nargin == 4
[v1, v2, v3, v4] = deal(varargin{:});
GNpk = v1; % HIGHT OF PEAK
GNnum = v2; % SIZE OF MASK
GNsd = v3; % STDEV OF SLOPE
GNres = v4; % RESOLUTION
else
warning('Too many inputs')
end
%% -- MASK SETUP
GNx0 = 0; % x-axis peak locations
GNy0 = 0; % y-axis peak locations
GNspr = ((GNnum-1)*GNres)/2;
a = .5/GNsd^2;
c = .5/GNsd^2;
[X, Y] = meshgrid((-GNspr):(GNres):(GNspr), (-GNspr):(GNres):(GNspr));
Z = GNpk*exp( - (a*(X-GNx0).^2 + c*(Y-GNy0).^2)) ;
Mask=Z;
doMASKfig = 1;
if doMASKfig
disp(Mask);
%----------------%
% FIGURE: 3D Gaussian Distribution
fh5 = figure(5); set(fh5,'OuterPosition',[200 200 500 300]);
%----------------%
subplot('Position',[.05 .55 .30 .40]);
ph5 = imagesc(Mask);
axis equal;
%set(gca,'XTick',[],'YTick',[])
subplot('Position',[.04 .08 .32 .42]);
ph6 = surf(X,Y,Z);
axis equal; shading interp; view(90,90);
subplot('Position',[.45 .05 .50 .90]);
ph7 = surf(X,Y,Z);
axis vis3d; shading interp;
view(-45,30);
xlabel('x-axis');ylabel('y-axis');zlabel('z-axis')
%-------------------------------%
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment