How to find the inclination of the back bone in the attached image?

Hello Everyone... Greetings...

I am working with some medical imaging problems.

I want to rotate the given image after finding the inclination angle,

So that backbone looks straight and performs further steps.

I have attached the image for the reference.

The image is high resolution with low quality.

Can anyone suggest, how to find the inclination of the backbone?

Thank you.

Answer 1

[Ryan Comeau]

Hello,

I am thinking the best way to do this will be analog. What i mean is just choose two sets of coordinates, compute their slope and apply the rotation yourself. It will be the least painful if you don't need to repeat this frequently.

Doing this autonomously will be quite a bit more difficult. There are ways to get it done, but guarenteeing good performance over your analog work will be difficult. Here is the way to do this analog:

image=imread('path/to/image');
%you choose 2 coordinate sets on the image
coord_1=[x1,y1];
coord_2=[x2,y2];
slope=(coord_1(2)-coord_2(2))/(coord_1(1)-coord_2(1)); %simple slope calculation
%now we convert this to degrees.
%assuming 90 is vertical and 0 is horizontal
slope_deg=atand(slope);
rotation_needed=90-slope_deg;
%perform the rotation
image=imrotate(image,rotation_needed);

Now to do this autonomouly will require some work. You will need to apply an image processing technique to segment the objects present, and identify the spine. You'll need to play around with the methods to see which one works the best. Here are some ideas that I have:

1. Use a hough transform,this algorithm identifies lines in an image(see: https://www.mathworks.com/help/images/ref/hough.html). Locate the longest line(which should be the spine, calculate its slope and perform the rotation). Try processing without the image pre processing ans see if you can detect the spine.

%code for the Hough transform
%this block gives the endpoints of the lines it locates. 
%try the Hough transform without the pre processing as well.
image=imread('path/to/image/')
image=rescale(image,0,1);%helps in binarization
BW_im=imbinarize(image);
BW_finel=bwskel(BW_im); %create image skeleton
[H,theta,rho] = hough(BW_final);
P = houghpeaks(H,45,'threshold',ceil(0.1*max(H(:))));
lines = houghlines(BW_im,theta,rho,P,'FillGap',5,'MinLength',7);
%lines is a struct, stick this in a for loop to extract
for i=1:length(lines)
    xy(i,:)=[lines(i).point1,lines(i).point2];
end
%parse for the longest one to identify the spine.

2. Use Haar like features

All of these methods except 4, will require you to parse the output and figure out the slope of the spine, then apply the rotation. A Haar like feature(just google these), can be used to identify objects. Here, we will construct one for the spine, and rotate it in intervals of 5 degrees(you can change this). The Haar like features will convolute the image, and whichever angle yields the highest value(from the convolution) will give a +-2.5 degrees on your spine angle. The following block of code explains how to use a Haar like feature and the logic to determine the rotation.

image=imread('path/to/image');
image=rescale(image,0,1); %helps with speed and other things
spine_average_width=30; %this is in pixels, change for your need(tune it)
spine_average_height=100;% in pixels change for your need(tune it)
spine_haar_feature(1:spine_average_height,1:2*spine_average_width)=zeros();
spine_haar_feature(1:spine_average_height,...
    spine_average_width/2:(3/2)*spine_average_width)=ones();
%we have created a binari image with a vertical bar. this is our Haar like feature
%you'll need to add a little routine to make sure the rotation is upwards
for i=1:18 %starting from vertical going horizontal
    spine_angle(i)=max(conv2(image,imrotate(spine_haar_feature,-i.*5)));
end
[~,index]=max(spine_angle(:));
angle=index*5; %because 5 degrees per iteration
image=imrotate(image,angle);

I hope this answer is not overbearing, i just like medical imaging problems and think they're fun to solve.

Hope this helps,

RC

[Mrutyunjaya Hiremath]

Hello Ryan Comeau,

Thanks for your time and answer.

I love the way of your coding and content presentation. Really it's nice.

I've worked out the same manual method (like, your first method).

close all;

[fileName folderName] = uigetfile({'*.jpg'}, 'Select an Image for Classification', '..\Dataset\Test');

imgTest = imread([folderName fileName]);
[~,name,~]= fileparts(fileName);
figure,imshow(imgTest);title('Query Image');


N = 5;
pts = readPoints(imgInput, N);

% fit straight line
x = pts(1,:); y = pts(2,:);
coeffs = polyfit(x, y, 1);
slope = coeffs(1);
intercept = coeffs(2);

% Get fitted values
fittedX = linspace(min(x), max(x), 200);
fittedY = polyval(coeffs, fittedX);
hold on;
plot(fittedX, fittedY, 'b-', 'LineWidth', 2);

% slope and angle of line
x1 = fittedX(1);
y1 = fittedY(1);
x2 = fittedX(end);
y2 = fittedY(end);
slope = (y2 - y1) ./ (x2 - x1); % just for better accuracy
angle = atand(slope);

% rotate image
if (angle < 0)
    angle = -angle;
    imgRotate = imrotate(imgInput, 90-angle, 'bicubic');
else
    imgRotate = imrotate(imgInput, -(90-angle), 'bicubic');
end
figure, imshow(imgRotate);

% get center part of image
X1 = size(imgRotate,1)*cosd(angle);
Y1 = size(imgRotate,2)*cosd(angle);
W = size(imgRotate,1) - 2*Y1; 
H = size(imgRotate,2) - 2*X1; 
imgCrop = imcrop(imgRotate, [X1, Y1, H, W]);
figure, imshow(imgCrop);

readPoints() is a function used in above code, from https://in.mathworks.com/matlabcentral/answers/118724-how-to-get-onclick-coordinate-pixel-value-and-location-from-an-image

Got good reuslts. like

But, the automotive method is good (the reason is, there are 100's of images :).

let me try your 3rd and 4th method.

Once again thank you, for the answer.

[Ryan Comeau]

I've updated the answer with some more Hough transform details.

Answer 2

[Image Analyst]

I disagree with the accepted answer. Perhaps he's never heard of the radon transform - the Nobel prize winning medical imaging algorithm that is at the heart of CT and MRI. You can do this automatically just by taking the radon transform and looking for the peak.

See my attached demo. Adapt as needed by applying the FAQ: FAQ How_can_I_process_a_sequence_of_files? for your hundreds of images.

[Ryan Comeau]

I mentioned the Hough transform (bullet point 1) which performs a very similar function to that of the radon transform. Here is a forum discussing the two: https://dsp.stackexchange.com/questions/470/whats-the-difference-between-the-hough-and-radon-transforms. I will fill in a more elaborate answer for the Hough transform.

[Image Analyst]

Thank you - that is a helpful comparison.

[Mrutyunjaya Hiremath]

Hello Ryan and Image Analyst,

Really it's a good discussion and explanation.

From our side, elderly people say, ‘we lose nothing, even when we are quarreling with good and knowledgeable people’. That’s true. J.

Learned a lot from this.

Thank you both of you.

[Image Analyst]

Girish, start a new question - let's not hijack this one. Post your image there and a definition of what direction you consider to be the "height".

Until then, see John D'Errico's outstanding File Exchange: https://www.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects