retinal fundus image branch issue

I am working on retinal fundus images. I have segmented the vascular network from the original image and have a black and white image as my vascular mask. From this, I would like to extract the main long vessels sans the small branches. I tried bwmorph to skeletonise the image, remove spurs, find the branchpoints and then to segment the vessels by removing the branchpoints but this creates false branches and breaks the larger vessels in to segments too as even the crossover points are being considered as branches. This leads to a final image with numerous tiny vessels whereas I want the main vessel tree to remain intact. I also tried bwareaopen to remove the small branches before skeletonising but that doesnt work either as the whole vascular network is being considered as one connected component.

I'd appreciate any ideas that can help me code so that the main vessel tree remains intact. I have looked in to various links too related to vessels, but they either deal with a single vessel at a time or just need to extract the longest vessel whereas I want my code to grab a good estimate of the tree which I can then segment in to seperate but proper vessels. Also I have more than 100 images to process so I am looking for some simple idea. Ofcourse there are papers available online that trace the whole tree but those are very complicated and again I already have the complete structure that I would like to simplify for further processing.

Pfa the vessel mask obtained and the final image. Thank you in advance.

<<

>>

[John BG]

Hi Mentor Geeek

this is John BG <mailto:[e-mail] [e-mail]>

You want to skeletonise and then trim, prune short branches.

How short should a short branch be to be pruned or remain in the picture?

For instance

.

Just trimming really short branches would remove the yellow ones, but it would also shorten the tips of the main vessels by the same amount.

When trimming longer branches such problem worsens, unless you want to consider defining in advance the length of the short branches to be removed, and then to sort out small branches as pegs or tips, and remove only pegs.

Spotting pegs means the script has to be able to tell when the short branch is 'end of line', or is just a peg off a main vessel.

Would you like to proceed this way?

[John BG]

bwmorph(B,'skel',Inf)

breaks the main vessels, rendering the result useless, correct?

[rabeeah]

Yes, it does. I am sorry I just saw your comments.

Answer 1

[John BG]

Hi again Mentor

this is John BG <mailto:[e-mail] [e-mail]>

The following script serves all short sections, let me call them 'pegs', ready ready for trimming.

There is a minor issue of command bwmorph that returns a few isolated white points. I have fixed it manually, in order to focus on the objective of compiling all and nothing else than pegs.

1.

Acquiring image

clear all;close all;clc

A=imread('001.jpg');     % make sure there is no thin (1 pixel wide) white frame around picture
B=imbinarize(A(:,:,1));
B2=bwmorph(B,'skel',Inf);             % skeletonize

2.

open point 1: bwmorph isolates a few points that should be connected: manual correction

B2(330,224)=1; 
B2(48,317)=1;B2(49,316)=1;B2(50,315)=1;

.

figure(1);imshow(B2);

3.

setting threshold for trimming, sections shorter that dth are considered for trimming, longer sections are ignored

dth=5;

4.

making sure picture frame is black

B2(1,:)=0;B2(end,:)=0;B2(:,1)=0;B2(:,end)=0;

5.

collecting locations of white pixels

[sz1_b,sz2_b]=size(B2)        % sz2 horizontal sz1 vertical         
[wx,wy,v]=find(B2);P=[wx wy];  % white points contained in P

PA={};                              % build cell container for white points
for k=1:1:length(P)
    PA=[PA;P(k,:)];
end

fc=zeros(length(P),1);    % vector with counters of how many white pixels around each white pixel
for k=1:1:length(P)
   p1=PA{k};    % read white point coordinates
   p1y=p1(1);p1x=p1(2);

     % mask 1
     p2=[p1y-1 p1x-1]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end        % top left, clock wise
     p2=[p1y-1 p1x]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end            % top
     p2=[p1y-1 p1x+1]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end       % top right
     p2=[p1y p1x+1]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end           % right
     p2=[p1y+1 p1x+1]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end       % down right
     p2=[p1y+1 p1x]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end           % down
     p2=[p1y+1 p1x-1]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end       % down left
     p2=[p1y p1x-1]; if B2(p2(1),p2(2)) p1=PA{k};PA{k}=[p1;p2];fc(k)=fc(k)+1; end            % left

end

6.

classifying points: whether tips sections or junctions

.

figure(2);stem(fc)          % comment: 1s are tips, 2s are sections, 3s and above are junctions

.

open point2: solving (1 green start on junction) thin and thick (more than 1 green start next to same junction location)

7.

compiling tip points

T=P(find(fc==1),:);                           
figure(1);hold all;plot(T(:,2),T(:,1),'r*');       % check

J=P;                       % compiling junction points
tp=find(fc==1);       % tips
sp=find(fc==2);       % sections
tps=(unique([tp;sp]));   % tips and sections
J(tps,:)=[];   % removing tip points and section points
figure(1);hold all;plot(J(:,2),J(:,1),'g*');       % check

.

8.

