Export Hounsfield Unit values as a table from each segmentation

Hi everyone,

I am a beginner at using 3D Slicer. Before that, I used to work with Mimics. I was wondering is it possible to extract/export the whole set of Hounsfield Unit values from a segmentation performed via threshold tool? I want to use them for calculating bone density. There is a possibility in Mimics that can export a 4 columns table (XYZ and HU value for each element that was created in Mimics after making the mask) and I am looking for the same function in 3D Slicer.

Thanks in advance,

Best,
Somayeh

Yes, of course, getting intensities of voxels in a segmented region is very simple, and the best is that you donā€™t even need to export the file, then load it into your own Python script for analysis and display, but you can do everything using a small Python code snippet within Slicer (copy-paste the code into the built-in Python console, use a Jupyter notebook, or create a small Python-scripted module with a custom GUI).

For example, to get all voxel values of the master volume corresponding to a segment:

# Get segmentation as labelmap volume node
labelmapVolumeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode")
slicer.modules.segmentations.logic().ExportVisibleSegmentsToLabelmapNode(segmentationNode, labelmapVolumeNode, masterVolumeNode)

# Extract all voxels of the segment as numpy array
volumeArray = slicer.util.arrayFromVolume(masterVolumeNode)
labelArray = slicer.util.arrayFromVolume(labelmapVolumeNode)
segmentVoxels = volumeArray[labelArray==labelValue]

See a complete example (that includes generating sample data, extracting voxel values of a segment, computing histogram, and plotting the results) here.

Dear Andras,

Thanks for your answer.

I did not work with Python, so I am a very beginner in this area.
I tried to copy-paste the code that you sent. I renamed the ā€œmasterVolumeNodeā€, ā€œsegmentationNodeā€ and ā€œlabelmapVolumeNodeā€ but it did not work. I was wondering could you please edit the code based on my modelā€™s subhierarchy?

Thanks in advance,

Best,
Somayeh
image

You can set segmentationNode like this:

segmentationNode = getNode('P3-L4')

Thank you for your answer.
I tried this code and it ran without any error, but I did not see any results:
segmentationNode = getNode(ā€˜P3-L4ā€™)
labelmapVolumeNode = getNode(ā€˜P3-L4ā€™)
masterVolumeNode = getNode(ā€˜P3-wholeā€™)

# Get segmentation as labelmap volume node
labelmapVolumeNode = slicer.mrmlScene.AddNewNodeByClass(ā€œvtkMRMLLabelMapVolumeNodeā€)
slicer.modules.segmentations.logic().ExportVisibleSegmentsToLabelmapNode(segmentationNode, labelmapVolumeNode, masterVolumeNode)

# Extract all voxels of the segment as numpy array
volumeArray = slicer.util.arrayFromVolume(masterVolumeNode)
labelArray = slicer.util.arrayFromVolume(labelmapVolumeNode)
labelValue = getNode(ā€˜LabelMapVolumeā€™)
*sgmentVoxels = volumeArray[labelArray==labelValue]

After that, I added two parts to the code (based on the webpage which you sent it):
# Compute histogram
import numpy as np
histogram = np.histogram(segmentVoxels, bins=50)

# Plot histogram
################################################

And I just got this table that does not make sense and can not be true:

Thanks in advance,

Best,
Somayeh

There are a few small issues in your script, for example:

  • remove labelmapVolumeNode = getNode('P3-L4') - you set the variable a few lines below
  • *sgmentVoxels should be segmentVoxels

Always put code between triple-ticks:

```python
your code is here
just another line
```

Otherwise, the code looks OK. I donā€™t see anything wrong in the output either. Why the tables do not make sense for you? What did you expect to see in the tables?

I have the same issue and followed your post for a few days. Today I tried your code but the same as your result. Also I was prompted that I can try this plugin. Thanks to the experts who wrote it. Still canā€™t get the whole set of Hounsfield Unit values from a segmentation. It was just enough for me


Iā€™ve tested this example script and it works perfectly for both latest Slicer Stable Release and Slicer Preview Release.

All you need to do is to replace the first part (that generates input data) with setting your data as masterVolumeNode, segmentationNode, labelValue.

Thanks for your answer.
The reason that I said the table does not make sense:
As I mentioned before, the table which I need should contain 4 columns, like this:

And considering the CT-scans that I used (0.3 mm) I expected around 250,000 data for this sample (something like the data that was exported for similar Vertebrae from Mimics).

Since you have the volume and labels as numpy arrays, you already have all these tables. You just need to specify what you need.

For example, to get coordinates and HU at all segment voxel positions:

coordinates = np.where(labelArray==labelValue)
hu = volumeArray[coordinates]

You can put together an array of i, j, k coordinates and HU values like this:

coordinatesWithHU = np.zeros([len(coordinates[0]), 4])
coordinatesWithHU[:,0:3] = np.array(coordinates).T
coordinatesWithHU[:,3] = hu

# save the dataframe as a csv file
pip_install('pandas')  # needs to executed only once
import pandas as pd
pd.DataFrame(coordinatesWithHU).to_csv("c:/tmp/coordshu.csv")

Such representation would be extremely inefficient for any kind of processing. Iā€™m just curious, why do you consider representing your data in such table?

Thank you for your answer dear Andras.
I still can not get the table that I want, but I am working on it. :frowning:
I need to get density for each element in the sample, so I wrote a code in Matlab to obtain the density. One of the input data for that code is the Hu value based on CT images which I am currently struggling with it now.

The code snippet above creates the table that you described you need.

What kind of processing do you plan to do?

Iā€™m just curious, why do you consider using Matlab? In Python you have many more tools, you donā€™t need to deal with licensing hassles, and Python programming skills are highly valued (while Matlab programming skills are largely irrelevant in industry, and increasingly in academia, too).

Unfortunately I could not get the table which I wanted. In fact, I just decided to learn Python and I need time to write code with it.
I worked with MATLAB because I was familiar with its functions and some of its rules.
I would be appreciated if you could send me any tutorial on learning Python for beginners.

Best,
Somayeh

1 Like

Function to convert HU values (calibrated HU) to density (g/cm3):

def HU2Dens(HU):
    if HU>-1000 and HU<0:
        return HU/1000+1
    elif HU>0:
        return HU/1955+1
    else:
        return 0

And now you can use lassoan script replacing this line

coordinatesWithHU[:,3] = hu

by this one

coordinatesWithHU[:,3] = HU2Dens(hu)

Thank you so much.
Could you get the table via that script (x, y, z, and Hu)?
I want to get that table, but the scripts I used do not work. Could please send me your script?

Best,
Somayeh

A post was split to a new topic: Convert between RAS and numpy array indices