Replace unwanted voxels with nearest neighbor labeled values

Hello Matlab experts,

I am trying to remove unwanted electrode voxels from 3D labeled images by replacing the electrodes with nearest skin, bone or other labeled tissue voxels. I experimented with the function here: https://www.mathworks.com/matlabcentral/fileexchange/24723-nearest_neighbour_3d

But that didn't generate the result I expected. The function takes given points(electrodes) and their minimum distant nearest candidate points(other tissue voxels). It considers Euclidean distance.

%sz is the size of the 3D label=> lab
[cxx, cyy, czz] = ind2sub(sz, find(lab==10|lab==4)); % candidate points from either skin(10) or bones(4) or other tissues
[gxx, gyy, gzz] = ind2sub(sz, find(lab==20)); % given poinst from electrodes(20)
[nxx, nyy, nzz] = ind2sub(sz, compute_nearest_neighbour([cxx,cyy,czz], [gxx,gyy,gzz])); % nearest points
lab(gxx,gyy,gzz) = lab(nxx,nyy,nzz);

Any other ways?

Thanks in advance.

Answer 1

[darova]

Solution

lab = lab_superior;     % the data shared
elec = zeros(size(lab));
elec(lab==20)=1;        % select electrodes
BB = lab;
BB(BB==20) = 0;         % remove electrodes
nb = unique(BB);        % unique regions
nb(nb==0) = [];         % remove '0' from list of regions
[EE,ne] = bwlabeln(elec);  
II = BB*0;
for i = 1:ne
    E1 = EE == i;   
    E2 = imdilate(E1,strel('sphere', 4)); 
    for j = nb(:)'
        B1 = BB == j;                                       % select one region
        tmp = B1 & E2;                                      % compare electrode and region
        if any(tmp(:))                                      % if region and electrode are close
            II = II+j*E1;
            break
        end        
    end
end
A = II + BB;

Original After processing

Answer 2

[Image Analyst]

The position of x and y are swapped, which makes a difference if the volume is not equal lengths in the x and y direction. Remember ind2sub() returns (row, column) which is NOT (x, y) -- it's (y, x).

[banikr]

I see. So the code should be as below?

lab= lab_skin_done;
[czz, cyy, cxx] = ind2sub(sz, find(lab~=20&lab~=0)); 
[gzz, gyy, gxx] = ind2sub(sz, find(lab==20)); % given poinst from electrodes
[nzz, nyy, nxx] = ind2sub(sz, compute_nearest_neighbour([cxx,cyy,czz], [gxx,gyy,gzz])); % nearest points
lab(gxx,gyy,gzz) = lab(nxx,nyy,nzz);

Is that what you mean?

Input to the function is [xx, yy, zz] sequence right? Opposite order of ind2sub output.

Still I am getting points at unexpected locations in the 3D matrix

[darova]

I think it should be switched

[banikr]

@darova

You are right. Thanks for the clarification.

sz = size(lab);
[gzz, gxx, gyy] = ind2sub(sz, find(lab==20));

provides the correct voxel indices to the image 'lab'. See the red objects(voxels) in the image below.

Also candidate points are taken accordingly from labeled 4 or 10 voxels.

[czz, cxx, cyy] = ind2sub(sz, find(lab==10|lab==4));

but the function output is not giving me desired output nearest voxels.

[nzz, nxx, nyy] = ind2sub(sz, compute_nearest_neighbour([czz,cxx,cyy], [gzz,gxx, gyy]));

which is nowhere near g(given) values.

>> max([gzz, gxx, gyy])
ans =
   129   203   159
>> max([nzz, nxx, nyy])
ans =
   170   254    30

Also the maximum values of n values are far from the g ones.

I guess the function is not working...

Any help appreciated.

[darova]

Can you use bwlabeln to separate each region?

[banikr]

@darova,

How separating each region may help? In fact I wanted to replace the 4 object voxels with nearest neighboring labels.

bwlabeln labels the connected components. I have already the four object voxels labelled as 20.

[gzz, gxx, gyy] = ind2sub(sz, find(lab==20));

Can you explain a bit more?

[darova]

Can we start from the start? Can you explain more what you are trying to achieve?

Is this correct?

Are you trying to remove electrodes?

[banikr]

