Help with the bwtraceboundary function

hey

i am a totally new beginner in matlab, and i hope someone can help me with this code, what i want to do is finding the intersection point between the edge of the bone and the horizontal line, i tried to modify the code in the demo section and i got it to work for the left side of the bone, i want to do the same thing for the other side of the bone as well.

please how can i achieve that ?

( i attached my code, and the output image )

      %%%%%%%%LOADING THE IMAGE %%%%%%%%%

RGB = imread('1.jpg');

%%%%%%%%EXTRACT THE ROI %%%%%%%%%%%%%%%

% you can obtain the coordinates of the rectangular region using
% pixel information displayed by imtool
start_row = 1;
start_col = 1;

start_row2 = 1100;
start_col2 = 1;

cropRGB = RGB(start_row:115, start_col:1210, :);
cropRGB2 = RGB(start_row2:1265, start_col2:1210, :);

% Store (X,Y) offsets for later use; subtract 1 so that each offset will
% correspond to the last pixel before the region of interest
offsetX = start_col-1;
offsetY = start_row-1;

offsetX2 = start_col2-1;
offsetY2 = start_row2-1;

%%%%%THRESHOLD THE IMAGE %%%%%%%%%

I = rgb2gray(cropRGB);
II = rgb2gray(cropRGB2);

threshold = graythresh(I);
BW = im2bw(I,threshold);
%BW = ~BW;  % complement the image (objects of interest must be white)
imshow(BW)

threshold2 = graythresh(II);
BW2 = im2bw(II,threshold2);
%BW = ~BW;  % complement the image (objects of interest must be white)
imshow(BW)

%%%%%%%%%%INITIAL POINT OF BOUNDERY %%%%%%%%%%%%%%%

dim = size(BW);
dim2 = size(BW2);

% horizontal beam
col1 = 420;
row1 = min(find(BW(:,col1)));

col11 = 420;
row11 = min(find(BW2(:,col11)));

% angled beam
row2 = 12;
col2 = min(find(BW(row2,:)))

row22 = 60;
col22 = min(find(BW2(row22,:)))

%%%%%%%%%TRACING THE BOUNDERY %%%%%%%%%%%

boundary1 = bwtraceboundary(BW, [row1, col1], 'N', 8, 70);

boundary11 = bwtraceboundary(BW2, [row11, col11], 'N', 8, 70);

% set the search direction to counterclockwise, in order to trace downward.
boundary2 = bwtraceboundary(BW, [row2, col2], 'E', 8, 90,'counter');

boundary22 = bwtraceboundary(BW2, [row22, col22], 'E', 8, 90,'counter');

%imshow(RGB); hold on;

% apply offsets in order to draw in the original image
plot(offsetX+boundary1(:,2),offsetY+boundary1(:,1),'g','LineWidth',2);
plot(offsetX+boundary2(:,2),offsetY+boundary2(:,1),'g','LineWidth',2);

plot(offsetX2+boundary11(:,2),offsetY2+boundary11(:,1),'g','LineWidth',2);
plot(offsetX2+boundary22(:,2),offsetY2+boundary22(:,1),'g','LineWidth',2);

%%%%%%%%%%FIT LINES TO THE BOUNDERY %%%%%%%

ab1 = polyfit(boundary1(:,2), boundary1(:,1), 1);
ab2 = polyfit(boundary2(:,2), boundary2(:,1), 1);

ab11 = polyfit(boundary11(:,2), boundary11(:,1), 1);
ab22 = polyfit(boundary22(:,2), boundary22(:,1), 1);

%%%%%%%%%FIND THE ANGLE OF INTERSECTION %%%%%%

vect1 = [1 ab1(1)]; % create a vector based on the line equation
vect2 = [1 ab2(1)];
dp = dot(vect1, vect2);

vect11 = [1 ab11(1)]; % create a vector based on the line equation
vect22 = [1 ab22(1)];
dp2 = dot(vect11, vect22);

