# basic dependencies
import os
import sys
import subprocess
from glob import glob
import math
# 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 embedded_antsreg(source_image, target_image,
run_rigid=False,
rigid_iterations=1000,
run_affine=False,
affine_iterations=1000,
run_syn=True,
coarse_iterations=40,
medium_iterations=50, fine_iterations=40,
cost_function='MutualInformation',
interpolation='NearestNeighbor',
regularization='High',
convergence=1e-6,
mask_zero=False,
ignore_affine=False, ignore_header=False,
save_data=False, overwrite=False, output_dir=None,
file_name=None):
""" Embedded ANTS Registration
Runs the rigid and/or Symmetric Normalization (SyN) algorithm of ANTs and
formats the output deformations into voxel coordinate mappings as used in
CBSTools registration and transformation routines.
Parameters
----------
source_image: niimg
Image to register
target_image: niimg
Reference image to match
run_rigid: bool
Whether or not to run a rigid registration first (default is False)
rigid_iterations: float
Number of iterations in the rigid step (default is 1000)
run_affine: bool
Whether or not to run a affine registration first (default is False)
affine_iterations: float
Number of iterations in the affine step (default is 1000)
run_syn: bool
Whether or not to run a SyN registration (default is True)
coarse_iterations: float
Number of iterations at the coarse level (default is 40)
medium_iterations: float
Number of iterations at the medium level (default is 50)
fine_iterations: float
Number of iterations at the fine level (default is 40)
cost_function: {'CrossCorrelation', 'MutualInformation'}
Cost function for the registration (default is 'MutualInformation')
interpolation: {'NearestNeighbor', 'Linear'}
Interpolation for the registration result (default is 'NearestNeighbor')
regularization: {'Low', 'Medium', 'High'}
Regularization preset for the SyN deformation (default is 'Medium')
convergence: float
Threshold for convergence, can make the algorithm very slow
(default is convergence)
mask_zero: bool
Mask regions with zero value
(default is False)
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)
* transformed_source (niimg): Deformed source image (_ants-def)
* mapping (niimg): Coordinate mapping from source to target (_ants-map)
* inverse (niimg): Inverse coordinate mapping from target to source
(_ants-invmap)
Notes
----------
Port of the CBSTools Java module by Pierre-Louis Bazin. The main algorithm
is part of the ANTs software by Brian Avants and colleagues [1]_. The
interfacing with ANTs is performed through Nipype [2]_. Parameters have
been set to values commonly found in neuroimaging scripts online, but not
necessarily optimal.
References
----------
.. [1] Avants et al (2008), Symmetric diffeomorphic
image registration with cross-correlation: evaluating automated labeling
of elderly and neurodegenerative brain, Med Image Anal. 12(1):26-41
.. [2] Gorgolewski et al (2011) Nipype: a flexible, lightweight and
extensible neuroimaging data processing framework in python. Front
Neuroinform 5. doi:10.3389/fninf.2011.00013
"""
# just overloading the multi-channel version
return embedded_antsreg_multi([source_image], [target_image],
run_rigid, rigid_iterations, run_affine, affine_iterations,
run_syn, coarse_iterations, medium_iterations, fine_iterations,
cost_function, interpolation, regularization, convergence,
mask_zero, ignore_affine, ignore_header,
save_data, overwrite, output_dir, file_name)
[docs]def embedded_antsreg_2d(source_image, target_image,
run_rigid=False,
rigid_iterations=1000,
run_affine=False,
affine_iterations=1000,
run_syn=True,
coarse_iterations=40,
medium_iterations=50, fine_iterations=40,
scaling_factor=32,
cost_function='MutualInformation',
interpolation='NearestNeighbor',
regularization='High',
convergence=1e-6,
mask_zero=False,
ignore_affine=False, ignore_header=False,
save_data=False, overwrite=False, output_dir=None,
file_name=None):
""" Embedded ANTS Registration 2D
Runs the rigid and/or Symmetric Normalization (SyN) algorithm of ANTs and
formats the output deformations into voxel coordinate mappings as used in
CBSTools registration and transformation routines.
Parameters
----------
source_image: niimg
Image to register
target_image: niimg
Reference image to match
run_rigid: bool
Whether or not to run a rigid registration first (default is False)
rigid_iterations: float
Number of iterations in the rigid step (default is 1000)
run_affine: bool
Whether or not to run a affine registration first (default is False)
affine_iterations: float
Number of iterations in the affine step (default is 1000)
run_syn: bool
Whether or not to run a SyN registration (default is True)
coarse_iterations: float
Number of iterations at the coarse level (default is 40)
medium_iterations: float
Number of iterations at the medium level (default is 50)
fine_iterations: float
Number of iterations at the fine level (default is 40)
cost_function: {'CrossCorrelation', 'MutualInformation'}
Cost function for the registration (default is 'MutualInformation')
interpolation: {'NearestNeighbor', 'Linear'}
Interpolation for the registration result (default is 'NearestNeighbor')
convergence: flaot
Threshold for convergence, can make the algorithm very slow
(default is convergence)
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)
* transformed_source (niimg): Deformed source image (_ants-def)
* mapping (niimg): Coordinate mapping from source to target (_ants-map)
* inverse (niimg): Inverse coordinate mapping from target to source
(_ants-invmap)
Notes
----------
Port of the CBSTools Java module by Pierre-Louis Bazin. The main algorithm
is part of the ANTs software by Brian Avants and colleagues [1]_. Parameters
have been set to values commonly found in neuroimaging scripts online, but
not necessarily optimal.
References
----------
.. [1] Avants et al (2008), Symmetric diffeomorphic
image registration with cross-correlation: evaluating automated labeling
of elderly and neurodegenerative brain, Med Image Anal. 12(1):26-41
"""
print('\nEmbedded ANTs Registration 2D')
# check if ants is installed to raise sensible error
try:
subprocess.run('antsRegistration', stdout=subprocess.DEVNULL)
except FileNotFoundError:
sys.exit("\nCould not find command 'antsRegistration'. Make sure ANTs is"
" installed and can be accessed from the command line.")
try:
subprocess.run('antsApplyTransforms', stdout=subprocess.DEVNULL)
except FileNotFoundError:
sys.exit("\nCould not find command 'antsApplyTransforms'. Make sure ANTs"
" is installed and can be accessed from the command line.")
# make sure that saving related parameters are correct
# filenames needed for intermediate results
output_dir = _output_dir_4saving(output_dir, source_image)
transformed_source_file = os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='ants-def'))
mapping_file = os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='ants-map'))
inverse_mapping_file = os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='ants-invmap'))
if save_data:
if overwrite is False \
and os.path.isfile(transformed_source_file) \
and os.path.isfile(mapping_file) \
and os.path.isfile(inverse_mapping_file) :
print("skip computation (use existing results)")
output = {'transformed_source': transformed_source_file,
'mapping': mapping_file,
'inverse': inverse_mapping_file}
return output
# load and get dimensions and resolution from input images
source = load_volume(source_image)
src_affine = source.affine
src_header = source.header
nsx = source.header.get_data_shape()[X]
nsy = source.header.get_data_shape()[Y]
nsz = 1
rsx = source.header.get_zooms()[X]
rsy = source.header.get_zooms()[Y]
rsz = 1
target = load_volume(target_image)
trg_affine = target.affine
trg_header = target.header
ntx = target.header.get_data_shape()[X]
nty = target.header.get_data_shape()[Y]
ntz = 1
rtx = target.header.get_zooms()[X]
rty = target.header.get_zooms()[Y]
rtz = 1
# in case the affine transformations are not to be trusted: make them equal
if ignore_affine or ignore_header:
mx = np.argmax(np.abs(src_affine[0][0:3]))
my = np.argmax(np.abs(src_affine[1][0:3]))
mz = np.argmax(np.abs(src_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(src_affine[0][mx])
new_affine[1][my] = rsy*np.sign(src_affine[1][my])
new_affine[2][mz] = rsz*np.sign(src_affine[2][mz])
if (np.sign(src_affine[0][mx])<0):
new_affine[0][3] = rsx*nsx/2.0
else:
new_affine[0][3] = -rsx*nsx/2.0
if (np.sign(src_affine[1][my])<0):
new_affine[1][3] = rsy*nsy/2.0
else:
new_affine[1][3] = -rsy*nsy/2.0
if (np.sign(src_affine[2][mz])<0):
new_affine[2][3] = rsz*nsz/2.0
else:
new_affine[2][3] = -rsz*nsz/2.0
#if (np.sign(src_affine[0][mx])<0): new_affine[mx][3] = rsx*nsx
#if (np.sign(src_affine[1][my])<0): new_affine[my][3] = rsy*nsy
#if (np.sign(src_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
src_img = nb.Nifti1Image(source.get_data(), new_affine, source.header)
src_img.update_header()
src_img_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srcimg'))
save_volume(src_img_file, src_img)
source = load_volume(src_img_file)
src_affine = source.affine
src_header = source.header
# create generic affine aligned with the orientation for the target
mx = np.argmax(np.abs(trg_affine[0][0:3]))
my = np.argmax(np.abs(trg_affine[1][0:3]))
mz = np.argmax(np.abs(trg_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(trg_affine[0][mx])
new_affine[1][my] = rty*np.sign(trg_affine[1][my])
new_affine[2][mz] = rtz*np.sign(trg_affine[2][mz])
if (np.sign(trg_affine[0][mx])<0):
new_affine[0][3] = rtx*ntx/2.0
else:
new_affine[0][3] = -rtx*ntx/2.0
if (np.sign(trg_affine[1][my])<0):
new_affine[1][3] = rty*nty/2.0
else:
new_affine[1][3] = -rty*nty/2.0
if (np.sign(trg_affine[2][mz])<0):
new_affine[2][3] = rtz*ntz/2.0
else:
new_affine[2][3] = -rtz*ntz/2.0
#if (np.sign(trg_affine[0][mx])<0): new_affine[mx][3] = rtx*ntx
#if (np.sign(trg_affine[1][my])<0): new_affine[my][3] = rty*nty
#if (np.sign(trg_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
trg_img = nb.Nifti1Image(target.get_data(), new_affine, target.header)
trg_img.update_header()
trg_img_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_trgimg'))
save_volume(trg_img_file, trg_img)
target = load_volume(trg_img_file)
trg_affine = target.affine
trg_header = target.header
# build coordinate mapping matrices and save them to disk
src_coordX = np.zeros((nsx,nsy))
src_coordY = np.zeros((nsx,nsy))
trg_coordX = np.zeros((ntx,nty))
trg_coordY = np.zeros((ntx,nty))
for x in range(nsx):
for y in range(nsy):
src_coordX[x,y] = x
src_coordY[x,y] = y
src_mapX = nb.Nifti1Image(src_coordX, source.affine, source.header)
src_mapX_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordX'))
save_volume(src_mapX_file, src_mapX)
src_mapY = nb.Nifti1Image(src_coordY, source.affine, source.header)
src_mapY_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordY'))
save_volume(src_mapY_file, src_mapY)
for x in range(ntx):
for y in range(nty):
trg_coordX[x,y] = x
trg_coordY[x,y] = y
trg_mapX = nb.Nifti1Image(trg_coordX, target.affine, target.header)
trg_mapX_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_trgcoordX'))
save_volume(trg_mapX_file, trg_mapX)
trg_mapY = nb.Nifti1Image(trg_coordY, target.affine, target.header)
trg_mapY_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_trgcoordY'))
save_volume(trg_mapY_file, trg_mapY)
if mask_zero:
# create and save temporary masks
trg_mask_data = (target.get_data()!=0)
trg_mask = nb.Nifti1Image(trg_mask_data, target.affine, target.header)
trg_mask_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_trgmask'))
save_volume(trg_mask_file, trg_mask)
src_mask_data = (source.get_data()!=0)
src_mask = nb.Nifti1Image(src_mask_data, source.affine, source.header)
src_mask_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srcmask'))
save_volume(src_mask_file, src_mask)
# run the main ANTS software: here we directly build the command line call
reg = 'antsRegistration --collapse-output-transforms 1 --dimensionality 2' \
+' --initialize-transforms-per-stage 0 --interpolation Linear'
# add a prefix to avoid multiple names?
prefix = _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_syn')
prefix = os.path.basename(prefix)
prefix = prefix.split(".")[0]
reg = reg+' --output '+prefix
if mask_zero:
reg = reg+' --masks ['+trg_mask_file+', '+src_mask_file+']'
srcfile = source.get_filename()
trgfile = target.get_filename()
# figure out the number of scales, going with a factor of two
n_scales = math.ceil(math.log(scaling_factor)/math.log(2.0))
iter_rigid = str(rigid_iterations)
iter_affine = str(affine_iterations)
iter_syn = str(coarse_iterations)
smooth = str(scaling_factor)
shrink = str(scaling_factor)
for n in range(n_scales):
iter_rigid = iter_rigid+'x'+str(rigid_iterations)
iter_affine = iter_affine+'x'+str(affine_iterations)
if n<n_scales/2: iter_syn = iter_syn+'x'+str(coarse_iterations)
elif n<n_scales-1: iter_syn = iter_syn+'x'+str(medium_iterations)
else: iter_syn = iter_syn+'x'+str(fine_iterations)
smooth = smooth+'x'+str(scaling_factor/math.pow(2.0,n+1))
shrink = shrink+'x'+str(math.ceil(scaling_factor/math.pow(2.0,n+1)))
# set parameters for all the different types of transformations
if run_rigid is True:
reg = reg + ' --transform Rigid[0.1]'
if (cost_function=='CrossCorrelation'):
reg = reg + ' --metric CC['+trgfile+', '+srcfile \
+', '+'1.0, 5, Random, 0.3 ]'
else:
reg = reg + ' --metric MI['+trgfile+', '+srcfile \
+', '+'1.0, 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+iter_rigid+', '+str(convergence)+', 10 ]'
reg = reg + ' --smoothing-sigmas '+smooth
reg = reg + ' --shrink-factors '+shrink
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_affine is True:
reg = reg + ' --transform Affine[0.1]'
if (cost_function=='CrossCorrelation'):
reg = reg + ' --metric CC['+trgfile+', '+srcfile \
+', '+'1.0, 5, Random, 0.3 ]'
else:
reg = reg + ' --metric MI['+trgfile+', '+srcfile \
+', '+'1.0, 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+iter_affine+', '+str(convergence)+', 10 ]'
reg = reg + ' --smoothing-sigmas '+smooth
reg = reg + ' --shrink-factors '+shrink
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_syn is True:
if regularization == 'Low': syn_param = [0.2, 1.0, 0.0]
elif regularization == 'Medium': syn_param = [0.2, 3.0, 0.0]
elif regularization == 'High': syn_param = [0.2, 4.0, 3.0]
else: syn_param = [0.2, 3.0, 0.0]
reg = reg + ' --transform SyN'+str(syn_param)
if (cost_function=='CrossCorrelation'):
reg = reg + ' --metric CC['+trgfile+', '+srcfile \
+', '+'1.0, 5, Random, 0.3 ]'
else:
reg = reg + ' --metric MI['+trgfile+', '+srcfile \
+', '+'1.0, 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+iter_syn+', '+str(convergence)+', 5 ]'
reg = reg + ' --smoothing-sigmas '+smooth
reg = reg + ' --shrink-factors '+shrink
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_rigid is False and run_affine is False and run_syn is False:
reg = reg + ' --transform Rigid[0.1]'
reg = reg + ' --metric CC['+trgfile+', '+srcfile \
+', '+'1.0, 5, Random, 0.3 ]'
reg = reg + ' --convergence [ 0, 1.0, 2 ]'
reg = reg + ' --smoothing-sigmas 1.0'
reg = reg + ' --shrink-factors 1'
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
reg = reg + ' --write-composite-transform 0'
# run the ANTs command directly
print(reg)
try:
subprocess.check_output(reg, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+str(e.returncode)+')\n Output: '+str(e.output)
raise subprocess.CalledProcessError(msg)
# output file names
results = sorted(glob(prefix+'*'))
forward = []
flag = []
for res in results:
if res.endswith('GenericAffine.mat'):
forward.append(res)
flag.append(False)
elif res.endswith('Warp.nii.gz') and not res.endswith('InverseWarp.nii.gz'):
forward.append(res)
flag.append(False)
#print('forward transforms: '+str(forward))
inverse = []
linear = []
for res in results[::-1]:
if res.endswith('GenericAffine.mat'):
inverse.append(res)
linear.append(True)
elif res.endswith('InverseWarp.nii.gz'):
inverse.append(res)
linear.append(False)
#print('inverse transforms: '+str(inverse))
# Transforms the moving image
at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
at = at+' --input '+source.get_filename()
at = at+' --reference-image '+target.get_filename()
at = at+' --interpolation '+interpolation
for idx,transform in enumerate(forward):
if flag[idx]:
at = at+' --transform ['+transform+', 1]'
else:
at = at+' --transform ['+transform+', 0]'
at = at+' --output '+transformed_source_file
print(at)
try:
subprocess.check_output(at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# Create coordinate mappings
src_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
src_at = src_at+' --input '+src_mapX.get_filename()
src_at = src_at+' --reference-image '+target.get_filename()
src_at = src_at+' --interpolation Linear'
for idx,transform in enumerate(forward):
if flag[idx]:
src_at = src_at+' --transform ['+transform+', 1]'
else:
src_at = src_at+' --transform ['+transform+', 0]'
# src_at = src_at+' --output '+mapping_file
src_mapX_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordX_map'))
src_at = src_at+' --output '+src_mapX_trans
print(src_at)
try:
subprocess.check_output(src_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
src_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
src_at = src_at+' --input '+src_mapY.get_filename()
src_at = src_at+' --reference-image '+target.get_filename()
src_at = src_at+' --interpolation Linear'
for idx,transform in enumerate(forward):
if flag[idx]:
src_at = src_at+' --transform ['+transform+', 1]'
else:
src_at = src_at+' --transform ['+transform+', 0]'
# src_at = src_at+' --output '+mapping_file
src_mapY_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordY_map'))
src_at = src_at+' --output '+src_mapY_trans
print(src_at)
try:
subprocess.check_output(src_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# combine X,Y mappings
mapX = load_volume(src_mapX_trans).get_data()
mapY = load_volume(src_mapY_trans).get_data()
src_map = np.stack((mapX,mapY),axis=-1)
mapping = nb.Nifti1Image(src_map, target.affine, target.header)
save_volume(mapping_file, mapping)
trans_mapping = []
trg_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
trg_at = trg_at+' --input '+trg_mapX.get_filename()
trg_at = trg_at+' --reference-image '+source.get_filename()
trg_at = trg_at+' --interpolation Linear'
for idx,transform in enumerate(inverse):
if linear[idx]:
trg_at = trg_at+' --transform ['+transform+', 1]'
else:
trg_at = trg_at+' --transform ['+transform+', 0]'
# trg_at = trg_at+' --output '+inverse_mapping_file
trg_mapX_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordX_map'))
trg_at = trg_at+' --output '+trg_mapX_trans
print(trg_at)
try:
subprocess.check_output(trg_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
trg_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
trg_at = trg_at+' --input '+trg_mapY.get_filename()
trg_at = trg_at+' --reference-image '+source.get_filename()
trg_at = trg_at+' --interpolation Linear'
for idx,transform in enumerate(inverse):
if linear[idx]:
trg_at = trg_at+' --transform ['+transform+', 1]'
else:
trg_at = trg_at+' --transform ['+transform+', 0]'
# trg_at = trg_at+' --output '+inverse_mapping_file
trg_mapY_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordY_map'))
trg_at = trg_at+' --output '+trg_mapY_trans
print(trg_at)
try:
subprocess.check_output(trg_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# combine X,Y mappings
mapX = load_volume(trg_mapX_trans).get_data()
mapY = load_volume(trg_mapY_trans).get_data()
trg_map = np.stack((mapX,mapY),axis=-1)
inverse_mapping = nb.Nifti1Image(trg_map, source.affine, source.header)
save_volume(inverse_mapping_file, inverse_mapping)
# pad coordinate mapping outside the image? hopefully not needed...
# clean-up intermediate files
if os.path.exists(src_mapX_file): os.remove(src_mapX_file)
if os.path.exists(src_mapY_file): os.remove(src_mapY_file)
if os.path.exists(trg_mapX_file): os.remove(trg_mapX_file)
if os.path.exists(trg_mapY_file): os.remove(trg_mapY_file)
if os.path.exists(src_mapX_trans): os.remove(src_mapX_trans)
if os.path.exists(src_mapY_trans): os.remove(src_mapY_trans)
if os.path.exists(trg_mapX_trans): os.remove(trg_mapX_trans)
if os.path.exists(trg_mapY_trans): os.remove(trg_mapY_trans)
if ignore_affine or ignore_header:
if os.path.exists(src_img_file): os.remove(src_img_file)
if os.path.exists(trg_img_file): os.remove(trg_img_file)
for name in forward:
if os.path.exists(name): os.remove(name)
for name in inverse:
if os.path.exists(name): os.remove(name)
if not save_data:
# collect saved outputs
output = {'transformed_source': load_volume(transformed_source_file),
'mapping': load_volume(mapping_file),
'inverse': load_volume(inverse_mapping_file)}
# remove output files if *not* saved
if os.path.exists(transformed_source_file): os.remove(transformed_source_file)
if os.path.exists(mapping_file): os.remove(mapping_file)
if os.path.exists(inverse_mapping_file): os.remove(inverse_mapping_file)
return output
else:
# collect saved outputs
output = {'transformed_source': transformed_source_file,
'mapping': mapping_file,
'inverse': inverse_mapping_file}
return output
def embedded_antsreg_2d_multi(source_images, target_images,
run_rigid=False,
rigid_iterations=1000,
run_affine=False,
affine_iterations=1000,
run_syn=True,
coarse_iterations=40,
medium_iterations=50, fine_iterations=40,
scaling_factor=32,
cost_function='MutualInformation',
interpolation='NearestNeighbor',
regularization='High',
convergence=1e-6,
mask_zero=False,
ignore_affine=False, ignore_orient=False, ignore_res=False,
save_data=False, overwrite=False, output_dir=None,
file_name=None):
""" Embedded ANTS Registration 2D Multi-contrasts
Runs the rigid and/or Symmetric Normalization (SyN) algorithm of ANTs and
formats the output deformations into voxel coordinate mappings as used in
CBSTools registration and transformation routines. Uses all input contrasts
with equal weights.
Parameters
----------
source_images: [niimg]
Images to register
target_images: [niimg]
Reference images to match
run_rigid: bool
Whether or not to run a rigid registration first (default is False)
rigid_iterations: float
Number of iterations in the rigid step (default is 1000)
run_affine: bool
Whether or not to run a affine registration first (default is False)
affine_iterations: float
Number of iterations in the affine step (default is 1000)
run_syn: bool
Whether or not to run a SyN registration (default is True)
coarse_iterations: float
Number of iterations at the coarse level (default is 40)
medium_iterations: float
Number of iterations at the medium level (default is 50)
fine_iterations: float
Number of iterations at the fine level (default is 40)
cost_function: {'CrossCorrelation', 'MutualInformation'}
Cost function for the registration (default is 'MutualInformation')
interpolation: {'NearestNeighbor', 'Linear'}
Interpolation for the registration result (default is 'NearestNeighbor')
convergence: flaot
Threshold for convergence, can make the algorithm very slow
(default is convergence)
ignore_affine: bool
Ignore the affine matrix information extracted from the image header
(default is False)
ignore_orient: bool
Ignore the orientation information and affine matrix information
extracted from the image header (default is False)
ignore_res: bool
Ignore the resolution 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)
* transformed_source (niimg): Deformed source image (_ants-def)
* mapping (niimg): Coordinate mapping from source to target (_ants-map)
* inverse (niimg): Inverse coordinate mapping from target to source
(_ants-invmap)
Notes
----------
Port of the CBSTools Java module by Pierre-Louis Bazin. The main algorithm
is part of the ANTs software by Brian Avants and colleagues [1]_. Parameters
have been set to values commonly found in neuroimaging scripts online, but
not necessarily optimal.
References
----------
.. [1] Avants et al (2008), Symmetric diffeomorphic
image registration with cross-correlation: evaluating automated labeling
of elderly and neurodegenerative brain, Med Image Anal. 12(1):26-41
"""
print('\nEmbedded ANTs Registration 2D Multi-contrasts')
# check if ants is installed to raise sensible error
try:
subprocess.run('antsRegistration', stdout=subprocess.DEVNULL)
except FileNotFoundError:
sys.exit("\nCould not find command 'antsRegistration'. Make sure ANTs is"
" installed and can be accessed from the command line.")
try:
subprocess.run('antsApplyTransforms', stdout=subprocess.DEVNULL)
except FileNotFoundError:
sys.exit("\nCould not find command 'antsApplyTransforms'. Make sure ANTs"
" is installed and can be accessed from the command line.")
# make sure that saving related parameters are correct
# filenames needed for intermediate results
output_dir = _output_dir_4saving(output_dir, source_images[0])
transformed_source_files = []
for idx,source_image in enumerate(source_images):
transformed_source_files.append(os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='ants-def'+str(idx))))
mapping_file = os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='ants-map'))
inverse_mapping_file = os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='ants-invmap'))
if save_data:
if overwrite is False \
and os.path.isfile(mapping_file) \
and os.path.isfile(inverse_mapping_file) :
missing = False
for trans_file in transformed_source_files:
if not os.path.isfile(trans_file):
missing = True
if not missing:
print("skip computation (use existing results)")
transformed = []
for trans_file in transformed_source_files:
transformed.append(trans_file)
output = {'transformed_sources': transformed,
'transformed_source': transformed[0],
'mapping': mapping_file,
'inverse': inverse_mapping_file}
return output
# load and get dimensions and resolution from input images
sources = []
targets = []
src_img_files = []
trg_img_files = []
for idx,img in enumerate(source_images):
source = load_volume(source_images[idx])
src_affine = source.affine
src_header = source.header
nsx = source.header.get_data_shape()[X]
nsy = source.header.get_data_shape()[Y]
nsz = 1
rsx = source.header.get_zooms()[X]
rsy = source.header.get_zooms()[Y]
rsz = 1.0
orig_src_aff = source.affine
orig_src_hdr = source.header
target = load_volume(target_images[idx])
trg_affine = target.affine
trg_header = target.header
ntx = target.header.get_data_shape()[X]
nty = target.header.get_data_shape()[Y]
ntz = 1
rtx = target.header.get_zooms()[X]
rty = target.header.get_zooms()[Y]
rtz = 1.0
orig_trg_aff = target.affine
orig_trg_hdr = target.header
# in case the affine transformations are not to be trusted: make them equal
if ignore_affine or ignore_orient or ignore_res:
mx = np.argmax(np.abs(src_affine[0][0:3]))
my = np.argmax(np.abs(src_affine[1][0:3]))
mz = np.argmax(np.abs(src_affine[2][0:3]))
new_affine = np.zeros((4,4))
if ignore_res:
new_affine[0][:] = src_affine[0][:]/rsx
new_affine[1][:] = src_affine[1][:]/rsy
new_affine[2][:] = src_affine[2][:]/rsz
rsx = 1.0
rsy = 1.0
rsz = 1.0
if ignore_orient:
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
elif ignore_affine:
new_affine[0][mx] = rsx*np.sign(src_affine[0][mx])
new_affine[1][my] = rsy*np.sign(src_affine[1][my])
new_affine[2][mz] = rsz*np.sign(src_affine[2][mz])
if (np.sign(src_affine[0][mx])<0):
new_affine[0][3] = rsx*nsx/2.0
else:
new_affine[0][3] = -rsx*nsx/2.0
if (np.sign(src_affine[1][my])<0):
new_affine[1][3] = rsy*nsy/2.0
else:
new_affine[1][3] = -rsy*nsy/2.0
if (np.sign(src_affine[2][mz])<0):
new_affine[2][3] = rsz*nsz/2.0
else:
new_affine[2][3] = -rsz*nsz/2.0
#if (np.sign(src_affine[0][mx])<0): new_affine[mx][3] = rsx*nsx
#if (np.sign(src_affine[1][my])<0): new_affine[my][3] = rsy*nsy
#if (np.sign(src_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
src_img = nb.Nifti1Image(source.get_data(), new_affine, source.header)
src_img.update_header()
src_img_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_srcimg'+str(idx)))
save_volume(src_img_file, src_img)
source = load_volume(src_img_file)
src_affine = source.affine
src_header = source.header
src_img_files.append(src_img_file)
# create generic affine aligned with the orientation for the target
mx = np.argmax(np.abs(trg_affine[0][0:3]))
my = np.argmax(np.abs(trg_affine[1][0:3]))
mz = np.argmax(np.abs(trg_affine[2][0:3]))
new_affine = np.zeros((4,4))
if ignore_res:
new_affine[0][:] = trg_affine[0][:]/rtx
new_affine[1][:] = trg_affine[1][:]/rty
new_affine[2][:] = trg_affine[2][:]/rtz
rtx = 1.0
rty = 1.0
rtz = 1.0
if ignore_orient:
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
elif ignore_affine:
new_affine[0][mx] = rtx*np.sign(trg_affine[0][mx])
new_affine[1][my] = rty*np.sign(trg_affine[1][my])
new_affine[2][mz] = rtz*np.sign(trg_affine[2][mz])
if (np.sign(trg_affine[0][mx])<0):
new_affine[0][3] = rtx*ntx/2.0
else:
new_affine[0][3] = -rtx*ntx/2.0
if (np.sign(trg_affine[1][my])<0):
new_affine[1][3] = rty*nty/2.0
else:
new_affine[1][3] = -rty*nty/2.0
if (np.sign(trg_affine[2][mz])<0):
new_affine[2][3] = rtz*ntz/2.0
else:
new_affine[2][3] = -rtz*ntz/2.0
#if (np.sign(trg_affine[0][mx])<0): new_affine[mx][3] = rtx*ntx
#if (np.sign(trg_affine[1][my])<0): new_affine[my][3] = rty*nty
#if (np.sign(trg_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
trg_img = nb.Nifti1Image(target.get_data(), new_affine, target.header)
trg_img.update_header()
trg_img_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_trgimg'+str(idx)))
save_volume(trg_img_file, trg_img)
target = load_volume(trg_img_file)
trg_affine = target.affine
trg_header = target.header
trg_img_files.append(trg_img_file)
sources.append(source)
targets.append(target)
# build coordinate mapping matrices and save them to disk
src_coordX = np.zeros((nsx,nsy))
src_coordY = np.zeros((nsx,nsy))
trg_coordX = np.zeros((ntx,nty))
trg_coordY = np.zeros((ntx,nty))
for x in range(nsx):
for y in range(nsy):
src_coordX[x,y] = x
src_coordY[x,y] = y
src_mapX = nb.Nifti1Image(src_coordX, source.affine, source.header)
src_mapX_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_srccoordX'))
save_volume(src_mapX_file, src_mapX)
src_mapY = nb.Nifti1Image(src_coordY, source.affine, source.header)
src_mapY_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_srccoordY'))
save_volume(src_mapY_file, src_mapY)
for x in range(ntx):
for y in range(nty):
trg_coordX[x,y] = x
trg_coordY[x,y] = y
trg_mapX = nb.Nifti1Image(trg_coordX, target.affine, target.header)
trg_mapX_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_trgcoordX'))
save_volume(trg_mapX_file, trg_mapX)
trg_mapY = nb.Nifti1Image(trg_coordY, target.affine, target.header)
trg_mapY_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_trgcoordY'))
save_volume(trg_mapY_file, trg_mapY)
if mask_zero:
# create and save temporary masks
target = targets[0]
trg_mask_data = (target.get_data()!=0)
trg_mask = nb.Nifti1Image(trg_mask_data, target.affine, target.header)
trg_mask_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_trgmask'))
save_volume(trg_mask_file, trg_mask)
source = sources[0]
src_mask_data = (source.get_data()!=0)
src_mask = nb.Nifti1Image(src_mask_data, source.affine, source.header)
src_mask_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_srcmask'))
save_volume(src_mask_file, src_mask)
# run the main ANTS software: here we directly build the command line call
reg = 'antsRegistration --collapse-output-transforms 1 --dimensionality 2' \
+' --initialize-transforms-per-stage 0 --interpolation Linear'
# add a prefix to avoid multiple names?
prefix = _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_syn')
prefix = os.path.basename(prefix)
prefix = prefix.split(".")[0]
reg = reg+' --output '+prefix
if mask_zero:
reg = reg+' --masks ['+trg_mask_file+', '+src_mask_file+']'
srcfiles = []
trgfiles = []
for idx,img in enumerate(sources):
print("registering "+sources[idx].get_filename()+"\n to "+targets[idx].get_filename())
srcfiles.append(sources[idx].get_filename())
trgfiles.append(targets[idx].get_filename())
weight = 1.0/len(srcfiles)
# figure out the number of scales, going with a factor of two
n_scales = math.ceil(math.log(scaling_factor)/math.log(2.0))
iter_rigid = str(rigid_iterations)
iter_affine = str(affine_iterations)
iter_syn = str(coarse_iterations)
smooth = str(scaling_factor)
shrink = str(scaling_factor)
for n in range(n_scales):
iter_rigid = iter_rigid+'x'+str(rigid_iterations)
iter_affine = iter_affine+'x'+str(affine_iterations)
if n<n_scales/2: iter_syn = iter_syn+'x'+str(coarse_iterations)
elif n<n_scales-1: iter_syn = iter_syn+'x'+str(medium_iterations)
else: iter_syn = iter_syn+'x'+str(fine_iterations)
smooth = smooth+'x'+str(scaling_factor/math.pow(2.0,n+1))
shrink = shrink+'x'+str(math.ceil(scaling_factor/math.pow(2.0,n+1)))
# set parameters for all the different types of transformations
if run_rigid is True:
reg = reg + ' --transform Rigid[0.1]'
if (cost_function=='CrossCorrelation'):
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
else:
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric MI['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+iter_rigid+', '+str(convergence)+', 10 ]'
reg = reg + ' --smoothing-sigmas '+smooth
reg = reg + ' --shrink-factors '+shrink
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_affine is True:
reg = reg + ' --transform Affine[0.1]'
if (cost_function=='CrossCorrelation'):
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
else:
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric MI['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+iter_affine+', '+str(convergence)+', 10 ]'
reg = reg + ' --smoothing-sigmas '+smooth
reg = reg + ' --shrink-factors '+shrink
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_syn is True:
if regularization == 'Low': syn_param = [0.2, 1.0, 0.0]
elif regularization == 'Medium': syn_param = [0.2, 3.0, 0.0]
elif regularization == 'High': syn_param = [0.2, 4.0, 3.0]
else: syn_param = [0.2, 3.0, 0.0]
reg = reg + ' --transform SyN'+str(syn_param)
if (cost_function=='CrossCorrelation'):
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
else:
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric MI['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+iter_syn+', '+str(convergence)+', 5 ]'
reg = reg + ' --smoothing-sigmas '+smooth
reg = reg + ' --shrink-factors '+shrink
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_rigid is False and run_affine is False and run_syn is False:
reg = reg + ' --transform Rigid[0.1]'
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
reg = reg + ' --convergence [ 0, 1.0, 2 ]'
reg = reg + ' --smoothing-sigmas 1.0'
reg = reg + ' --shrink-factors 1'
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
reg = reg + ' --write-composite-transform 0'
# run the ANTs command directly
print(reg)
try:
subprocess.check_output(reg, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+str(e.returncode)+')\n Output: '+str(e.output)
raise subprocess.CalledProcessError(msg)
# output file names
results = sorted(glob(prefix+'*'))
forward = []
flag = []
for res in results:
if res.endswith('GenericAffine.mat'):
forward.append(res)
flag.append(False)
elif res.endswith('Warp.nii.gz') and not res.endswith('InverseWarp.nii.gz'):
forward.append(res)
flag.append(False)
#print('forward transforms: '+str(forward))
inverse = []
linear = []
for res in results[::-1]:
if res.endswith('GenericAffine.mat'):
inverse.append(res)
linear.append(True)
elif res.endswith('InverseWarp.nii.gz'):
inverse.append(res)
linear.append(False)
#print('inverse transforms: '+str(inverse))
# Transforms the moving image
for idx,source in enumerate(sources):
at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
at = at+' --input '+sources[idx].get_filename()
at = at+' --reference-image '+targets[idx].get_filename()
at = at+' --interpolation '+interpolation
for idx2,transform in enumerate(forward):
if flag[idx2]:
at = at+' --transform ['+transform+', 1]'
else:
at = at+' --transform ['+transform+', 0]'
at = at+' --output '+transformed_source_files[idx]
print(at)
try:
subprocess.check_output(at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# Create coordinate mappings
src_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
src_at = src_at+' --input '+src_mapX.get_filename()
src_at = src_at+' --reference-image '+target.get_filename()
src_at = src_at+' --interpolation Linear'
for idx,transform in enumerate(forward):
if flag[idx]:
src_at = src_at+' --transform ['+transform+', 1]'
else:
src_at = src_at+' --transform ['+transform+', 0]'
# src_at = src_at+' --output '+mapping_file
src_mapX_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordX_map'))
src_at = src_at+' --output '+src_mapX_trans
print(src_at)
try:
subprocess.check_output(src_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
src_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
src_at = src_at+' --input '+src_mapY.get_filename()
src_at = src_at+' --reference-image '+target.get_filename()
src_at = src_at+' --interpolation Linear'
for idx,transform in enumerate(forward):
if flag[idx]:
src_at = src_at+' --transform ['+transform+', 1]'
else:
src_at = src_at+' --transform ['+transform+', 0]'
# src_at = src_at+' --output '+mapping_file
src_mapY_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordY_map'))
src_at = src_at+' --output '+src_mapY_trans
print(src_at)
try:
subprocess.check_output(src_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# combine X,Y mappings
mapX = load_volume(src_mapX_trans).get_data()
mapY = load_volume(src_mapY_trans).get_data()
src_map = np.stack((mapX,mapY),axis=-1)
mapping = nb.Nifti1Image(src_map, target.affine, target.header)
save_volume(mapping_file, mapping)
trans_mapping = []
trg_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
trg_at = trg_at+' --input '+trg_mapX.get_filename()
trg_at = trg_at+' --reference-image '+source.get_filename()
trg_at = trg_at+' --interpolation Linear'
for idx,transform in enumerate(inverse):
if linear[idx]:
trg_at = trg_at+' --transform ['+transform+', 1]'
else:
trg_at = trg_at+' --transform ['+transform+', 0]'
# trg_at = trg_at+' --output '+inverse_mapping_file
trg_mapX_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordX_map'))
trg_at = trg_at+' --output '+trg_mapX_trans
print(trg_at)
try:
subprocess.check_output(trg_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
trg_at = 'antsApplyTransforms --dimensionality 2 --input-image-type 0'
trg_at = trg_at+' --input '+trg_mapY.get_filename()
trg_at = trg_at+' --reference-image '+source.get_filename()
trg_at = trg_at+' --interpolation Linear'
for idx,transform in enumerate(inverse):
if linear[idx]:
trg_at = trg_at+' --transform ['+transform+', 1]'
else:
trg_at = trg_at+' --transform ['+transform+', 0]'
# trg_at = trg_at+' --output '+inverse_mapping_file
trg_mapY_trans = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='tmp_srccoordY_map'))
trg_at = trg_at+' --output '+trg_mapY_trans
print(trg_at)
try:
subprocess.check_output(trg_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# combine X,Y mappings
mapX = load_volume(trg_mapX_trans).get_data()
mapY = load_volume(trg_mapY_trans).get_data()
trg_map = np.stack((mapX,mapY),axis=-1)
inverse_mapping = nb.Nifti1Image(trg_map, source.affine, source.header)
save_volume(inverse_mapping_file, inverse_mapping)
# pad coordinate mapping outside the image? hopefully not needed...
# clean-up intermediate files
if os.path.exists(src_mapX_file): os.remove(src_mapX_file)
if os.path.exists(src_mapY_file): os.remove(src_mapY_file)
if os.path.exists(trg_mapX_file): os.remove(trg_mapX_file)
if os.path.exists(trg_mapY_file): os.remove(trg_mapY_file)
if os.path.exists(src_mapX_trans): os.remove(src_mapX_trans)
if os.path.exists(src_mapY_trans): os.remove(src_mapY_trans)
if os.path.exists(trg_mapX_trans): os.remove(trg_mapX_trans)
if os.path.exists(trg_mapY_trans): os.remove(trg_mapY_trans)
if ignore_affine or ignore_orient or ignore_res:
for src_img_file in src_img_files:
if os.path.exists(src_img_file): os.remove(src_img_file)
for trg_img_file in trg_img_files:
if os.path.exists(trg_img_file): os.remove(trg_img_file)
for name in forward:
if os.path.exists(name): os.remove(name)
for name in inverse:
if os.path.exists(name): os.remove(name)
# if ignoring header and/or affine, must paste back the correct headers
if ignore_affine or ignore_orient or ignore_res:
mapping = load_volume(mapping_file)
save_volume(mapping_file, nb.Nifti1Image(mapping.get_data(), orig_trg_aff, orig_trg_hdr))
inverse = load_volume(inverse_mapping_file)
save_volume(inverse_mapping_file, nb.Nifti1Image(inverse.get_data(), orig_src_aff, orig_src_hdr))
for trans_file in transformed_source_files:
trans = load_volume(trans_file)
save_volume(trans_file, nb.Nifti1Image(trans.get_data(), orig_trg_aff, orig_trg_hdr))
if not save_data:
# collect saved outputs
transformed = []
for trans_file in transformed_source_files:
transformed.append(load_volume(trans_file))
output = {'transformed_sources': transformed,
'transformed_source': transformed[0],
'mapping': load_volume(mapping_file),
'inverse': load_volume(inverse_mapping_file)}
# remove output files if *not* saved
for idx,trans_image in enumerate(transformed_source_files):
if os.path.exists(trans_image): os.remove(trans_image)
if os.path.exists(mapping_file): os.remove(mapping_file)
if os.path.exists(inverse_mapping_file): os.remove(inverse_mapping_file)
return output
else:
# collect saved outputs
transformed = []
for trans_file in transformed_source_files:
transformed.append(trans_file)
output = {'transformed_sources': transformed,
'transformed_source': transformed[0],
'mapping': mapping_file,
'inverse': inverse_mapping_file}
return output
[docs]def embedded_antsreg_multi(source_images, target_images,
run_rigid=True,
rigid_iterations=1000,
run_affine=False,
affine_iterations=1000,
run_syn=True,
coarse_iterations=40,
medium_iterations=50, fine_iterations=40,
cost_function='MutualInformation',
interpolation='NearestNeighbor',
regularization='High',
convergence=1e-6,
mask_zero=False,
ignore_affine=False, ignore_header=False,
save_data=False, overwrite=False, output_dir=None,
file_name=None):
""" Embedded ANTS Registration Multi-contrasts
Runs the rigid and/or Symmetric Normalization (SyN) algorithm of ANTs and
formats the output deformations into voxel coordinate mappings as used in
CBSTools registration and transformation routines. Uses all input contrasts
with equal weights.
Parameters
----------
source_images: [niimg]
Image list to register
target_images: [niimg]
Reference image list to match
run_rigid: bool
Whether or not to run a rigid registration first (default is False)
rigid_iterations: float
Number of iterations in the rigid step (default is 1000)
run_affine: bool
Whether or not to run a affine registration first (default is False)
affine_iterations: float
Number of iterations in the affine step (default is 1000)
run_syn: bool
Whether or not to run a SyN registration (default is True)
coarse_iterations: float
Number of iterations at the coarse level (default is 40)
medium_iterations: float
Number of iterations at the medium level (default is 50)
fine_iterations: float
Number of iterations at the fine level (default is 40)
cost_function: {'CrossCorrelation', 'MutualInformation'}
Cost function for the registration (default is 'MutualInformation')
interpolation: {'NearestNeighbor', 'Linear'}
Interpolation for the registration result (default is 'NearestNeighbor')
regularization: {'Low', 'Medium', 'High'}
Regularization preset for the SyN deformation (default is 'Medium')
convergence: float
Threshold for convergence, can make the algorithm very slow (default is convergence)
mask_zero: bool
Mask regions with zero value
(default is False)
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)
* transformed_source ([niimg]): Deformed source image list (_ants_def0,1,...)
* mapping (niimg): Coordinate mapping from source to target (_ants_map)
* inverse (niimg): Inverse coordinate mapping from target to source (_ants_invmap)
Notes
----------
Port of the CBSTools Java module by Pierre-Louis Bazin. The main algorithm
is part of the ANTs software by Brian Avants and colleagues [1]_. Parameters
have been set to values commonly found in neuroimaging scripts online, but
not necessarily optimal.
References
----------
.. [1] Avants et al (2008), Symmetric diffeomorphic
image registration with cross-correlation: evaluating automated labeling
of elderly and neurodegenerative brain, Med Image Anal. 12(1):26-41
"""
print('\nEmbedded ANTs Registration Multi-contrasts')
# check if ants is installed to raise sensible error
try:
subprocess.run('antsRegistration', stdout=subprocess.DEVNULL)
except FileNotFoundError:
sys.exit("\nCould not find command 'antsRegistration'. Make sure ANTs is"
" installed and can be accessed from the command line.")
try:
subprocess.run('antsApplyTransforms', stdout=subprocess.DEVNULL)
except FileNotFoundError:
sys.exit("\nCould not find command 'antsApplyTransforms'. Make sure ANTs"
" is installed and can be accessed from the command line.")
# make sure that saving related parameters are correct
# output files needed for intermediate results
output_dir = _output_dir_4saving(output_dir, source_images[0])
transformed_source_files = []
for idx,source_image in enumerate(source_images):
transformed_source_files.append(os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_image,
suffix='ants-def'+str(idx))))
mapping_file = os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='ants-map'))
inverse_mapping_file = os.path.join(output_dir,
_fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='ants-invmap'))
if save_data:
if overwrite is False \
and os.path.isfile(mapping_file) \
and os.path.isfile(inverse_mapping_file) :
missing = False
for trans_file in transformed_source_files:
if not os.path.isfile(trans_file):
missing = True
if not missing:
print("skip computation (use existing results)")
transformed = []
for trans_file in transformed_source_files:
transformed.append(trans_file)
output = {'transformed_sources': transformed,
'transformed_source': transformed[0],
'mapping': mapping_file,
'inverse': inverse_mapping_file}
return output
# load and get dimensions and resolution from input images
sources = []
targets = []
for idx,img in enumerate(source_images):
source = load_volume(source_images[idx])
src_affine = source.affine
src_header = source.header
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]
orig_src_aff = source.affine
orig_src_hdr = source.header
target = load_volume(target_images[idx])
trg_affine = target.affine
trg_header = target.header
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]
orig_trg_aff = target.affine
orig_trg_hdr = target.header
# 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
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:
#mx = np.argmax(np.abs(src_affine[0][0:3]))
#my = np.argmax(np.abs(src_affine[1][0:3]))
#mz = np.argmax(np.abs(src_affine[2][0:3]))
#new_affine[0][mx] = rsx*np.sign(src_affine[0][mx])
#new_affine[1][my] = rsy*np.sign(src_affine[1][my])
#new_affine[2][mz] = rsz*np.sign(src_affine[2][mz])
#if (np.sign(src_affine[0][mx])<0):
# new_affine[0][3] = rsx*nsx/2.0
#else:
# new_affine[0][3] = -rsx*nsx/2.0
#
#if (np.sign(src_affine[1][my])<0):
# new_affine[1][3] = rsy*nsy/2.0
#else:
# new_affine[1][3] = -rsy*nsy/2.0
#
#if (np.sign(src_affine[2][mz])<0):
# new_affine[2][3] = rsz*nsz/2.0
#else:
# new_affine[2][3] = -rsz*nsz/2.0
mx = np.argmax(np.abs([src_affine[0][0],src_affine[1][0],src_affine[2][0]]))
my = np.argmax(np.abs([src_affine[0][1],src_affine[1][1],src_affine[2][1]]))
mz = np.argmax(np.abs([src_affine[0][2],src_affine[1][2],src_affine[2][2]]))
new_affine[mx][0] = rsx*np.sign(src_affine[mx][0])
new_affine[my][1] = rsy*np.sign(src_affine[my][1])
new_affine[mz][2] = rsz*np.sign(src_affine[mz][2])
if (np.sign(src_affine[mx][0])<0):
new_affine[mx][3] = rsx*nsx/2.0
else:
new_affine[mx][3] = -rsx*nsx/2.0
if (np.sign(src_affine[my][1])<0):
new_affine[my][3] = rsy*nsy/2.0
else:
new_affine[my][3] = -rsy*nsy/2.0
if (np.sign(src_affine[mz][2])<0):
new_affine[mz][3] = rsz*nsz/2.0
else:
new_affine[mz][3] = -rsz*nsz/2.0
#if (np.sign(src_affine[0][mx])<0): new_affine[mx][3] = rsx*nsx
#if (np.sign(src_affine[1][my])<0): new_affine[my][3] = rsy*nsy
#if (np.sign(src_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
src_img = nb.Nifti1Image(source.get_data(), new_affine, source.header)
src_img.update_header()
src_img_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_srcimg'+str(idx)))
save_volume(src_img_file, src_img)
source = load_volume(src_img_file)
src_affine = source.affine
src_header = source.header
# create generic affine aligned with the orientation for the target
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:
#mx = np.argmax(np.abs(trg_affine[0][0:3]))
#my = np.argmax(np.abs(trg_affine[1][0:3]))
#mz = np.argmax(np.abs(trg_affine[2][0:3]))
#new_affine[0][mx] = rtx*np.sign(trg_affine[0][mx])
#new_affine[1][my] = rty*np.sign(trg_affine[1][my])
#new_affine[2][mz] = rtz*np.sign(trg_affine[2][mz])
#if (np.sign(trg_affine[0][mx])<0):
# new_affine[0][3] = rtx*ntx/2.0
#else:
# new_affine[0][3] = -rtx*ntx/2.0
#
#if (np.sign(trg_affine[1][my])<0):
# new_affine[1][3] = rty*nty/2.0
#else:
# new_affine[1][3] = -rty*nty/2.0
#
#if (np.sign(trg_affine[2][mz])<0):
# new_affine[2][3] = rtz*ntz/2.0
#else:
# new_affine[2][3] = -rtz*ntz/2.0
mx = np.argmax(np.abs([trg_affine[0][0],trg_affine[1][0],trg_affine[2][0]]))
my = np.argmax(np.abs([trg_affine[0][1],trg_affine[1][1],trg_affine[2][1]]))
mz = np.argmax(np.abs([trg_affine[0][2],trg_affine[1][2],trg_affine[2][2]]))
#print('mx: '+str(mx)+', my: '+str(my)+', mz: '+str(mz))
#print('rx: '+str(rtx)+', ry: '+str(rty)+', rz: '+str(rtz))
new_affine[mx][0] = rtx*np.sign(trg_affine[mx][0])
new_affine[my][1] = rty*np.sign(trg_affine[my][1])
new_affine[mz][2] = rtz*np.sign(trg_affine[mz][2])
if (np.sign(trg_affine[mx][0])<0):
new_affine[mx][3] = rtx*ntx/2.0
else:
new_affine[mx][3] = -rtx*ntx/2.0
if (np.sign(trg_affine[my][1])<0):
new_affine[my][3] = rty*nty/2.0
else:
new_affine[my][3] = -rty*nty/2.0
if (np.sign(trg_affine[mz][2])<0):
new_affine[mz][3] = rtz*ntz/2.0
else:
new_affine[mz][3] = -rtz*ntz/2.0
#if (np.sign(trg_affine[0][mx])<0): new_affine[mx][3] = rtx*ntx
#if (np.sign(trg_affine[1][my])<0): new_affine[my][3] = rty*nty
#if (np.sign(trg_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
#print("\nbefore: "+str(trg_affine))
#print("\nafter: "+str(new_affine))
trg_img = nb.Nifti1Image(target.get_data(), new_affine, target.header)
trg_img.update_header()
trg_img_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_trgimg'+str(idx)))
save_volume(trg_img_file, trg_img)
target = load_volume(trg_img_file)
trg_affine = target.affine
trg_header = target.header
sources.append(source)
targets.append(target)
# build coordinate mapping matrices and save them to disk
src_coord = np.zeros((nsx,nsy,nsz,3))
trg_coord = np.zeros((ntx,nty,ntz,3))
for x in range(nsx):
for y in range(nsy):
for z in range(nsz):
src_coord[x,y,z,X] = x
src_coord[x,y,z,Y] = y
src_coord[x,y,z,Z] = z
src_map = nb.Nifti1Image(src_coord, source.affine, source.header)
src_map_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_srccoord'))
save_volume(src_map_file, src_map)
for x in range(ntx):
for y in range(nty):
for z in range(ntz):
trg_coord[x,y,z,X] = x
trg_coord[x,y,z,Y] = y
trg_coord[x,y,z,Z] = z
trg_map = nb.Nifti1Image(trg_coord, target.affine, target.header)
trg_map_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_trgcoord'))
save_volume(trg_map_file, trg_map)
if mask_zero:
# create and save temporary masks
target = targets[0]
trg_mask_data = (target.get_data()!=0)
trg_mask = nb.Nifti1Image(trg_mask_data, target.affine, target.header)
trg_mask_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_trgmask'))
save_volume(trg_mask_file, trg_mask)
source = sources[0]
src_mask_data = (source.get_data()!=0)
src_mask = nb.Nifti1Image(src_mask_data, source.affine, source.header)
src_mask_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_srcmask'))
save_volume(src_mask_file, src_mask)
# run the main ANTS software: here we directly build the command line call
reg = 'antsRegistration --collapse-output-transforms 1 --dimensionality 3' \
+' --initialize-transforms-per-stage 0 --interpolation Linear'
# add a prefix to avoid multiple names?
prefix = _fname_4saving(module=__name__,file_name=file_name,
rootfile=source_images[0],
suffix='tmp_syn')
prefix = os.path.basename(prefix)
prefix = prefix.split(".")[0]
#reg.inputs.output_transform_prefix = prefix
reg = reg+' --output '+prefix
if mask_zero:
reg = reg+' --masks ['+trg_mask_file+', '+src_mask_file+']'
srcfiles = []
trgfiles = []
for idx,img in enumerate(sources):
print("registering "+sources[idx].get_filename()+"\n to "+targets[idx].get_filename())
srcfiles.append(sources[idx].get_filename())
trgfiles.append(targets[idx].get_filename())
weight = 1.0/len(srcfiles)
# set parameters for all the different types of transformations
if run_rigid is True:
reg = reg + ' --transform Rigid[0.1]'
if (cost_function=='CrossCorrelation'):
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
else:
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric MI['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+str(rigid_iterations)+'x' \
+str(rigid_iterations)+'x'+str(rigid_iterations) \
+', '+str(convergence)+', 10 ]'
reg = reg + ' --smoothing-sigmas 3.0x2.0x1.0'
reg = reg + ' --shrink-factors 4x2x1'
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_affine is True:
reg = reg + ' --transform Affine[0.1]'
if (cost_function=='CrossCorrelation'):
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
else:
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric MI['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+str(affine_iterations)+'x' \
+str(affine_iterations)+'x'+str(affine_iterations) \
+', '+str(convergence)+', 10 ]'
reg = reg + ' --smoothing-sigmas 3.0x2.0x1.0'
reg = reg + ' --shrink-factors 4x2x1'
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_syn is True:
if regularization == 'Low': syn_param = [0.2, 1.0, 0.0]
elif regularization == 'Medium': syn_param = [0.2, 3.0, 0.0]
elif regularization == 'High': syn_param = [0.2, 4.0, 3.0]
else: syn_param = [0.2, 3.0, 0.0]
reg = reg + ' --transform SyN'+str(syn_param)
if (cost_function=='CrossCorrelation'):
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
else:
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric MI['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 32, Random, 0.3 ]'
reg = reg + ' --convergence ['+str(coarse_iterations)+'x' \
+str(coarse_iterations)+'x'+str(medium_iterations)+'x' \
+str(fine_iterations)+', '+str(convergence)+', 5 ]'
reg = reg + ' --smoothing-sigmas 2.0x1.0x0.5x0.0'
reg = reg + ' --shrink-factors 8x4x2x1'
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
if run_rigid is False and run_affine is False and run_syn is False:
reg = reg + ' --transform Rigid[0.1]'
for idx,img in enumerate(srcfiles):
reg = reg + ' --metric CC['+trgfiles[idx]+', '+srcfiles[idx] \
+', '+'{:.3f}'.format(weight)+', 5, Random, 0.3 ]'
reg = reg + ' --convergence [ 0x0x0, 1.0, 2 ]'
reg = reg + ' --smoothing-sigmas 3.0x2.0x1.0'
reg = reg + ' --shrink-factors 4x2x1'
reg = reg + ' --use-histogram-matching 0'
reg = reg + ' --winsorize-image-intensities [ 0.001, 0.999 ]'
reg = reg + ' --write-composite-transform 0'
# run the ANTs command directly
print(reg)
try:
subprocess.check_output(reg, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+str(e.returncode)+')\n Output: '+str(e.output)
raise subprocess.CalledProcessError(msg)
# output file names
results = sorted(glob(prefix+'*'))
forward = []
flag = []
for res in results:
if res.endswith('GenericAffine.mat'):
forward.append(res)
flag.append(False)
elif res.endswith('Warp.nii.gz') and not res.endswith('InverseWarp.nii.gz'):
forward.append(res)
flag.append(False)
#print('forward transforms: '+str(forward))
inverse = []
linear = []
for res in results[::-1]:
if res.endswith('GenericAffine.mat'):
inverse.append(res)
linear.append(True)
elif res.endswith('InverseWarp.nii.gz'):
inverse.append(res)
linear.append(False)
#print('inverse transforms: '+str(inverse))
# Transforms the moving image
for idx,source in enumerate(sources):
at = 'antsApplyTransforms --dimensionality 3 --input-image-type 0'
at = at+' --input '+sources[idx].get_filename()
at = at+' --reference-image '+targets[idx].get_filename()
at = at+' --interpolation '+interpolation
for idx2,transform in enumerate(forward):
if flag[idx2]:
at = at+' --transform ['+transform+', 1]'
else:
at = at+' --transform ['+transform+', 0]'
at = at+' --output '+transformed_source_files[idx]
print(at)
try:
subprocess.check_output(at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# Create coordinate mappings
src_at = 'antsApplyTransforms --dimensionality 3 --input-image-type 3'
src_at = src_at+' --input '+src_map.get_filename()
src_at = src_at+' --reference-image '+target.get_filename()
src_at = src_at+' --interpolation Linear'
for idx,transform in enumerate(forward):
if flag[idx]:
src_at = src_at+' --transform ['+transform+', 1]'
else:
src_at = src_at+' --transform ['+transform+', 0]'
src_at = src_at+' --output '+mapping_file
print(src_at)
try:
subprocess.check_output(src_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
trans_mapping = []
trg_at = 'antsApplyTransforms --dimensionality 3 --input-image-type 3'
trg_at = trg_at+' --input '+trg_map.get_filename()
trg_at = trg_at+' --reference-image '+source.get_filename()
trg_at = trg_at+' --interpolation Linear'
for idx,transform in enumerate(inverse):
if linear[idx]:
trg_at = trg_at+' --transform ['+transform+', 1]'
else:
trg_at = trg_at+' --transform ['+transform+', 0]'
trg_at = trg_at+' --output '+inverse_mapping_file
print(trg_at)
try:
subprocess.check_output(trg_at, shell=True, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
msg = 'execution failed (error code '+e.returncode+')\n Output: '+e.output
raise subprocess.CalledProcessError(msg)
# pad coordinate mapping outside the image? hopefully not needed...
# clean-up intermediate files
if os.path.exists(src_map_file): os.remove(src_map_file)
if os.path.exists(trg_map_file): os.remove(trg_map_file)
if ignore_affine or ignore_header:
if os.path.exists(src_img_file): os.remove(src_img_file)
if os.path.exists(trg_img_file): os.remove(trg_img_file)
for name in forward:
if os.path.exists(name): os.remove(name)
for name in inverse:
if os.path.exists(name): os.remove(name)
# if ignoring header and/or affine, must paste back the correct headers
if ignore_affine or ignore_header:
mapping = load_volume(mapping_file)
save_volume(mapping_file, nb.Nifti1Image(mapping.get_data(), orig_trg_aff, orig_trg_hdr))
inverse = load_volume(inverse_mapping_file)
save_volume(inverse_mapping_file, nb.Nifti1Image(inverse.get_data(), orig_src_aff, orig_src_hdr))
for trans_file in transformed_source_files:
trans = load_volume(trans_file)
save_volume(trans_file, nb.Nifti1Image(trans.get_data(), orig_trg_aff, orig_trg_hdr))
if not save_data:
# collect saved outputs
transformed = []
for trans_file in transformed_source_files:
transformed.append(load_volume(trans_file))
output = {'transformed_sources': transformed,
'transformed_source': transformed[0],
'mapping': load_volume(mapping_file),
'inverse': load_volume(inverse_mapping_file)}
# remove output files if *not* saved
for idx,trans_image in enumerate(transformed_source_files):
if os.path.exists(trans_image): os.remove(trans_image)
if os.path.exists(mapping_file): os.remove(mapping_file)
if os.path.exists(inverse_mapping_file): os.remove(inverse_mapping_file)
return output
else:
# collect saved outputs
transformed = []
for trans_file in transformed_source_files:
transformed.append(trans_file)
output = {'transformed_sources': transformed,
'transformed_source': transformed[0],
'mapping': mapping_file,
'inverse': inverse_mapping_file}
return output