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

[map sepTime]=VBSregionMerge(appTimeImg)
function [map sepTime]=VBSregionMerge(appTimeImg)

% delete small regions
appTimeImg=appTimeImg.*CutSmallRegions(appTimeImg);
% load vesselClustering.mat;
cc=AppTime2ConnComp(uint8(appTimeImg));
appTimeImg=double(appTimeImg);
labelMap=labelmatrix(cc);%---------------------------
fprintf('Initial Number of Segments: %d \n',cc.NumObjects);



[maxInArtery, maxInVein]=FindLocalMax(appTimeImg);

lt=zeros(1,cc.NumObjects);% lt: label time;
for n=1:cc.NumObjects;
    lt(n)=appTimeImg(cc.PixelIdxList{n}(1));
end;

clusters.img=zeros(cc.ImageSize);
clusters.num=1;
%prevNumObj=-1000;
n=0;
tic;
fprintf('Merge regions...\n')
while true;
    n=n+1;
    
    % check histogram
    sepTime=CheckHistoHasZero(lt,maxInArtery,maxInVein);
    if sepTime(1)~=false;
        msgbox('STOP')
        break;
    end
    %     if prevNumObj==cc.NumObjects || cc.NumObjects < 1 ;
    %         break;
    %     end
    
    %     prevNumObj=cc.NumObjects;
    
    [areaVec index]=sort(struct2array(regionprops(cc,'area')));
    
    % ------ for debug
        if mod(n,1000)==0
            if n>25000; keyboard; end
            toc
            % figure('name',sprintf('iteration: %d',n));imagesc(CCLT2appTimeMap(cc,index,lt));drawnow;
            % drawnow;
            disp(n);
            tic;
        end
    % ------ end debug
    
    tmpbw=false(cc.ImageSize);
    tmpbw(cc.PixelIdxList{index(1)})=true;
    tmpbw=xor(bwmorph(tmpbw,'dilate'),tmpbw) & (labelMap~=0);
    edgeLabel= labelMap(tmpbw(:));
    
    if isempty(edgeLabel)
        % surrounded by backgound -> defined as cluster
        [cc labelMap clusters lt]=DefineAsCluster(cc,clusters,index(1),lt);
        continue; % next step
    end
    
    maxEdgeLabel = max(edgeLabel);
    minEdgeLabel = min(edgeLabel);
    
    % surrounded by singlelabel
    if minEdgeLabel==maxEdgeLabel
        edgeTimes  =appTimeImg(tmpbw(:));
        if abs(lt(index(1)) - edgeTimes(1))<4
            [cc labelMap lt]=MergeRegions(cc,lt,minEdgeLabel,index(1));
        else
            [cc labelMap clusters lt]=DefineAsCluster(cc,clusters,index(1),lt);
        end
        continue; % next step
    end
    
    % surrounded by multiple labels ----------------------
    % Edge has same appearrance time -> merge
    currentTime=appTimeImg(cc.PixelIdxList{index(1)}(1));
    edgeTimes  =appTimeImg(tmpbw(:));
    doesEdgeHasCurrentTime=(edgeTimes==currentTime);
    if any(doesEdgeHasCurrentTime)
        dilateLabel=edgeLabel(doesEdgeHasCurrentTime);
        [cc labelMap lt]=MergeRegions(cc,lt,dilateLabel(1),index(1));
        continue; % next step
    end
    %-----------------------------------------------------
    timeDiff = abs(edgeTimes-currentTime);
    minDiff  = min(timeDiff);
    if minDiff < 4 % difference is 1 - 3
        minIndex = (timeDiff==minDiff);
        uniqueMinTime= unique(edgeTimes(minIndex));
        if numel(uniqueMinTime)==1
            % appearance time around current region has unique
            % neighbouring level
            dilateLabel=edgeLabel(minIndex);
        else
            % clustering by frequency of neighbouring AppTimes in Edge
            umtIndex1=(edgeTimes==uniqueMinTime(1)); umtIndex2=(edgeTimes==uniqueMinTime(2));
            switch WhichIsBigger(sum(umtIndex1),sum(umtIndex2))
                case 1; dilateLabel = edgeLabel(umtIndex1);
                case 2; dilateLabel = edgeLabel(umtIndex2);
                case 0
                    % iterative clustering by frequency of neighbouring AppTimes
                    % in image
                    switch IterativeSelectionByTimeCountsInImage(uniqueMinTime,lt,cc);
                        case 1; dilateLabel = edgeLabel(umtIndex1);
                        case 2; dilateLabel = edgeLabel(umtIndex2);
                    end
            end
        end
        [cc labelMap lt]=MergeRegions(cc,lt,dilateLabel(1),index(1));
        continue; % next step
    end
    
    %     % regional mode filter (current pixel )
    %     modeTimeInEdge=mode(edgeTimes);
    %     % !!! WARNING !!!: mode([1 1 2 2]) returns 1; thus artery region will
    %     % be broader.
    %     dilateLabel=edgeLabel(edgeTimes==modeTimeInEdge);
    %     [cc labelMap lt]=MergeRegions(cc,lt,dilateLabel(1),index(1));
    
    [cc labelMap clusters lt]=DefineAsCluster(cc,clusters,index(1),lt);