When considering iteration, successive capture and trimming of all short tips, it's going to cut all end-of-line branches that that as commented earlier on, it shouldn't happen.

P2=P;
P2(tp,:)=[];  % trimming tips P2 = P - T
figure(10);hold all;plot(P2(:,2),P2(:,1),'b.');

.

open point 3: there are a few tips that have not been captured

.

.

% 004

So far

P; % white points
T; % all tips. A few tips missed, mask 1 can be refined to catch all tips
J; % junction points

9.

approach 1: J points, T points, for each T point find nearest J point

JT=zeros(length(T),2);   % nearest junction point to each tip point
for k=1:1:length(T)
    pt1=T(k,:);
    d=[];
    for s=1:1:length(J)
        pt2=J(s,:);
        d=[d;((pt1(1)-pt2(1))^2+(pt1(2)-pt2(2))^2)^.5]; % euclidean distance
    end
    dmin=find(d==min(d));
    JT(k,:)=J(dmin(1),:);
end

%  check
close all;figure(1);imshow(B2);hold all
for k=1:1:length(T)
    figure(1); plot([T(k,2) JT(k,2)],[T(k,1) JT(k,1)],'y')
end

comment: too many tips jump to nearest (euclidean) junction that happens to be on a different branch

% 005

10.

approach 2: move pointer from tip to nearest junction and measure distance

d=zeros(length(T),1); % vector containing distances from each tip to nearest junction

TR={}      % log for branches to trim

for k=1:1:length(T)
    pn1=T(k,:);          % read tip
    TR{k,1}=pn1;        % copy tip into short branch log

      pn1x=pn1(2);;pn1y=pn1(1);

     % mask 2
     pr1=[pn1y-1 pn1x-1]; if B2(pr1(1),pr1(2)) TR{k,2}=pr1;pn2=pr1;d(k)=d(k)+1; end        % top left, clock wise
     pr2=[pn1y-1 pn1x]; if B2(pr2(1),pr2(2)) TR{k,2}=pr2;pn2=pr2;d(k)=d(k)+1; end            % top
     pr3=[pn1y-1 pn1x+1]; if B2(pr3(1),pr3(2)) TR{k,2}=pr3;pn2=pr3;d(k)=d(k)+1; end       % top right
     pr4=[pn1y pn1x+1]; if B2(pr4(1),pr4(2)) TR{k,2}=pr4;pn2=pr4;d(k)=d(k)+1; end           % right
     pr5=[pn1y+1 pn1x+1]; if B2(pr5(1),pr5(2)) TR{k,2}=pr5;pn2=pr5;d(k)=d(k)+1; end       % down right
     pr6=[pn1y+1 pn1x]; if B2(pr6(1),pr6(2)) TR{k,2}=pr6;pn2=pr6;d(k)=d(k)+1; end           % down
     pr7=[pn1y+1 pn1x-1]; if B2(pr7(1),pr7(2)) TR{k,2}=pr7;pn2=pr7;d(k)=d(k)+1; end       % down left
     pr8=[pn1y pn1x-1]; if B2(pr8(1),pr8(2)) TR{k,2}=pr8;pn2=pr8;d(k)=d(k)+1; end            % left
    s=3;  % index for spear head pixel
     while ~sum(and(J(:,1)==pn2(1),J(:,2)==pn2(2)))   % while next point not junction, 2D and a(and(a(:,1)==b(1),a(:,2)==b(2)),:)
             r=[pn2(1)-1 pn2(2)-1;
                  pn2(1)-1 pn2(2);
                  pn2(1)-1 pn2(2)+1;
                  pn2(1) pn2(2)+1; 
                  pn2(1)+1 pn2(2)+1;
                  pn2(1)+1 pn2(2);
                  pn2(1)+1 pn2(2)-1;
                  pn2(1) pn2(2)-1]

             r(and(r(:,1)==pn1(1),r(:,2)==pn1(2)),:)=[];  % remove previous point to avoid endless back forward on same couple pixels

             for n=1:1:length(r)
                 p03=r(n,:);
                 if B2(p03(1),p03(2)) 
                     TR{k,s}=p03;d(k)=d(k)+1; 
                 end 
             end
             pn3=TR{k,s}
             pn1=pn2   % shift 1 pixel towards junction
             pn2=pn3
             s=s+1;

     end

end

distances tips to nearest junctions, counting pixels 'following the pipes'

figure(11);stem(d)

.

close all;figure(1);imshow(B2);hold all    % visualize sections tips to nearest junction only, any length
figure(1);hold all;plot(T(:,2),T(:,1),'r*');    
[sz1tr sz2tr]=size(TR)
for k=1:1:sz1tr
    L=[]
    for s=1:1:sz2tr
         L=[L;TR{k,s}];
   end
   figure(1);plot(L(:,2),L(:,1),'g.')    
end

TR2=TR;           % compiling tip to nearest junction sections
for k=1:1:sz1tr
    L=[]
    for s=1:1:sz2tr
         L=[L;TR2{k,s}]
    end
    if length(L)>dth
        for p=1:1:sz2tr
                TR2{k,p}=[];
        end
    end