Yes, the white region is electrode(lab==20) in a 3D segmented image(170x256x256). I want to replace those electrode voxels with bones(lab==4) or skin(lab==10) or any other tissue labels.

The background is black with lab==0 so I want to avoid that.

[darova]

Simply use this

lab(lab==20) = 10;

[banikr]

I wish it was that simple.

some electrodes are embedded inside the skull and there are electrodes objects inside deep brain. So can not wrap up with one label in all the voxels.

[darova]

Can you explain more? What electrodes you want to fill/remove and how are you detecting them

[banikr]

The 3D label image has 10 different labels. Tissues and electrode.

>> unique(lab)
ans =
     0
     1
     2
     3
     4
     8
    10
    20
    56

20 for electrode and 4 for bones and 10 for skin and others.

The white voxels in the above axial slice are the electrodes.

[gzz, gxx, gyy] = ind2sub(sz, find(lab==20)); % given points from electrodes
dem_lab = zeros(size(lab));
for i=1:numel(gxx)
    dem_lab(gzz(i),gxx(i),gyy(i))=1; 
end
volumeViewer(dem_lab)

This chunk of code gives the following only electrode voxels.

or simply:

lab(lab~=20)=0;

will provide the same result too.

I want to replace the electrodes with neighboring non-electrode labels. I was thinking about nearest neighbor algorithm. I can hard code it but that will be both memory and time consuming. Trying to find simple or efficient matlab function.

I have segmented the electrodes from CT images which has individual intensity for metals but that is another story.

Please let me know if that is not clear.

Thanks

[darova]

You can dilate each electrode and try to merge them with each area/region separately. If there are commond nodes - close region

if any(A1 & A2)
    disp('electrode is close to region')
end

See this idea

[banikr]

@darova

I think there has been a misunderstanding. The idea you've shared does not GET RID OF the electrode voxels, does it?

Can you clarify how dilation of morphological operations can help? I want to replace the electrodes with nearest labeled voxels not with background or dilated electrode region.

Thanks.

[darova]

Try this example

[Image Analyst]

I've been sitting out for a few reasons. (1) darova seems to be helping you fine. (2) it looks like it might take more time that I can afford to devote, especially if I have to invent data. And mostly, (3) I have no data to work with. If you had uploaded a small 3-D image in a .mat file that illustrates the need, I might have done more on this. Any way you could upload a small chunk of your volume, say a 500x500x300 labeled subvolume that has all 3 or 4 classes in it?

[darova]

I can handle it.

[Image Analyst]

Yes, but it would be easier if you had data. And start your own answer so you can get credit for it if he Accepts or Votes for your Answer. Let's promote you from a "Rising Star" to an MVP.

[banikr]

@Image Analyst

Yeah, I usually upload the .mat file but the medical image files have some confidentiality restraints. But I will share small part of it that relates to the problem. Please see the link below:

Google Drive File link

[banikr]

Hi @darova,

Could you explain why did you take the inverted binary?

I1 = ~im2bw(I);

In my case, I is an image with 10 different labels.

lab= lab_superior;
elec = zeros(size(lab));
elec(lab==20)=1;
% [BB,nb] = bwlabeln((lab~=20&lab~=0)&~elec);    % label regions?
[BB,nb] = bwlabeln((lab~=0)&~elec);                  % region without background and without electrodes?
[EE,ne] = bwlabeln(elec);                        % label all electrode objects. 
II = BB*0;
for i = 1:ne
    E1 = EE == i;                                           % select one electrode
    E2 = imdilate(E1,strel('sphere', 4)); 
    for j = 1:nb
        B1 = BB == j;                                       % select one region
        % imshowpair(E2,E1|B1)                              % 3D volumes won't work
        % pause(1)
        tmp = B1 & E2;                                      % compare electrode and region
        if any(tmp(:))                                      % if region and electrode are close
            II = II + j*(B1 | E2);                          % merge region and electrode
            % imshow(II)                                    % for 3D? 
            % pause(1)
            break
        end        
    end
end

In the end, I get II with 3 unique labels. I changed the dilation considering 3D volume. And not clear about this line here:

II = II + j*(B1 | E2); 

[darova]

• Could you explain why did you take the inverted binary?
• I1 = ~im2bw(I);

It's for bwselect and bwlabel. They want regions to be 'ones' not 'zeros'

