Need of new voxel calculation in DICOM example

[Jugaris89]

Hi,

first of all thanks for this framework, it is very useful.

I had a question regarding the notebook for the DICOM example under "./examples/LCTSC_DICOMInterface.ipynb"

what is the need of the "Find out good voxel spacing" part? I see, the "x", "y" and "z" calculated there are not used in the rest of the notebook?

and also, when the resampling is done on next step, where do these value come from? "sf_resample = Resampling((3.22, 1.62, 1.62))"

thanks in advance

[MLRadfys]

Hi,

I think the main purpose is to show the relationship between patch size and resampled volume. Like Dominik wrote in the example:

"A good patch shape <-> median image ratio is 1/8 up to 1/64."

So its just to verify that your ratio is feasible.

3.22, 1.62, 162 is just one choice for the voxel size. You could also resample the volumes to be isotropic (1,1,1), which is often seen in different papers.

Cheers,

M

[muellerdo]

Hey @Jugaris89,

what is the need of the "Find out good voxel spacing" part? I see, the "x", "y" and "z" calculated there are not used in the rest of the notebook?

Voxel spacings are normally inhomogeneous across samples. This makes it quite hard for the model to learn a common pattern if the sizes between voxels are not uniform.

The general idea is that you normalize these voxel spacing to a uniform ratio across all samples.

As @MichaelLempart correctly stated, many papers normalizes the voxel spacings to be isotropic (1,1,1), which is totally fine and sufficient.

However, in my experience, you will gain a performance increase if you watch out that your patch size takes up a big enough part of the median image shape (after the voxel space normalization) of your samples. I recommend a patch <-> median image ratio of about 1/8 up to 1/64, in order to ensure that the patch contains a large part of the image, but still in a high enough resolution.

Be aware, that the voxel space normalization directly affect the image size and therefore your image resolution. The task here is to find a balance between:

• How much information (area of the image) is available in my patch vs
• How good is the quality (resolution) of the information in my patch

Cheers, Dominik

[Jugaris89]

thank you for your responses.

If I may, I have a couple of short questions in which I hope you could help me too please:

• In the visualization of the predicted classes (utils/visualizer.py) you consider only 2 classes to visualize: I understand this is because in your example you are segmenting both lungs but if you would be segmenting more than 2 organs, you would have to extend these lines for as many classes as you are segmenting, right?

• I am adpating your framework to segment the pancreas by labeling it in the same LCTSC dataset you use for lungs. I am complementing this dataset, with own images and the ones from TCIA (the cancer imaging archive), where I labeled the pancreas too. I am following a similar approach as you propose for the lungs, but changed the patch-size and voxel size. However, when I include the "resample" subfunction to the set of subfunctions, I get always the following error (see at the end of my message). That does not happen when I train only with LCTSC images. Have you already experienced this error before? Could you give me a hint on why this may be happening? It happens only with "resample", i can include "clipping" and "normalize", and no error comes.

Thanks in advance

error I am getting:

File "pancreas_seg.py", line 96, in draw_figures=True, run_detailed_evaluation=True) File "/data/MIScnn/miscnn/evaluation/split_validation.py", line 69, in split_validation iterations=iterations, callbacks=callbacks) File "/data/MIScnn/miscnn/neural_network/model.py", line 205, in evaluate max_queue_size=self.batch_queue_size) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/training.py", line 66, in _method_wrapper return method(self, *args, *kwargs) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/training.py", line 872, in fit return_dict=True) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/training.py", line 66, in _method_wrapper return method(self, args, kwargs) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/training.py", line 1057, in evaluate model=self) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/data_adapter.py", line 1112, in init model=model) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/data_adapter.py", line 908, in init kwargs) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/data_adapter.py", line 772, in init peek, x = self._peek_and_restore(x) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/tensorflow/python/keras/engine/data_adapter.py", line 912, in _peek_and_restore return x[0], x File "/data/MIScnn/miscnn/neural_network/data_generator.py", line 63, in getitem else : batch = self.generate_batch(idx) File "/data/MIScnn/miscnn/neural_network/data_generator.py", line 147, in generate_batch self.validation) File "/data/MIScnn/miscnn/processing/preprocessor.py", line 123, in run sf.preprocessing(sample, training=training) File "/data/MIScnn/miscnn/processing/subfunctions/resampling.py", line 71, in preprocessing order=3, order_seg=1, cval_seg=0) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/batchgenerators/augmentations/spatial_transformations.py", line 70, in augment_resize sample_data = resize_multichannel_image(sample_data, target_size_here, order) File "/home/ubuntu/anaconda3/envs/tensorflow2_p36/lib/python3.6/site-packages/batchgenerators/augmentations/utils.py", line 617, in resize_multichannel_image result = np.zeros(new_shp, dtype=multichannel_image.dtype) ValueError: negative dimensions are not allowed

