Image extraction, segmentation, resampling, 2D, MRI cardiac images

[tommydino93]

Hi everyone!
Since I am about to begin a radiomic study ex novo, I just wanted to ask some questions so as to avoid annoying future problems. We will extract images and masks of interest next week from cardiac MRI with the intent of performing a future binary classification.

1. Should I necessarily extract images in .nrrd format or the .dcm (DICOM) format is fine and easy to convert into .nrrd?

2. Should the mask simply be an image with ones in the ROI and zeros outside of it?

3. Since we will probably work with single slices (we will select the most significant 2D slice for each patient), how should I set my third dimension since pyradiomics wants volumes as input? Moreover, which are the features that I will be able to extract? (meaning...are there some features that are only applicable in 3D?)

4. This is a strange question: one of the two ground truths of the final binary classification is myocarditis and this has the problem that it often doesn't appear as a single concentrated region in the image. As a matter of fact, since it is an inflammatory phenomenon, it often appears as sparse sub regions on the image. How should we deal with masks in this case? Should we just select the biggest region of inflammation? Should we somehow average the inflamed regions?

Thank you very much in advance,
Tommaso

[JoostJM]

1. Should I necessarily extract images in .nrrd format or the .dcm (DICOM) format is fine and easy to convert into .nrrd?

Yes, see also #361

2. Should the mask simply be an image with ones in the ROI and zeros outside of it?

Yes, although other values than 1 are allowed (you'd then need to specify that value in the label parameter, though). Importantly, the direction, origin and size must match to the image. If this is not the case, but the region of interest occupies a physical space contained within the image bounds, you can still use it by setting correctMask to True.

3. Since we will probably work with single slices (we will select the most significant 2D slice for each patient), how should I set my third dimension since pyradiomics wants volumes as input? Moreover, which are the features that I will be able to extract? (meaning...are there some features that are only applicable in 3D?)

Your third dimension is allowed to have size 1, it's just that it has to be present. All features can be extracted from 2D slices. However, take care when applying Wavelet or LoG filters, these assume true 3D input.

4. This is a strange question: one of the two ground truths of the final binary classification is myocarditis and this has the problem that it often doesn't appear as a single concentrated region in the image. As a matter of fact, since it is an inflammatory phenomenon, it often appears as sparse sub regions on the image. How should we deal with masks in this case? Should we just select the biggest region of inflammation? Should we somehow average the inflamed regions?

It's not a strange question. Quite the opposite, it's a very important, but difficult question. I don't have a direct answer, but you can try a large segmentation (containing all areas) or an average of multiple sub regions.

[tommydino93]

Thank you very much again @JoostJM! I'll keep in mind the mask origin/direction/size issue.

As for question 4, I will speak with my supervisor (the doctor in charge of the research project) on Monday and we'll try to figure out what makes more sense from a clinical point of view.

I was thinking that by averaging all the sub-regions we could lose the gray level dependencies, and it would also be difficult because from what I understood the sub-regions often have different dimensions and shapes. Anyway, I'll keep you updated!

Thanks again for your support!

[tommydino93]

Hi @JoostJM!
At the end, as far as myocarditis is concerned, I think we'll just select the widest inflamed region. In the meantime, since I don't have the cardiac images yet, I tried to load in PyCharm a simple brain image and corresponding mask that I downloaded from the local PACS just as a trial.

I have stored, in the same directory (All_Data_nrrd), two files which are: brain_image.nrrd and brain_label.nrrd. Both of them have size (512,512,1) and label is of course the mask. I checked the dimension through this.

Then on PyCharm I just wrote the following code trying to emulate the "brain1" example:

from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
import SimpleITK as sitk
import six
from radiomics import firstorder, getTestCase, glcm, glrlm, glszm, imageoperations, shape, getImageTypes, \
    getFeatureClasses, getParameterValidationFiles


imageName, maskName = getTestCase('brain', dataDirectory= "C:\\Users\\tommaso\\Desktop\\Lavoro\\Zurich_Fellowship\\Programmi\\All_Data_nrrd")

and when I run it I get the following warning:

Testcase "brain" not recognized!

What am I doing wrong? Thank you very much!

[JoostJM]

@tommydino93, that function finds the pyradiomics test cases. Try brain1.

[tommydino93]

@JoostJM I managed to run everything with the brain1 example.....but what I need now for MY brain image and for my future cardiac images will be to load images (and masks) from one of my own directories and not from the pyradiomics data test cases.

Which function could I use to do that? I thought that getTestCase could be used, just modifying the "dataDirectory" parameter..

[JoostJM]

@tommydino93, The getTestCase Function cannot be used for such a case, as it also contains code to download the test data from github when it's not found (which would be quite difficult for your data).

Easiest to use your own data is like this:

import os
datadir = "C:\\Users\\tommaso\\Desktop\\Lavoro\\Zurich_Fellowship\\Programmi\\All_Data_nrrd"
imageName = os.path.join(datadir, 'brain_image.nrrd')
maskName = os.path.join(datadir, 'brain_label.nrrd')

This sets the last two variables to the path of your data, which is similar to what getTestCase does for the pyradiomics test data (with the addition that it downloads the data if not available). It does not load or preprocess the data in any way.

[tommydino93]

Ok, thank you very much! And how about the settings (e.g. binWidth, interpolator, resampledPixelSpacing, voxelArrayShift)? Is there a rule of thumb to choose them?

P.S: would it be a problem if I commented in a fork some of the code lines? I am trying to figure it out on my own just following the code (for instance in pyradiomics/radiomics/first order) but coming from Matlab (they spoiled me too much, I admit that) I am having trouble following even the single inputs and outputs for every functions (e.s. input datatype, output datatype, parameters datatype)

[JoostJM]

@tommydino93, for that, check out the example settings in the repository. Besides that the only advice I can give you is to just see what works.

P.S: would it be a problem if I commented in a fork some of the code lines? I am trying to figure it out on my own just following the code (for instance in pyradiomics/radiomics/first order) but coming from Matlab (they spoiled me too much, I admit that) I am having trouble following even the single inputs and outputs for every functions (e.s. input datatype, output datatype, parameters datatype)

Sure no problem.

[tommydino93]

Hi @JoostJM! I managed to run this code for my brain image, following the brain1 example:

from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
import SimpleITK as sitk
import six
from radiomics import firstorder, getTestCase, glcm, glrlm, glszm, imageoperations, shape, getImageTypes, \
    getFeatureClasses, getParameterValidationFiles
import os

datadir = "C:\\Users\\tommaso\\Desktop\\Lavoro\\Zurich_Fellowship\\Programmi\\All_Data_nrrd"
imageName = os.path.join(datadir, 'brain_image.nrrd')
maskName = os.path.join(datadir, 'brain_label.nrrd')

image = sitk.ReadImage(imageName)
mask = sitk.ReadImage(maskName)

applyLog = False
applyWavelet = False

#settings for the feature calculation. Currently resampling is disabled.
settings = {'binWidth':25, 'interpolator': sitk.sitkBSpline, 'resampledPixelSpacing':None}

if settings['interpolator'] is not None and settings['resampledPixelSpacing'] is not None:
  image, mask = imageoperations.resampleImage(image,mask,settings['resampledPixelSpacing'],settings['interpolator'])
else:
  bb,correctedMask = imageoperations.checkMask(image,mask) #bb stands for bounding box
  if correctedMask is not None: #if the mask has been correctly resampled to the input image geometry
    mask = correctedMask
  image,mask = imageoperations.cropToTumorMask(image,mask,bb)

#let's now show the first order feature calculations
firstOrderFeatures = firstorder.RadiomicsFirstOrder(image, mask, **settings) #la notazione ** serve per passare i valori del dict settings
firstOrderFeatures.enableAllFeatures()

print('Calculating first order features...')
firstOrderFeatures.calculateFeatures()
print('Done')

print('Calculated first order features: ')
for (key, val) in six.iteritems(firstOrderFeatures.featureValues):
  print('  ',key, ':', val)

Everything works fine, meaning that the program extracts all numerical first order features. What I wanted to ask you is:

1. Is there a way to check if the checkMask function worked properly? For instance, could I visualize (maybe with pyplot) the segmented mask?

2. Same thing for cropToTumorMask...how can I visualize what has been assigned to image and mask in the line image,mask = imageoperations.cropToTumorMask(image,mask,bb) and how can I be sure that it's correct?

3. When in the checkMask help you say "check if the label is present in the mask"...what do you mean by "label"? Do you mean the ones of the ROI? As I told you my mask is composed by ones in the ROI and zeros outside of it.

[JoostJM]

1. The corrected mask is usually None (That means no correction applied). The only case where it is not None is when it has to be resampled due to image geometry mismatch between image and mask (see also parameter correctMask). If the checks fail, (None, None) is returned.

2. image,mask = imageoperations.cropToTumorMask(image,mask,bb) returns the cropped image and mask as SimpleITK Images. You can store them as follows:

sitk.WriteImage(Image, "image.nrrd")  # Don't forget the extension!
sitk.WriteImage(Mask, "mask.nrrd", True)  # True here specifies it can apply compression

This will store it as "image.nrrd" and "mask.nrrd" in your current working directory.

3. Yes. Sometimes a labelmap contains multiple labels (e.g. multiple lesions). A "label" here means just that certain pixels have value 1, others value 2, etc. The check here is intended to ensure that the desired label (specified in the customization, default 1) is actually present in the mask, i.e. that there is at least 1 voxel in the mask that has the same value as label.

[tommydino93]

Hi again @JoostJM!
Thank you very much for the information. If it's not a problem, I also wanted to ask you something about the settings:

1. For instance, what it the meaning of "resampledPixelSpacing"? I get that it has to be a list of 3 floats (x,y,z) and that it represents the size of the voxel but....which is the order of magnitude? which is the unit of measurement? (e.g. what does 4,4,4 would mean?)

2. Same thing for interpolator and padDistance...how can I link them so a reasoning similar to this (see the part where it says "That is, for each 5 pixels in the original image, the interpolated image has 6 pixels")

3. What does padDistance = 5 mean? Sorry if the question is trivial, but I don't understand why would I need to add so many values outside my bounding box

4. More subtle question: is resampling/pixel spacing step mandatory? I have several MRI images all extracted from the short-axis view but we just selected, case by case, the most significant slice (where the lesion was most evident). What does the literature say about resampling?

5. Trivial question: when I compute the shape features, Flatness and LeastAxis are equal to zero. Is this because I only have 2D images and not 3D volumes? Which leads to the more serious question....do you have a reference (or could you briefly explain to me) the correlation between those features and the eigenvalues? I don't understand what is the difference between major, minor and least axes?

Thanks a lot in advance for your great work!
Best,
Tommaso

[JoostJM]

For instance, what it the meaning of "resampledPixelSpacing"? I get that it has to be a list of 3 floats (x,y,z) and that it represents the size of the voxel but....which is the order of magnitude? which is the unit of measurement? (e.g. what does 4,4,4 would mean?)

Spacing in mm

Same thing for interpolator and padDistance...how can I link them so a reasoning similar to this (see the part where it says "That is, for each 5 pixels in the original image, the interpolated image has 6 pixels")

interpolator specifies the method to use to calculate new pixel values (which order of funtion)
To reduce computation time, both for resampling and subsequent feature extraction, only the part of the image that is of interest is resampled. This area is determined by the bounding box of the mask, with an additional padding (specified in padDistance) to ensure correct application of some filters (which can require pixel intensities outside of the mask)

What does padDistance = 5 mean? Sorry if the question is trivial, but I don't understand why would I need to add so many values outside my bounding box

padDistance is the number of pixels to pad in the new image space (i.e. if you resample to (3, 3, 3) and pad 5, the size of the cropped region will increase by 2 * 5 * 3 = 30 mm in all directions, compared to bounding box size). As mentioned above, this is necessary for some filters (i.e. LoG, Wavelet and LBP2D/LBP3D)

More subtle question: is resampling/pixel spacing step mandatory? I have several MRI images all extracted from the short-axis view but we just selected, case by case, the most significant slice (where the lesion was most evident). What does the literature say about resampling?

Depends on how you do your extraction and how your dataset looks like. 2 main rules:

1. If you extract features in 3D (not your case as you apply single slice) you need to ensure that either the voxels are isotropic, or you take the different distances-to-neighbor into account (see weightingNorm).
   To ensure isotropic voxels, you can utilize resampling
2. If you have varying voxelsizes (in your case, in-plane voxel sizes) you resampling is mandatory, as many features are depedent on the number of voxels in the ROI and therefore, by extension, on the voxelsize. You need to harmonize these to reduce the confounding by acquisition protocol.

Additionally, resampling can also be used to focus on more coarse structures (when you use large resampled sizes), which of course describes a different texture than features focussing on fine structures.

Trivial question: when I compute the shape features, Flatness and LeastAxis are equal to zero. Is this because I only have 2D images and not 3D volumes? Which leads to the more serious question....do you have a reference (or could you briefly explain to me) the correlation between those features and the eigenvalues? I don't understand what is the difference between major, minor and least axes?

Yes they are 0 because your ROI is flat. Look at the documentation of these features. Lambda means eigenvalue here, i.e. major axis is 4 * sqrt (largest eigenvalue). If I'm correct, these axis features are the lengths of the enclosing ellipsoid of the ROI.
Elongation and Flatness are ratios between these axes. NB. both Elongation and Flatness are the inverse of the more common definition, this is to prevent division by 0.

[tommydino93]

Ok thanks a lot!!

If you have varying voxelsizes (in your case, in-plane voxel sizes) you resampling is mandatory

How do I know if they are varying?...I mean...Which is the easiest way to see ImagePositionPatient, ImageOrientationPatient, PixelSpacing, etc.?

[JoostJM]

Run a simple extraction for your data (e.g. with just firstorder features). PyRadiomics includes the original spacing in the output by default (general_info_Spacing)

[tommydino93]

Yeah, that's what I thought...but I'm not getting them. My code is:

from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
import SimpleITK as sitk
import six
from radiomics import firstorder, getTestCase, glcm, glrlm, glszm, imageoperations, shape, getImageTypes, \
    getFeatureClasses, getParameterValidationFiles, featureextractor, generalinfo
import os

datadir = "C:\\Users\\tommaso\\Desktop\\Lavoro\\Zurich_Fellowship\\Programmi\\All_Data_nrrd"
imageName = os.path.join(datadir, 'brain_image.nrrd')
maskName = os.path.join(datadir, 'brain_label.nrrd')

image = sitk.ReadImage(imageName)
mask = sitk.ReadImage(maskName)

applyLog = False
applyWavelet = False

#settings for the feature calculation. Currently resampling is disabled.
settings = {'binWidth':25, 'interpolator': sitk.sitkBSpline, 'resampledPixelSpacing':None}

if settings['interpolator'] is not None and settings['resampledPixelSpacing'] is not None:
  image, mask = imageoperations.resampleImage(image,mask,settings['resampledPixelSpacing'],settings['interpolator'])
else:
  bb,correctedMask = imageoperations.checkMask(image,mask) #bb stands for bounding box
  if correctedMask is not None: #if the mask has been correctly resampled to the input image geometry
    mask = correctedMask

  cropped_image,cropped_mask = imageoperations.cropToTumorMask(image,mask,bb)

#let's now store the image in the directory just to check whether if the cropping was successful
#sitk.WriteImage(cropped_image, "image_cropped.nrrd") #FUNZIONA, è effettivamente l'immagine croppata e quadrata

#let's now show the first order feature calculations
firstOrderFeatures = firstorder.RadiomicsFirstOrder(cropped_image, cropped_mask, **settings) #la notazione ** serve per passare i valori del dict settings
firstOrderFeatures.enableAllFeatures()
firstOrderFeatures.calculateFeatures()

print('Calculated first order features: ')
for (key, val) in six.iteritems(firstOrderFeatures.featureValues):
  print('  ',key, ':', val)

and as output I get:

C:\Users\tommaso\PycharmProjects\Prova\venv\Scripts\python.exe D:/Getting_Started_PyCharm/Pyradiomics/Prova_Github.py
Calculated first order features: 
   10Percentile : 1215.0
   90Percentile : 5548.7
   Energy : 44379355767.0
   Entropy : 7.47192856194
   InterquartileRange : 3496.75
   Kurtosis : 1.57285375225
   Maximum : 6564.0
   MeanAbsoluteDeviation : 1569.92729285
   Mean : 3668.53786192
   Median : 4239.0
   Minimum : 3.0
   Range : 6561.0
   RobustMeanAbsoluteDeviation : 1348.28049435
   RootMeanSquared : 4058.74438726
   Skewness : -0.339285695793
   TotalEnergy : 44379355767.0
   Uniformity : 0.00730788923556
   Variance : 3015235.95681

Process finished with exit code 0

What am I missing?
Thanks again in advance!

[JoostJM]

@tommydino93 You are missing that data because you are implementing the featureclasses directly, which is not advised. Use the featureextractor interface, it's a lot easier to use (follow the helloRadiomics.py example)

[tommydino93]

@JoostJM Alright! I've done that, thanks.

It appears that my pixel spacing is now (1.0, 1.0, 1.0) which is weird because in my original .dcm image it was different. I figured that maybe during the conversion from DICOM to nrrd something went wrong. I followed this guide for the conversion. Here are my two images opened in ImageJ:

as you can see the .nrrd image (on the right) has lost the pixel spacing (highlighted) and is now in microns and not in mm. I also unticked the "Compress" option while doing the conversion, but the problem persists.

Is there an alternative way for converting and preserving dimensions/pixel spacing?

Thanks a lot again!

[JoostJM]

You can try this, it's a script/package I built for conversion in Python. It will scan the folder and create NRRD (or NIFTII). If there is something wrong with the image, it will give you a warning what is happening.

The spacing (1, 1, 1) is usually a sign that something went wrong during conversion.

[JoostJM]

Alternatively, you can also try to see if Slicer found some issue by checking advanced in the DICOM browser and pressing the buttom Examine. After you've loaded your volume, you can also check the characteristics of that volume in the volumes model.

[tommydino93]

Thanks again @JoostJM !

I am having problems with PyCharm..maybe you can help me. I think it's a problem of configuration between PyCharm and GitHub. I typed this simple code:

import logging
import os
import pydicom
import SimpleITK as sitk
import tqdm

Source = "C:\\Users\\tommaso\\Desktop\\Lavoro\\Zurich_Fellowship\\Programmi\\Images_dcm"
Destination = "C:\\Users\\tommaso\\Desktop\\Lavoro\\Zurich_Fellowship\\Programmi\\Images_nrrd"
walk_folder(Source,Destination)

and I get as error:

Traceback (most recent call last):
  File "D:/Getting_Started_PyCharm/Pyradiomics/Prova_Conversion_to_nrrd.py", line 9, in <module>
    walk_folder(Source,Destination)
NameError: name 'walk_folder' is not defined

1. Did I call the function properly?
2. Why walk_folder is not recognized? I cloned the GitHub HTTPS in PyCharm as explained here

[JoostJM]

@tommydino93, You'll have to import the module before you can use it:

from nrrdify import walk_folder

That being said, Nrrdify also has a commandline interface. To use it, it is most easy to install it.
from the root of the package (i.e. the folder containing setup.py)

> python setup.py install
> nrrdify --help

This will print the help message detailing the usage (with optional configuration). Once it is installed, it should be usable throughout your computer. E.g.

C:\Users\tommaso\Desktop\Lavoro\Zurich_Fellowship\Programmi> nrrdify Images_dcm -o Images_nrrd

This will convert all dicoms in the folder above into nrrd and store them in Images_nrrd

[tommydino93]

Wow thanks a lot!
I managed to do it both ways. Nevertheless, the .nrrd image created is still a 512x512 microns and not a 260x260 mm as the original DICOM one :(

Could it be a problem just with the opening in ImageJ?

[JoostJM]

@tommydino93, could be. Can you share an anonymised version of your DICOMs and NRRD?

[tommydino93]

@JoostJM, I managed to anonymize only the DICOM image (but not the .nrrd) with DicomCleaner. I hope it is sufficient for you to try it out :) Here it's the MEGA link to the image:

https://mega.nz/#!gPhHRayB

Thanks a lot again!

[tommydino93]

Sorry, that link needs a key! This already has the key:

https://mega.nz/#!gPhHRayB!JTen715bznaNqNxIzNKez_Wutb0M5RIZv1vS0_tg9BQ

[JoostJM]

As far as I can see, it appears to be caused by ITK, as it also shows the incorrect spacing of (1, 1, 1) in ITKsnap. When I looked at the tags, spacing was defined as [5.0781200E-001\5.0781200E-001], it is possible ITK is failing on the E-001, but I'm not sure. Also, I noticed the SOP class is secondaryCaptureImageStorage (usually the case for key images), do you have access to the original series? Can you try to convert that one?

As to why imageJ gives you 512 microns, I do not know. Maybe a bug of some kind? the NRRD image says the spacing is 1, 1, 1, and does not encode the Field of View size (as it can easiliy be obtained by multiplying the spacing * the matrix size.

On potential explanation is that it just does not know how to cope with the fact that there is no unit encoded in NRRD and it sort of assumes its microns (in which case 1 x 512 is indeed 512 microns).

Here are also some pointers to other tools you could try. @fedorov is currently also busy with the DICOM to NIFTII/NRRD conversion problem, this is from a project he started at NAMIC project week in Boston last January:

dcm2niix
https://github.com/rordenlab/dcm2niix

dicom2nifti
https://github.com/icometrix/dicom2nifti

dcmstack
https://github.com/moloney/dcmstack

vtk-dicom/dicomtools/dicomtonifti
https://github.com/dgobbi/vtk-dicom

FreeSurfer/mri_convert
https://surfer.nmr.mgh.harvard.edu/fswiki/mri_convert

Plastimatch/convert
http://plastimatch.org/plastimatch.html

mriconvert and mcverter
https://lcni.uoregon.edu/downloads/mriconvert/mriconvert-and-mcverter

[tommydino93]

Hi @JoostJM and @fedorov !
Unfortunately I don't have access to the original images, because images are not exportable from IMPAX (the software they use here as PACS). I had to load the images from IMPAX to Siemens SyngoVia (a dedicated Siemens software for visualizing images) and then I extracted images from there.

Anyway, I will try the "dicom2nifti" way, hoping to solve the issue. I'll keep you updated!

Thanks again!

[tommydino93]

@fedorov , I have been trying this simple code:

import dicom2nifti

dicom_directory = "C:\\Users\\tommaso\\Desktop\\Prova_Loading"
output_folder = "C:\\Users\\tommaso\\Desktop\\Outut_Images_Nifti"

dicom2nifti.convert_directory(dicom_directory, output_folder, compression=False, reorient=False)

but no images appear in my output directory "Outut_Images_Nifti".

Could this be because all my dicom images are uncompressed? @JoostJM Could this be the problem also for the other issues I had? I mean...Do I need to extract the images as COMPRESSED rather than UNCOMPRESSED?

thanks a lot guys in advance!

[JoostJM]

@tommydino93, as far as I know, uncompressed or compressed should not matter, if anything, uncompressed should be easier. I also used this program, for the example you sent, it is able to correctly read in the DICOM and you can re-save in NIFTII (which is also accepted by PyRadiomics).

[tommydino93]

@JoostJM , @fedorov ! Thanks a lot for advising me Mango. Actually, I have great news: while doing some trials between Mango and ImageJ, I discovered that ImageJ itself allows to convert DICOM to nrrd and the pixel spacing is preserved!! We can just simply do:

File>Save As>Nrrd

And this works both for the image and for the mask. I don't know why I didn't think about it before.

Anyway, before converting all images, I just tried to load one image and one mask simply typing:

datadir = "C:\\Users\\tommaso\\Desktop\\Prova_PyCharm"
imageName = os.path.join(datadir, '87951608.nrrd')
maskName = os.path.join(datadir, 'mask.nrrd')

image = sitk.ReadImage(imageName)
mask = sitk.ReadImage(maskName)

applyLog = False
applyWavelet = False

#settings for the feature calculation. Currently resampling is disabled.
settings = {'binWidth':25, 'interpolator': sitk.sitkBSpline, 'resampledPixelSpacing':None}

if settings['interpolator'] is not None and settings['resampledPixelSpacing'] is not None:
  image, mask = imageoperations.resampleImage(image,mask,settings['resampledPixelSpacing'],settings['interpolator'])
else:
  bb,correctedMask = imageoperations.checkMask(image,mask) #bb stands for bounding box
  if correctedMask is not None: #if the mask has been correctly resampled to the input image geometry
    mask = correctedMask

  cropped_image,cropped_mask = imageoperations.cropToTumorMask(image,mask,bb)

#let's now store the image in the directory just to check whether if the cropping was successful
#sitk.WriteImage(cropped_image, "image_cropped.nrrd") #FUNZIONA, è effettivamente l'immagine croppata e quadrata

# Initialize feature extractor
extractor = featureextractor.RadiomicsFeaturesExtractor(**settings)
#enable all features
extractor.enableAllFeatures()

featureVector = extractor.execute(cropped_image, cropped_mask)

for featureName in featureVector.keys():
  print("Computed %s: %s" % (featureName, featureVector[featureName]))

and then I tried to run the shape feature program but it throws me this error:

C:\Users\tommaso\PycharmProjects\Prova\venv\Scripts\python.exe D:/Getting_Started_PyCharm/Pyradiomics/Feature_Extraction_New.py
Traceback (most recent call last):
  File "D:/Getting_Started_PyCharm/Pyradiomics/Feature_Extraction_New.py", line 14, in <module>
    image = sitk.ReadImage(imageName)
  File "C:\Users\tommaso\PycharmProjects\Prova\venv\lib\site-packages\SimpleITK\SimpleITK.py", line 8332, in ReadImage
    return _SimpleITK.ReadImage(*args)
RuntimeError: Exception thrown in SimpleITK ReadImage: C:\d\VS14-Win64-pkg\SimpleITK-build\ITK\Modules\IO\NRRD\src\itkNrrdImageIO.cxx:249:
itk::ERROR: NrrdImageIO(000001E932882740): ReadImageInformation: Error reading C:\Users\tommaso\Desktop\Prova_PyCharm\87951608.nrrd:
[nrrd] nrrdLoad: trouble reading "C:\Users\tommaso\Desktop\Prova_PyCharm\87951608.nrrd"
[nrrd] nrrdRead: trouble
[nrrd] _nrrdRead: trouble reading NRRD file
[nrrd] _nrrdFormatNRRD_read: trouble parsing space directions info |(0.89843749,0,0) (0,0.89843749,0)|
[nrrd] _nrrdReadNrrdParse_space_directions: trouble getting space vector 1 of 2
[nrrd] _nrrdSpaceVectorParse: space dimension is 2, but seem to have 3 coefficients

From what I understand, it's a spacing problem, which leads us back to my initial question: "how could I set the resampling among all my images, since they have different spacing"? E.g. should I choose the average pixel spacing among all images? Or maybe the smallest one and apply it to all other images?

Thanks a lot again, I really appreciate your help!

[JoostJM]

255 is accepted, but you'll need to specify it in the label parameter

[tommydino93]

Ok thanks a lot, I'll keep that in mind!
In the following mega links, the nrrd image and mask that I've extracted through ImageJ

Image: https://mega.nz/#!oLQVFYZZ!PNs-pOWlnXWM9Dau5kwrwYWnsdys4QZjka8uZxAh7ek
Mask: https://mega.nz/#!wDY0gTDI!88ZEMhDpHNkJZDfYE3o9DoRthtGi6BpxsWIMtKmU-14

In the ImageJ menu: Image>Show Info it says (like in the error that PyCharm gave me) that the space direction is (0.89843749,0,0) (0,0.89843749,0) but I still don't understand why this is a problem for ITK

Thanks in advance!

[fedorov]

@fedorov , I have been trying this simple code:

You don't need to write any code to use dicom2nifti, it provides the command line tool:

$ dicom2nifti --help                                                                                                                                                                     2.3.6
usage: dicom2nifti [-h] [-G] [-M] [-C] [-R] input_directory output_directory

dicom2nifti, convert dicom files into nifti format.

Most of the tools that @JoostJM referenced above in #360 (comment) provide command line interface. Using a GUI-based interface works when you need to convert just a few datasets, but it does not scale.

If I were you, I would try several command line tools (dicom2nifti, plastimatch, and dcm2niix are good ones), and if any of those do not work, it would be good to document this and ideally submit a de-identified dataset, if possible at all, so that the issue can be fixed in the converter.

[tommydino93]

Thanks a lot @fedorov ! You're right, I should definitely become more agile with the command line interface. It's way faster.

Anyhow, in the end I managed to convert the images from DICOM to nrrd with ImageJ, but, as I was explaining JoostJM in my previous comment, now I have a different problem with ITK which is unable to read my nrrd images, probably because of the space direction.

Any suggestion?
Thanks again for your amazing work!

[fedorov]

Can you open the NRRD file in the text editor, and copy paste the header here?

[tommydino93]

Yes, sure. Here it is:

header.txt

[fedorov]

I am not an NRRD expert, but it looks to me like an invalid header. Indeed, it looks inconsistent to have 3 components for the direction vectors, but only 2 vectors, and only 2 dimensions of the image.

Is your dataset volumetric, or is it a single slice? If it is supposed to be volumetric, I would really recommend you go back and try other converters. It maye well be that ImageJ NRRD writing may have issues. You could also ask ImageJ support.

If you want to just try something, you can maybe manually modify this line

space directions: (0.89843749,0,0) (0,0.89843749,0)

to this:

space directions: (0.89843749,0) (0,0.89843749)

Maybe it will fix the ITK issue ...

[tommydino93]

The dataset is not volumetric. We just selected a single slice for each patient (the slice where the inflamed region was most evident)...so what should I expect for the third dimension in the header?

[fedorov]

you can also try this header, but it's just a guess for a hack...

type: uint16
encoding: raw
endian: big
dimension: 3
sizes: 512 512 1
space dimension: 3
space directions: (0.89843749,0,0) (0,0.89843749,0) (0,0,1)
space units: "mm" "mm" "mm"

[JoostJM]

@tommydino93, @fedorov, I agree that the header looks weird. Additionally, changing the directions to tuples of 2 will mean that PyRadiomics does not accept it (as it only accepts volumetric sets, even if only 1 slice).

I agree with @fedorov that it is best to try to use some different tools, but if you want to manually edit the header, do so according to @fedorov's last post.

[fedorov]

@tommydino93 if you will go back to try other converters, you might also consider converting the 3d dataset. You do not need to have just one image slice if you want to extract features from just that slice. You can just define your label in one slice, but keep the image volumetric, so that if/when you decide to look at the analysis beyond single slice, you do not need to go back and deal with the data conversion again.

[tommydino93]

Thank you both @JoostJM @fedorov !
For now I think I'll stick with the single slices just because I'm very late with the project.
I am trying to run dicom2nifti from command line and I'm having problems calling the help. Sorry guys if this is a very trivial question but I'm a newbie in command line programming. What I've done so far is:

1. I have unzipped from GitHub the dicom2nifti directory to D:\Getting_Started_PyCharm\dicom2nifti-master
2. From cmd, I've navigated to that directory with: D: and then cd Getting_Started_PyCharm\dicom2nifti-master
3. Once inside the directory I've typed python setup.py install
4. Then I've typed pip install dicom2nifti
5. Then I've typed dicom2nifti -h but this throws me the warning dicom2nifti is not recognized as an internal or external command

I just wanted to understand what the various [-h] [-G] [-r] [-o RESAMPLE_ORDER] [-p RESAMPLE_PADDING] [-M] [-C] [-R] mean in the line:

dicom2nifti [-h] [-G] [-r] [-o RESAMPLE_ORDER] [-p RESAMPLE_PADDING] [-M] [-C] [-R] input_directory output_directory

Thanks in advance!

[tommydino93]

HI @JoostJM ,
I think I found a work around my problem: I went back to the conversion from DICOM to nrrd with 3D slicer and I noticed, also asking in the 3D slicer forum, that 3D slicer automatically reshapes the pixel spacing to (1,1,1) for all images, which basically solves my problem. I mean....the important thing is that all images have the same pixel spacing right? In this way features extracted from them are comparable. Let me know, whenever you have time, if this reasoning makes sense :)

Furthermore, I was also trying to extract features from the LoG filtered image through the feature extractor. I tried to use this code from one of the examples:

# Initialize feature extractor
extractor = featureextractor.RadiomicsFeaturesExtractor(**settings)
#enable all features
extractor.enableAllFeatures()

featureVector = extractor.execute(cropped_image, cropped_mask)

extractor.enableAllImageTypes() #enable all possible filters (e.g. LoG, Wavelet, Logarithm, etc.)

for decompositionImage, decompositionName, inputKwargs in imageoperations.getWaveletImage(cropped_image, cropped_mask):
    waveletFirstOrderFeaturs = firstorder.RadiomicsFirstOrder(decompositionImage, cropped_mask, **inputKwargs)
    waveletFirstOrderFeaturs.enableAllFeatures()
    waveletFirstOrderFeaturs.calculateFeatures()
    print('Calculated firstorder features with wavelet ', decompositionName)
    for (key, val) in six.iteritems(waveletFirstOrderFeaturs.featureValues):
      waveletFeatureName = '%s_%s' % (str(decompositionName), key)
      print('  ', waveletFeatureName, ':', val)

but I get this error:

File "D:/Getting_Started_PyCharm/Pyradiomics/Feature_Extraction.py", line 43, in <module>
    for decompositionImage, decompositionName, inputKwargs in imageoperations.getWaveletImage(cropped_image, cropped_mask):
TypeError: getWaveletImage() takes 1 positional argument but 2 were given

Shouldn't the function take 2 parameters as input? Why it's only expecting 1 positional argument?

Thanks a lot in advance!

[QiChen2014]

Hi @JoostJM , I come again. Is there a way to check if the image and mask is under the same slice? If it's possible, how can I check it? Thanks a lot!

[QiChen2014]

I have another question: when I debug the code to extract features, I see the boundingBox is

but after excute the below code segment, it changed. I'm puzzled, could you tell me what has done inner the program? Look forward to your reply and thank you again.

[JoostJM]

@QiChen2014

Is there a way to check if the image and mask is under the same slice? If it's possible, how can I check it

This is done automatically by PyRadiomics. If you want to know more, check out this function (checkmask())

I debug the code to extract features, I see the boundingBox is but after excute the below code segment, it changed

The difference here is how the bounding box is defined. I use 2 separate functions for that. Internally, the bounding box is handled as a tuple specifying the lower and upper bounds. However, the bounding box that is returned as part of the provenance info is just the lower bounds as the first 3 elements. The last 3 elements there specify the size of the bounding box. E.g. for the first dimension, x, the lower bound is 213 (first element in both), the upper bound is 278 (2nd element in the internal bounding box) and the size is 278 - 213 + 1 = 66 (4th element in the provenance-returned bounding box)

[QiChen2014]

I understand, thank you @JoostJM , I am very grateful!

[tommydino93]

Hi everyone! With the help of @fedorov and @JoostJM, part of the resampling problem was solved in this 3D slicer forum question.

Going back to PyRadiomics now: let's say I want a resampledPixelSpacing = [0.625, 0.625, 0]
Following the helloFeatureClass.py example, I wrote this on my PyCharm code:

settings = {
  'binWidth': 25,
  'interpolator': sitk.sitkBSpline,
  'resampledPixelSpacing': [0.625, 0.625, 0]
}

if settings['interpolator'] is not None and settings['resampledPixelSpacing'] is not None:
    image, mask = imageoperations.resampleImage(image, mask, settings['resampledPixelSpacing'], settings['interpolator'])
else:
    bb, correctedMask = imageoperations.checkMask(image, mask)  # bb stands for bounding box
    if correctedMask is not None:  # if the mask has been correctly resampled to the input image geometry
        mask = correctedMask

    cropped_image, cropped_mask = imageoperations.cropToTumorMask(image, mask, bb)

My question is:
"since the else is not computed, checkMask and cropToTumorMask are skipped...Isn't this a problem?"

[tommydino93]

I have answered myself. By looking at the helloResampling.py it is explained that resampleImage already resamples and crops onto bounding box defined by the label!

[tommydino93]

Hi @JoostJM!
Sorry for re-opening this after some time. I just had a quick question. It's about when you said:

"If you extract features in 3D you need to ensure that either the voxels are isotropic, or you take the different distances-to-neighbor into account (see weightingNorm). To ensure isotropic voxels, you can utilize resampling"

my questions are:

1. what is the problem of extracting 3D features from the stack of images resampled only along xy (and not along the z axis)?

2. Also, what does the resampling along the z axis do? I am having difficulties visualizing/understanding it.

3. would the 3D feature extraction be wrong if using non isotropic voxels or just not optimal?

Thanks again a lot for your support :)
Tommaso

