Prev: L moment for Gruped/Classed data
Next: GLMfit run time
From: astro mmi on 4 Mar 2010 12:45
Hi everyone,

I have written the following code for boundary detection. It is using the image processing toolbox for angle measurement. I am trying to measure the angle between two beams/cylindrical surfaces(vertebra image). Iam not getting the detected boundary as a line. I understand that the problem is with defining the initial points. Pls help me identify where I am making an error with this regard. Thanx in advance.
clc;
clear all;
close all;

i=imread('spine.jpg');
figure,imshow(i)
title('original image')

%imtool
start_row = 1;
start_col = 196;
cropRGB = i(start_row:445, start_col:334, :);
imshow(cropRGB)
title('cropped image')

j= imrotate(cropRGB,-60,'nearest','loose');
figure,imshow(j)
title('rotated image')
hold on
h = imshow(cropRGB, gray(256));
set(h, 'AlphaData', 0.9)
%This figure is saved as rotangle.jpg

offsetX = start_col-1;
offsetY = start_row-1;

rangle=imread('rotangle.jpg');
im =adapthisteq(rgb2gray(rangle));
figure,imshow(im)

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

dim = size(BW);
% vertical spine
col1 =192;
row1 = min(find(BW(:,col1)));

% angled spine
col2 = 300;
row2 = min(find(BW(:,col2)));

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

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

figure,imshow(im); 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);

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

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

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

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

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]

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

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

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

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

text(inter_x-10, inter_y+20, interString,...
'Color','y','FontSize',14,'FontWeight','bold');

 | 
Pages: 1
Prev: L moment for Gruped/Classed data
Next: GLMfit run time