[muellerdo]

Hey @Jugaris89,

sorry for the late response.

In the visualization of the predicted classes (utils/visualizer.py) you consider only 2 classes to visualize: imagen I understand this is because in your example you are segmenting both lungs but if you would be segmenting more than 2 organs, you would have to extend these lines for as many classes as you are segmenting, right?

You are right. The provided functionality for visualization is, sadly, still highly experimental and for a long time now on my rework agenda. It is only functional for 3D data with one channel (grayscale, HU, ...) and, as you correctly noticed, for 3 classes (including background). I would recommend copy&pasting the code and just modifying the classes to your needs.
If you implement a more dynamic function for visualization, you are more than welcome to contribute your code! :+1:

I am adpating your framework to segment the pancreas by labeling it in the same LCTSC dataset you use for lungs. I am complementing this dataset, with own images and the ones from TCIA (the cancer imaging archive), where I labeled the pancreas too. I am following a similar approach as you propose for the lungs, but changed the patch-size and voxel size. However, when I include the "resample" subfunction to the set of subfunctions, I get always the following error (see at the end of my message). That does not happen when I train only with LCTSC images. Have you already experienced this error before? Could you give me a hint on why this may be happening? It happens only with "resample", i can include "clipping" and "normalize", and no error comes.

The error arises from wrong/not logical voxel spacings/slice thickness values in the DICOM files.

Would it be possible that you run the following code and share the output with us:

# Initialize your Data IO class on your data set
data_io = ...

# Obtain sample list
sample_list = data_io.get_indiceslist()

# Print out voxel spacings
for index in sample_list:
   sample = data_io.sample_loader(index)
   print(sample.index, sample.details)

What I suspect is that your new data have voxel spacings encoded in negative values. I already noted negative voxel spacings in NIfTI data, but haven't saw ones in DICOM data. Therefore, the DICOM interface still misses a line in which the spacings will be transformed to their absolute values.

If you can confirm my supposition, I will add the absolute transformation of voxel spacings in the DICOM interface and release a new version of MIScnn. Then you can just update your version and continue working.

Cheers, Dominik

EDITED: Fixed a little bug in the code section. The variable sample_list contains the indices and not the sample objects -> added a line to load the samples given the index.

[Jugaris89]

Hi Dominik,

thanks for your response. Here the voxel spacing after including your lines of code:

(tensorflow2_p36) ubuntu@ip-...:/data/MIScnn$ python lungs_seg.py Found structures in sample : [["1", 'SpinalCord'], ["2", 'Lung_R'], ["3", 'Lung_L'], ["4", 'Heart'], ["5", 'Esophagus'], ["6", 'pancreas']] {'pancreas': 1, 'Esophagus': 2} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([2. , 1.171875, 1.171875], dtype=float32)} {'spacing': array([3. , 1.269531, 1.269531], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2.5 , 0.976562, 0.976562], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)} {'spacing': array([2. , 1.171875, 1.171875], dtype=float32)} {'spacing': array([3. , 0.9765625, 0.9765625], dtype=float32)}

In the end, I reached the same conclusion as you, that the before commented issue was due to a unproper configuration of the voxel or slice thikness. The negative voxel spacing comes when I include DICOM files from this public dataset: https://wiki.cancerimagingarchive.net/display/Public When I use the dataset you use in the example, I think it does not come up.

Regarding the visualization, I am adding line for up to 4 classes (which are the ones I need).

May I ask you your opinion on another thing, please: I am having problems finding the pancreas. I labeled it in the dataset you use for the lungs (and the labeling is correct, because I visualize correctly). I have configured the proper HU factor in the clipping subfunction but still in the prediction, it finds the esophagus or other parts, but not the pancreas. Do you think this is because the pancreas is present in very few slices from the CT scan? Since it is a small organ, do you think I shall make small voxels?

thanks again for your help

[muellerdo]

I added the absolute transformation of voxel spacing in the DICOM IO interface. Could you update to the latest MIScnn version 1.0.4 and try it again.

pip install miscnn --upgrade

Please keep me updated if it works or not, now.

May I ask you your opinion on another thing, please: I am having problems finding the pancreas. I labeled it in the dataset you use for the lungs (and the labeling is correct, because I visualize correctly). I have configured the proper HU factor in the clipping subfunction but still in the prediction, it finds the esophagus or other parts, but not the pancreas. Do you think this is because the pancreas is present in very few slices from the CT scan? Since it is a small organ, do you think I shall make small voxels?

What dice score do you get for the pancreas? What min & max values are you currently using for clipping?
And what resampling ratio are you using to normalize voxel spacing?

If your data set has a similar voxel spacing as the LCTSC data, I would definitely stay with the resampling to 3.22x1.62x1.62 and a patch shape of 80x160x160. Maybe, even increase the ratio between patch <-> resampled image a little bit more. For clipping I would start using a HU range of -100 to +250 for pancreas identification, maybe even decrease it to -50 to +200.

It is hard to give advice out of the blue, because I personally like to play around a bit with my dataset to see what works and what not. The size of the pancreas is definitely a challenge in detection, but the fact that it is hopefully present in all patients with about the same shape and location, it should, in my experience, be able to get a high performance.

How did you generate the segmentation mask for labelling the pancreas? Per hand-drawn or with a VOI & threshold or some other way? Be aware that most of the time, the annotated segmentation masks are refined/smoothed via some kind of anti-aliasing technique. I'm not an expert in annotation methods, but Nicholas Heller et al. wrote an interesting paper on what they did for creating their Kidney Tumor data set: https://arxiv.org/abs/1904.00445
Maybe this paper can give you some interesting insights.

Hope that I was able to help you.

Cheers, Dominik

[Jugaris89]

Hi Dominik,

thanks for the info. Regarding your questions, I answer here: What dice score do you get for the pancreas? --> It starts by 0.3... and decreases until <0.1 during training What min & max values are you currently using for clipping? --> min:-50 max:30 And what resampling ratio are you using to normalize voxel spacing? --> I tried a few, (1,1,1), smaller and the one you suggested. None provides a trustworthy prediction after a few epochs

I labeled the pancreas by hand on the same dataset your are using with the tool "Slice3D" which is opensource and from it I exported the DICOM RT-STRUCT file.

The I use the function "split_validation(sample_list, model, epochs=40, iterations=30, draw_figures=True, run_detailed_evaluation=True)" to train it

I will take a look at the paper you suggested, in case I am doing sth wrong.

Thanks again

[muellerdo]

No problem.

What dice score do you get for the pancreas? --> It starts by 0.3... and decreases until <0.1 during training

This means that the model doesn't find it at all. Mhm.

Try out to increasing the maximum clipping value to at least +100 or +150. -> min:-50 max:+150 I have the feeling that CNN models performing the best when the clipping is not too strict, but also not too vague.

The resampling heavily depends on your dataset, gpu hardware/patch or batch size and the actual task. Personally, I always recommend using a median image ratio to patch size between 1/8 - 1/64. The voxel spacing of 3.22x1.62x1.62 and 80x160x160 patch size should work quite well on this dataset.

I labeled the pancreas by hand on the same dataset your are using with the tool "Slice3D" which is opensource and from it I exported the DICOM RT-STRUCT file. ... I will take a look at the paper you suggested, in case I am doing sth wrong.

Theoretically, it should work out of the box. Could you try another run with the increased clipping values and have a look on the dice score again, then?

[Jugaris89]

Hi,

thanks for your comment. I already tried what you suggested (I am making several trials playing a bit with everything xD). Here my configuration and below the dice score during training:

sf_normalize = Normalization(mode="z-score") sf_resample = Resampling((3.22, 1.62, 1.62)) sf_clipping = Clipping(min=-79, max=64) subfunctions = [sf_clipping, sf_normalize, sf_resample]

from miscnn import Preprocessor pp = Preprocessor(data_io, data_aug=None, batch_size=2, subfunctions=subfunctions, prepare_subfunctions=False, prepare_batches=False, analysis="patchwise-crop", patch_shape=(80, 160, 160)) ... Epoch 15/20 30/30 [==============================] - 340s 11s/step - loss: 1.6120 - dice_soft: 0.4637 - dice_crossentropy: 5.8534 - val_loss: 2.0236 - val_dice_soft: 0.3129 - val_dice_crossentropy: 5.9576 Epoch 16/20 30/30 [==============================] - 343s 11s/step - loss: 1.5663 - dice_soft: 0.4784 - dice_crossentropy: 5.7650 - val_loss: 1.8221 - val_dice_soft: 0.3879 - val_dice_crossentropy: 5.8166 Epoch 17/20 30/30 [==============================] - 377s 13s/step - loss: 1.5679 - dice_soft: 0.4775 - dice_crossentropy: 5.7276 - val_loss: 1.8259 - val_dice_soft: 0.3839 - val_dice_crossentropy: 6.0892 Epoch 18/20 30/30 [==============================] - 381s 13s/step - loss: 1.5291 - dice_soft: 0.4904 - dice_crossentropy: 5.6957 - val_loss: 1.7839 - val_dice_soft: 0.3912 - val_dice_crossentropy: 4.9475 Epoch 19/20 30/30 [==============================] - 346s 12s/step - loss: 1.5105 - dice_soft: 0.4956 - dice_crossentropy: 5.6046 - val_loss: 1.7488 - val_dice_soft: 0.4072 - val_dice_crossentropy: 5.7244 Epoch 20/20 30/30 [==============================] - 393s 13s/step - loss: 1.4639 - dice_soft: 0.5115 - dice_crossentropy: 5.6339 - val_loss: 1.7103 - val_dice_soft: 0.4108 - val_dice_crossentropy: 4.1493

The results are a bit better but still not good. For example, in this case I am looking for pancreas and Esophagus. The Esophagus is found good but the pancreas not...so I am still struggling with the proper configuration for it :S

[muellerdo]

Hey @Jugaris89,

sorry for the late reply again...

How is it going with your pancreas segmentation? Did you make progress?

Have you e.g. tried further increasing the maximum clipping range to e.g. +150? Also try to activate data augmentation like just removing data_aug parameter from the Preprocessor() in order to run data augmentation with default settings. I would also highly recommend to run prepare_subfunctions=True in the Preprocessor. With active data augmentation, it will then perform on-the-fly data augmentation which could significantly improve your performance.

Cheers, Dominik

[Jugaris89]

Hi @muellerdo ,

thaks for your help and interest.

I am still looking for the issue here. I am not able to do a proper segmentation of the pancreas still. I tried with and without data augmentation and changing the clipping range.

I am segmenting the pancreas and the esophagus both, just to see the effect of my changes in both, and the esophagus is very good segmented but not the pancreas.

I will try now setting prepare_subfunction as True, as you suggest. I tried it in the past but it consumed a lot of memory and that's why I set it to False.

regarding the clipping range, I understand it is linked to the Hounsfield factor, right? I set it to the proper one (in my opinion) for the pancreas, which is quite different than the one for lungs, but still I do not think it is finding the pancreas at all during training... I will keep on trying things in this regard

thanks!

[Jugaris89]

Hi @muellerdo

I think we have gone a bit further with our problem but we still encounter the issue with the "ValueError: negative dimensions are not allowed" when loading some datasets (not all).

We have updated the miscnn version as you suggested, but it still happens.

Could you help us please and tell us where could we manually apply the patch inside miscnn to avoid this error?

Thanks in advance

[muellerdo]

Hey @Jugaris89,

mhm, could you please provide the full error log for this issue using the updated miscnn version? As well as your code how you applied the resampling?

Cheers, Dominik