How to find the upper and lower curve from this binary image?

I first segmented the OCT scan using thresholding. After thresholding, I want to extract the upper and lower curves as marked with red and red line on the right image.

[Walter Roberson]

Does it need to follow the bright lines, or could it instead be the outside edges?

Answer 1

[Mathieu NOE]

hi

I really do not consider myself as an image processing guy but I wanted to try it

so here we go .... see result below based on image provided (a bit cropped to have only the "valid" area, see attachment)

A = imread('image.png');
A = double(A);
A = flipud(A);
[m,n] = size(A);
[y,x] = find(A>0.5);

[x3,y3,x4,y4] = top_bottom_boundary(x,y);

%% top line
yfit3 = parabola_fit(x3,y3);
% eliminate outliers and redo fit on cleaned data
loops = 3;
threshold = -1;
for ck = 1:loops 
    ind = find((y3-yfit3)<-1);
    x3(ind) = [];
    y3(ind) = [];
    yfit3 = parabola_fit(x3,y3);
end

%% bottom line
yfit4 = parabola_fit(x4,y4);
% eliminate outliers and redo fit on cleaned data
loops = 3;
threshold = 3;
for ck = 1:loops 
    ind = find((y4-yfit4)>threshold);
    x4(ind) = [];
    y4(ind) = [];
    yfit4 = parabola_fit(x4,y4);
end


figure(1)
imagesc(A);
set(gca,'YDir','normal');
colormap('gray');
colorbar('vert');
hold on
plot(x,y, '.', x3, yfit3, '-r', x4, yfit4, '-g','linewidth',5)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function yfit = parabola_fit(x,y)
    % PARABOLA FIT
    m = length(x);
    %Set up the appropriate matrix A to find the best-fit parabola of the form y=C+Dx+Ex^2.  The 
    %first column of A will contain all 1's, using the ones() command.  The second column of A 
    %contains x values that are stored in X.  The third column of A contains the squared x values
    %that are stored in X.  Elementwise multiplication of X by itself, using .* operator, will 
    %produce the desired values for the third column.  
    A = [ones(m,1) x x.*x];
    A_transposeA = A.' * A;
    A_transposeY = A.' * y;
    %backslash operation to solve the overdetermined system.
    Soln2 = A_transposeA\A_transposeY;
    %x values to use for plotting the best-fit parabola.  
    yfit = Soln2(1) + Soln2(2)*x + Soln2(3)*x.*x; %
    % 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x3,y3,x4,y4] = top_bottom_boundary(x,y)

% based on FEX : https://fr.mathworks.com/matlabcentral/answers/299796-tight-boundary-around-a-set-of-points


    %split data into classes and find max/min for each class
    class_label = unique(x);
    upper_boundary = zeros(size(class_label));
    lower_boundary = zeros(size(class_label));
    for idx = 1:numel(class_label)
        class = y(x == class_label(idx));
        upper_boundary(idx) = max(class);
        lower_boundary(idx) = min(class);
    end

%     left_boundary = y(x == class_label(1));
%     right_boundary = y(x == class_label(end));
% 
%     % left_boundary
%     x1 = class_label(1)*ones(size(left_boundary));
%     y1 = left_boundary;
% 
%     % right_boundary
%     x2 = class_label(end)*ones(size(right_boundary));
%     y2 = right_boundary;

    % top boundary
    x3 = class_label;
    y3 = upper_boundary;

    % bottom boundary
    x4 = class_label;
    y4 = lower_boundary;

%     x_out = [x1(:);x2(:);x3(:);x4(:);];
%     y_out = [y1(:);y2(:);y3(:);y4(:);];

end

[J. Alex Lee]

Hi, nice approach, I too as a non-expert was interested in this problem and tried a few things that did not work out, so it was nice to see this solution.

Would "discretize" make the "top_bottom_boundary" function any better?

Also, I like how the tasks are broken up into functions, so if the investigator is specifically interested in a circle the call to parabola_fit can easily be replaced by circle fitting, which I'm sure can be found on FEX all over the place.

[Mathieu NOE]

hello @Ounkhin

Problem solved ?