Source code for nighres.registration.simple_align

# basic dependencies
import os
import sys

# main dependencies: numpy, nibabel
import numpy as np
import nibabel as nb

# nighresjava and nighres functions
import nighresjava
from ..io import load_volume, save_volume
from ..utils import _output_dir_4saving, _fname_4saving, \
                    _check_topology_lut_dir

# convenience labels
X=0
Y=1
Z=2
T=3

[docs]def simple_align(source_image, target_image, copy_header=False, align_center=False, rescale=False, data_type='intensity', ignore_affine=False, ignore_header=False, save_data=False, overwrite=False, output_dir=None, file_name=None): """ Simple alignment routines Simple routines to align image headers (image data is unchanged) Parameters ---------- source_image: niimg Image to align target_image: niimg Reference image to match copy_header: bool To copy the target header to the source (default is False) align_center: bool To align the source center of mass to the target (default is False) rescale: bool To rescale the source to the volume of the target (default is False) data_type: {'intensity','nonzero','boundingbox'} The type of datato consider for alignment (default is 'intensity') ignore_affine: bool Ignore the affine matrix information extracted from the image header (default is False) ignore_header: bool Ignore the orientation information and affine matrix information extracted from the image header (default is False) save_data: bool Save output data to file (default is False) overwrite: bool Overwrite existing results (default is False) output_dir: str, optional Path to desired output directory, will be created if it doesn't exist file_name: str, optional Desired base name for output files with file extension (suffixes will be added) Returns ---------- dict Dictionary collecting outputs under the following keys (suffix of output files in brackets) * result (niimg): Aligne source image (_al-img) Notes ---------- This uses the ANTs/ITK conventions with regard to Nifti headers, for better or for worse. Note that Nibabel conventions, of course, are different. """ print('\nSimple align') # make sure that saving related parameters are correct output_dir = _output_dir_4saving(output_dir, source_image) # needed for intermediate results if save_data: result_file = os.path.join(output_dir, _fname_4saving(file_name=file_name, rootfile=source_image, suffix='al-img')) if overwrite is False \ and os.path.isfile(result_file) : print("skip computation (use existing results)") output = {'result': load_volume(result_file)} return output # load and get dimensions and resolution from input images source = load_volume(source_image) nsx = source.header.get_data_shape()[X] nsy = source.header.get_data_shape()[Y] nsz = source.header.get_data_shape()[Z] rsx = source.header.get_zooms()[X] rsy = source.header.get_zooms()[Y] rsz = source.header.get_zooms()[Z] target = load_volume(target_image) ntx = target.header.get_data_shape()[X] nty = target.header.get_data_shape()[Y] ntz = target.header.get_data_shape()[Z] rtx = target.header.get_zooms()[X] rty = target.header.get_zooms()[Y] rtz = target.header.get_zooms()[Z] # in case the affine transformations are not to be trusted: make them equal if ignore_affine or ignore_header: # create generic affine aligned with the orientation for the source mx = np.argmax(np.abs(source.affine[0][0:3])) my = np.argmax(np.abs(source.affine[1][0:3])) mz = np.argmax(np.abs(source.affine[2][0:3])) new_affine = np.zeros((4,4)) if ignore_header: new_affine[0][0] = rsx new_affine[1][1] = rsy new_affine[2][2] = rsz new_affine[0][3] = -rsx*nsx/2.0 new_affine[1][3] = -rsy*nsy/2.0 new_affine[2][3] = -rsz*nsz/2.0 else: new_affine[0][mx] = rsx*np.sign(source.affine[0][mx]) new_affine[1][my] = rsy*np.sign(source.affine[1][my]) new_affine[2][mz] = rsz*np.sign(source.affine[2][mz]) if (np.sign(source.affine[0][mx])<0): new_affine[0][3] = rsx*nsx/2.0 else: new_affine[0][3] = -rsx*nsx/2.0 if (np.sign(source.affine[1][my])<0): new_affine[1][3] = rsy*nsy/2.0 else: new_affine[1][3] = -rsy*nsy/2.0 if (np.sign(source.affine[2][mz])<0): new_affine[2][3] = rsz*nsz/2.0 else: new_affine[2][3] = -rsz*nsz/2.0 #if (np.sign(source.affine[0][mx])<0): new_affine[mx][3] = rsx*nsx #if (np.sign(source.affine[1][my])<0): new_affine[my][3] = rsy*nsy #if (np.sign(source.affine[2][mz])<0): new_affine[mz][3] = rsz*nsz #new_affine[0][3] = nsx/2.0 #new_affine[1][3] = nsy/2.0 #new_affine[2][3] = nsz/2.0 new_affine[3][3] = 1.0 source = nb.Nifti1Image(source.get_data(), new_affine, source.header) source.update_header() # create generic affine aligned with the orientation for the target mx = np.argmax(np.abs(target.affine[0][0:3])) my = np.argmax(np.abs(target.affine[1][0:3])) mz = np.argmax(np.abs(target.affine[2][0:3])) new_affine = np.zeros((4,4)) if ignore_header: new_affine[0][0] = rtx new_affine[1][1] = rty new_affine[2][2] = rtz new_affine[0][3] = -rtx*ntx/2.0 new_affine[1][3] = -rty*nty/2.0 new_affine[2][3] = -rtz*ntz/2.0 else: new_affine[0][mx] = rtx*np.sign(target.affine[0][mx]) new_affine[1][my] = rty*np.sign(target.affine[1][my]) new_affine[2][mz] = rtz*np.sign(target.affine[2][mz]) if (np.sign(target.affine[0][mx])<0): new_affine[0][3] = rtx*ntx/2.0 else: new_affine[0][3] = -rtx*ntx/2.0 if (np.sign(target.affine[1][my])<0): new_affine[1][3] = rty*nty/2.0 else: new_affine[1][3] = -rty*nty/2.0 if (np.sign(target.affine[2][mz])<0): new_affine[2][3] = rtz*ntz/2.0 else: new_affine[2][3] = -rtz*ntz/2.0 #if (np.sign(target.affine[0][mx])<0): new_affine[mx][3] = rtx*ntx #if (np.sign(target.affine[1][my])<0): new_affine[my][3] = rty*nty #if (np.sign(target.affine[2][mz])<0): new_affine[mz][3] = rtz*ntz #new_affine[0][3] = ntx/2.0 #new_affine[1][3] = nty/2.0 #new_affine[2][3] = ntz/2.0 new_affine[3][3] = 1.0 target = nb.Nifti1Image(target.get_data(), new_affine, target.header) target.update_header() # compute the various options if copy_header: result = nb.Nifti1Image(source.get_data(), target.affine, target.header) result.update_header() else: if align_center: # compute source center etc etc src_center = np.zeros(4) trg_center = np.zeros(4) if data_type is 'intensity': for x in range(nsx): for y in range(nsy): for z in range(nsz): src_center[X] += x*source.get_data()[x,y,z] src_center[Y] += y*source.get_data()[x,y,z] src_center[Z] += z*source.get_data()[x,y,z] src_center[T] += source.get_data()[x,y,z] for x in range(ntx): for y in range(nty): for z in range(ntz): trg_center[X] += x*target.get_data()[x,y,z] trg_center[Y] += y*target.get_data()[x,y,z] trg_center[Z] += z*target.get_data()[x,y,z] trg_center[T] += target.get_data()[x,y,z] elif data_type is 'nonzero': for x in range(nsx): for y in range(nsy): for z in range(nsz): if source.get_data()[x,y,z]>0: src_center[X] += x src_center[Y] += y src_center[Z] += z src_center[T] += 1 for x in range(ntx): for y in range(nty): for z in range(ntz): if target.get_data()[x,y,z]>0: trg_center[X] += x trg_center[Y] += y trg_center[Z] += z trg_center[T] += 1 elif data_typ is 'boundingbox': src_center[X] = nsx/2 src_center[Y] = nsy/2 src_center[Z] = nsz/2 src_center[T] = 1 trg_center[X] = nsx/2 trg_center[Y] = nsy/2 trg_center[Z] = nsz/2 trg_center[T] = 1 src_center[X] /= src_center[T] src_center[Y] /= src_center[T] src_center[Z] /= src_center[T] trg_center[X] /= trg_center[T] trg_center[Y] /= trg_center[T] trg_center[Z] /= trg_center[T] source.affine[X,T] = target.affine[X,T] - rsx*src_center[X] \ + ntx*trg_center[X] source.affine[Y,T] = target.affine[Y,T] - rsy*src_center[Y] \ + nty*trg_center[Y] source.affine[Z,T] = target.affine[Z,T] - rsz*src_center[Z] \ + ntz*trg_center[Z] if rescale: src_size=0 trg_size=0 if data_type is 'intensity': for x in range(nsx): for y in range(nsy): for z in range(nsz): src_size += source.get_data()[x,y,z] for x in range(ntx): for y in range(nty): for z in range(ntz): trg_size += target.get_data()[x,y,z] elif data_type is 'nonzero': for x in range(nsx): for y in range(nsy): for z in range(nsz): if source.get_data()[x,y,z]>0: src_size += 1 for x in range(ntx): for y in range(nty): for z in range(ntz): if target.get_data()[x,y,z]>0: trg_size += 1 elif data_typ is 'boundingbox': src_size = nsx*nsy*nsz trg_size = ntx*nty*ntz source.affine = trg_size/src_size*source.affine source.affine[T,T] = 1 result = nb.Nifti1Image(source.get_data(), source.affine, source.header) result.update_header() outputs = {'result': result} if save_data: save_volume(result_file, result) return outputs