Help with the bwtraceboundary function

3 views (last 30 days)
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);

Accepted Answer

Image Analyst
Image Analyst on 9 Dec 2012
Edited: Image Analyst on 9 Dec 2012
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');
  21 Comments
Image Analyst
Image Analyst on 11 Dec 2012
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
mounim on 12 Dec 2012
Edited: mounim on 12 Dec 2012
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 )

Sign in to comment.

More Answers (1)

mounim
mounim on 14 Dec 2012
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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!