# Image model use cases

In which we lay out the various things that users and developers may want to do to images. See also Resampling use cases

## Taking a mean over a 4D image

We could do this much more simply than below, this is just an example of reducing over a particular axis:

    # take mean of 4D image
    from glob import glob
    import numpy as np
    import nipy as ni

    fname = 'some4d.nii'

    img_list = ni.load_list(fname, axis=3)
    vol0 = img_list[0]
    arr = vol0.array[:]
    for vol in img_list[1:]:
       arr += vol.array
    mean_img = ni.Image(arr, vol0.coordmap)
    ni.save(mean_img, 'mean_some4d.nii')

## Taking mean over series of 3D images

Just to show how this works with a list of images:

    # take mean of some PCA volumes
    fnames = glob('some3d*.nii')
    vol0 = ni.load(fnames[0])
    arr = vol0.array[:]
    for fname in fnames[1:]:
        vol = ni.load(fname)
        arr += vol.array
    mean_img = ni.Image(arr, vol0.coordmap)
    ni.save(mean_img, 'mean_some3ds.nii')

## Simple motion correction

This is an example of how the 4D -> list of 3D interface works:

    # motion correction
    img_list = ni.load_list(fname, axis=3)
    reggie = ni.interfaces.fsl.Register(tol=0.1)
    vol0 = img_list[0]
    mocod = [] # unresliced
    rmocod = [] # resliced
    for vol in img_list[1:]:
        rcoord_map = reggie.run(moving=vol, fixed=vol0)
        cmap = ni.ref.compose(rcoord_map, vol.coordmap)
        mocovol = ni.Image(vol.array, cmap)
        # But...
        try:
           a_vol = ni.Image(vol.array, rcoord_map)
        except CoordmapError, msg
           assert msg == 'need coordmap with voxel input'
        mocod.append(mocovol)
        rmocovol = ni.reslice(mocovol, vol0)
        rmocod.append(rmocovol)
    rmocod_img = ni.list_to_image(rmocovol)
    ni.save(rmocod_img, 'rsome4d.nii')
    try:
        mocod_img = ni.list_to_image(mocovol)
    except ImageListError:
        print 'That is what I thought; the transforms were not the same'

## Slice timing

Here putting 3D image into an image list, and back into a 4D image / array:

    # slice timing
    img_list = ni.load_list(fname, axis=2)
    slicetimer = ni.interfaces.fsl.SliceTime(algorithm='linear')
    vol0 = img_list[0]
    try:
       vol0.timestamp
    except AttributeError:
       print 'we do not have a timestamp'
    try:
       vol0.slicetimes
    except AttributeError:
       print 'we do not have slicetimes'
    try:
       st_list = slicetimer.run(img)
    except SliceTimeError, msg:
       assert msg == 'no timestamp for volume'
    TR = 2.0
    slicetime = 0.15
    sliceaxis = 2
    nslices = vol0.array.shape[sliceaxis]
    slicetimes = np.range(nslices) * slicetime
    timestamps = range(len(img_list)) * TR
    # Either the images are in a simple list
    for i, img in enumerate(img_list):
       img.timestamp = timestamps[i]
       img.slicetimes = slicetimes
       img.axis['slice'] = sliceaxis # note setting of voxel axis meaning
    # if the sliceaxes do not match, error when run
    img_list[0].axis['slice'] = 1
    try:
       st_list = slicetimer.run(img)
    except SliceTimeError, msg:
       assert msg == 'images do not have the same sliceaxes']
    # Or - with ImageList object
    img_list.timestamps = timestamps
    img_list.slicetimes = slicetimes
    img_list.axis['slice'] = sliceaxis
    # Either way, we run and save
    st_list = slicetimer.run(img)
    ni.save(ni.list_to_image(st_img), 'stsome4d.nii')

## Creating an image given data and affine

Showing how we would like the image creation API to look:

    # making an image from an affine
    data = img.array
    affine = np.eye(4)
    scanner_img = ni.Image(data, ni.ref.voxel2scanner(affine))
    mni_img = ni.Image(data, ni.ref.voxel2mni(affine))

## Coregistration / normalization

Demonstrating coordinate maps and non-linear resampling:

    # coregistration and normalization
    anat_img = ni.load_image('anatomical.nii')
    func_img = ni.load_image('epi4d.nii')
    template = ni.load_image('mni152T1.nii')

    # coreg
    coreger = ni.interfaces.fsl.flirt(tol=0.2)
    coreg_cmap = coreger.run(fixed=func_img, moving=anat_img)
    c_anat_img = ni.Image(anat_img.data, coreg_cmap.compose_with(anat_img.cmap))

    # calculate normalization parameters
    template_cmap = template.coordmap
    template_dims = template.data.shape
    c_anat_cmap = c_anat_img.coordmap
    normalizer = ni.interfaces.fsl.fnirt(param=3)
    norm_cmap = normalizer.run(moving=template, fixed=c_anat_img)

    # resample anatomical using calculated coordinate map
    full_cmap = norm_cmap.composed_with(template_cmap)
    w_anat_data = img.resliced_to_grid(full_cmap, template_dims)
    w_anat_img = ni.Image(w_anat_data, template.coordmap)

    # resample functionals with calculated coordinate map
    w_func_list = []
    for img in ni.image_list(func_img, axis=3):
      w_img_data = img.resliced_to_grid(full_cmap, template_dims)
      w_func_list.append(ni.Image(w_img_data, template_cmap))
    ni.save(ni.list_to_image(w_func_list), 'stsome4d.nii')