Creating Volume from Numpy

Support
1

Hi all,

I have slicer version 4.7.0-2017-07-07 r26146 installed.

I am trying to convert a numpy 3d array into a new volume. I tried the following as listed in the nightly documentation:

resultVolumeNode = getNode(‘fixed’)
resultVtkArray = vtk.util.numpy_support.numpy_to_vtk(a, deep=True, array_type=vtk.VTK_SHORT)
resultVolumeNode.GetImageData().SetScalars(resultVtkArray)

I’m running into a couple of issues.

  1. numpy_to_vtk only seems to support 1 or 2d arrays. What is your recommended approach for 3d arrays?

  2. SetScalars isn’t an available function on the (vtkCommonDataModelPython.vtkImageData) instance.

If anyone has any suggestions on how to tackle these issues, I would greatly appreciate it!

Hailey

2

Hi Hailey -

If you are creating a new volume from scratch you need to create and initialize the vtkImageData as well as the vtkMRML nodes.

Here’s a fully worked out example in this postNRRD method which starts from a nrrd file in memory but it should give you the info you need to adapt to your use case.

Here are the key steps:

-Steve

3

I’ve added better numpy array adapter methods and added working examples to the wiki:

https://www.slicer.org/wiki/Documentation/Nightly/Developers/Python_scripting#Accessing_Volume_data_as_numpy_array

The new methods will be available in tomorrow’s nightly build. If you don’t want to wait till then then you can get the updated util.py from here:

2 Likes
4

A post was split to a new topic: Get each slice of a volume as numpy array

5 
import numpy as np
import math

def some_func(x, y, z):
  return 0.01*x*x + 0.01*y*y + 0.01*z*z

a = np.fromfunction(some_func,(512,512,512))

volumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLScalarVolumeNode')
volumeNode.CreateDefaultDisplayNodes()
updateVolumeFromArray(volumeNode, a)
setSliceViewerLayers(background=volumeNode)

I went created a empty volume with a 512 size, but I get a “Thread” map( the upper python code)

6

Do you have any question related to this?

7

What I am doing is improve the resolution for “import model as segment”
Firstly, I need create a volume by python script…
Secondly, I will porform crop volume, and then the model could be imported as a segment…Finally, the resolution for 3d reconstruction will be higher
So, would you like to gvie me some advice to such a script? Thanks

2018-09-09_172133
2018-09-09_172133758×486 108 KB

8

See the script repository for an example how to create a new volume: https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Create_a_new_volume

In recent Slicer nightly builds, there is a feature introduced exactly for this, as explained in this post: Improve the resolution of segment for "import a model"

9

This is not thread safe, I have tried to process several volumes independently at the same time but it does not work. The linear code works good but when trying to use some parallelization technique the calls to slicer.mrmlScene.AddNewNodeByClass(“vtkMRMLScalarVolumeNode”, nodeName) breaks the flow with silly errors such as: module slicer does not contain member mrmlScene. even the call to slicer.util.getNode(‘MRHead’) which is read only breaks the flow.
Sadly one has to do this volume by volume, I am dealing with 15 volumes here, so it would be nice some way to parallelize the code.
I use joblib that can be installed easily.

10

Yes, of course it is not thread-safe. Thread safety has enormous overhead in complexity of development and also has significant performance impact, so no software component is ever developed thread-safe unless there is a strong reason to do so. You always need to add your own thread-safe synchronization and communication mechanisms around general-purpose code to make it thread-safe.

11

I dont need synchronization/communication in my problem, I just want to modify one volume with some calculations, then the next, and so on. Therefore there are no interactions between them. Basically this is parallelization, not multi-threading. Slicer should be able to handle this, I mean several calls to its mrmlscene and util.
Minimum example

def calculate(name):
    node = slicer.util.getNode(name)
    volumeNode = slicer.vtkSlicerVolumesLogic().CloneVolume(slicer.mrmlScene,node,name.replace('ACSC', 'DOSE'),True)
    print('Node '+name.replace('ACSC', 'DOSE')+' created...')

nodes = slicer.util.getNodesByClass('vtkMRMLScalarVolumeNode')
nodenames = [node.GetName() for node in nodes]
nodenames = [name for name in nodenames if 'ACSC' in name]

for name in nodenames: # this works fine
    calculate(name)

from joblib import Parallel, delayed  # install joblib first using internal pip
Parrallel(n_jobs=4)(delayed(calculate)(name) for name in nodenames) # This launches the above error
12

Multi-threading is parallelization. Read up a bit more on this here.

Operations that may seem independent may often interfere with each other. For example, many Slicer modules observe node modification events, so any node code that triggers node modifications must be run from the main thread (cannot be run from worker threads).

You actually do need them and use them if you use joblib. Joblib can hide some of the low-level details and choose reasonable defaults, so you may get performance improvement and correct results without thinking for some very simple cases. However, in general, you cannot just blindly accept defaults expect that any functions that you put in calculate will work correctly.

I would suggest to use only a very restricted set of operations within calculate. For example, only use simple native python and numpy calls (and carefully study joblib and numpy documentation about what operations and parallelization options are permitted). Do all MRML method calls before/after the parallel processing.

2 Likes
13

Hello,

I have slicer version 4.10.2 r28257 installed on an iMac

I am trying to convert two volume nodes into a numpy 3d arrays, do some processing to them, and then convert the numpy arrays back into a volume. I tried the method as listed in the nightly documentation that Andras posted:

import numpy as np
import math
a = arrayFromVolume(input1Volume)
b = arrayFromVolume(input2Volume)

def some_func(x, y):
return y(x-y)

a = np.fromfunction(some_func,(30,20,15))

volumeNode = slicer.mrmlScene.AddNewNodeByClass(‘vtkMRMLScalarVolumeNode’)
volumeNode.CreateDefaultDisplayNodes()
updateVolumeFromArray(volumeNode, a)
setSliceViewerLayers(background=volumeNode)

But for that, I got the error that “arrayFromVolume” global name not defined.

Then, I put slicer.util before what I took to be the slicer related commands, like so:

import numpy as np
import slicer.util
import math
a = slicer.util.arrayFromVolume(input1Volume)
b = slicer.util.arrayFromVolume(input2Volume)

def some_func(x, y):
return y(x-y)

a = np.fromfunction(some_func,(30,20,15))

volumeNode = slicer.mrmlScene.AddNewNodeByClass(‘vtkMRMLScalarVolumeNode’)
volumeNode.CreateDefaultDisplayNodes()
slicer.util.updateVolumeFromArray(volumeNode, a)
setSliceViewerLayers(background=volumeNode)

After that, though, I got the error that ‘only integer scalar arrays can be converted to a scalar index’

The volumes are from DICOM files, but I (think?) I have them as mrml volume nodes in my script, since I used the Qt GUI to treat the DICOM files I select as inputs as volume nodes.
I tried the methods listed on this page, but I they either didn’t work or I implemented them incorrectly. Could you please advise me on how to fix this problem?

Thanks,
Timothy

14

You correctly retrieve the volume from Slicer and you get the error message in the part where you manipulate the arrays using numpy functions. I’ve Google for the error message and found this solution: python - numpy array TypeError: only integer scalar arrays can be converted to a scalar index - Stack Overflow

1 Like
15

Hi Andras,
I don’t think that solution applies here, since I’m not trying to simply make a three dimensional array from two arrays, and I tried the alternative solution of flattening the arrays, but I’m not sure exactly how that is supposed to help; and I got the same error message when I put it in my script.

import numpy as np
import slicer.util
import math
a = slicer.util.arrayFromVolume(input1Volume)
b = slicer.util.arrayFromVolume(input2Volume)

a.reshape((1, -1))
b.reshape((1, -1))

def subtract_mult(x, y):
  return y(x-y)

c = np.fromfunction(some_func,(a,b))

volumeNode = slicer.mrmlScene.AddNewNodeByClass(‘vtkMRMLScalarVolumeNode’)
volumeNode.CreateDefaultDisplayNodes()
slicer.util.updateVolumeFromArray(volumeNode, c)
setSliceViewerLayers(background=volumeNode)

Do I have to somehow transform the arrays obtained with slicer.util into numpy arrays? I thought they already were numpy arrays, at least it seemed that way from the documentation on the wiki.
Thank you for your help,
Timothy

16

slicer.util.array and slicer.util.arrayFromVolume already return numpy array.

You can Subtract an array from the other by calling d = a - b. If the geometry of the two volumes are different then you need to resample one to have the same geometry as the other using one of the volume resampling modules.

1 Like
17

Hi Andras,
The volumes have the same geometry. One is a Dixon fat image, and the other is a Dixon water image; both acquired simultaneously and they must have the same geometry. I assumed they did return numpy arrays, but I am still having the same error with ‘only integer scalar arrays can be converted to a scalar index’. I think this error means that the arrays I am getting from calling

slicer.util.arrayFromVolume(volumeNode)

must somehow not be integer scalar arrays. Would you happen to know if this is the case, or if

slicer.util.updateVolumeFromArray(volume2Node,a)

does not accept the numpy arrays generated from the volume nodes?

Am I oversimplifying this? Because from the above thread I can see that the process of creating a volume array seems much more complicated than it does in the documentation.

Thanks again,
Timothy Pinkhassik

18

This works well for me:

import SampleData
[input1Volume, input2Volume] = SampleData.SampleDataLogic().downloadDentalSurgery()

import slicer.util
a = slicer.util.arrayFromVolume(input1Volume)
b = slicer.util.arrayFromVolume(input2Volume)

c = b-a

volumeNode = slicer.modules.volumes.logic().CloneVolume(input1Volume, "Difference")
slicer.util.updateVolumeFromArray(volumeNode, c)
setSliceViewerLayers(background=volumeNode)
1 Like
19

Hi Andras,
Thank you so much! That worked perfectly.
-Timothy

20

Hi Andras,
I get new volumes but its empty. It has same dimensions as the one I have and its created but its all black.

I am using fuzzy classification to get two classes. whenI create volume its empty. Ithink I need to set scalars or is there anything I am doing wrong. Please help the code is below.

import numpy as np
import skfuzzy as fuzz
# Compute the thresholded output volume using the Threshold Scalar Volume CLI module
“”“cliParams = {
‘InputVolume’: inputVolume.GetID(),
‘OutputVolume’: outputVolume.GetID(),
‘ThresholdValue’ : imageThreshold,
‘ThresholdType’ : ‘Below’ if invert else ‘Above’
}
cliNode = slicer.cli.run(slicer.modules.thresholdscalarvolume, None, cliParams, wait_for_completion=True)”""
mriImage = slicer.util.getNode(inputVolume.GetID())
brainImg = slicer.util.arrayFromVolume(mriImage)

maskImage = slicer.util.getNode(outputVolume.GetID())
brainMask = slicer.util.arrayFromVolume(maskImage)

voxelIntensities = brainImg[brainMask > 0]
print(voxelIntensities)

ncenters = 2
cntr, u, u0, d, jm, p, fpc = fuzz.cluster.cmeans(voxelIntensities.reshape(1, voxelIntensities.size), ncenters, 2, error=0.005, maxiter=1000, init=None)
brain_ventricle_membership = np.zeros(brainImg.shape)
brain_ventricle_membership[brainMask > 0] = u[0]

brain_parenchima_membership = np.ones(brainImg.shape)
brain_parenchima_membership[brainMask > 0] = u[1]

volumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLScalarVolumeNode')
volumeNode.CreateDefaultDisplayNodes()
slicer.util.updateVolumeFromArray(volumeNode, brain_ventricle_membership)