% compute vector lengths
length1 = sqrt(sum(vect1.^2));
length2 = sqrt(sum(vect2.^2));

% compute vector lengths
length11 = sqrt(sum(vect11.^2));
length22 = sqrt(sum(vect22.^2));

% obtain the larger angle of intersection in degrees
angle = 180-acos(dp/(length1*length2))*180/pi

angle1 = 180-acos(dp2/(length11*length22))*180/pi

%%%%%%%%FIND POINT OF INTERSECTION

intersection = [1 ,-ab1(1); 1, -ab2(1)] \ [ab1(2); ab2(2)];
% apply offsets in order to compute the location in the original,
% i.e. not cropped, image.
intersection = intersection + [offsetY; offsetX]

intersection1 = [1 ,-ab11(1); 1, -ab22(1)] \ [ab11(2); ab22(2)];
% apply offsets in order to compute the location in the original,
% i.e. not cropped, image.
intersection1 = intersection1 + [offsetY2; offsetX2]

%%%%%%%%%%%%PLOT THE RESULTS %%%%%%%%%

inter_x = intersection(2);
inter_y = intersection(1);

inter_x1 = intersection1(2);
inter_y1 = intersection1(1);

% draw an "X" at the point of intersection
plot(inter_x,inter_y,'rx','LineWidth',2);

text(inter_x-100, inter_y-30, [sprintf('%1.3f',angle),'{\circ}'],...
     'Color','r','FontSize',14,'FontWeight','bold');

%interString = sprintf('(%2.1f,%2.1f)', inter_x, inter_y);

plot(inter_x1,inter_y1,'rx','LineWidth',2);

text(inter_x1-100, inter_y1-30, [sprintf('%1.3f',angle1),'{\circ}'],...
     'Color','r','FontSize',14,'FontWeight','bold');

%interString = sprintf('(%2.1f,%2.1f)', inter_x1, inter_y1);

Answer 1

[Image Analyst]

Don't use bwtraceboundary() - it's almost never necessary. I think you're making this way way too complicated. Just scan down the left and right columns to find the line, then find the first and last point on the row where it's white. First convert the image to binary.

rightColumn = binaryImage(:, end); % Get right column
topLine = find(rightColumn, 1, 'first'); % Find row where line is
% Extract the row above
oneRow = binaryImage(topLine-1, :);
% Find first and last places - the bone edges
leftEdge1 = find(oneRow, 1, 'first');
rightEdge1 = find(oneRow, 1, 'last');

bottomLine = find(rightColumn, 1, 'last'); % Find row where line is
% Extract the row above
oneRow = binaryImage(bottomLine -3, :); % 3 (or whatever) lines up.
% Find first and last places - the bone edges
leftEdge2 = find(oneRow, 1, 'first');
rightEdge2 = find(oneRow, 1, 'last');

[mounim]

thanks for your answer, but i do have a question though, i could see at the documentation that the find(x) function returns nonzero elements of array X, so if i understood it right should

leftEdge1 = x value for the left intersection point and topLine = y value for the left and right intersection point

but when plotting it like this :

plot(leftEdge1,topLine,'rx','LineWidth',2);
plot(rightEdge1,topLine,'rx','LineWidth',2);

plot(leftEdge2,bottomLine,'rx','LineWidth',2);
plot(rightEdge2,bottomLine,'rx','LineWidth',2);

Only the first plot is plotted, what i am doing wrong ?

and another question, let's say i have 4 horizontal lines ( my goal is finding the axis of the bone )

topLine = find(rightColumn, 1, 'first');

will return the first line and 'last' will return the last one, what about the lines in between ?

[Image Analyst]

Here is the complete code. Note that I had to go up an extra line because the line is not straight - perhaps to to jepg artifacts on your uploaded image. Also, I just found the intersections with the top of the thick horizontal lines but it is a straightforward adaptation to have it do it again for underneath the thick horizontal lines.

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.
format longg;
format compact;
fontSize = 20;

% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Users\Mounim\Documents';
baseFileName = 'bone.jpg';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- 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 in the search path folders.', 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)
% Convert to grayscale
if numberOfColorBands > 1
  grayImage = grayImage(:, :, 2);
end
% Crop away white boundary
grayImage = grayImage(20:705, 56:743);
% Binarize the image
binaryImage = grayImage > 128;
% 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]);
% Give a name to the title bar.
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off') 

rightColumn = binaryImage(:, end); % Get right column
topLine = find(rightColumn, 1, 'first'); % Find row where line is
% Find the first black row above it.  Note there is a half row 
% on top due to jpg artifact that may not be present in the original image.
topLine = topLine - 2;
% Extract the row above
oneRow = binaryImage(topLine, :);
% Find first and last places - the bone edges
leftEdge1 = find(oneRow, 1, 'first');
rightEdge1 = find(oneRow, 1, 'last');
hold on;
plot(leftEdge1, topLine, 'ro');
plot(rightEdge1, topLine, 'ro');

bottomLine = find(rightColumn, 1, 'last'); % Find row where line is
% Extract the row 6 rows above because the line is 4 rows thick
% but there is an extra half row do to jpg artifacts.
bottomLine = bottomLine - 6; % Actually the black line above the bottom thick line.
oneRow = binaryImage(bottomLine, :);
% Find first and last places - the bone edges
leftEdge2 = find(oneRow, 1, 'first');
rightEdge2 = find(oneRow, 1, 'last');
plot(leftEdge2, bottomLine, 'ro');
plot(rightEdge2, bottomLine, 'ro');

[mounim]

Thanks a lot for your great help really appreciate it, one last question.

now i added multiple lines ( 2 at the top, and to at the bottom )

the intersections with the first line and second line work just great, when i try to mark the intersection with the 2 next lines, it doesnt work, here what i tried :

topLine = find(rightColumn, 1, 'first');
topLine = topLine - 1;
nextLine = topline(2) - 1
.
.
plot(leftEdge1, topLine, 'ro');
plot(rightEdge1, topLine, 'ro');

plot(leftEdge1, nextLine, 'ro');
plot(rightEdge1, nextLine, 'ro');

how should i do it propely, when doing that, i still can only show one row.

[Image Analyst]

Question: is that thick line always the same width (4 pixels)?

[mounim]

no now i changed it ... it is only one pixel

[Image Analyst]

topline is not a 2 element array - it's a single number.

nextLine = topline + 1;

[mounim]

hmm .. then i need to understand, because :

topLine = 92 142 472 877 1268

and nextline = 93 143 473 878 1269

that's why i wrote nextLine = topline(2) - 1 before

the next line = the next white horizontal line ( i have 4 ), sorry

[Image Analyst]

Then you changed something without telling me. Because this line:

topLine = find(rightColumn, 1, 'first');

will return only ONE value. See the second argument to find() - the 1? That means return only one value, not 5 as you have it.

[mounim]

i am sorry, because i am trying to understand the code and test it, before coming here and asking questions, but yeah you are right

i changed it to this :

topLine = find(rightColumn);

because i wanted to get the y-coordinates for every white line, which did, but when i tried

nextLine = topline(2) - 1
plot(leftEdge1, nextLine, 'ro');

it does not work.

[Image Analyst]

ALright, let's settle on some known characteristics.

1) Will you have a fixed number of lines going across the image (like 4 or 1), or will it be variable and unknown?

2) Will there always be a half-line in that strip going across the image, or was that just there (in the one I downloaded) because your image got converted to jpeg upon uploading to the web site, and it's just a jpeg artifact and not really there in your original image? If so, can you upload a perfect PNG image with no jpeg artifacts?

[mounim]

Ok .. thanks for your patience, what i basically want to do is to find the angle between the top bone (femur) and the bottom bone(tibia).

