Splitting a binary image into 2 parts

HI everybody! I have a binary image matrix uploaded as A.mat:

White pixels correspond to 1, black regions correspond to 0. I would like to split the curve into left and right parts, something like the following:

and

Initially, I thought of finding the indices of midponts of elements having 1 at their positions, at each row. Using the indices, I can substitute ones on the left (and right) of the midpoints (at each row) using relational operators. Maybe there might be a better way to do this, but this method comes to my mind right now.

So I start from the bottommost row, and I keep on moving upwards. As evident from the A.mat figure (first figure), a row can have more than 2 cells having values one. To choose the appropriate cells for calculating midpoints at row "r", I compare the number of elements having one on left and right side of midpoint at row "r+1". I am attaching the code below:

clc;
clearvars;
load('A.mat');
[m,n] = size(A);
midpoint = zeros(m,n);
midpt_idx = zeros(m,1);
%%
A_sep = A;

for r = m:-1:1
    k = find(A(r,:));
    if size(k,2)==2
        midpt_idx(r) = ceil((k(1)+k(2))*0.5);
    elseif sum(A(r,1:midpt_idx(r+1))) <  sum(A(r,midpt_idx(r+1)+1:end))
        midpt_idx(r) = ceil((k(1)+k(2))/2);
    elseif sum(A(r,1:midpt_idx(r+1))) >  sum(A(r,midpt_idx(r+1)+1:end))
        midpt_idx(r) = ceil((k(end-1)+k(end))/2);
    elseif sum(A(r,1:midpt_idx(r+1))) ==  sum(A(r,midpt_idx(r+1)+1:end))
        midpt_idx(r) = ceil((k(size(k,2)/2)+k((size(k,2)/2)+1))*0.5);
    end
    A_sep(r,midpt_idx(r)) = 2; % Midpoints have value 2 in their cells, I change this to 1 for checking solution by using imshow(A_sep)
end

An error is dispayed when execution comes at row 64. There is only one cell having 1 in row 64. For now, this error can be ignored (I think!).

This code runs correctly for some parts, but fails at other parts:

If you run the code along with A.mat file, starting from row 416, it gives correct results till row 278. At row 277, the code fails because there are a larger number of points adjacent to each other on right side, compared to left side of the midpoint index at row 279.

So, basically my method fails in these situations. Can you please help me with this? Maybe an advice for my code, or some other method which is better than what I am doing right now.

Thanks!

Answer 1

[Matt J]

There is a small break in the curve which I assumed was supposed to indicate the desired dividing point between the right and left sides.

Under that assumption, I propose the following:

load A

[I,J]=find(A);
B=[I,J];

[b,j]=sortrows(B,[+1,-2]);
j=j(end);

C=regionprops(A,'PixelList');

for i=1:numel(C)
   if ismember(flip(B(j,:)),C(i).PixelList,'rows')
      C=fliplr(C(i).PixelList); break 
   end
end

Left=accumarray(C,1,size(A)); %the resultant images
Right=A&~Left;

[Ankit Sahay]

This is excellent! Thank you so much. I will try to understand how your code works. Thanks a lot!

[Ankit Sahay]

Hi Matt! I used your code for my data, and it worked perfectly except for a few cases. I am attaching a screenshot of a case where things go wrong:

The middle binary image is the original image, which I want to separate into 2 parts. I have uploaded the binary matrix for the middle binary image (flame_edge_right.mat).

As you can see, the separation point is not correct. maybe I should have been more careful in formulating my original question, I would like the separation point to be the cell at the highest point, which, in the data I have uploaded, is [X,Y] = [283,60] or [X,Y] = [284,60]. I am writing 'or' because it is possible that the highest row having cell value 1 can have multiple columns having 1 in their cells, in which case, it is completely fine to take any point.

To solve this problem, I started going through your code. I was able to understand till

C=regionprops(A,'PixelList');

command, where the edge is divided into 7 parts. If we colormap the 7 parts individually, we get the following results:

Unfortunately, I could not figure out the working of the loop that you have used in your code.