end

close all;figure(1);imshow(B2);hold all      % visualizing tip to nearest junction short sections only
figure(1);hold all;plot(T(:,2),T(:,1),'r*');  
for k=1:1:sz1tr
    L=[]
    for s=1:1:sz2tr
         L=[L;TR2{k,s}];
    end
    if ~isempty(L)
        figure(1);plot(L(:,2),L(:,1),'g.')   
    end
end

% 008

.

Dinner is ready, who wants to trim, go for the green.

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

open point 4: end-of-line twin short pegs, with this basic script, both get trimmed. Further refinement required to, when 2 or more short pegs right at the end of line, cur all but one

% 009

open point 5: there is just one T shaped double tip not caught yet

[Jan]

The iterative growing of arrays wastes resources. To find the smallest distance, you can omit to calculate the SQRT() of the sum of squares, because the minimum of X (with positive values) is at the same index as the minimum of SQRT(X). The part:

JT=zeros(length(T),2);
for k=1:1:length(T)
  pt1=T(k,:);
  d=[];
  for s=1:1:length(J)
      pt2=J(s,:);
      d=[d;((pt1(1)-pt2(1))^2+(pt1(2)-pt2(2))^2)^.5];
  end
  dmin=find(d==min(d));
  JT(k,:)=J(dmin(1),:);
end

can be simplified to:

JT = zeros(size(T, 1), 2);
for k = 1:size(T, 1)
  [~, dmin] = min((T(k,1)-J(:,1)).^2 + (T(k,2)-J(:,2)).^2);
  JT(k,:)   = J(dmin,:);
end

Or

PA={};
for k=1:1:length(P)
  PA=[PA;P(k,:)];
end

by:

PA = num2cell(P, 2);

And

 L = [];
 for s=1:1:sz2tr 
   L = [L;TR2{k,s}]
 end

by

L = cat(1, TR2{k, :});

Avoiding of an iterative growing of arrays does not change the nature of this suggested solution, but it is essential for a clean and efficient programming. Therefore it is suggested every day in this forum. Even if the job is not time-critical for the OP mentor geeek, it helps to save energy and it is worth to mention this good programming practice again.

[John BG]

Hi Jan

thanks for the constructive comments and possible improvements.

1.

the Euclidean distance calculation lines of script would be improved if using your suggestion. However, for the sake of saving uploaded images, I did not show on picture the comment that given the nature of this network, many tips end up connected to junctions which are not their nearest junction, following the white pixels.

As mentioned, the Euclidean distance doesn't help in this case, but I have already made a quick to make sure that in next question requiring Euclidean distance calculations I am considering first your suggestion.

2.

In

PA={};
for k=1:1:length(P)
  PA=[PA;P(k,:)];
end

out of many ways to append points to an already list of points, I am adding on the tail of PA, a point contained in P. Still do not see how your suggestion

PA = num2cell(P, 2);

really appends a point of P 2xN doubles on PA, a cell containing lines of points, but each line may have variable length.

3.

regarding replacing

L = [];
 for s=1:1:sz2tr 
   L = [L;TR2{k,s}]
 end

with

L = cat(1, TR2{k, :});

Yes it may work.

Jan, I bear your suggestions in high regard, and always learn whenever reading your answers and questions, but right now the core of the problem is not syntactic, but like Image Analyst points out, the tools applied so far either don't work, or perform excessive trimming, highlighted in point 8 of my answer.

This ocular background question may help doctors at work, so it's primordial to remove all short pegs but not the tips of long sections.

As important as it is in any coding, right now I focus on meeting the question requirements, putting for a but latter the syntactic optimization that you undoubtedly have already mastered, thanks for sharing.

My script gets all short pegs, allowing further processing damaging long sections, thus answering the question.

Having said this, I have listed 4 open points that make the difference between a MATLAB forum correct answer and the sketch for a robust application, like

• the skeleton still returns a few isolated pixels, shouldn't,
• improving the masks to discern T-shapped end of lines or
• what to do where there are end of lines with more than 2 short pegs

More robust script should cut all but one of those identical or almost identical pegs right at the end of a line. But there's now need for Mentor Geeek feedback.

I expect Mentor to comment on what way such processing is going to be used regarding my guess that, once spotting correctly all small pegs as I have just did may be enough for then perform a visual inspection that the trained eye can focus on the pegs, and then with perhaps a sigle mouse click on a given peg, the script would remove it, rather than going for a fully automated removal, that may not be required.

Appreciating time and attention, awaiting Mentor feedback

regards

John BG

<mailto:[e-mail] [e-mail]>

[Jan]

To 2.: "[...] appends a point of P 2xN doubles on PA, a cell containing lines of points, but each line may have variable length." Compare it: These codes create the same output PA:

P = rand(1e4, 2);
tic;
PA1 = {};                   % Your code
for k = 1:1:length(P)
   PA1 = [PA1;P(k,:)];
end
toc

tic;                       
PA2 = cell(size(P, 1), 1);  % Pre-allocate
for k = 1:size(P, 1)
   PA2{k} = P(k, :);
end
toc

tic;                       % Nicer and slower:
PA3 = num2cell(P, 2);
toc

isequal(PA1, PA2, PA3)     % 1: Same result

% Elapsed time is 1.228125 seconds.
% Elapsed time is 0.027210 seconds.
% Elapsed time is 0.052116 seconds.

Such a small modification of the code can give a speedup of factor 45 for relatively small inputs. I could not resist to mention this common pitfall of omitted pre-allocation, although it is a side note only and in the problem of the OP, P is tiny. Let's see, what the OP comments about the real problem.

Answer 2

[John BG]

Hi Rabeeah, Mentor Geek

this is John BG <mailto:[e-mail] [e-mail]>

would this result be acceptable

1.

.

2.

.

this fine smoothing of just the small pegs is achieved with the initial script modified the following way:

clear all;close all;clc

dth=5;  % trimming threshold, all pegs closer or at 3 distance to nearest junction are trimmed

A=imread('001.jpg');     % load image
[sz1A sz2A sz3A]=size(A)  %sz2A is x, horizontal. sz1A is y, vertical
A1=A(:,:,1);
hfg1=figure(1);imshow(A1)
% hfg1.ToolBar='none';
hfg1.MenuBar='none';hfg1.DockControls='off';hfg1.Name='start image';hfg1.Position=[500 200 2*sz2A 2*sz1A];

% hfg1.OuterPosition=[500 200 2*sz2A 2*sz1A];
% hfg1.InnerPosition=[500 200 2*sz2A 2*sz1A];

B=imbinarize(A1);
B0=B;   % reference
figure(25);imshow(B0)
hfg25.MenuBar='none';hfg25.Name='skeletonized image';hfg25.Position=[510 210 2*sz2A 2*sz1A];

for s=1:1:10
B2=bwmorph(B,'skel',Inf);             % skeletonize
hfg2=figure(2);imshow(B2);
hfg2.MenuBar='none';hfg2.Name='skeletonized image';hfg2.Position=[510 210 2*sz2A 2*sz1A];
% 001
% hfg2.MenuBar='none';
% hfg2.ToolBar='none';hfg1.DockControls='off';

% bwmorph has no error: a few short sections of the skeleton are disconnected. 
% Manual correction for following code that requires junctions to exist along all sections
% 001-2 001-3
% Pc=[330 224;48 317;49 316;50 315];
% B2(330,224)=1;B2(48,317)=1;B2(49,316)=1;B2(50,315)=1;
% hfg3=figure(3);imshow(B2);
% hfg3.MenuBar='none';hfg3.Name='bwmorph isolated sections correction';hfg3.Position=[520 220 2*sz2A 2*sz1A];
% figure(3);hold all;plot(Pc(:,2),Pc(:,1),'r*')

% B2=B3  % iteration point for simple trimming, cuts too much

B2(1,:)=0;B2(end,:)=0;B2(:,1)=0;B2(:,end)=0; % make sure 1 pixel wide picture frame is black

[sz1_b,sz2_b]=size(B2)        % sz2 horizontal sz1 vertical         
[wy,wx,v]=find(B2);P=[wy wx];  % locations of white points, now contained in P

PA={};                              % build cell container to log how many white points around each white point
for k=1:1:length(P)             % allocating each white point location at the head of each PA line
    PA=[PA;P(k,:)];
end

fc=zeros(length(P),1);    % vector with counters of how many white pixels around each white pixel
for k=1:1:length(P)
   p1=PA{k};    % read white point coordinates
   p1y=p1(1);p1x=p1(2);

     % mask 1.1 capturing straight tips
     p2=[p1y-1 p1x-1]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end        % top left, clock wise
     p2=[p1y-1 p1x]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end            % top
     p2=[p1y-1 p1x+1]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end       % top right
     p2=[p1y p1x+1]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end           % right
     p2=[p1y+1 p1x+1]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end       % down right
     p2=[p1y+1 p1x]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end           % down
     p2=[p1y+1 p1x-1]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end       % down left
     p2=[p1y p1x-1]; if B2(p2(1),p2(2)) p0=PA{k};PA{k}=[p0;p1];fc(k)=fc(k)+1; end            % left

end

PA2={};                              % build cell container for bent tip points
for k=1:1:length(P)             
    PA2=[PA2;P(k,:)];
end
fc2=zeros(length(P),1);        % meter for bent tip points