The most accurate way to do it - i think - is to :

• draw two horizontal lines in the top of the bone
• find the intersection points between the bone and the horisontal lines
• find the center of each segment ( the top segment and the second )
• draw a line between the two centers - now i've got the axis of the first bone
• repeat the same procedure for bottom bone, and finde the second axis
• simple finde the angle between the 2 axis

i need to do that for 6000 pictures, so i just wanted to segment my question to get some help.

i will upload a picture shoing my idea, and another PNG picture

[Image Analyst]

OK - so you can just ignore practically everything we did up until now. The horizontal lines are not really there and don't need to be found at the edges and then scoot along them until they hit the bone. So can we say this: the upper bone will have a strong signal at the outer edge and no other white pixels outside of that, and that you will have two bones on the bottom?

[mounim]

OK then ... that's right, the horizontal line are not there ( i just add them my self ). No other white pixels will be outside the edge of the bone ( i will. at the bottom, yes i will have two bones, but i am basically only interested on getting the axis of the the thicker bone in the bottom ( and of course the one in the top ).

[Image Analyst]

What I would do is to look in the range from row1 to row2. Take one row and use find() to find the first and last white pixel. Then take the average to get the middle column. Get all those middle column values in the range and then pass them into polyfit to get the line that fits the centerline in that range.

Then look in the next range, from row3 to row4. You might want to extract those rows into a new subimage and then skeletonize. Then take each row and use find() without using 'first' and 'last' so that you should get 4 values this time. The first two are the left bone and the next two are the right bone. Then, like before, find the centerline and call polyfit to get the average centerline. Can you code that up - because I don't have time to do it for you tonight?

[mounim]

Thanks for your directions .. to be honest, it was difficult for me to do it the way you discribed it, but i did it another way and it works!

here is the code :

row = 50;
row2 = 100;

oneRow = binaryImage(row, :); 
oneRow2 = binaryImage(row2, :); 

leftEdge1 = find(oneRow, 1, 'first'); 
rightEdge1 = find(oneRow, 1, 'last');

leftEdge2 = find(oneRow2, 1, 'first'); 
rightEdge2 = find(oneRow2, 1, 'last');

Center1=[((leftEdge1+rightEdge1)/2) ((leftEdge2+rightEdge2)/2)];
Center2=[row row2];

hold on;

plot(leftEdge1, row, 'rx','LineWidth',2);
plot(rightEdge1, row, 'rx','LineWidth',2);

plot(leftEdge2, row2, 'rx','LineWidth',2);
plot(rightEdge2, row2, 'rx','LineWidth',2);

plot([Center1(1), Center1(2)],[Center2(1), Center2(2)],'g','LineWidth',2);

%%%%BUTTOM BONE %%%%

row3 = 1200;
row4 = 1100;

oneRow3 = binaryImage(row3, :); 
oneRow4 = binaryImage(row4, :); 

leftEdge3 = find(oneRow3); 
rightEdge3 = find(oneRow3);

leftEdge4 = find(oneRow4); 
rightEdge4 = find(oneRow4);

Center3=[((leftEdge3(5)+rightEdge3(6))/2) ((leftEdge4(5)+rightEdge4(9))/2)];
Center4=[row3 row4];

plot(leftEdge3(5), row3, 'rx','LineWidth',2);
plot(rightEdge3(6), row3, 'rx','LineWidth',2);

plot(leftEdge4(5), row4, 'rx','LineWidth',2);
plot(rightEdge4(9), row4, 'rx','LineWidth',2);

plot([Center3(1), Center3(2)],[Center4(1), Center4(2)],'g','LineWidth',2);

- the axis are drawn as vectors, but i couldnt find out how to elongate them so they are lines, because i have to calculate the angle between those 2 lines. - and i guess i have to clean my pictures much more, so i only end with the border and absolutly no other pixels outside.

[Image Analyst]