Can you please help me solve my issue? I would like my separation point to be the cell having 1 and positioned at the highest row (or possibly, the extreme right cell in that particular row, if multiple cells have 1 in the highest row). Please let me know if I need to make any clarification in my question. If the images are not visible properly, I can upload them again with higher resolution.

Thanks!

[Matt J]

I used your code for my data, and it worked perfectly except for a few cases. I am attaching a screenshot of a case where things go wrong:

I'm not sure how it could have "worked perfectly", if you say the dividing point is supposed to be the "highest point". My code definitely doesn't give that, even in the first data set that you presented to us.

The reason my code is confounded is that your curves are not fully connected. They have multiple breaks which, as you have now updated us, do not represent the desired breaking point between "left" and "right". If you can find a way to seal the breaks between the connected pixels (except for the breaks where you want the left/right division to occur), then my code should work fine for you.

If you cannot seal the breaks, then your problem is simply ill-defined. There is no way to unambiguously define "left" and "right", if your curve has 3 or more disconnected parts. Imagine, for example, if you had the simple 10x10 image below. To which side do the ones in the 5th column belong? I could argue they could be assigned either to the left or right section.

A =

     0     0     0     0     1     0     0     0     0     0
     0     0     0     0     1     0     0     0     0     0
     0     0     1     0     1     0     1     1     0     0
     0     0     1     0     1     0     0     1     0     0
     0     0     1     0     0     0     0     1     0     0
     0     0     1     0     0     0     0     1     0     0
     0     0     1     0     0     0     0     1     0     0
     0     0     1     0     0     0     0     1     0     0
     0     0     1     0     0     0     0     1     0     0
     0     0     1     0     0     0     0     1     0     0

[Ankit Sahay]

Yeah, I get it now. I did not anticipate such a situation when I started the analysis. It is only after analysing a few more data did I come across a situation where the code failed. Let me reformulate the question such that no ambiguities arise I will take help from the 10X10 matrix that you have provided. Thanks for your input.

Answer 2

[Image Analyst]

Why not simply scan your image column by column with find() finding the topmost and bottommost row?

[Ankit Sahay]

Yes I can definitely do that. But how should I proceed after that? What should I do to ensure that the separation point lies on a cell having 1 in the topmost row?

[Image Analyst]

Let's say your blob looks like a capital J. Where would you want the boundary to start?

[Ankit Sahay]

I would like my boundary to start at the midpoints of the horizontal portions of letter J. I am attaching an image for your reference.

The middle image is the letter J. The cells colored yellow have value 1, the cells colored blue have value 0. I have used imagesc() instead of imshow(); only for illustration purposes.

The left and the right images are the left and right portions of the blob, respectively. I hope this is what you meant in your response when you talk about a "blob like a capital J". Please let me know if I have misunderstood you or if I need to make any clarifications.

P.S. Maybe there is a slight discrepancy in fixing the separation point when many adjacent cells having value 1 are present in topmost row, like the situation for blob J in the image above. In that case, we can start the boundary at the middle point if the number of adjacent cells having 1 in the topmost row are odd, or at one of the middle points (no such preference) if the number of adjacent cells (having 1 in the topmost row) are even.

[Image Analyst]

What is they are in different segments, like a U? Start at the middle of the left most segment in the top row?

Why do you want this anyway? What does it matter? And how did you get the edge image in the first place? I'm wondering if we even need it. Lots of people call edge() when there is no reason to and they could segment the blob perfectly well, and better, using thresholding. Then if they need the oundary, they can use bwboundary. Using edge() causes lots of problems, like disconnected edges and then you'll have to do "edge linking" if you want a closed figure (demo attached).

Let's assume you could split this into upper/lower or left/right parts. What would you then do next with that information?

[Matt J]

I would like my boundary to start at the midpoints of the horizontal portions of letter J.

There are 3 horizontal portions in the J. How did you know to exclude the horizontal portion along the bottom?

[Ankit Sahay]

Yes, this is a valid question. I am really sorry, I should have posed my question properly. I now understand the contradictions that my argument has. I will try to reformulate my question such that there are no more contradictions. Please give me some time to work on this. Thanks a lot for your time!

