Coding volume segmentation with multiple thresholds (CT data)

I'm working with CT data to visualise sediment structures in 3D. I'd like to apply several masks at different CT number thresholds to create a multi-colour overlay. For example, voxels with CT numbers between 250 HU and 299 HU will have a different colour to those between 300-349 HU, 350-399 HU, etc. (This is to visualise spatial relationships between sediments of different density).

Using the Volume Segmenter app, I can run Filter & Threshold, manually adding a new label (e.g. HU250, HU300...) and colour for each threshold. I can also save these as a 3D categorical to be read back into Volume Segmenter. But with ~18 thresholds this isn't practical to do manually on several datasets. I'd like to program an equivalent to generate a set of masks for each dataset, and display the result through volshow.

I can iterate through different CT number thresholds, using a function to create a binary mask (a 3D matrix) at each. However, I'm not sure how to link these masks to labels and create a categorical that can be read back into Volume Segmenter.

It is possible to bypass masks entirely by using the normalised CT number data as OverlayData in volshow, but I would like to adjust opacity and the thresholds of what densities are shown on the CT scan (which can be done in the Volume Segmenter app), and programming an interface to do that seems far more complicated.

I would be grateful for advice either on creating the categorical set of labels and masks, or on alternative methods.

Answer 1

[Matt J]

I can iterate through different CT number thresholds, using a function to create a binary mask (a 3D matrix) at each. However, I'm not sure how to link these masks to labels and create a categorical

If you have the thresholds already, you can just use,

[___] = discretize(___,'categorical')

[EL]

Thank you for the suggestion! That does look useful, but I may be completely misunderstanding how to use it.

The creation of the masks is something like

for n=1:18 %iterating through each threshold 
    mask=imbinarize(CTdata,n*50) %creating the mask (thresholds are 50 HU apart).
end

but this loop is no use if I can't save the data. To keep everything together I could do something like

mask=zeroes(18, datadim1, datadim2,datadim3) %creating a variable to contain 18 labels and their 3D masks
for n=1:18 %iterating through each threshold 
    mask(n)=imbinarize(CTdata,n*50) %creating the mask (thresholds are 50 HU apart).
end

which seems clunky. So far I've avoided using 4D matrices because of memory issues - I think there must be a more efficient approach that I'm missing? Or is this completely the wrong track? I'm not sure how to distribute the data into threshold categories using discretize if the data isn't put into a single structure to start with.

[Matt J]

I don't think you should be making masks at all. Since you know the thresholds, you can go straight to the label map. There is nothing you can do with masks that you can't do with a label map.

CTdata=randi(200,3,3,2)

CTdata = 
CTdata(:,:,1) =

   120   122    78
    71   181   130
   113     6   136


CTdata(:,:,2) =

    36   102   136
   116    74   169
    68   158   122

edges=(0:4)*50;
labelmap=discretize(CTdata,edges,'categorical',{'CSF', 'Adipose','Muscle','Bone'})

labelmap = 3×3×2 categorical array
labelmap(:,:,1) = 

     Muscle       Muscle      Adipose 
     Adipose      Bone        Muscle  
     Muscle       CSF         Muscle  


labelmap(:,:,2) = 

     CSF          Muscle       Muscle 
     Muscle       Adipose      Bone   
     Adipose      Bone         Muscle 

[EL]

This is perfect, thank you! It is the concept of a label map that I was missing, and that simplifies this hugely - now if I set the colours in Volume Segmenter I should be able to open any dataset with its corresponding label map and generate a colour-coded volume.