OK. Not as robust as the method I recommended since it's more sensitive to the particular lines you picked, but whatever..... Perhaps someday you'll have a greater need to learn and use polyfit().

[mounim]

All right, i tried to follow your advice and do exactly what you suggest ( i hope). the code is working, But it not showing the centerline ( polyfit), what i am doing wrong ?

% for loop to define the range of interest for the bone 1
for i=50:300

   oneRow = cell(50,1);
   leftEdge = cell(50,1);
   rightEdge = cell(50,1);

   oneRow{i} = binaryImage(i, :);

   leftEdge{i} = find(oneRow{i}, 1, 'first');
   rightEdge{i} = find(oneRow{i}, 1, 'last');

   % Write the Polyfit Function
   P = polyfit(((rightEdge{i} + leftEdge{i})/2), i, 1); 
   r = polyval(P, ((rightEdge{i} + leftEdge{i})/2));

   hold on;

   %plot the left and right edge
   plot(leftEdge{i}, i, 'rx','LineWidth',2);
   plot(rightEdge{i}, i, 'rx','LineWidth',2);

   % CENTER POINTS 
   plot(((rightEdge{i} + leftEdge{i})/2), i, 'rx','LineWidth',2);

   % Plot the Polyfit
   plot(((rightEdge{i} + leftEdge{i})/2), r);

 end

and one more question, can i calculate the angle between the polyfit line for the topBone, and the one for the bottomBone ?

[Image Analyst]

Not quite right. First of all oneRow, leftEdge, and rightEdge don't need to be cell arrays - just regular numbers is fine. Don't use i (the imaginary variable as a loop index. Use k. Third, you need to compute

midpoints(k) = (leftEdge + rightEdge)/2;

Then don't call polyfit() until that loop has exited and all the midpoints (the midpoint on each row) have been calculated. Then do

coeffs = polyfit(midpoints, 1);

The first coeff, coeffs(1), is the slope. And of course the angle in degrees is just atand(slope).

[mounim]

i am kind of confused, that's the closest i got :

for k=50:100

oneRow{k} = binaryImage(k, :);

leftEdge{k} = find(oneRow{k}, 1, 'first');
rightEdge{k} = find(oneRow{k}, 1, 'last');

hold on;

midpoints(k) = (leftEdge{k} + rightEdge{k})/2;

plot(leftEdge{k}, k, 'rx','LineWidth',2);
plot(rightEdge{k}, k, 'rx','LineWidth',2);

end
Y=[50:100];
coeffs = polyfit(midpoints,Y,1);

i am sorry i dont think i am getting how to fit the polyfit in. i tried looking at the documentation. but i dont know.

[Image Analyst]

That looks about right except that you're still using cell arrays for oneRow and leftEdge and rightEdge when you don't have to, but it doesn't hurt anything other than making it more coplicated. Just get rid of the {k} from all of them. Other than that, what's the problem? It looks like it should work once you flip Y and midpoints, since Y is your "x" or independent axis in this case.

coeffs = polyfit(Y, midpoints, 1);

[mounim]

Thanks again, well i just cannot draw that line, what i am doing wrong ?. Here is my code :

for k=50:100

oneRow = binaryImage(k, :);

leftEdge = find(oneRow, 1, 'first');
rightEdge = find(oneRow, 1, 'last');

hold on;

midpoints(k) = (leftEdge + rightEdge)/2;

plot(leftEdge, k, 'rx','LineWidth',2);
plot(rightEdge, k, 'rx','LineWidth',2);

end

X=[1:100];
coeffs = polyfit(X,midpoints, 1);
Y = coeffs(1) .* X + coeffs(2);
plot(X,Y, '-','LineWidth',2)

And the output in matlab ( i used a white rectangle, to be sure that things are clean )

Answer 2

[mounim]

It is okay, i figured that out !! ... i had to remove the zeros - from the array in my for-loop. so now it is working.

thanks for your help