diff --git a/ciftify/bin/ciftify_falff.py b/ciftify/bin/ciftify_falff.py index 018656c..b061953 100755 --- a/ciftify/bin/ciftify_falff.py +++ b/ciftify/bin/ciftify_falff.py @@ -21,24 +21,23 @@ -h,--help Print help """ - +import logging import os -import numpy as np + +from docopt import docopt import nibabel as nib +import numpy as np from scipy.fftpack import fft -from docopt import docopt -from ciftify.utils import run, TempDir -from ciftify.meants import NibInput + import ciftify.config -import logging -import sys +from ciftify.utils import run, TempDir config_path = os.path.join(os.path.dirname(ciftify.config.find_ciftify_global()), 'bin', "logging.conf") logging.config.fileConfig(config_path, disable_existing_loggers=False) logger = logging.getLogger(os.path.basename(__file__)) + def main(): - '''''' arguments = docopt(__doc__) funcfile = arguments[''] outputname = arguments[''] @@ -59,45 +58,53 @@ def main(): with ciftify.utils.TempDir() as tmpdir: - # IF INPUT IS A NIFTI FILE - # Sets input funcfile equal to inputfile - func = NibInput(funcfile) - - - # IF INPUT IS CIFTI FILE - # Convert cifti input file to nifti input file - if func.type == "cifti": - inputfile = convert_cifti_to_nifti(func.path, tmpdir) - if maskfile: - maskinput = convert_cifti_to_nifti(maskfile, tmpdir) - else: - maskinput = None - elif func.type == "nifti": - inputfile = func.path - if maskfile: - maskinput = maskfile - else: - maskinput = None - else: - logger.critical("Could not read as nifti or cifti file") - sys.exit(1) - - - falff_nifti_output = calc_nifti(inputfile, maskinput, min_low_freq, max_low_freq, min_total_freq, max_total_freq, tmpdir, calc_alff) + func, inputfile, maskinput, repetition_time = get_inputs( + funcfile, maskfile, tmpdir) + + falff_nifti_output = calc_nifti( + inputfile, + maskinput, + min_low_freq, + max_low_freq, + min_total_freq, + max_total_freq, + tmpdir, + calc_alff, + repetition_time + ) # Convert nifti output file to cifti output file - if func.type == "cifti": + if isinstance(func, nib.Cifti2Image): convert_nifti_to_cifti(falff_nifti_output, funcfile, outputname) # IF INPUT IS NIFTI FILE # If funcfile was not cifti file, save as nifti file to outputname - if func.type == "nifti": + if isinstance(func, nib.Nifti1Image): run("mv {} {}".format(falff_nifti_output, outputname)) -# If input is cifti - convert to fake nifti (fake_nifti_input) -# Convert to nifti +def get_inputs(funcfile, maskfile, tmpdir): + """Ensure funcfile and maskfile are niftis and get repetition time. + """ + func = nib.load(funcfile) + + if isinstance(func, nib.Cifti2Image): + inputfile = convert_cifti_to_nifti(funcfile, tmpdir) + mask = convert_cifti_to_nifti(maskfile, tmpdir) if maskfile else None + repetition_time = func.header.get_axis(0).step + elif isinstance(func, nib.Nifti1Image): + inputfile = funcfile + mask = maskfile + repetition_time = func.header['pixdim'][4] + else: + raise IOError("Could not read as nifti or cifti file") + + return func, inputfile, mask, repetition_time + + def convert_cifti_to_nifti(funcfile, tmpdir): + """ If input is cifti - convert to fake nifti. + """ fake_nifti_input = os.path.join(tmpdir, 'input_fake.nii.gz') run('wb_command -cifti-convert -to-nifti {} {} '.format(funcfile, fake_nifti_input)) return fake_nifti_input @@ -108,7 +115,8 @@ def convert_nifti_to_cifti(falff_nifti_output, funcfile, outputname): run('wb_command -cifti-convert -from-nifti {} {} {} -reset-scalars'.format(falff_nifti_output, funcfile, outputname)) -def calc_nifti(inputfile, maskfile, min_low_freq, max_low_freq, min_total_freq, max_total_freq, tmpdir, calc_alff): +def calc_nifti(inputfile, maskfile, min_low_freq, max_low_freq, min_total_freq, + max_total_freq, tmpdir, calc_alff, repetition_time): ''' calculates falff from nifti input and retruns nifti output @@ -141,7 +149,15 @@ def calc_nifti(inputfile, maskfile, min_low_freq, max_low_freq, min_total_freq, # Loop through x,y,z indices, send to calculate_falff function for x,y,z in zip(indx,indy,indz): - falff_vol[x,y,z] = calculate_falff(func_data[x,y,z,:], min_low_freq, max_low_freq, min_total_freq, max_total_freq, calc_alff) + falff_vol[x,y,z] = calculate_falff( + func_data[x, y, z, :], + min_low_freq, + max_low_freq, + min_total_freq, + max_total_freq, + calc_alff, + repetition_time + ) # Save falff values to fake nifti output temp file output_3D = nib.Nifti1Image(falff_vol, affine) @@ -151,27 +167,31 @@ def calc_nifti(inputfile, maskfile, min_low_freq, max_low_freq, min_total_freq, return falff_nifti_output # CALCULATES FALFF -def calculate_falff(timeseries, min_low_freq, max_low_freq, min_total_freq, max_total_freq, calc_alff): - ''' this will calculate falff from a timeseries''' +def calculate_falff(timeseries, min_low_freq, max_low_freq, min_total_freq, + max_total_freq, calc_alff, repetition_time): + """This will calculate falff from a timeseries""" n = len(timeseries) - time = (np.arange(n))*2 # Takes fast Fourier transform of timeseries fft_timeseries = fft(timeseries) - # Calculates frequency scale - freq_scale = np.fft.fftfreq(n, 1/1) - + freq_scale = np.fft.fftfreq(n, repetition_time) # Calculates power of fft mag = (abs(fft_timeseries))**0.5 # Finds low frequency range (0.01-0.08) and total frequency range (0.0-0.25) - low_ind = np.where((float(min_low_freq) <= freq_scale) & (freq_scale <= float(max_low_freq))) - total_ind = np.where((float(min_total_freq) <= freq_scale) & (freq_scale <= float(max_total_freq))) + low_ind = np.where( + (float(min_low_freq) <= freq_scale) + & (freq_scale <= float(max_low_freq)) + ) + total_ind = np.where( + (float(min_total_freq) <= freq_scale[1:]) + & (freq_scale[1:] <= float(max_total_freq)) + ) # Indexes power to low frequency index, total frequency range low_power = mag[low_ind] - total_power = mag[total_ind] + total_power = mag[1:][total_ind] # Calculates sum of lower power and total power low_pow_sum = np.sum(low_power) total_pow_sum = np.sum(total_power)