for k=1:1:length(P)
     p1=PA2{k};   
     p1y=p1(1);p1x=p1(2);
     r1=[p1y-1 p1x-1;        % top left, clock wise
          p1y-1 p1x;           % top
          p1y-1 p1x+1;     % top right
          p1y p1x+1;          % right
          p1y+1 p1x+1;        % down right
          p1y+1 p1x;      % down
          p1y+1 p1x-1;    % down left
          p1y p1x-1];          % left
     rm=zeros(1,length(r1));
     for k2=1:1:length(r1)
         rm(k2)=B2(r1(k2,1),r1(k2,2));
     end
     rm2=find(rm);
     if numel(rm2)==2 && (abs(rm2(1)-rm2(2))==1 || (abs(rm2(1)-rm2(2))==length(r1)-1))
         p0=PA2{k};PA2{k}=[p0;p1];
         fc2(k)=fc2(k)+1;
     end
end

figure(10);stem(fc)  % 002        % 1s are tips, 2s are sections, 3s and above are junctions
figure(12);stem(fc2)         % fc2 marks P points that are bent tips 

% pending 1: merge thin and thick junctions difference

T=P(find(fc==1),:);                            % compiling straight tip points
TL=P(find(fc2==1),:);                            % compiling bent tip points
% T=[T;TL]; % including bent tips in list of all tip points, then the while loop measuring distance 
% to nearest junction requires a more complicationed logical condition 

figure(4);imshow(B2);hold all;plot(T(:,2),T(:,1),'r*');       % check
figure(4);plot(TL(:,2),TL(:,1),'Marker','*','Color','c','LineStyle','none');  
% 201

figure(25);hold all;plot(T(:,2),T(:,1),'r*');       % check
figure(25);plot(TL(:,2),TL(:,1),'Marker','*','Color','c','LineStyle','none');  

for k=1:1:length(T)
    B(T(k,1),T(k,2))=0;
end

end

figure(26);imshow(B)
hfg26.MenuBar='none';hfg26.Name=' ';hfg26.Position=[510 210 2*sz2A 2*sz1A];

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

[rabeeah]

Hi john, I think theres some problem with my internet so I emailed you the image on the id that appears here on the forum.I apologize for the late replies.

[rabeeah]

Wow. bless you john bg. So if I increase dth, it could eliminate some more of the shorter vessels? The short ones attached to the large vessels? My question would be nore clear if you could please check your email wheb u have time. Thank you so much for all the help.

[Walter Roberson]

Akadne Noah comments to John BG:

thanks so much, the code works perfectly well.More power to your elbow

Answer 3

[Image Analyst]

To extract the largest blob only use bwareafilt():

mainTreeOnly = bwareafilt(binaryImage, 1);

[rabeeah]

Thanks a lot image analyst.we Have been stuck on this issue for a while now. Will try and get back to you.

[rabeeah]

@ Image analyst, This doesnt seem to work.it gives me the whole mask as it is. I think the reason being that the whole of it is being considered as one connected component. I'd be grateful for any other ideas. Thank you

[Image Analyst]

Well it looks like they're all connected so it's all the main tree. Though there is one small tree in the lower right that is separate.

[rabeeah]

True that. But I want some idea to eliminate the small branches u can see on the large vessels, so that when I skeletonise it, the large vessels are not broken into small segments leading to "numerous false vessels" as u can see in the colored image. Or is there some way I can fill it up again so that it looks like the original?

[rabeeah]

And really appreciate your prompt responses. Thanks !!

[Image Analyst]

Post an image of what you'd want. For every vessel there it's not clear which you want to keep and which you want to reject. Do you have any criteria that you can explain mathematically to make that distinction?

[John BG]

IA, the image that Mentor Geeek wants to process is the second image in the question.

The faulty result is the first image in the question, the one with the coloured branches.

[Image Analyst]

Yes, I knew that. But you can't just skeletonize and prune or else you'll get the problems you stated above in your comment. I've talked with The Mathworks about getting code to extract the longest path ("spine") of a skeleton but so far it's not in there yet. And anyway, mentor would have lots of such long paths. It's further complicated by the structure not being a pure tree-like structure, but often a lattice-like structure with loops or crossing branches. So it's not an easy, straightforward problem.

The literature is filled with papers on how to extract the vasculature from retinal images so maybe he should look for a method that gives a better initial segmentation here: http://www.visionbib.com/bibliography/contentsmedical.html#Medical%20Applications,%20CAT,%20MRI,%20Ultrasound,%20Heart%20Models,%20Brain%20Models

Another thing: I'm not sure why he wants to split the vessels into short segments between branchpoints. However it might help. If you skeletonized, then identified endpoints with bwmorph(binaryImage, 'endpoints'), then found branchpoints with bwmorph(binaryImage, 'branchpoints'), then you could label the split up image and use find() and bwselect() (or alternatively bwlabel() and ismember) to identify and remove those branches that touch endpoints. Of course, it will remove the tip of the longest branch also, so you might have to then find the shortest path from the starting point (which in itself is difficult to identify) to each of the endpoints, and restore the branch which was part of the longest of the shortest paths. See Steve's blog: http://blogs.mathworks.com/steve/2011/11/01/exploring-shortest-paths-part-1/

