Code covered by the BSD License  

Highlights from
Vessel branch segmentation

image thumbnail

Vessel branch segmentation

by

 

Segment the vessel branches from dynamic image of fluorescent microscopy

[appTime, mask]=VBSvesselMask(vimg)
function [appTime, mask]=VBSvesselMask(vimg)

tic
orgImg=vimg;
%% temporal smoothing 
% injection time was 2.23 ml/min = 2.715 ul/frame
% if blood flow is almost constant, 15 ul is 5 frame
% so window size was defined 5
windowSize=5;
vimg=MedfiltTime(vimg,windowSize);
vimg=single(vimg);
vimg=smooth3(vimg,'box',[1 1 windowSize]);
vimgSize= size(vimg);
toc
%% calculate vessel mask

% when the maximum intensity at each pixel was used as the criterion, noise
% pixel was also identified as vessel pixel. So, the pixels which maximum 
% intensity was less then mean + 2 sigma were eliminated.
meanImg=mean(single(orgImg(:,:,windowSize+1:30)),  3);
stdImg = std(single(orgImg(:,:,windowSize+1:30)),0,3);
maxImg=max(vimg,[],3);
meanStdMask=(maxImg>(meanImg+2*stdImg)); %clear meanImg stdImg
% maxImg     = maxImg.*meanStdMask;
vimg=vimg.*repmat(meanStdMask, [1,1,vimgSize(3)]);
toc
% extract more than half max. 
% vimg was normarized to 0, then maximum intensity at each pixel was
% reculculated.

vimg   = vimg-repmat(min(vimg(:,:,windowSize+1:50),[],3),[1,1,vimgSize(3)]);
maxImg = max(vimg,[],3);
binVimg= (vimg >repmat(maxImg.*0.5,[1,1,vimgSize(3)]));
% eliminate isolated voxels continuous half maximum less than window size 
binVimg= TemporalAreaOpen(binVimg,windowSize);
binVimg= TemporalSumFilter(vimg.*binVimg);
binVimg= SpatialAreaOpen(binVimg);
binVimg= TemporalAreaOpen(binVimg,windowSize);
toc

[~,appTime]=min(~binVimg,[],3);
mask=appTime>1;

function vfoo=TemporalAreaOpen(foo, cutSize)
cutSize=cutSize+1;
fooSize = size(foo);
vfoo=ipermute(reshape(bwareaopen(reshape(permute(foo,[3 1 2]), [],1),cutSize),[fooSize(3) fooSize(1) fooSize(2)]),[3 1 2]);

function foo=SpatialAreaOpen(foo)
fooSize = size(foo);
for n=1:fooSize(3)
     vfoo=medfilt2(foo(:,:,n),[3,1]);
     hfoo=medfilt2(foo(:,:,n),[1,3]);
     foo(:,:,n)= bwareaopen(vfoo.*hfoo,2,4);
end

function fimg=MedfiltTime(vimg,windowSize)
padSize=floor(windowSize*0.5);
vimgVec=zeros([size(vimg), windowSize] , class(vimg));
vimg=padarray(vimg,[0 0 padSize],'symmetric');
for n=1:windowSize
    vimgVec(:,:,:,n)=vimg(:,:,n:end-windowSize+n);
end

vimgVec=sort(vimgVec,4);
fimg=squeeze(vimgVec(:,:,:,padSize+1));


function vimg=TemporalSumFilter(vimg)
% Temporal Window by sum of each frame
vimgSize= size(vimg);
timg=sum(sum(vimg,1),2);
timg=timg-min(timg(:));
timg=timg./max(timg(:));
tmask=(timg>=0.01);
vimg=vimg.*repmat(tmask,[vimgSize(1:2), 1]);
return;

Contact us