Commit 8283ebee authored by Antonio Politi's avatar Antonio Politi
Browse files

Upload of Julius Hossein Segmentation code

To do is to remove tiffread dependency and use bioformat for matlab
parent e768fa31
function [ threshold ] = Otsu_2D(image)
[xSize, ySize] = size(image);
maxVal = max(max(image));
histN = round(maxVal)+1;
histogram = zeros(1,histN);
threshold=0;
w0=0.0;
w1=0.0;
m0=0.0;
m1=0.0;
N=0;
p=0;
sum=0.0;
mean=0.0;
var_bet_class=0.0;
var_max=0.0;
mu_k=0.0;
for i = 1:xSize
for j = 1:ySize
hIndex = round(image(i,j))+1;
histogram(hIndex) = histogram(hIndex)+1;
end
end
% histogram
% pause;
for i=1:histN
N=N+histogram(i);
sum=sum+histogram(i)*i;
end
mean=sum/N;
for i=1:histN
p=histogram(i)/N;
%cumulative for class 1 and class 2
w0=w0+p;
w1=1-w0;
%mean for class 1 and class 2
mu_k=mu_k+i*p;
if(w0==0)
m0=0;
else
m0=(mu_k/w0);
end
if(1-w0==0)
m1=0;
else
m1=(mean-mu_k)/(1-w0);
end
var_bet_class = w0*(m0-mean)*(m0-mean)+w1*(m1-mean)*(m1-mean);
if var_bet_class>var_max
var_max=var_bet_class;
threshold=i;
end
end
end
function [ threshold ] = Otsu_2D_Img(image, masked)
[xSize, ySize] = size(image);
maxVal = max(max(image));
histN = round(maxVal)+1;
histogram = zeros(1,histN);
threshold=0;
w0=0.0;
w1=0.0;
m0=0.0;
m1=0.0;
N=0;
p=0;
sum=0.0;
mean=0.0;
var_bet_class=0.0;
var_max=0.0;
mu_k=0.0;
for i = 1:xSize
for j = 1:ySize
hIndex = round(image(i,j))+1;
histogram(hIndex) = histogram(hIndex)+1;
end
end
if masked == 1
histogram(1) = 0;
end
for i=1:histN
N=N+histogram(i);
sum=sum+histogram(i)*i;
end
mean=sum/N;
for i=1:histN
p=histogram(i)/N;
%cumulative for class 1 and class 2
w0=w0+p;
w1=1-w0;
%mean for class 1 and class 2
mu_k=mu_k+i*p;
if(w0==0)
m0=0;
else
m0=(mu_k/w0);
end
if(1-w0==0)
m1=0;
else
m1=(mean-mu_k)/(1-w0);
end
var_bet_class = w0*(m0-mean)*(m0-mean)+w1*(m1-mean)*(m1-mean);
if var_bet_class>var_max
var_max=var_bet_class;
threshold=i;
end
end
end
function [threshold] = Otsu_3D_Hist(histogram,cut_thresh, masked)
%This function receives histogram as input
%It computes a threshold without considering the values in the histogram
%greater than cut_threshold
histN = cut_thresh;
threshold=1;
w0=0.0;
w1=0.0;
m0=0.0;
m1=0.0;
N=0;
p=0;
sum=0.0;
mean=0.0;
var_bet_class=0.0;
var_max=0.0;
mu_k=0.0;
if masked ==1
histogram(1) = 0;
end
for i=1:histN
N=N+histogram(i);
sum=sum+histogram(i)*i;
end
mean=sum/N;
for i=1:histN
p=histogram(i)/N;
%cumulative for class 1 and class 2
w0=w0+p;
w1=1-w0;
%mean for class 1 and class 2
mu_k=mu_k+i*p;
if(w0==0)
m0=0;
else
m0=(mu_k/w0);
end
if(1-w0==0)
m1=0;
else
m1=(mean-mu_k)/(1-w0);
end
var_bet_class = w0*(m0-mean)*(m0-mean)+w1*(m1-mean)*(m1-mean);
if var_bet_class>var_max
var_max=var_bet_class;
threshold=i;
end
end
end
function [threshold, histogram] = Otsu_3D_Img(image, masked)
[xSize, ySize, zSize] = size(image);
maxVal = max(max(max(image)));
histN = round(maxVal)+1;
histogram = zeros(1,histN);
threshold=1;
w0=0.0;
w1=0.0;
m0=0.0;
m1=0.0;
N=0;
p=0;
sum=0.0;
mean=0.0;
var_bet_class=0.0;
var_max=0.0;
mu_k=0.0;
for i = 1:xSize
for j = 1:ySize
for k = 1:zSize
hIndex = round(image(i,j,k))+1;
histogram(hIndex) = histogram(hIndex)+1;
end
end
end
if masked == 1
histogram(1) = 0;
end
for i=1:histN
N=N+histogram(i);
sum=sum+histogram(i)*i;
end
mean=sum/N;
for i=1:histN
p=histogram(i)/N;
%cumulative for class 1 and class 2
w0=w0+p;
w1=1-w0;
%mean for class 1 and class 2
mu_k=mu_k+i*p;
if(w0==0)
m0=0;
else
m0=(mu_k/w0);
end
if(1-w0==0)
m1=0;
else
m1=(mean-mu_k)/(1-w0);
end
var_bet_class = w0*(m0-mean)*(m0-mean)+w1*(m1-mean)*(m1-mean);
if var_bet_class>var_max
var_max=var_bet_class;
threshold=i;
end
end
end
% Version data: 2015-11-11
% Author: Julius Hossain, EMBL, Heidelberg, Germany: julius.hossain@embl.de
% modified 21.04.2016 Antonio Politi introduce channel parameter
%This scripts search for folders containing high resoluation images and
%call the function that segments chromosomes and save output in a mat file
clc; clear all; close all;
rootDir = 'Z:\AntonioP_elltier1\CelllinesImaging\160413\mEGFP-Nup107\';
fn = dir([rootDir 'HR5.lsm']);
channel = 3;
for i=1:length(fn)
try
fn_SegmentChromosome_ST(rootDir, fn(i).name, channel);
catch
disp(['Could not process the cell: ' fn(i).name]);
end
end
% rootDir = 'U:\Antonio_151214\toSegment\Ki67_1336_meta\';
% fn = dir([rootDir '*.lsm']);
%
% for i=1:length(fn)
% try
% fn_SegmentChromosome_ST(rootDir, fn(i).name);
% catch
% disp(['Could not process the cell: ' fn(i).name]);
% end
% end
\ No newline at end of file
function D=bwdistsc(bw,aspect)
% D=BWDISTSC(BW,ASPECT)
% BWDISTSC computes Euclidean distance transform of the binary 3D image
% BW. For each pixel in BW, the distance transform assignes a number
% that is the distance between that pixel and the nearest nonzero pixel
% of BW. BW may be a single 2D image, 3D array or a cell array of
% 2D slices. ASPECT is 3-component vector defining aspect ratio in
% the dataset BW. If ASPECT is not specified, isotropic aspect
% ratio [1 1 1] is assumed.
%
% BWDISTSC uses fast optimized scanning algorithm and cell-arrays to
% represent internal data, so that it is less demanding to physical
% memory. In many cases BWDISTSC is actually faster than MATLAB's
% optimized kd-tree algorithm used for Euclidean distance
% transform in 3D.
%
% BWDISTSC tries to use MATLAB bwdist for 2D scans if possible, which
% is significantly faster. Otherwise BWDISTSC uses internal algorithm
% to perform 2D scans.
%
% Yuriy Mishchenko JFRC HHMI Chklovskii Lab JUL 2007
% This code is free for use or modifications, just please give credit
% where appropriate. And if you modify code or fix bugs, please drop
% me a message at gmyuriy@hotmail.com.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Scan algorithms below use following Lema: %
% LEMA: let F(X,z) be lower envelope of a family of parabola: %
% F(X,z)=min_{i} [G_i(X)+(z-k_i)^2]; %
% and let H_k(X,z)=A(X)+(z-k)^2 be a parabola. %
% Then for H_k(X,z)==F(X,z) at each X there exist at most %
% two solutions k1<k2 such that H_k12(X,z)=F(X,z), and %
% H_k(X,z)<F(X,z) is restricted to at most k1<k2. %
% Here X is any-dimensional coordinate. %
% %
% Thus, simply scan away from any z such that H_k(X,z)<F(X,z) %
% in either direction as long as H_k(X,z)<F(X,z) and update %
% F(X,z). Note that need to properly choose starting point; %
% starting point is any z such that H_k(X,z)<F(X,z); z==k is %
% usually, but not always the starting point!!! %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parse inputs
if(nargin<2 || isempty(aspect)) aspect=[1 1 1]; end
% determine geometry of data
if(iscell(bw)) shape=[size(bw{1}),length(bw)]; else shape=size(bw); end
% fix to handle 1D & 2D images
if(length(shape)<3) shape(length(shape)+1:3)=1; end
if(length(aspect)<3) aspect(length(aspect)+1:3)=1; end
% allocate space
D=cell(1,shape(3)); for k=1:shape(3) D{k}=zeros(shape(1:2)); end
%%%%%%%%%%%%% scan along XY %%%%%%%%%%%%%%%%
for k=1:shape(3)
if(iscell(bw)) bwXY=bw{k}; else bwXY=bw(:,:,k); end
% initialize arrays
DXY=zeros(shape(1:2));
D1=zeros(shape(1:2));
DK=zeros(shape(1:2));
% if can, use 2D bwdist from image processing toolbox
if(exist('bwdist') && aspect(1)==aspect(2))
D1=aspect(1)^2*bwdist(bwXY).^2;
else % if not, use full XY-scan
%%%%%%%%%%%%%%% X-SCAN %%%%%%%%%%%%%%%
% reference nearest bwXY "on"-pixel in x direction downward:
% scan bottow-up, copy x-reference from previous row unless
% there is bwXY "on"-pixel in that point in current row
xlower=repmat(Inf,shape(1:2));
xlower(1,find(bwXY(1,:)))=1; % fill in first row
for i=2:shape(1)
xlower(i,:)=xlower(i-1,:); % copy previous row
xlower(i,find(bwXY(i,:)))=i;% unless there is pixel
end
% reference nearest bwXY "on"-pixel in x direction upward:
xupper=repmat(Inf,shape(1:2));
xupper(end,find(bwXY(end,:)))=shape(1);
for i=shape(1)-1:-1:1
xupper(i,:)=xupper(i+1,:);
xupper(i,find(bwXY(i,:)))=i;
end
% find points for which distance needs to be updated
idx=find(~bwXY); [x,y]=ind2sub(shape(1:2),idx);
% set distance as the shortest to upward or to downward
DXY(idx)=aspect(1)^2*min((x-xlower(idx)).^2,(x-xupper(idx)).^2);
%%%%%%%%%%%%%%% Y-SCAN %%%%%%%%%%%%%%%
% this will be the envelop
% envelop is initialized at Inf to ensure single scan direction,
% otherwise may end up in infinite loop when trying to find
% starting point
D1=repmat(Inf,shape(1:2));
% these will be the references to parabolas defining the envelop
DK=repmat(Inf,shape(1:2));
% starting points
i0=zeros(shape(1),1);
% convenience x-coords array
x=(1:shape(1))';
for i=1:shape(2)
% need to select starting point for each X:
% * starting point should be below current envelop
% * i0==i is not necessarily a starting point
% * there is at most one starting point
% * there may be no starting point
% i0 is the starting points for each X: i0(X) is the first
% y-index such that parabola from line i is below the envelop
% first guess is the current y-line
i0(:)=i;
% some auxiliary datasets
d0=DXY(:,i);
% L0 indicates for which X starting point had been fixed
L0=isinf(d0) | (d0==0);
while(~isempty(find(~L0,1)))
% reference starting points in DXY
idx=sub2ind(shape(1:2),x(~L0),i0(~L0));
% reduce out trivial points (DXY==0)
L=(DXY(idx)==0);
L0(~L0)=L;
idx=idx(~L);
if(isempty(idx)) continue; end
% these are current best parabolas for starting points
ik=DK(idx);
% these are new values from parabola from line #i
dtmp=d0(~L0)+aspect(2)^2*(i0(~L0)-i).^2;
% these starting points are OK - below the envelop
L=D1(idx)>dtmp; D1(idx(L))=dtmp(L);
% points which are still above the envelop but ik==i0,
% will not get any better, so fix them as well
L=L | (ik==i0(~L0));
% all other points are not OK, need new starting point:
% starting point should be at least below parabola
% beating us at current choice of i0
% solve quadratic equation to find where this happens
ik=(ik-i);
di=(D1(idx(~L))-dtmp(~L))./ik(~L)/2/aspect(2)^2;
% should select next highest index to the equality
di=fix(di)+sign(di);
% the new starting points
idx=find(~L0);
i0(idx(~L))=i0(idx(~L))+di;
% update L0 to indicate which points we've fixed
L0(~L0)=L; L0(idx(~L))=(di==0);
% points that went out of boundaries can't get better;
% fix them as well
idx=idx(~L);
idx=idx((i0(idx)<1) | (i0(idx)>shape(2)));
i0(idx)=i;
L0(idx)=1;
end
% reduce out trivial points DXY(idx)<DXY(:,i)
idx=sub2ind(shape(1:2),x,i0);
L=(DXY(idx)>0) | (i0==i);
idx=idx(L);
% these will keep track along which X should
% keep updating distances
map_lower=L;
map_upper=L;
idx_lower=idx;
idx_upper=idx;
% set trivial pixels D==0 in line #i:
% this has to be done b/s we manually discarded them from L0
D1(d0==0,i)=0;
% scan from starting points for each X,i0 in increments of 1
di=0; % distance from current y-line
eols=2; % end-of-line-scan flag
totlen=prod(shape(1:2));
while(eols)
eols=2;
di=di+1;
% select X which can be updated for di<0;
% i.e. X which had been below envelop all way till now
if(~isempty(idx_lower))
% shift y by -1
idx_lower=idx_lower-shape(1);
% prevent index dropping below 1st
L=(idx_lower>=1);
map_lower(map_lower)=L;
idx_lower=idx_lower(L);
if(~isempty(idx_lower))
dtmp=d0(map_lower)+...
aspect(2)^2*(i0(map_lower)-di-i).^2;
% these pixels are to be updated with i0-di
L=D1(idx_lower)>dtmp & DXY(idx_lower)>0;
map_lower(map_lower)=L;
idx_lower=idx_lower(L);
D1(idx_lower)=dtmp(L);
DK(idx_lower)=i;
end
else % if this is empty, get ready to quit
eols=eols-1;
end
% select X which can be updated for di>0;
% i.e. X which had been below envelop all way till now
if(~isempty(idx_upper))
% shift y by +1
idx_upper=idx_upper+shape(1);
% prevent index from going over array limits
L=(idx_upper<=totlen);
map_upper(map_upper)=L;
idx_upper=idx_upper(L);
if(~isempty(idx_upper))
dtmp=d0(map_upper)+...
aspect(2)^2*(i0(map_upper)+di-i).^2;
% check which pixels are to be updated with i0+di
L=D1(idx_upper)>dtmp & DXY(idx_upper)>0;
map_upper(map_upper)=L;
idx_upper=idx_upper(L);
D1(idx_upper)=dtmp(L);
DK(idx_upper)=i;
end
else % if this is empty, get ready to quit
eols=eols-1;
end
end
end
end
D{k}=D1;
end
%%%%%%%%%%%%% scan along Z %%%%%%%%%%%%%%%%
% this will be the envelop:
% envelop has to be initialized at Inf to ensure single direction of scan,
% otherwise may end up in infinite loop when trying to find starting point
D1=cell(size(D));
for k=1:shape(3) D1{k}=repmat(Inf,shape(1:2)); end
% these will be the references to parabolas defining the envelop
DK=cell(size(D));
for k=1:shape(3) DK{k}=repmat(Inf,shape(1:2)); end
% start building the envelope
for k=1:shape(3)
% need to select starting point for each X:
% * starting point should be below current envelop
% * k0==k is not necessarily a starting point
% * there may be no starting point
% k0 is the starting points for each XY: k0(XY) is the first
% z-index such that parabola from line k is below the envelop
% initial starting point guess is current slice
k0=repmat(k,shape(1:2));
% L0 indicates which starting points had been fixed
L0=isinf(D{k}) | (D{k}==0);
idxtot=find(~L0);
while(~isempty(idxtot))