[rabeeah]

@Image analyst. Will post an image of what I would like to keep. Also will try the method u have mentioned and will get back to you. I really hope it helps. I am actually trying to determine the tortuousity of the vessels and for the code I took help from one of matlabs retinal image processing videos online but I am not getting good results and so I needed help from the forum.

Answer 4

[rabeeah]

This is what I would like to keep from the entire image. Some more vessels here and there are ok. But I want these ones to remain intact. And since I have a lot of images,and this is one od the preprocessing steps I would like the system to be fully automated. I am not sure why the image isnt getting uploaded.

[Jan]

UIploading an image required three steps: 1. Hit the "Image" button, 2. select the image from your disk, 3. Hit the "Insert Image" button.

[John BG]

Hi Rabeeah

this is John BG <mailto:[e-mail] [e-mail]>

The supplied image is enough for me to complete a working solution, but I need to know how do you want to proceed regarding the following points:

1.

So far the shorted branches are now available, on the skeleton. Is it ok for you to work on the skeleton, or would you like the skeleton with tips and short highlighted branches copied on the original image?

2.

There are many end of lines with 2 or more identical or quite similar short branches. In such cases, do you want to cut them all, or cut all but one.

If cut all but one, which one is the one you want to keep, the one that is aligned (same or similar slope) with the branch connecting to nearest junction?

3.

There's also need to make sure that the skeleton does not have isolated branches, and so far, there are a couple of isolated branches that I fixed manually.

In order to automate it, would you please be so kind to comment whether the way I manually connected the isolated branches is ok?

I tried to connect the isolated branches to the nearest branch that had similar slope.

4.

Do you want the trimmed short branches still visible, for instance in red colour? or would you rather remove the trimmed short branches and replace with background colour?

5.

a few really short branches not caught as tips.

Have a look at picture next to open point 3in my answer.

Is it ok to ignore the few pegs that get away?

[John BG]

Hi Rabeaah and Mentor

this is John BG <mailto:[e-mail] [e-mail]>

the example shown so far has a threshold of 5 pixels (variable dth), green pixels show short branches of 5 or shorter length while the red dots are branch tips.

.

Would you please clarify how you want the following specific cases of short ending trimmed:

.

.

On this image, in yellow, cases 1 2 and 3

1.

three leaf clovers

do you cut all short branches?

or do you want to keep the one that has same or similar slope compared to main branch?

2.

branches ending with 2 bent pegs

how do you decide which peg to chop?

or since the length of the main vessel is only reduced by 5 pixels, would it be ok to just cot both?

3.

longer than 5 pixels short branches that still may look like short branches.

Please tell whether boxed examples with number 3 on the above picture can be considered main vessels, or you still want them removed?

As soon as your decision criteria is known on these points, the script can further evolve towards the desired automated removal.

John BG

<mailto:[e-mail] [e-mail]>

Answer 5

[John BG]

Hi again

this is John BG <mailto:[e-mail] [e-mail]>

ok, let's recap

1.

WRONG START IMAGE

The image supplied in the question is not the outcome of human eye fundus scans.

It probably is the result of Mentor Geeek or Rabeaah initial effort.

The processing should start with images like the ones obtained in

https://www.google.com/search?site=&tbm=isch&source=hp&biw=987&bih=697&q=retinal+fundus&oq=retinal+fundus&gs_l=img.3..0l6j0i5i30k1j0i8i30k1l3.1545.3355.0.6748.14.14.0.0.0.0.121.1063.11j3.14.0....0...1.1.64.img..0.14.1057.vb_qeQrfyg8#imgrc=_FxY4EpMtFJboM:

example

.

2.

NOT ENOUGH PIXELS

The stair effect shown in the star image initially supplied in the question does not have enough resolution.

.

3.

FILTERING WITH BIFURCATION REFERENCE

BLACK & WHITE RESULTS UTAH UNIVERSITY REFERENCE

The University of Utah published this short note:

Retinal Fundus Images, Ground Truth of Vascular Bifurcations and Crossovers

http://webvision.med.utah.edu/2011/10/retinal-fundus-images-ground-truth-of-vascular-bifurcations-and-crossovers/

.

Where it's mentioned that a database of interest has been compiled, paying attention to bifurcations and crossovers detection. It points at

http://www.cs.rug.nl/~imaging/databases/retina_database/retinalfeatures_database.html

The paper

Detection of Retinal Vascular Bifurcations by Trainable V4-Like Filters

that explains how such B&W results have been obtained was written by G. Azzopardi N. Petkov

http://www.cs.rug.nl/~imaging/databases/retina_database/CAIP2011.pdf

This paper shows how a particular junction has been chosen, and then using Gabor filters, more or less alike junctions are obtained. Interesting that it's mentioned the difference of count whether using filter without rotation or with rotation.

So, Azzopardi Petrov paper doesn't scan for tortuosity, but for similar junctions to an initial selection.

4.

MATHEMATICAL TORTUOSITY

