Source code for nighres.intensity.t2s_optimal_combination

import numpy as np
import nibabel as nb
import os
import sys
import nighresjava
from import load_volume, save_volume
from ..utils import _output_dir_4saving, _fname_4saving, \
                    _check_topology_lut_dir, _check_available_memory

[docs]def t2s_optimal_combination(image_list, te_list, depth=None, save_data=False, overwrite=False, output_dir=None, file_name=None): """ T2*-Optimal Combination Estimate T2*/R2* and use it to combine multi-echo data Parameters ---------- image_list: [niimg] List of 4D input images, one per echo time te_list: [float] List of input echo times (TE) depth: [int] List of echo depth to keep for input time points save_data: bool, optional Save output data to file (default is False) overwrite: bool, optional 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) * combined (niimg): Optimally combined 4D image * t2s (niimg): Map of estimated T2* times (_qt2fit-t2s) * r2s (niimg): Map of estimated R2* relaxation rate (_qt2fit-r2s) * s0 (niimg): Estimated PD weighted image at TE=0 (_qt2fit-s0) * residuals (niimg): Estimated residuals between input and estimated echoes (_qt2fit-err) Notes ---------- Original Java module by Pierre-Louis Bazin. """ print('\nT2* Optimal Combination') # make sure that saving related parameters are correct if save_data: output_dir = _output_dir_4saving(output_dir, image_list[0]) comb_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name, rootfile=image_list[0], suffix='qt2scomb-combined')) t2s_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name, rootfile=image_list[0], suffix='qt2scomb-t2s')) r2s_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name, rootfile=image_list[0], suffix='qt2scomb-r2s')) s0_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name, rootfile=image_list[0], suffix='qt2scomb-s0')) err_file = os.path.join(output_dir, _fname_4saving(module=__name__,file_name=file_name, rootfile=image_list[0], suffix='qt2scomb-err')) if overwrite is False \ and os.path.isfile(comb_file) \ and os.path.isfile(t2s_file) \ and os.path.isfile(r2s_file) \ and os.path.isfile(s0_file) \ and os.path.isfile(err_file) : output = {'combined': comb_file, 't2s': t2s_file, 'r2s': r2s_file, 's0': s0_file, 'residuals': err_file} return output # start virtual machine, if not already running try: mem = _check_available_memory() nighresjava.initVM(initialheap=mem['init'], maxheap=mem['max']) except ValueError: pass # create algorithm instance qt2scomb = nighresjava.T2sOptimalCombination() # set algorithm parameters qt2scomb.setNumberOfEchoes(len(image_list)) # load first image and use it to set dimensions and resolution img = load_volume(image_list[0]) data = img.get_data() #data = data[0:10,0:10,0:10] affine = img.affine header = img.header resolution = [x.item() for x in header.get_zooms()] dimensions = data.shape dim3d = (dimensions[0], dimensions[1], dimensions[2]) if len(dimensions) is 3: qt2scomb.setDimensions(dimensions[0], dimensions[1], dimensions[2]) else: qt2scomb.setDimensions(dimensions[0], dimensions[1], dimensions[2], dimensions[3]) qt2scomb.setResolutions(resolution[0], resolution[1], resolution[2]) # input images # important: set image number before adding images for idx, image in enumerate(image_list): #print('\nloading ('+str(idx)+'): '+image) data = load_volume(image).get_data() #data = data[0:10,0:10,0:10] qt2scomb.setEchoImageAt(idx, nighresjava.JArray('float')( (data.flatten('F')).astype(float))) qt2scomb.setEchoTimeAt(idx, te_list[idx]) if depth is not None: qt2scomb.setImageEchoDepth(nighresjava.JArray('int')(depth)) # execute the algorithm try: qt2scomb.execute() except: # if the Java module fails, reraise the error it throws print("\n The underlying Java code did not execute cleanly: ") print(sys.exc_info()[0]) raise return # reshape output to what nibabel likes comb_data = np.reshape(np.array(qt2scomb.getCombinedImage(), dtype=np.float32), dimensions, 'F') t2s_data = np.reshape(np.array(qt2scomb.getT2sImage(), dtype=np.float32), dim3d, 'F') r2s_data = np.reshape(np.array(qt2scomb.getR2sImage(), dtype=np.float32), dim3d, 'F') s0_data = np.reshape(np.array(qt2scomb.getS0Image(), dtype=np.float32), dim3d, 'F') err_data = np.reshape(np.array(qt2scomb.getResidualImage(), dtype=np.float32), dim3d, 'F') # adapt header max for each image so that correct max is displayed # and create nifiti objects header['cal_min'] = np.nanmin(comb_data) header['cal_max'] = np.nanmax(comb_data) comb = nb.Nifti1Image(comb_data, affine, header) header['cal_min'] = np.nanmin(t2s_data) header['cal_max'] = np.nanmax(t2s_data) t2s = nb.Nifti1Image(t2s_data, affine, header) header['cal_min'] = np.nanmin(r2s_data) header['cal_max'] = np.nanmax(r2s_data) r2s = nb.Nifti1Image(r2s_data, affine, header) header['cal_min'] = np.nanmin(s0_data) header['cal_max'] = np.nanmax(s0_data) s0 = nb.Nifti1Image(s0_data, affine, header) header['cal_min'] = np.nanmin(err_data) header['cal_max'] = np.nanmax(err_data) err = nb.Nifti1Image(err_data, affine, header) if save_data: save_volume(comb_file, comb) save_volume(t2s_file, t2s) save_volume(r2s_file, r2s) save_volume(s0_file, s0) save_volume(err_file, err) return {'combined': comb_file, 't2s': t2s_file, 'r2s': r2s_file, 's0': s0_file, 'residuals': err_file} else: return {'combined': comb, 't2s': t2s, 'r2s': r2s, 's0': s0, 'residuals': err}