• And not clear about this line here:
• II = II + j*(B1 | E2);

You are right, i think better be

II = II + j*E1;         % fill not dilated electrode with the same value as region

Did it work? Do you have more questions?

[banikr]

Hi @darova,

lab = lab_superior;     % the data shared
elec = zeros(size(lab));
elec(lab==20)=1;
I1 = lab;
[BB,nb] = bwlabeln(I1&~elec);  
[EE,ne] = bwlabeln(elec);  
II = BB*0;
for i = 1:ne
    E1 = EE == i;   
    E2 = imdilate(E1,strel('sphere', 4)); 
    for j = 1:nb
        B1 = BB == j;                                       % select one region
        tmp = B1 & E2;                                      % compare electrode and region
        if any(tmp(:))                                      % if region and electrode are close
            II = II+j*E1;
            break
        end        
    end
end

Executing this code gives 'II' with 0 and 1 binary values. The same as elec(with same 475) electrode voxels.

So it didn't work.

Another reason for nearest neighbor is some electrodes share locality with different labels.

So if we dilute the electrodes and assign them common region, will it take the whole electrode as single neighboring label?

I have shared the data in previous comment, were you able to take a look at that?

[darova]

give me permission

[banikr]

Yes, check now.

[darova]

do you smaller file? it's too large 164Mb

[banikr]

Check this .nii.gz format of the same file:

https://drive.google.com/open?id=1K2GDbycS45DG6ZBw2ecGgGIBq_K2wvIY

[darova]

An error occurs

Error using load
Unknown text on line number 1 of ASCII file lab_superior.nii
"\". 

What format is it? .nii

[banikr]

Yes it is saved in nifti file format.

You can access the file as below:

lab = niftiread('lab_superior.nii.gz');
% volumeViewer(lab) % works
>> unique(lab)
ans =
     0
     1
     2
     3
     4
     8
    10
    20

[darova]

It's because all electrodes merge with BB==1

try

A = BB + II;            % eletrodes and brain together
volumeViewer(A)

[banikr]

Did you run with the code I mentioned here in Previous comment ?

So I include the following:

A = BB + II; 

after that, it gives me an image with 28 separate regions. How do I get original label image from that with one less label(20 for electrodes).

>> numel(unique(BB)) == numel(unique(A))
ans =
  logical
   1

[darova]

• Did you run with the code I mentioned here in Previous comment ?

Of course

• How do I get original label image from that with one less label(20 for electrodes).

You didn't tell that you want original labeled image. Just replace BB with your original image (without electrodes)

[banikr]

Ahhh, I mentioned couple of times that I want the electrode voxels replaced with neighboring non-electrode labeled voxels.

If we can do that we have original labeled image without electrodes. If there was a misunderstanding sorry for that.

And could you demonstrate with the codes since you have the data too, so far words are very confusing so I post codes which are easy to understand... for example which BB to replace with which lab (do 'without electrodes' in original image means zeroing out lab==20 voxels?) and which line in the code?

[darova]

Try this

lab = lab_superior;     % the data shared
elec = zeros(size(lab));
elec(lab==20)=1;        % select electrodes
BB = lab;
BB(BB==20) = 0;         % remove electrodes
nb = unique(BB);        % unique regions
[EE,ne] = bwlabeln(elec);  
II = BB*0;
for i = 1:ne
    E1 = EE == i;   
    E2 = imdilate(E1,strel('sphere', 4)); 
    for j = nb(:)'
        B1 = BB == j;                                       % select one region
        tmp = B1 & E2;                                      % compare electrode and region
        if any(tmp(:))                                      % if region and electrode are close
            II = II+j*E1;
            break
        end        
    end
end
A = II + BB;

[banikr]

Hi @darova,

If A ought to be the final result then it removed the electrodes but didn't replace them with neighboring tissues.

lab with electrodes(white voxels)

Thanks for the support, I think I will try some other applications.

[darova]

Im sorry, forgot about '0'. The script merged electrodes with background (zero value)

nb = unique(BB);            % unique regions
nb(nb==0) = [];             % remove '0' from list of regions

Original

After processing

[banikr]

@darova,

I think it worked.

The electrode voxels were replaced with mostly skin(brown) voxels.

Thanks for the support!

[darova]

Can you accept the answer?