How can I create an x,y curve from an edge?

Question

I wanna explain my problem with MATLAB. My objective is to obtain a free liquid surface in a tank from a photo and to do this I have used this algorithm:

A = 'C:\foto\img3.bmp';
B = imread(A, 'bmp');
figure(1), imshow(B);
C = rgb2gray(B);
level = graythresh(C);
bw = im2bw(C,level);
bw = bwareaopen(bw, 50);
figure, imshow(bw);
BW1 = edge(bw,'canny');
figure(2), imshow(BW1);
imwrite(BW1, 'C:\foto\im1_edge.bmp', 'bmp')

Now I have the surface edge in white, with the background in black. Next I have also detected the position of only white pixels:

I= imread('C:\foto\img3_edge.bmp');
imshow(I);
[r c] =size(I);
for j=1:c
  for i=1:r
   if(I(i,j)==1)                  
   [i j]
   end
  end
end

At this point, how can I report (with a macro, automatically possibly) each couple of coordinates on a cartesian (x,y) plane? My objective is to obtain from the edge, so reconstructed, a function of the type "y = f(x)". I have tried with another edge and have modified with paint the image to delete all useless pixels, the example is this one:

with the code:

 I = im2bw(I);

it returns me an error "Warning: The input image is already binary." Next using the code:

 [r c] = find(I), output = [r c];
 plot(r,c,'.')

I obtained this one:

Moreover, when I try to insert r as xdata and c as ydata in cftool, I obtain the same problem and when I use "Interpolant" fitting, it returns an error. Why?

Answer 1

I don't know what type of fit you want to apply to your data. The following code will fit a 3rd order polynomial on your data.

I = imread('20jh1g2.jpg');
I = im2bw(I);
imshow(I);

[r, c] = find(I);
figure;
plot(c,r,'.');
hold on;
f = fit(c, r, 'poly3');
plot((min(c):max(c)),f(min(c):max(c)), 'red', 'LineWidth', 3);

Which will produce:

The rotation can be explained by how the axis are defined. In your image the Y axis goes from 0 at the top to 374 at the bottom. You can transform the result of your fit back into a binary image with the following code;

x = (min(c):max(c))';
y = round(f(x));
I = zeros(size(I));
I(y +((x-1)*size(I,1))) = 1;
figure
imshow(I);

Which will produce:

The result of the fit, f is stored in a cfit object. You can evaluate this function by feeding it values for x, as shown above. The function coefficients can be found by printing the fields of the cfit object in the command window;

f = 

     Linear model Poly3:
     f(x) = p1*x^3 + p2*x^2 + p3*x + p4
     Coefficients (with 95% confidence bounds):
       p1 =  -6.252e-06  (-6.542e-06, -5.963e-06)
       p2 =    0.001753  (0.001588, 0.001918)
       p3 =     -0.3667  (-0.3932, -0.3401)
       p4 =       290.4  (289.3, 291.6)

To flip the function inside your reference frame and compute the centroid you can use this;

I = imread('20jh1g2.jpg');
I = im2bw(I);

[r, c] = find(I);
r = -r + size(I,1);
f = polyfit(c, r, 3);
plot((min(c):max(c)),polyval(f,(min(c):max(c))), 'red', 'LineWidth', 3);
hold on;

xf = [f 0];
fx2 = sym2poly(poly2sym(f)^2);
centroid = 1/polyval(polyint(f),size(I,2)) * [polyval(polyint(xf),size(I,2)) 1/2 * polyval(polyint(fx2),size(I,2))];
plot(centroid(1),centroid(2),'X');

Which will produce: