Create a path given a binary image?

Hi everyone

I'm trying to automatically measure the distance between 2 layers in a given image.

But I can't figure out how to do it! What I want is to automatically recognize each layer and then measure the distance between them in certain region, as shown in the follow image:

Assuming that layers are not parametrizable forms, I can't use the Houg transformation. Once they are not a completely bright path (I mean, once they have discontinuity points), I can't simply look for points where bright is highest.

Any idea?

Answer 1

[Image Analyst]

First I'd get a mask for each "line". Call imdilate() to thicken each line. Then call bwareafilt(binaryImage, 2) to make sure you have only two blobs. Then call bwlabel to label them.

labeledImage = bwlabel(binaryImage);

Then get each one as a mask

mask1 = labeledImage == 1;
mask2 = labeledImage == 2;

Now mask the image and threshold

line1Image = grayImage;
line1Image(~mask1) = 0; % Erase outside the mask.

Next threshold and get the mean y (row) value for each column

for col = 1 : columns
    thisColumn = binaryImage(:, col);
    if max(thisColumn) > 0 % If there's something there in that column...
        y(col) = mean(find(thisColumn));
        x(xcol) = col;
    end
end

Then you have the mean y of the line, which is hopefully at the center. But because it might not be, and you want something smooth and continuous, I'd pass x and y into sgolayfilt() to smooth out the path. Then repeat for the second mask and subtract the y values to get the delta Y between the two "lines". And of course you have to watch out for columns where only one line is present instead of both (I kind of showed you how to do that already).

[Jose Pedro Dias Oliveira]

Hi Image Analyst!

Many thanks for your help!

I understood everything except the loop. Are y and x a vector? What is xcol? Can you please explain?

Bests, Pedro.

[Image Analyst]

What we need to do is to get two arrays, one for x and one for y, that describe where the center point of that messy, blurry line is. So for each x, or column, I get the x value (which is just the column number itself), and I get the y value (which is the average of the y values of all pixels in the line for that column that are above the threshold level). Does that explain it better?

[Jose Pedro Dias Oliveira]

I think I understand now! I'll try to code, then I post a feedback. Thank you.

[Jose Pedro Dias Oliveira]

Well, I think I left something behind, just don't know what.

se = strel('line',1,1);
aa = imdilate(BinaryImage,se);
bb = bwareafilt(aa, 3);

labeledImage = bwlabel(bb);
mask1 = labeledImage == 1;
mask2 = labeledImage == 2;

line1Image = BinaryImage;
line1Image(~mask1) = 0; % Erase outside the mask.

line2Image = BinaryImage;
line2Image(~mask2) = 0; % Erase outside the mask.

for col = 1 : 926
    thisColumn = BinaryImage(:, col);
    if max(thisColumn) > 0
        y(col) = mean(find(thisColumn));
        x(col) = col;
    end
end

where:

thisColumn = BinaryImage(:, col);
Is it in BinaryImage or line Image?

And... where:

    y(col) = mean(find(thisColumn));

I think I miss something... Didn't I?

At last, x and y are 2 vectors, ok... Y contains the mean value of each column, and x is the column index... It does not makes sense for me... I should have 2 y vector and 2 x vectors, corresponding to each of the layers. Or didn't I?

EDIT: where is the code that works:

BinaryImage = BinImgCornea;

se = strel('line',5,3);
aa = imdilate(BinaryImage,se);
%imshow(aa)
bb = bwareafilt(aa, 3);

labeledImage = bwlabel(bb);
mask1 = labeledImage == 1;
mask2 = labeledImage == 2;

line1Image = BinaryImage;
line1Image(~mask1) = 0; % Erase outside the mask.

line2Image = BinaryImage;
line2Image(~mask2) = 0; % Erase outside the mask.

%
% for col = 1 : length(BinaryImage)
%     thisColumn = BinaryImage(:, col);
%     if max(thisColumn) > 0 % If there's something there in that column...
%         y(col) = mean(find(thisColumn==1));
%         x(col) = col;
%     end
% end

y1 = [];
x1 = [];
for col1 = 1 : length(line1Image)
    thisColumn = line1Image(:, col1);
    if max(thisColumn) > 0 % If there's something there in that column...
        y1 = [y1 mean(find(thisColumn))];
        x1 = [x1 col1];
    end
end

x2=[];y2=[];
for col2 = 1 : length(line2Image)
    thisColumn = line2Image(:, col2);
    if max(thisColumn) > 0 % If there's something there in that column...
        y2 = [y2 mean(find(thisColumn))];
        x2 = [x2 col2];
    end
end

surface = zeros(size(line1Image));

for ii = 1:length(x1)
    surface(round(y1(ii)),x1(ii)) = 1;
end

for jj = 1:length(x2)
    surface(round(y2(jj)),x2(jj)) = 1;
end

[Image Analyst]

I can try it later this evening if you post the top half of the image. Just the two lines, not any screenshot or image with red and green annotations.