[JoostJM]

@tommydino93 No problem!

1. Computationally, there is no problem. However, unless weighting by distance (weightingNorm) or force2D is used, PyRadiomics considers neighbors in-plane and out-of-plane to be comparable. If, however, z is e.g. 5 mm and xy are 1 mm, neighbours in z are 5 times further away, and therefore not exactly comparable. This then violates the comparable-neighbors assumption, upon which the feature formulas are based. There are 3 options in PyRadiomics for dealing with this:
   1. Enable resampling to ensure isotropic voxels (see below).
   2. Enable weighting, i.e. taking the distance into account when building the feature matrix.
   3. Enable force2D, ignoring ofsets out-of-plane, thereby only considering neighbors that are in-plane and therefore comparable.
2. With resampling, you define new gridpoints in your volume, and then calculate the values at those points using interpolation (different interpolators exist, such as linear, nearest neighbor and b-spline). You can best visualize this as your volume being a continuous space, and the values of voxels being the value of discreet points in that space. Resampling means you just 'sample' the volume at different discreet locations.
3. This is more a philosophical question. I think I'd go for "not optimal", as there will still be valuable information extracted. It mostly reduces the rotation invariance of the extracted features.

[tommydino93]

@JoostJM Thanks a lot for the immediate answer and detailed explanation!!

So, suppose I wanted to go for option:

i. Enable resampling to ensure isotropic voxels

should I just set the three values of resampledPixelSpacing to the same number like below?

        # settings for the feature calculation.
        settings = {
          'binWidth': 25,
          'interpolator': sitk.sitkLinear,
          'resampledPixelSpacing': [0.8, 0.8, 0.8]  # in [mm]
        }

Also, completely off-topic question:
is there, to your knowledge, a gold standard on if/how many filters to apply to the original images for extracting features? I mean, is it standard, optional, rare?

Thanks again :)

[JoostJM]

@tommydino93 Yes, that enables resampling to isotropic voxels. On what size to use, I usually advise to take a compromise with deleting information in-plane and 'creating' it out-of plane. E.g. in the case of x,y = 1mm, z = 5 mm, I usually go for 2 or 3mm isotropic. Always keep in mind that larger voxels are not necessariliy bad, as small voxels are more sensitive to noise. On the other hand, if your lesions are very small, large voxels may result in ROIs with only few number of voxels, which is also unstable.

as to your off-topic question. No gold standard. There is not even a standard (yet) on which to apply and how. Generally some are usually applied and most good radiomics packages do implement at least some.

[tommydino93]

Thank you very much again! I will open a new quick issue since it's a completely different topic.

[harun-cu]

@JoostJM @fedorov
Hi,
Is there any way to convert .tiff image to .nifti. My original images are in DICOM but after creating mask using ImageJ it saved as .tiff.
Would you suggest me anyway?

Thanks

[JoostJM]

The easiest way to do it in python is to load both image and mask in SimpleITK, then set the image as the reference volume for the mask, the save the mask as niftii. One question though, is your image 2d or 3d? I believe tiff will only save in 2d, so you'd need to specify which slice it is.
Additionally, a good thing to do is to enter a vector index from the mask if it has a vector pixel type (I.e. color image). If it appears white, most likely all color channels have the value aet to 255, so the index doesn't matter.