Create a gist now

Instantly share code, notes, and snippets.

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