Measure distances in binary image along specified directions

Hello,

In the attached image (binary_im_example), I would like to, for every black pixel, measure the distance to the nearest white pixel, along specified directions.

The second attached image (binary_dist_example) shows an example of what I mean by this, for a particular pixel. Along 8 lines (vertical up and down, horizontal left and right and along lines at 45 degrees to these) I would like to measure the distance from every black pixel to the nearest white pixel.

Overall, my goal is to find an "average" distance for each black pixel, to the nearest white pixel (by averaging the distance along these 8 lines). For this purpose the distance transform is not suitable, because this would only give the minimum distance.

If there are alternative ways to measure an average distance in this way, it would also be interesting to hear about them!

Cheers

[KSSV]

There are many pixles which are nearest, why did you choose only those pixels?

[Angus Coyne]

In the attached image (Binary_dist_example), I have found the nearest white pixels to the indicated black pixel along 8 specified lines. These specified lines are vertically up and down, horizontally left and right, and the 4 lines at 45 degrees to these.

So I think the answer to your question is: those 8 pixels were chosen because they lie along the specified lines. Closer pixels may exist, but they do not lie along the specified lines.

The reason I am looking along specified lines, is that for every black pixel, I would like a measure of the "average" distance to a white pixel. By average, I mean an average over all possible directions.

[DGM]

In this scheme, is the distance between two diagonal neighbors 1 or sqrt(2)?

[Angus Coyne]

I would like to use Euclidian distances, so the distance between diagonal neighbours is sqrt(2).

[Jan]

What does "average over all possible directions" exactly mean? There is an infinite number of directions. Maybe you want to use the directions between a pixel and all pixels on the border of the image? But does this satisfy your needs for a pixel on the border? It has valid directions over 180° only, while center pixels have 360°. So the pixels at the center might be weighted higher than the pixels near to the borders.

Start with defining the wanted procedure mathematically. With an exact definition, the implementation in code is not compicated anymore.

[DGM]

The definition seems fairly clear to me. OP wants a distance transform that only considers pixels along the four cardinal directions and four diagonals WRT any background pixel, so it's like a restricted case of a regular euclidean distance xform.

I'm not coming up with remotely elegant solutions at the moment.

[Angus Coyne]

In this case, "average over all possible directions" means an average over 8 specified directions - vertically up and down, horizontally left and right, and along the 4 lines at 45 degrees to these. This is of course not actually all possible directions, sorry if this was unclear!

Pixels near or on the border indeed need special consideration. For the approach I have in mind (using the lines along 8 specified directions) I would ideally like to ignore distances which are the result of a line (along one of the 8 specified directions) intersecting the border, and only take the average of distances which are the result of a line (along one of the 8 specified directions) intersecting a white pixel.

Ultimately, I would like to use this procedure on several different binary images, and visually compare these average distances between the images (e.g. by colouring the black pixels according to their average distance). For these purposes, having artefacts related to pixels at or near the border would not be a serious issue.

[DGM]

I was going to say something to the effect of "I haven't come up with an elegant solution, but I'm sure something better than a naive approach would exist", but given the conditional handling of cases where no directional neighbor exists, I'm not so sure about that anymore.

I'm probably going to have to tap out anyway, considering I can't even stay online for more than 30 seconds at a time right now.

[Matt J]

Certain pixels have no nearest white pixel along some of the 8 directions. For example, the black pixels along the lower edge of the image can have at most one vertical neighbor. How are those cases to be handled?

[Angus Coyne]

For my purposes overall, it wouldn't matter if there were artefacts related to pixels at the border, so these needn’t be handled very "elegantly".

They could be handled, for example, by discounting distances which result from an intersection with the edge of an image, rather than an intersection with a white pixel. So for pixels along the lower edge, at least 3 of the 8 directions would be discounted when calculating the "average" distance to a white pixel (south, south-west and south-east would be discounted).

Answer 1

[Matt J]

Well, getting the horizontal and vertical distances are pretty easy (see below). For the diagonal distances, I would apply imrotate(__,'loose'), similar to what I've done with rot90 below. There is some interpolation and post-cropping that will be involved with this, which will introduce approximation errors, but I think it's probably a worthwhile compromise.

load BWImage.mat

distances=cat(3,verticalDistances(BW),...                          %vertical
                rot90(  verticalDistances( rot90(BW,1) )  ,-1 ) ); %horizontal