The basic definition of Tortuosity may be

T=L/C

L: arch length

C: Chord length

This seems pretty simple for a single arch.

For this definition a straight line has T=1 and a circle T=0.

But for an eye fundus, with several vessels, and professionals usually focusing on the main vessels, putting aside the smaller vessels, and deciding that on experience and human observations, then one understands why there is not a unique definition of mathematical tortuosity

4.1.- E. Grisan, M. Foracchia, A. Ruggeri: A novel method for automatic evaluation of retinal vessel tortuosity. Proceedings of the 25th Annual International Conference of the IEEE EMBS, Cancun, Mexico, 2003

4.2.- M. Mächler: Very smooth nonparametric curve estimation by penalizing change of curvature, Technical Report 71, ETH Zurich, May 1993

For such complicated tracking, Mächler proposes to use as measure of roughness = tortuosity

d/dx(log(T))=d(T)/d(x)/T

instead of T itself.

4.3.- In road building, when designing alternative roads between places, one of the parameters comparing different possible roads is the degree of curvature or tortuosity of each option.

For this definition a straight line has T=Inf and a circle T=R where R is the radius of the circle.

In such context, the curvature of a bent is the Radius of the arch, for each arch.

For a linear series of curves, the total curvature could be defined as the average, or the [min max] interval or whatever custom definition one may consider suitable.

However, the eye fundus has many junctions, and it 's not as simple as just 'following the road' of the main vessels, ignoring the minor ones, as it's common practice.

4.4.- Patasius, M.; Marozas, V.; Lukosevicius, A.; Jegelevicius, D.. Evaluation of tortuosity of eye blood vessels using the integral of square of derivative of curvature // EMBEC'05: proceedings of the 3rd IFMBE European Medical and Biological Engineering Conference, November 20–25, 2005, Prague. - ISSN 1727-1983. - Prague. - 2005, Vol. 11, p. [1-4]

uses

For this definition, both a straight section and a circle have T=0

4.5.- Some may choose to use the Fractal Dimension of a figure as a measure of the tortuosity.

However, as commented by

• Iannaccone, Khokha (1996). Fractal Geometry in Biological Systems. ISBN 978-0-8493-7636-8.
• Vicsek, Tamás (2001). Fluctuations and scaling in biology. Oxford [Oxfordshire]: Oxford University Press. ISBN 0-19-850790-9.

Fractal dimensions .. do not uniquely define patterns.

https://en.wikipedia.org/wiki/Fractal_dimension

Please find in comment 1 a list of references regarding usage of fractal dimensions, but also questioning the validity of such tool in image analysis.

5.-

In

Comparative study of retinal vessel segmentation methods on a new publicly available database

M. Niemeijer, J.J. Staal, B. van Ginneken, M. Loog, M.D. Abramoff

abstract: http://www.isi.uu.nl/Research/Publications/publicationview/id=866.html

compares 5 retinal vessel segmentation methods, 4 automated, and 1 requires supervision.

But there's no actual reference to a mathematical quantification of the amount of tortuosity that a retinal image has.

6.-

In the Simon Fraser University Image Analysis Laboratory

http://mial.fas.sfu.ca/researchProject.php@s=393.html

there is an application to precisely analyse retinal vessels, attached zip livevessel.zip

I haven't tried it but perhaps it does what you ask for.

7.-

In same MIAL website, the MATLAB - ITK is available

http://mial.fas.sfu.ca/researchProject.php@s=308.html

also attached here

matitk_win64.zip

The retinal processing and obtaining of the B&W start images, for further processing has been traditionally done in other languages than MATLAB.

Now this interface claims to harness with MATLAB the C++ libraries used to obtain major vessels.

Please let us know if you decide to use it.

MATITK revision IJ_66_MATITK_revision_02.pdf and source code IJ_66_src.zip available here

http://www.insight-journal.org/browse/publication/66

.

So

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

Comment 1:

Pierre Soille and Jean-F. Rivest (1996). "On the Validity of Fractal Dimension Measurements in Image Analysis" (PDF). Journal of Visual Communication and Image Representation. 7 (3): 217–229. ISSN 1047-3203. doi:10.1006/jvci.1996.0020.

Liu, Jing Z.; Zhang, Lu D.; Yue, Guang H. (2003). "Fractal Dimension in Human Cerebellum Measured by Magnetic Resonance Imaging". Biophysical Journal. 85 (6): 4041–4046. Bibcode:2003BpJ....85.4041L. PMC 1303704Freely accessible. PMID 14645092. doi:10.1016/S0006-3495(03)74817-6.

Smith, T. G.; Lange, G. D.; Marks, W. B. (1996). "Fractal methods and results in cellular morphology — dimensions, lacunarity and multifractals". Journal of Neuroscience Methods. 69 (2): 123–136. PMID 8946315. doi:10.1016/S0165-0270(96)00080-5.

