Asked by mounim
on 9 Dec 2012

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 by Image Analyst
on 9 Dec 2012

Edited by Image Analyst
on 9 Dec 2012

Accepted answer

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');

Show 18 older comments

mounim
on 11 Dec 2012

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
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
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 )

Opportunities for recent engineering grads.

## 0 Comments