0

I am trying to get a Linear best fit line through certain points in an image but am unable to get a line that would seem logical. I am not sure if I am doing something wrong here, or if that really is the best fit line in this case.

This is the original image:

Original Image

This is the result after selecting an area of interest, filtering the edge, and then plotting the best fit line through the pixels of that edge.

Results

I expected the best fit line (blue line in the results) to be almost vertical. However, as shown in the image, it is tilted at an angle while the red points are all close to a vertical and not scattered around the image. How can I correctly detect the vertical line?

Code:

image = imread('https://i.stack.imgur.com/rQ1SQ.png');
image = rgb2gray(image(85:270, 165:210, :));
filter = [ones(1, 6);zeros(1,6);-ones(1,6)];
filter = filter.';
filtered =imfilter(image,filter);
regions = regionprops(bwlabel(filtered,8), 'Area', 'PixelList');
[maxarea, maxindex] = max([regions.Area]);
p = polyfit(regions(maxindex).PixelList(:,1), regions(maxindex).PixelList(:,2), 1);
x = 1:size(image,2);
y = polyval(p,x);
subplot(1,3,1), imshow(image, 'border','tight');
hold on;
plot(x,y,'b');
title('Area of Interest in Original Image')
subplot(1,3,2), imshow(filtered);
hold on;
plot(x,y,'b');
title('Vertical Edge after Filtering')
subplot(1,3,3), imshow(fim);
hold on;
plot(x,y,'b');
plot(regions(maxindex).PixelList(:,1), regions(maxindex).PixelList(:,2), 'r.');
title('Image with points plotted')

Code Reference: Shai's answer on a related question.

Adriaan
  • 17,741
  • 7
  • 42
  • 75
Usama
  • 1
  • 1

2 Answers2

0

First a few comments. First you clearly have a vertical line here, but polyfit isn't fitting a vertical line. Why? because you are using a least squares fit on points with a high aspect ratio (4 px wide, 168 px high). One way to work around it is to rotate your data, fit a line, then rotate it back. 

% rotate imroi points by 45 degrees
X = regions(maxindex).PixelList(:,1);
Y = regions(maxindex).PixelList(:,2);
N = size(X,1);
a = 45; %rotation angle
A = [cosd(a), sind(a); -sind(a), cosd(a)];
Z = (A * [X,Y]')'; %rotate points

figure;
subplot(1,2,1)
% plot rotated points
plot(Z(:,1), Z(:,2), 'r.')
hold on;
XRot = ones(N, 2); XRot(:,1) = Z(:,1);
b = XRot \ Z(:,2);
z = [linspace(min(Z(:,1)),max(Z(:,1)),N)', ones(N,1)];
plot(z(:,1), z*b, 'b-');
axis ij
title('Rotated imroi points')

subplot(1,2,2)
% plot un-rotated points
a = -45; %rotation angle
A = [cosd(a), sind(a); -sind(a), cosd(a)];
imshow(filtered);
hold on;
plot(X, Y, 'r.')
W = (A * [z(:,1) ,z * b]')';
plot(W(:,1), W(:,2), 'b-')
title('Image with points plotted')

I'd say in general your best bet for this kind of work is a computer vision package like vl_feat, which has a nice MATLAB (MEX) wrapper.

peterfranciscook
  • 311
  • 2
  • 7
  • Thank you for your answer @peter. You're absolutely right about the high aspect ratio being the reason for this particular problem, and your answer seems to solve that issue. While digging around, I found another approach which I think is easier to implement and uses an Orthogonal Least Squares approach instead of the Vertical Least Squares approach. I will definitely look into the computer vision package you have suggested here. – Usama Aug 05 '19 at 20:47
0

As pointed out by Peter in his answer, the issue in this problem is that the aspect ratio of the data points is very high: they are almost a vertical line already. The default Least Squares Method for a best fit line does not work well in this case. This can be resolved by using Orthogonal Least Squares method and to better understand that concept, have a look here which I found from Andrey's answer to a similar problem.

Luckily, there is a package developed by F.Carr which makes it extremely easy to implement the Orthogonal Linear Regression method.

The results by using this method are shown here.

This seems to have solved the problem for now although I haven't tried it on a large number of cases yet.

Usama
  • 1
  • 1