Table of Contents
AIMS is a C++ library, but has python language bindings: PyAIMS. This means that the C++ classes and functions can be used from python. This has many advantages compared to pure C++:
- Writing python scripts and programs is much easier and faster than C++: there is no fastidious and long compilation step.
- Scripts are more flexible, can be modified on-the-fly, etc
- It can be used interactively in a python interactive shell.
- As pyaims is actually C++ code called from python, it is still fast to execute complex algorithms. There is obviously an overhead to call C++ from python, but once in the C++ layer, it is C++ execution speed.
The data for the examples in this section can be downloaded here: ftp://ftp.cea.fr/pub/dsv/anatomist/data/demo_data.zip. To use the examples directly, users should go to the directory where this archive was uncompressed, and then run ipython from this directory. A cleaner alternative, especially if no write access is allowed on this data directory, is to make a symbolic link to the "data_for_anatomist" subdirectory:
cd $HOME mkdir bvcourse cd bvcourse ln -s <path_to_data>/data_for_anatomist . ipython
More information about the API is available here.
In python, the aimsdata library is available as the soma.aims module.
import soma.aims # the module is actually soma.aims: vol = soma.aims.Volume( 100, 100, 100, dtype='int16' )or:
from soma import aims # the module is available as aims (not soma.aims): vol = aims.Volume( 100, 100, 100, dtype='int16' ) # in the following, we will be using this form because it is shorter.
Reading operations are accessed via a single read() function, and writing through a single write() function. The read() function reads any object from a given file name, in any supported file format, and returns it:
from soma import aims obj = aims.read( 'data_for_anatomist/subject01/subject01.nii' ) print obj obj2 = aims.read( 'data_for_anatomist/subject01/Audio-Video_T_map.nii' ) print obj2 obj3 = aims.read( 'data_for_anatomist/subject01/subject01_Lhemi.mesh' ) print obj3The returned object can have various types according to what is found in the disk file(s).
Writing is just as easy. The file name extension generally determines the output format. An object read from a given format can be re-written in any other supported format, provided the format can actually store the object type.
from soma import aims obj2 = aims.read( 'data_for_anatomist/subject01/Audio-Video_T_map.nii' ) aims.write( obj2, 'Audio-Video_T_map.ima' ) obj3 = aims.read( 'data_for_anatomist/subject01/subject01_Lhemi.mesh' ) aims.write( obj3, 'subject01_Lhemi.gii' )
Volumes are array-like containers of voxels, plus a set of additional information kept in a header structure. In AIMS, the header structure is generic and extensible, and does not depend on a specific file format. Voxels may have various types, so a specific type of volume should be used for a specific type of voxel. The type of voxel has a code that is used to suffix the Volume type: Volume_S16 for signed 16-bit ints, Volume_U32 for unsigned 32-bit ints, Volume_DOUBLE for 32-bit floats, Volume_RGBA for RGBA colors, etc.
# create a 3D volume of signed 16-bit ints, of size 192x256x128 vol = aims.Volume( 192, 256, 128, dtype='int16' ) # fill it with zeros vol.fill(0) # set value 12 at voxel ( 100, 100, 60 ) vol.setValue( 12, 100, 100, 60 ) # get value at the same position x = vol.value( 100, 100, 60 ) # set the voxels size vol.header()[ 'voxel_size' ] = [ 0.9, 0.9, 1.2, 1. ] print vol.header()
Whole volume operations:
# multiplication, addition etc vol *= 2 vol2 = vol * 3 + 12 vol /= 2 vol3 = vol2 - vol - 12 vol4 = vol2 * vol / 6Voxel-wise operations:
# fill the volume with the distance to voxel ( 100, 100, 60 )
vs = vol.header()[ 'voxel_size' ]
pos0 = ( 100 * vs[0], 100 * vs[1], 60 * vs[2] ) # in millimeters
for z in xrange( vol.getSizeZ() ):
for y in xrange( vol.getSizeY() ):
for x in xrange( vol.getSizeX() ):
# get current position in an aims.Point3df structure, in mm
p = aims.Point3df( x * vs[0], y * vs[1], z * vs[2] )
# get relative position to pos0, in voxels
p-= pos0
# distance: norm of vector p
dist = p.norm()
# set it into the volume
vol.setValue( dist, x, y, z )
# save the volume
aims.write( vol, 'distance.nii' )
Now look at the "distance.nii" volume in Anatomist.
data_for_anatomist/subject01/Audio-Video_T_map.nii and thresholds it so as to keep values above 3.
from soma import aims
vol = aims.read( 'data_for_anatomist/subject01/Audio-Video_T_map.nii' )
for z in xrange( vol.getSizeZ() ):
for y in xrange( vol.getSizeY() ):
for x in xrange( vol.getSizeX() ):
if vol.value( x, y, z ) < 3.:
vol.setValue( 0, x, y, z )
aims.write( vol, 'Audio-Video_T_thresholded.nii' )
data_for_anatomist/subject01/subject01.nii and keeps one voxel out of two in every direction.
from soma import aims
vol = aims.read( 'data_for_anatomist/subject01/subject01.nii' )
# allocate a new volume with half dimensions
vol2 = aims.Volume( vol.getSizeX() / 2, vol.getSizeY() / 2, vol.getSizeZ() / 2, dtype='DOUBLE' )
# set the voxel size to twice it was in vol
vs = vol.header()[ 'voxel_size' ]
vs2 = [ x * 2 for x in vs ]
vol2.header()[ 'voxel_size' ] = vs2
for z in xrange( vol2.getSizeZ() ):
for y in xrange( vol2.getSizeY() ):
for x in xrange( vol2.getSizeX() ):
vol2.setValue( vol.value( x*2, y*2, z*2 ), x, y, z )
aims.write( vol2, 'resampled.nii' )
The first thing that comes to mind when running these examples, is that they are slow. Indeed, python is an interpreted language and loops in any interpreted language are slow. In addition, accessing individually each voxel of the volume has the overhead of python/C++ bindings communications. The conclusion is that that kind of example is probably a bit too low-level, and should be done, when possible, by compiled libraries or specialized array-handling libraries. This is the role of numpy.
Accessing numpy arrays to AIMS volume voxels is supported:
import numpy vol.fill( 0 ) arr = numpy.array( vol, copy=False ) # set value 100 in a whole sub-volume arr[60:120, 60:120, 40:80] = 100 # note that arr is a shared view to the volume contents, # modifications will also affect the volume aims.write( vol, "cube.nii" )
Now we can re-write the thresholding example using numpy:
from soma import aims vol = aims.read( 'data_for_anatomist/subject01/Audio-Video_T_map.nii' ) arr = numpy.array( vol, copy=False ) arr[ numpy.where( arr < 3. ) ] = 0. aims.write( vol, 'Audio-Video_T_thresholded2.nii' )Here,
arr < 3. returns a boolean array with the same size as arr, and numpy.where() returns arrays of coordinates where the specified contition is true.
The distance example, using numpy, would like the following:
from soma import aims import numpy vol = aims.Volume( 192, 256, 128, 'S16' ) vol.header()[ 'voxel_size' ] = [ 0.9, 0.9, 1.2, 1. ] vs = vol.header()[ 'voxel_size' ] pos0 = ( 100 * vs[0], 100 * vs[1], 60 * vs[2] ) # in millimeters arr = numpy.array( vol, copy=False ) # build arrays of coordinates for x, y, z x, y, z = numpy.ogrid[ 0.:vol.getSizeX(), 0.:vol.getSizeY(), 0.:vol.getSizeZ() ] # get coords in millimeters x *= vs[0] y *= vs[1] z *= vs[2] # relative to pos0 x -= pos0[0] y -= pos0[1] z -= pos0[2] # get norm, using numpy arrays broadcasting arr[:,:,:,0] = numpy.sqrt( x**2+y**2+z**2 ) # and save result aims.write( vol, 'distance2.nii' )This example appears a bit more tricky, since we must build the coordinates arrays, but is way faster to execute, because all loops within the code are executed in compiled routines in numpy. One interesting thing to note is that this code is using the famous "array broadcasting" feature of numpy, where arrays of heterogeneous sizes can be combined, and the "missing" dimensions are extended.
To make a deep-copy of a volume, use the copy constructor:
vol2 = aims.Volume( vol ) vol2.setValue( 12, 100, 100, 60 ) # now vol and vol2 have different values print 'vol.value( 100, 100, 60 ):', vol.value( 100, 100, 60 ) print 'vol2.value( 100, 100, 60 ):', vol2.value( 100, 100, 60 )If you need to build another, different volume, with the same structure and size, don't forget to copy the header part:
vol2 = aims.Volume( vol.getSizeX(), vol.getSizeY(), vol.getSizeZ(), vol.getSizeT(), 'FLOAT' ) vol2.header().update( vol.header() )Important information can reside in the header, like voxel size, or coordinates systems and geometric transformations to other coordinates systems, so it is really very important to carry this information with duplicated or derived volumes.
You can also build a volume from a numpy array:
arr = numpy.array( numpy.diag( xrange( 40 ) ), dtype=numpy.float32 ).reshape( 40, 40, 1 ) \
+ numpy.array( xrange( 20 ), dtype=numpy.float32 ).reshape( 1, 1, 20 )
# WARNING: the array must be in Fortran ordering for AIMS, at leat at the moment
# whereas the numpy addition always returns a C-ordered array
arr = numpy.array( arr, order='F' )
arr[ 10, 12, 3 ] = 25
vol = aims.Volume( arr )
print 'vol.value( 10, 12, 3 ):', vol.value( 10, 12, 3 )
# data are shared with arr
vol.setValue( 35, 10, 15, 2 )
print 'arr[ 10, 15, 2 ]:', arr[ 10, 15, 2 ]
arr[ 12, 15, 1 ] = 44
print 'vol.value( 12, 15, 1 ):', vol.value( 12, 15, 1 )
4D volumes work just like 3D volumes. Actually all volumes are 4D in AIMS, but the last dimension is commonly of size 1. In value() and setValue() methods, only the first dimension is mandatory, others are optional and default to 0, but up to 4 coordinates may be used. In the same way, the constructor takes up to 4 dimension parameters:
from soma import aims # create a 4D volume of signed 16-bit ints, of size 30x30x30x4 vol = aims.Volume( 30, 30, 30, 4, 'S16' ) # fill it with zeros vol.fill(0) # set value 12 at voxel ( 10, 10, 20, 2 ) vol.setValue( 12, 10, 10, 20, 2 ) # get value at the same position x = vol.value( 10, 10, 20, 2 ) # set the voxels size vol.header()[ 'voxel_size' ] = [ 0.9, 0.9, 1.2, 1. ] print vol.header()Similarly, 1D or 2D volumes may be used exactly the same way.
For historical reasons, another set of classes may also represent volumes. These classes are the older API in AIMS, and tend to be obsolete. But as they were used in many many routines and programs, they have still not been eradicated. Many C++ routines build volumes and actually return those older classes, so we could not really hide them, and they also have python bindings. These classes are aims.AimsData_<type>. Converting from and to aims.Volume_ classes is rather simple since the newer Volume classes are used internally in the AimsData API.
from soma import aims # create a 4D volume of signed 16-bit ints, of size 30x30x30x4 vol = aims.Volume( 30, 30, 30, 4, 'S16' ) vol.header()[ 'voxel_size' ] = [ 0.9, 0.9, 1.2, 1. ] advol = aims.AimsData( vol ) # vol and advol share the same header and voxel data vol.setValue( 12, 10, 10, 20, 2 ) print 'advol.value( 10, 10, 20, 2 ):', advol.value( 10, 10, 20, 2 ) advol.setValue( 44, 12, 12, 24, 1 ) print 'vol.value( 12, 12, 24, 1 ):', vol.value( 12, 12, 24, 1 )And, in the other direction:
from soma import aims # create a 4D volume of signed 16-bit ints, of size 30x30x30x4 advol = aims.AimsData( 30, 30, 30, 4, 'S16' ) advol.header()[ 'voxel_size' ] = [ 0.9, 0.9, 1.2, 1. ] vol = advol.volume() # vol and advol share the same header and voxel data vol.setValue( 12, 10, 10, 20, 2 ) print 'advol.value( 10, 10, 20, 2 ):', advol.value( 10, 10, 20, 2 ) advol.setValue( 44, 12, 12, 24, 1 ) print 'vol.value( 12, 12, 24, 1 ):', vol.value( 12, 12, 24, 1 )
AimsData has a bit richer API, since it includes minor processing functions that have been removed from the newer Volume for the sake of API simplicity and minimalism.
# minimum / maximum print 'min:', advol.minimum(), 'at', advol.minIndex() print 'max:', advol.maximum(), 'at', advol.maxIndex() # clone copy advol2 = advol.clone() advol2.setValue( 12, 4, 8, 11, 3 ) # now advol and advol2 have different values print 'advol.value( 4, 8, 11, 3 ):', advol.value( 4, 8, 11, 3 ) print 'advol2.value( 4, 8, 11, 3 ):', advol2.value( 4, 8, 11, 3 ) # Border handling # Border width is th 5th parameter of AimsData constructor advol = aims.AimsData( 192, 256, 128, 1, 2, 'S16' ) advol.header()[ 'voxel_size' ] = [ 0.9, 0.9, 1.2, 1. ] advol.fill( 0 ) advol.setValue( 15, 100, 100, 60 ) vol = advol.volume() # then vol is 4 voxels wider in each direction, and shifted: print 'vol.value( 100, 100, 60 ):', vol.value( 100, 100, 60 ) # ... it is 0, not 15... print 'vol.value( 102, 102, 62 ):', vol.value( 102, 102, 62 ) # here we get 15 # some algorithms require this border to exist, otherwise fail or crash... from soma import aimsalgo aimsalgo.AimsDistanceFrontPropagation( advol, 0, -1, 3, 3, 3, 10, 10 ) aims.write( advol, 'distance3.nii' )
A surfacic mesh represents a surface, as a set of small polygons (generally triangles, but sometimes quads). It has two main components: a vector of vertices (each vertex is a 3D point, with coordinates in millimeters), and a vector of polygons: each polygon is defined by the vertices it links (3 for a triangle). It also optionally has normals (unit vectors). In our mesh structures, there is one normal for each vertex.
from soma import aims mesh = aims.read( 'data_for_anatomist/subject01/subject01_Lhemi.mesh' ) vert = mesh.vertex() print 'vertices:', len( vert ) poly = mesh.polygon() print 'polygons:', len( poly ) norm = mesh.normal() print 'normals:', len( norm )To build a mesh, we can instantiate an object of type
aims.AimsTimeSurface_<n>, with n being the number of vertices by polygon. Then we can add vertices, normals and polygons to the mesh:
# build a flying saucer mesh from soma import aims import numpy mesh = aims.AimsTimeSurface( 3 ) # a mesh has a header mesh.header()[ 'toto' ] = 'a message in the header' vert = mesh.vertex() poly = mesh.polygon() x = numpy.cos( numpy.ogrid[ 0.:20 ] * numpy.pi / 10. ) * 100 y = numpy.sin( numpy.ogrid[ 0.:20 ] * numpy.pi / 10. ) * 100 z = numpy.zeros( 20 ) c = numpy.vstack( ( x, y, z ) ).transpose() vert.assign( [ aims.Point3df( 0., 0., -40. ), aims.Point3df( 0., 0., 40. ) ] + [ aims.Point3df( x ) for x in c ] ) pol = numpy.vstack( ( numpy.zeros( 20, dtype=numpy.int32 ), numpy.ogrid[ 3:23 ], numpy.ogrid[ 2:22 ] ) ).transpose() pol[ 19, 1 ] = 2 pol2 = numpy.vstack( ( numpy.ogrid[ 2:22 ], numpy.ogrid[ 3:23 ], numpy.ones( 20, dtype=numpy.int32 ) ) ).transpose() pol2[19,1] = 2 poly.assign( [ aims.AimsVector(x,'U32',3) for x in numpy.vstack( ( pol, pol2 ) ) ] ) # write result aims.write( mesh, 'saucer.mesh' ) # automatically calculate normals mesh.updateNormals()
# slightly inflate a mesh from soma import aims import numpy mesh = aims.read( 'data_for_anatomist/subject01/subject01_Lwhite.mesh' ) vert = mesh.vertex() varr = numpy.array( vert ) norm = numpy.array( mesh.normal() ) varr += norm * 2 # push vertices 2mm away along normal vert.assign( [ aims.Point3df(x) for x in varr ] ) mesh.updateNormals() aims.write( mesh, 'subject01_Lwhite_semiinflated.mesh' )Now look at both meshes in Anatomist...
Alternatively, without numpy, we could have written the code like this:
from soma import aims mesh = aims.read( 'data_for_anatomist/subject01/subject01_Lwhite.mesh' ) vert = mesh.vertex() norm = mesh.normal() for v, n in zip( vert, norm ): v += n * 2 mesh.updateNormals() aims.write( mesh, 'subject01_Lwhite_semiinflated.mesh' )
In AIMS, meshes are actually time-indexed dictionaries of meshes. This way a deforming mesh can be stored in the same object. To copy a timestep to aonother, use the following:
from soma import aims mesh = aims.read( 'data_for_anatomist/subject01/subject01_Lwhite.mesh' ) # mesh.vertex() is equivalent to mesh.vectex( 0 ) mesh.vertex( 1 ).assign( mesh.vertex( 0 ) ) # same for normals and polygons mesh.normal( 1 ).assign( mesh.normal( 0 ) ) mesh.polygon( 1 ).assign( mesh.polygon( 0 ) ) print 'number of time steps:', mesh.size()
from soma import aims import numpy mesh = aims.read( 'data_for_anatomist/subject01/subject01_Lwhite.mesh' ) vert = mesh.vertex() varr = numpy.array( vert ) norm = numpy.array( mesh.normal() ) for i in xrange( 1, 10 ): mesh.normal( i ).assign( mesh.normal() ) mesh.polygon( i ).assign( mesh.polygon() ) varr += norm * 0.5 mesh.vertex( i ).assign( [ aims.Point3df(x) for x in varr ] ) mesh.updateNormals() aims.write( mesh, 'subject01_Lwhite_semiinflated_time.mesh' )
A texture is merely a vector of values, each of them is assigned to a mesh vertex, with a one-to-one mapping, in the same order. A texture is also a time-texture.
from soma import aims tex = aims.TimeTexture( 'FLOAT' ) t = tex[0] # time index, inserts on-the-fly t.reserve( 10 ) # pre-allocates memory for i in xrange( 10 ): t.append( i / 10. )
data_for_anatomist/subject01/subject01.niifrom soma import aims
mesh = aims.read( 'subject01_Lwhite_semiinflated_time.mesh' )
vol = aims.read( 'data_for_anatomist/subject01/subject01.nii' )
tex = aims.TimeTexture( 'FLOAT' )
vs = vol.header()[ 'voxel_size' ]
for i in xrange( mesh.size() ):
t = tex[i]
vert = mesh.vertex( i )
t.reserve( len( vert ) )
for p in vert:
t.append( vol.value( *[ int( round(x/y) ) for x,y in zip( p, vs ) ] ) )
aims.write( tex, 'subject01_Lwhite_semiinflated_texture.tex' )
Now look at the texture on the mesh (inflated or not) in Anatomist. Compare it to a 3D fusion between the mesh and the MRI volume.
Bonus: We can do the same for functional data. But in this case we may have a spatial transformation to apply between anatomical data and functional data (which may have been normalized, or acquired in a different referential).
from soma import aims
import numpy
mesh = aims.read( 'subject01_Lwhite_semiinflated_time.mesh' )
vol = aims.read( 'data_for_anatomist/subject01/Audio-Video_T_map.nii' )
# get header info from anatomical volume
f = aims.Finder()
f.check( 'data_for_anatomist/subject01/subject01.nii' )
anathdr = f.header()
# get functional -> MNI transformation
m1 = aims.AffineTransformation3d( vol.header()[ 'transformations' ][1] )
# get anat -> MNI transformation
m2 = aims.AffineTransformation3d( anathdr[ 'transformations' ][1] )
# make anat -> functional transformation
anat2func = m1.inverse() * m2
# include functional voxel size to get to voxel coordinates
vs = vol.header()[ 'voxel_size' ]
mvs = aims.AffineTransformation3d( numpy.diag( vs[:3] + [ 1. ] ) )
anat2func = mvs.inverse() * anat2func
# now go as in the previous program
tex = aims.TimeTexture( 'FLOAT' )
for i in xrange( mesh.size() ):
t = tex[i]
vert = mesh.vertex( i )
t.reserve( len( vert ) )
for p in vert:
t.append( vol.value( *[ int(round(x)) for x in anat2func.transform( p ) ] ) )
aims.write( tex, 'subject01_Lwhite_semiinflated_audio_video.tex' )
See how the functional data on the mesh changes across the depth of the cortex. this demonstrates the need to have a proper projection of functional data before dealing with surfacic functional processing.
"Buckets" are voxels lists. They are typically used to represent ROIs. A BucketMap is a list of Buckets. Each Bucket contains a list of voxels coordinates.
from soma import aims
bck_map=aims.read( 'data_for_anatomist/roi/basal_ganglia.data/roi_Bucket.bck' )
print 'Bucket map: ', bck_map
print 'Nb buckets: ', bck_map.size()
for i in xrange(bck_map.size()):
b=bck_map[i]
print "Bucket ", i, ", nb voxels: ", b.size()
if b.keys():
print " Coordinates of the first voxel: ", b.keys()[0].list()
Graphs are data structures that may contain various elements. They can represent sets of smaller structures, and also relations between such structures. The main usage we have for them is to represent ROIs sets, sulci, or fiber bundles.
A graph contains:
- properties of any type, like a volume or mesh header.
- nodes (also called vertices), which represent structured elements (a ROI, a sulcus part, etc), which in turn can store properties, and geometrical elements: buckets, meshes...
- optionally, relations, which link nodes and can also contain properties and geometrical elements.
Properties are stored in a dictionary-like way. They can hold almost anything, but a restricted set of types can be saved and loaded. It is exactly the same thing as headers found in volumes, meshes, textures or buckets.
from soma import aims graph = aims.read( 'data_for_anatomist/roi/basal_ganglia.arg' ) print graph print 'properties:', graph.keys() for p, v in graph.iteritems(): print p, ':', v graph[ 'gudule' ] = [ 12, 'a comment' ]Note: Only properties declared in a "syntax" file may be saved and re-loaded. Other properties are just not saved.
Vertices (or nodes) can be accessed via the vertices() method. Each vertex is also a dictionary-like properties set.
for v in graph.vertices(): print v['name']
To insert a new vertex, the addVertex() method should be used:
v = graph.addVertex( 'roi' ) print v v[ 'name' ] = 'new ROI'
An edge, or relation, links nodes together. Up to now we have always used binary, unoriented, edges. They can be added using the addEdge() method. Edges are also dictionary-like properties sets.
v2 = graph.vertices().list()[1] e = graph.addEdge( v, v2, 'roi_link' ) print graph.edges() # get vertices linked by this edge print e.vertices()
Setting meshes or buckets in vertices properties is OK internally, but for saving and loading, additional consistancy must be ensured and internal tables update is required. Then, use the aims.GraphManip.storeAims function:
mesh = aims.read( 'data_for_anatomist/subject01/subject01_Lwhite.mesh' ) # store mesh in the 'roi' property of vertex v of graph graph aims.GraphManip.storeAims( graph, v, 'roi', mesh )
There are other examples for pyaims referenced here.