[Ankit Sahay]

To Image Analyst: My actual data, similar to what I had included in the original post, is actually a flame front of a combustor. I need to separate the blob with the 2 sub-parts starting from the topmost cell, and further continue my work on the separated sub-parts. Thanks to Matt, I now understand the contradictions that my question poses to itself. Please give me some time to work on this and reformulate the question so that no further contradictions arise. Thank you! I will look into the demo that you have attached to see if it can help me.

[Image Analyst]

Again, Matt and I need to see your original gray scale image of your flame front to see if an edge detection is even the right approach.

[Ankit Sahay]

Yeah sure.

Set340_00H2_PLIF50.mat is the original .mat file that I have. I don't have a tiff image as my collaborators have shared the mat files with us.

I extract the relevant data from the .mat file, use median filtering, perform some morphological operations and some editing to extract the relevant region. Please note that the following section of code wasn't originally written by me. It has been used by my collegues for the past few years, I think.

load('Set340_00H2_PLIF50.mat');
x = Xtt.Ax;
y = Xtt.Ay;
rowcon = [44:463];
colcon = [15:762];

x1 = x(colcon);
y1 = y(rowcon);
ms = 3; % used for median filtering

Xtt_AW_act = Xtt.Aw;
Ivec = Xtt_AW_act';
Xtt_AW_nor = (Ivec./max(Ivec, [], 'all'))';
Xtt_AW_nor_red = Xtt_AW_nor(rowcon,colcon);
Kmedian = medfilt2(Xtt_AW_nor_red,[ms,ms]); % median filtering, removes
                                            % salt-and-pepper noise
ttmean = mean(mean(adaptthresh(Kmedian)));
ImbL1 = imbinarize(Kmedian,adaptthresh(Kmedian)-(0.1*ttmean)); 
ImbL2 = ImbL1;
ImbL2(125:end,300:440)=1;
ImbL2(1:end,end:end)=1;
ImbL2(1:end,1:10)=1;
ImbL2(end-20:end,1:end)=1;

BW=imfill(ImbL2,'holes');
dil_ima = imdilate(BW,strel('disk',3));
FlamBouSmo = bwmorph(dil_ima,'majority');

fbf= bwmorph(FlamBouSmo,'remove','Inf');
fb= bwmorph(fbf,'thin','Inf');

fb(:,1)=0;
fb(:,end)=0;
fb(1,:)=0;
fb(end,:)=0;

props = regionprops(fb, 'Area');
propsmat = cell2mat( struct2cell( props));
aremaxord = sort(propsmat, 'descend');
fb2 = bwareaopen(fb, aremaxord(2));

fb2(1:4,:) = [];
FlaBou = flipud(fb2);

figure(1)
imshow(Xtt_AW_nor_red,[0 0.2]);
axis xy

figure(2)
imshowpair(Xtt_AW_nor_red,fb2,'falsecolor');
axis xy

This gives me the following 2 figures:

I now perform some more operations according to my requirements. The relevant code is attached as flame_edge_code.m. The code is quite rusty and inefficient right now, as I am still working on it and the project is in its nascent stages. Out of the many results that I obtain from flame_edge_code.m, I get the following (from the very last line of flame_edge_code.m code):

My original post in the present question starts from this point onwards. I need to separate the 2 flame fronts into 4 parts. I now understand that my question is quite ambiguous, and I am working on reframing the problem correctly.

I fully understand that this post is quite huge and maybe not very clear. Please let me know if I need to clarify anything. Again, I cannot thank you guys enough for giving so much time to my doubts! I really hope that someday I can be of help to others too.

[Image Analyst]

No it doesn't. It's not a complete program. It gives this:

>> flame_edge_code

Unrecognized function or variable 'fb2'.

Error in flame_edge_code (line 2)

fb2_l = fb2(:,1:end/2);

Please delay us (and you) no longer and give us a program that runs.

[Ankit Sahay]

I apologize profusely for this mistake. I should have been more careful. I have attached the complete program along with the .mat file that the programs needs to run. Again, I am really sorry for wasting your time.