Li, J.; Du, Q.; Sun, C. (2009). "An improved box-counting method for image fractal dimension estimation". Pattern Recognition. 42 (11): 2460–2469. doi:10.1016/j.patcog.2009.03.001.

Cheng, Qiuming (1997). "Multifractal Modeling and Lacunarity Analysis". Mathematical Geology. 29 (7): 919–932. doi:10.1023/A:1022355723781.

Tolle, C. R.; McJunkin, T. R.; Gorsich, D. J. (2003). "Suboptimal minimum cluster volume cover-based method for measuring fractal dimension". IEEE Transactions on Pattern Analysis and Machine Intelligence. 25: 32–41. doi:10.1109/TPAMI.2003.1159944.

Comment 2:

The references of G. Azzopardi N. Petkov paper are worth being checked to find the specifics of the coding

1. Ali, C., Hong, S., Turner, J., Tanenbaum, H., Roysam, B.: Rapid automated tracing and feature extraction from retinal fundus images using direct exploratory algorithms. IEEE Transactions on Information Technology in Biomedicine 3, 125–138 (1999)

2. Bhuiyan, A., Nath, B., Chua, J., Ramamohanarao, K.: Automatic detection of vascular bifurcations and crossovers from color retinal fundus images. In: Third International IEEE Conference on Signal-Image Technologies and Internet-Based System (SITIS), pp. 711–718 (2007)

3. Chanwimaluang, T., Guoliang, F.: An efficient blood vessel detection algorithm for retinal images using local entropy thresholding. In: Proceedings of the 2003 IEEE International Symposium on Circuits and Systems (Cat. No.03CH37430) 5

4. Chapman, N., Dell’omo, G., Sartini, M., Witt, N., Hughes, A., Thom, S., Pedrinelli, R.: Peripheral vascular disease is associated with abnormal arteriolar diameter relationships at bifurcations in the human retina. Clinical Science 103 (2002)

5. Eunhwa, J., Kyungho, H.: Automatic retinal vasculature structure tracing and vascular landmark extraction from human eye image. In: International Conference on Hybrid Information Technology, vol. 7 (2006)

6. Gheorghiu, E., Kingdom, F.: Multiplication in curvature processing. Journal of Vision 9 (2009)

7. Lowe, D.: Distinctive image features from scale-invariant keypoints. International Journal of Computer Vision 60, 91–110 (2004)

8. Martinez-Perez, M., Hughes, A., Stanton, A., Thom, S., Chapman, N., Bharath, A., Parker, K.: Retinal vascular tree morphology: A semi-automatic quantification. IEEE Transactions on Biomedical Engineering 49, 912–917 (2002)

9. Pasupathy, A., Connor, C.: Responses to contour features in macaque area v4. Journal of Neurophysiology 82, 2490–2502 (1999)

10. Pasupathy, A., Connor, C.: Shape representation in area v4: Position-specific tuning for boundary conformation. Journal of Neurophysiology 86, 2505–2519 (2001)

11. Pasupathy, A., Connor, C.: Population coding of shape in area v4. Nature Neuroscience 5, 1332–1338 (2002)

12. Patton, N., Aslam, T., MacGillivray, T., Deary, I., Dhillon, B., Eikelboom, R., Yogesan, K., Constable, I.: Retinal image analysis: Concepts, applications and potential. Progress in Retinal and Eye Research 25, 99–127 (2006) 13. Petkov, N.: Biologically motivated computationally intensive approaches to image pattern-recognition. Future Generation Computer Systems 11, 451–465 (1995)

14. Sherman, T.: On connecting large vessels to small - the meaning of murray law. Journal of General Physiology 78, 431–453 (1981)

15. Staal, J., Abramoff, M., Niemeijer, M., Viergever, M., van Ginneken, B.: Ridgebased vessel segmentation in color images of the retina. IEEE Transactions on Medical Imaging 23, 501–509 (2004)

16. Tsai, C., Stewart, C., Tanenbaum, H., Roysam, B.: Model-based method for improving the accuracy and repeatability of estimating vascular bifurcations and crossovers from retinal fundus images. IEEE Transactions on Information Technology in Biomedicine 8, 122–130 (2004)

17. Tso, M., Jampol, L.: Path-physiology of hypertensive retinopathy. Opthalmology 89 (1982)

[rabeeah]

1. Yes we started witht the same colored image but needed to convert it to bw to proceed further calculations.

2. I dont get this one. Not enough pixels is doing what?

3. So, we use thesw ground truth database as our comparative data set?

4. Yes theres a lot of research available on tortuosity but i'd like to achieve my desired results with the vessels first.

5. My matlab 2015a does not support these c++ libraries, so they dont run. Is there some way I can add it in matlab 2015a windows8 64 bits?

Really appreciate your time and effort. Thank you once again.

[Akande Noah]

Thanks so much Mr John, the code works perfectly well. Please what do i need to add to extract crossover points and bifurcation points from the blood vessels. I can use different colours to distinguish both. I will appreciate your prompt response. Thanks