Find degree of branching using circular neighbourhood method matlab

I have a skeleton image. as shown on following link.

http://www.flickr.com/photos/92388309@N03/8397759970/in/photostream/

I have detected branchpoints and endpoints and labelled the image. Now I am drawing circle on each branch and I am using following code to draw a circle on each branch.

mn=bwmorph(y,'branchpoints');
[row column] = find(mn);
branchPts    = [row column];
endImg    = bwmorph(y, 'endpoints');
[row column] = find(endImg);
endPts       = [row column];
figure;imshow(y);
hold on ; 
plot(branchPts(:,2),branchPts(:,1),'rx');
hold on; plot(endPts(:,2),endPts(:,1),'*');
% Labeling the Branches
branches = (y & ~mn); % set branch points to zero
figure; imshow(branches);
branchesLabeled = bwlabel( branches); % label connected components
vislabels(branchesLabeled)
% Calculation of Length of Branches and Circular Neighbourhood Method for
% for detection of normal and abnormal branches.
sts = regionprops( branchesLabeled,'Area', 'Perimeter','MajorAxisLength','Centroid' );
figure
imshow(branchesLabeled)
hold on
% extract properties
% Loop for circles
for i=1:size(sts)
r= sts(i).MajorAxisLength/2 ; %desired radius
centerx = sts(i).Centroid(1);  
centery = sts(i).Centroid(2);
th = 0:pi/50:2*pi;
xunit = r * cos(th) + centerx;
yunit = r * sin(th) + centery;
plot(xunit, yunit);
end

I have couple of problems with this code.

The circles are drawn on center line of skeleton(It is not a branch it is center line of skeleton). I want circles only on branches and not on skeleton center line.

http://www.flickr.com/photos/92388309@N03/8403454589/in/photostream

How to detect only branches in skeleton? I am storing branchpoints in mn variable in above code and to label it i am inverting those branchpoints with my skeleton image

branches = (y & ~mn);

Because of this my center line is breaking and code is considering non connected center line as branch and drawing circle on it. Is there a way to label the branch points and drawing circle on branches only?

[Kirti]

Can someone help me in above problem? It is urgent and I need solution badly.

[Kirti]

Can someone attend this question? I am still searching for solution?

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/29922-why-your-question-is-not-urgent-or-an-emergency

Answer 1

[Image Analyst]

I'm not really sure what you want to do. You have a skeleton and then identify branchpoints and break the skeleton apart into individual brand segments. Then you place a circle at the centroid of each branch which has a diameter that is approximately the length of the branch. It's doing what you told it to do. If you want circles only at the "branchpoints" only, then don't loop over all branches placing circles at the centroid of each branch, loop over branchpoints (not branches) and place circles at the branchpoint locations (not the branch centroids).

[Kirti]

I want circles on branches only. Not on skeleton center line which is indicated by arrows in following link.

http://www.flickr.com/photos/92388309@N03/8403454589/in/photostream

When I execute above code it gives me circle on center line too.

http://www.flickr.com/photos/92388309@N03/8403454589/lightbox/

I also want to know co-ordinates of each branch. By drawing circle on each branch I want to check for normal and abnormal branch in MRI image(based on length of branch) and also want to check order of branching(1st order,2nd order). What changes I need to do in above code to achieve this?

[Image Analyst]

I can't see - can you plot in a color that you can see like red or yellow and repost? A larger image would also help, the image is so tiny (at least on my hi-res wide screen monitor) - pick a different web site if necessary. Blue on black is very hard to see. Are you saying that you want circles only on the branches that are also considered endpoints, and not on branches that have connections on both ends?

[Kirti]

Sorry. My image with circle is ROI of my skeleton. I have plotted circles in red and indicated circles which are on center line with yellow arrows in following link. These are circles plotted on center line and not on the branch.

https://plus.google.com/photos/111799079095235591256/albums/5837842199082357233/5837842451151872290?hl=en_US "Are you saying that you want circles only on the branches that are also considered endpoints, and not on branches that have connections on both ends?" yes, that's what I want

[Image Analyst]

Then you need to use imreconstruct(). Have the marker image be the endpoints of the original skeleton, so it will have only the tips of branches. Then have the other image be the image with branchpoints removed. So only branches which also have an endpoint will be reconstructed.

[Kirti]

Thanks for the reply. Can you please explain it with code.

[Image Analyst]

clc;
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
fontSize = 20;

% Change the current folder to the folder of this m-file.
if(~isdeployed)
  cd(fileparts(which(mfilename)));
end
  
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

grayImage = peaks(256);
% Display the original image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Image for deriving a binary image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Position', get(0,'Screensize')); 

% Get a binary image.
binaryImage = grayImage > 3;
% Display the image.
subplot(2, 2, 2);
imshow(binaryImage, []);
title('The Starting Binary Image', 'FontSize', fontSize);

% Erode it to form a marker image.
% Erode it wnough to have only 2 blobs instead of 3
markerImage = imerode(binaryImage, ones(15));
% Display the image.
subplot(2, 2, 3);
imshow(markerImage, []);
title('Marker Image = Eroded Binary Image', 'FontSize', fontSize);

reconst = imreconstruct(markerImage, binaryImage);
% Display the reconstructed image.
subplot(2, 2, 4);
imshow(reconst, []);
title('Final Reconstructed Image', 'FontSize', fontSize);

[Image Analyst]

Second demo:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
fontSize = 20;

% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

% Read in a standard MATLAB gray scale demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'eight.tif';
fullFileName = fullfile(folder, baseFileName);
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % Didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% Display the original gray scale image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]); % Maximize figure.
set(gcf,'name','Image Analysis Demo','numbertitle','off') 

% Let's compute and display the histogram.
[pixelCount grayLevels] = imhist(grayImage);
subplot(2, 2, 2); 
bar(pixelCount);
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

maskImage = imfill(grayImage < 220, 'holes');
% Display the image.
subplot(2, 2, 2);
imshow(maskImage, []);
title('Mask Image (Original Image Thresholded)', 'FontSize', fontSize);

markerImage = false(size(maskImage));
markerImage(20:70, 225:230) = true;
% Display the image.
subplot(2, 2, 3);
imshow(markerImage, []);
title('Marker Image', 'FontSize', fontSize);

finalImage = imreconstruct(markerImage, maskImage);
% Display the image.
subplot(2, 2, 4);
imshow(finalImage, []);
title('Final Reconstructed Image', 'FontSize', fontSize);