end

map=CCLT2appTimeMap(cc,lt);
map=map+clusters.img;
clim=map(map(:)~=0);
figure;imagesc(map,[min(clim)-1, max(clim)+1]); colormap([0 0 0; flipud(jet(max(clim)-min(clim)+1)); 0 0 0 ])
return;

% function bwImg=CutSmallRegionsAndFillSmallHoles(appImg)
function bwImg=CutSmallRegions(appImg)
% return 
bwImg=appImg>0;
bwImg=bwareaopen(bwImg,8,8);


function torf=CheckHistoHasZero(lt,maxInArtery,maxInVein)
tmplt=unique(lt(lt>maxInArtery & lt<maxInVein));
if numel(tmplt)==(maxInVein-maxInArtery-1)
    torf=false;
else
    torf=setdiff(maxInArtery+1:maxInVein-1,tmplt);    
end

function cc=AppTime2ConnComp(appTimeImg)
bw2d=appTimeImg>0;
bwSize=size(bw2d);
labelMap=zeros(bwSize);
numObj=0;
for n= min(appTimeImg(bw2d(:))) : max(appTimeImg(bw2d(:)))
    currentFlameBin=appTimeImg==n;
    [tmpLabel tmpNum]=bwlabel(currentFlameBin,8);
    labelMap=labelMap+(currentFlameBin.*(tmpLabel+numObj));
    numObj  =numObj+tmpNum;
end

cc=label2conncomp(labelMap);

function cc=label2conncomp(labelImg)
cc.Connectivity=8;
cc.NumObjects=max(labelImg(:));
cc.ImageSize=size(labelImg);
cc.PixelIdxList=cell(1,cc.NumObjects);
index=1:prod(cc.ImageSize);
for n=1:cc.NumObjects
    cc.PixelIdxList{n}=index(labelImg(:)==n)';
end
return;

function uniqueTime=IterativeSelectionByTimeCountsInImage(uniqueMinTimes,lt,cc)
% return value: 1 currentpixel should be merged uniqueMinTimes(1)
%               2 currentpixel should be merged uniqueMinTimes(2)
val=WhichIsBigger(uniqueMinTimes(1),uniqueMinTimes(2));
if val==1
    earlier=uniqueMinTimes(2); later =uniqueMinTimes(1);
else
    earlier=uniqueMinTimes(1); later =uniqueMinTimes(2);
end

while true;
    if earlier<0;error('cannot find time.');end
    index=WhichIsBigger(CountSameAppTime(earlier,lt,cc),CountSameAppTime(later,lt,cc));
    switch index;
        case {1,2}; break
        case 0; earlier=earlier-1; later=later+1; continue;
    end
end

if val==index % uniqueMinTimes(2)
    uniqueTime=2;
else          % uniqueMinTimes(1)
    uniqueTime=1;
end