function D=verticalDistances(BW)

 [m,n]=size(BW);
 BWc=~BW;

 I=(1:m)';

 A=repmat(I,1,n);
 
 C=I.*BW;
 C(BWc)=nan;

 D(:,:,2)=(A-fillmissing(C,'previous',1)).*BWc; %upper distances
 D(:,:,1)=(fillmissing(C,'next',1)-A).*BWc;     %lower distances

end

[Angus Coyne]

Hi Matt,

Thank you very much for the suggestion. I've taken your idea (which works perfectly) and implemented it as shown below. Instead of rotating the image, to find the horizontal distances, I check distances along the rows rather than the columns.

Using the below I am able to exactly half-solve the problem, as I can find the distances in 4 out of 8 directions (vertically up and down, and horizontally left and right).

To find the distances along the 4 lines at 45 degrees to these I could rotate the image, as you suggest. As you mentioned this introduces some approximation errors. Also, rotating changes the dimensions of the image, meaning the 4 distances measured this way don't have the same indexes as the 4 distances we have measured so far, so if I want the average distance (along the 8 directions), I have to convert the indexes in the rotated image to their corresponding indexes in the unrotated image.

So, I would like to find the remaining 4 distances (along the diagonals) from the unrotated image, if possible. This means searching diagonally along the image, from each black pixel, for the nearest white pixel.

Could we somehow use linear indexing? Taking steps of ±(n-1) or ±(-n-1) from the index of each black pixel to find the nearest white pixel (where n is the number of rows in the image).

[m,n]=size(BW);
BWc=~BW;
I=(1:m)';
J=(1:n);
row_matrix=repmat(I,1,n);
col_matrix=repmat(J,m,1);
C=I.*BW;
R=J.*BW;
C(BWc)=nan;
R(BWc)=nan;
N=(row_matrix-fillmissing(C,'previous',1)).*BWc; %upper distances
S=(fillmissing(C,'next',1)-row_matrix).*BWc; %lower distances
W=(col_matrix-fillmissing(R,'previous',2)).*BWc; %left distance 
E=(fillmissing(R,'next',2)-col_matrix).*BWc; %right distances

[Matt J]

Perhaps as follows:

load BWImage.mat


sz=size(BW);
rmask=rasterMask(sz);
forw=@(z) forwRast(z,rmask);
back=@(z) invRast(z,rmask);

D1=westeast(BW);                      %west-east
D2=rot90( westeast(rot90(BW) )  ,-1); %north-south

D3=sqrt(2)*back(westeast(forw(BW))); %southwest-northeast
D4=sqrt(2)*fliplr(back(westeast(forw(fliplr(BW))))); %southeast-northwest


distances=cat(3, D1,D2,D3,D4);

imshow(mean(distances,3,'omitnan'),[])

function D=westeast(BW)

    [m,n]=size(BW);
    BWc=~BW;
    J=(1:n);

    col_matrix=repmat(J,m,1);
    R=J.*BW;
    R(BWc)=nan;
    W=(col_matrix-fillmissing(R,'previous',2)).*BWc; %left distance 
    E=(fillmissing(R,'next',2)-col_matrix).*BWc; %right distances

    D=cat(3,W,E);

end



function mask=rasterMask(sz)

 [m,n]=deal(sz(1),sz(2));

 rows=hankel(1:m,m+(0:n-1));
 cols=repmat(1:n,m,1);

 mask=logical( accumarray([rows(:),cols(:)],1,[max(rows(:)),n]));

end

function Arast=forwRast(A,mask)
 Arast=double(mask);
 Arast(mask)=A(:);
end

function A=invRast(Arast,mask)

 m=sum(mask(:,1));
 n=size(mask,2);

 A=nan(m,n);
 for i=1:size(Arast,3)
   tmp =Arast(:,:,i);
   tmp=reshape( tmp(mask),[m,n]);
   A(:,:,i)=tmp;
 end
end

Answer 2

[Image Analyst]

Is it going to be in 3-D? That measurement is called "star volume". There are papers on it.

[Angus Coyne]

For me this analysis is only in 2D, although I had not heard of star volume previously. After a quick Google I found it first mentioned in this paper: "Star volume of marrow space and trabeculae in iliac crest: Sampling procedure and correlation to star volume of first lumbar vertebra", Bone, 1990.

[Image Analyst]

Yes. It's commonly used in open cell foams, like trabecular bone, where you want a "cell size" but the concept of a "cell" is not so well defined, as opposed to a closed cell foam where all the cells are like air bubbles completely enclosed by the material.