function [maxInArtery, maxInVein]=FindLocalMax(appTimeImg)
nonZeroAppTime=appTimeImg(appTimeImg(:)>1);
histCentre=min(nonZeroAppTime):max(nonZeroAppTime);
histBins=hist(nonZeroAppTime,histCentre);
histCentre=[histCentre(1)-1, histCentre, histCentre(end)+1];
histBins  =[ 0 histBins 0];
%%%% Slow but more acculate
oldNum = inf;
while 1
    localMaxs =(histBins(2:end-1)>histBins(1:end-2)) & (histBins(2:end-1)>=histBins(3:end));
    localMaxBins  = histBins([false, localMaxs , false]);
    localMaxCentre= histCentre([false, localMaxs , false]);
    localMaxCentre= localMaxCentre(localMaxBins>(100));
    switch numel(localMaxCentre)
        case {0,1}
            error('Cannot find out two local maximums in histogram')
        case 2
            break
        otherwise
            if oldNum==numel(localMaxCentre);
                %%% DBG
                %tmpFgh=figure;
                %plot(histCentre, histBins, 'b.', localMaxCentre, localMaxBins, 'ro');
                %%% DBG
                break
            else
                oldNum=numel(localMaxCentre);
                histBins=medfilt1(histBins);
            end
    end
end
maxInArtery = localMaxCentre(1);
maxInVein   = localMaxCentre(end);
%%%% Fast but sometimes fault
%             localMins =(histBins(2:end-1)<histBins(1:end-2)) & (histBins(2:end-1)<histBins(3:end));
%             localMaxMinBins=    histBins([false, localMaxs | localMins, false]);
%             localMaxMinCentre=histCentre([false, localMaxs | localMins, false]);
%             localMaxMinBins=    localMaxMinBins(localMaxMinBins>(100));
%             localMaxMinCentre=  localMaxMinCentre(localMaxMinBins>(100));
%             [deepVal deepInd]=max(localMaxMinBins(1:end-2) + localMaxMinBins(3:end) - 2*localMaxMinBins(2:end-1));
%             [tmpVal tmpInd]=max(localMaxMinBins(1:deepInd));     maxInArtery=localMaxMinCentre(tmpInd(1));
%             [tmpVal tmpInd]=max(localMaxMinBins(deepInd+2:end)); maxInVein  =localMaxMinCentre(tmpInd(1)+deepInd+1);
% error('cannot find max freq time in artery and vein side');
% error('cannot find max freq time in artery and vein side');
return;


function num=CountSameAppTime(t,lt,cc)
index = 1:numel(lt);
sameTimeIndex = index(lt==t);
num=0;
for n=1:numel(sameTimeIndex)
    num=num+numel(cc.PixelIdxList{sameTimeIndex(n)});
end

function bigger=WhichIsBigger(a1, a2)
% 0: a1==a2
% 1: a1 >a2
% 2: a1 <a2
if     a1>a2
    bigger=1;
elseif a1<a2
    bigger=2;
else
    bigger=0;    
end

function map=CCLT2appTimeMap(cc,lt)
map=zeros(cc.ImageSize);
for n=1:cc.NumObjects
    map(cc.PixelIdxList{n})=lt(n);
end

% function indexupdateConnComp(cc)
function [cc labelMap lt]=MergeRegions(cc,lt,dilateLabel,index)
cc.PixelIdxList{dilateLabel} = vertcat(cc.PixelIdxList{dilateLabel},cc.PixelIdxList{index});
[cc lt]= DeleteMergedRegion(cc,lt,index);
labelMap=labelmatrix(cc);


function [cc labelMap clusters lt]=DefineAsCluster(cc,clusters,index,lt)
clusters.img(cc.PixelIdxList{index})=lt(index);
[cc lt]= DeleteMergedRegion(cc,lt,index);
labelMap=labelmatrix(cc);
clusters.num=clusters.num+1;


function [cc lt]= DeleteMergedRegion(cc,lt,index)
% delete integrated region from PixelIdxList
switch index
    case 1
        cc.PixelIdxList= cc.PixelIdxList(2:end);
        lt=lt(2:end);
    case cc.NumObjects
        cc.PixelIdxList= cc.PixelIdxList(1:end-1);
        lt=lt(1:end-1);
    otherwise
        cc.PixelIdxList=[cc.PixelIdxList(1:index-1) cc.PixelIdxList(index+1:end)];
        lt = [lt(1:index-1) lt(index+1:end)];
end
cc.NumObjects=cc.NumObjects-1;

Contact us