faster voxel loading

tfce_mediation/misc_scripts/tm_massunivariatemodels.py

@@ -174,7 +174,7 @@ def gram_schmidt_orthonorm(X, columns=True):
 	else:
 		Q, _ = np.linalg.qr(X.T)
 		Q = Q.T
-	return Q
+	return -Q
 
 def full_glm_results(endog_arr, exog_vars, return_resids = False, only_tvals = False, PCA_whiten = False, ZCA_whiten = False,  orthogonalize = True, orthogNear = False, orthog_GramSchmidt = False):
 	if np.mean(exog_vars[:,0])!=1:
@@ -184,10 +184,12 @@ def full_glm_results(endog_arr, exog_vars, return_resids = False, only_tvals = F
 
 	if orthogonalize:
 		exog_vars = sm.add_constant(orthog_columns(exog_vars[:,1:]))
-	if orthogNear:
-		exog_vars = sm.add_constant(orthog_columns(ortho_neareast[:,1:]))
-	if orthog_GramSchmidt: # for when order matters AKA type 2 sum of squares
-		exog_vars = sm.add_constant(orthog_columns(gram_schmidt_orthonorm[:,1:]))
+	elif orthogNear:
+		exog_vars = sm.add_constant(ortho_neareast(exog_vars[:,1:]))
+	elif orthog_GramSchmidt: # for when order matters AKA type 2 sum of squares
+		exog_vars = sm.add_constant(gram_schmidt_orthonorm(exog_vars[:,1:]))
+	else:
+		pass
 
 	invXX = np.linalg.inv(np.dot(exog_vars.T, exog_vars))
 
@@ -303,10 +305,8 @@ def run_permutations_med(endog_arr, exog_vars, medtype, leftvar, rightvar, num_p
 				for m in mixed_blocks:
 					rotate_groups.append(index_groups[uniq_groups == unique_blocks[m]])
 				index_groups = np.array(rotate_groups).flatten()
-			nx = exog_vars[index_groups]
 		else:
 			index_groups = np.random.permutation(list(range(n)))
-			nx = exog_vars[index_groups]
 
 		if medtype == 'I':
 			EXOG_A = sm.add_constant(np.column_stack((leftvar, strip_ones(exog_vars))))
@@ -456,7 +456,7 @@ def run(opts):
 			inteaction_vars = int_terms.split("*")
 			for scale_var in inteaction_vars:
 				var_temp = scalevar(pdCSV[scale_var])
-				var_tempname = '%s_p' % scale_var
+				var_tempname = '%s_std' % scale_var
 				if var_tempname in opts.exogenousvariables:
 					pass
 				else:
@@ -466,10 +466,10 @@ def run(opts):
 			inteaction_vars = int_terms.split("*")
 			for i, scale_var in enumerate(inteaction_vars):
 				if i == 0:
-					int_temp = pdCSV['%s_p' % scale_var]
+					int_temp = pdCSV['%s_std' % scale_var]
 					int_tempname = '%s' % scale_var
 				else:
-					int_temp = int_temp * pdCSV['%s_p' % scale_var]
+					int_temp = int_temp * pdCSV['%s_std' % scale_var]
 					int_tempname = int_tempname + '.X.' + scale_var
 			if int_tempname in opts.exogenousvariables:
 				pass

tfce_mediation/pyfunc.py

@@ -23,13 +23,13 @@
 import math
 import sys
 import struct
+import uuid
 from scipy.stats import linregress
 import matplotlib.pyplot as plt
 import matplotlib.colors as colors
 import matplotlib.patches as mpatches
 from time import time
 
-
 from tfce_mediation.cynumstats import calc_beta_se, resid_covars, cy_lin_lstsqr_mat
 
 # Creation of adjacencty sets for TFCE connectivity
@@ -1103,4 +1103,36 @@ def erode_3D_image(img_data, erode_iter = 2, do_binary_opening = False, do_binar
 	img_data *= erode_mask
 	return img_data, erode_mask
 
+# load voxel neuorimages, checks for file size to conserve RAM, return nibabel img or masked data
+def import_voxel_neuroimage(image_path, mask_index = None):
+	if not os.path.exists(image_path):
+		print('Error: %s not found' % image_path)
+		quit()
+	base, file_ext = os.path.splitext(image_path)
+	if file_ext == '.gz':
+		file_ext = os.path.splitext(base)[1]
+		if file_ext == '.nii':
+			if os.path.getsize(image_path) < 50000000:
+				image = nib.load(image_path)
+			else:
+				tempname = str(uuid.uuid4().hex) + '.nii'
+				os.system("zcat %s > %s" % (image_path,tempname))
+				image = nib.load(tempname)
+				os.system("rm %s" % tempname)
+		else:
+			print('Error: filetype for %s is not supported' % image_path)
+	elif file_ext == '.nii':
+		image = nib.load(image_path)
+	elif file_ext == '.mnc':
+		image = nib.load(image_path)
+	else:
+		print('Error filetype for %s is not supported' % img_all_name)
+		quit()
+	print("Imported:\t%s" % image_path)
+	if mask_index is not None:
+		image_data = image.get_data()
+		image_data = image_data[mask_index]
+		return image_data
+	else:
+		return image
 

tfce_mediation/tmanalysis/STEP_0_load_volumes.py

@@ -19,7 +19,10 @@
 import os
 import numpy as np
 import nibabel as nib
-import argparse 
+import argparse
+
+from tfce_mediation.pyfunc import import_voxel_neuroimage
+
 
 DESCRIPTION = "Initial step to load in a 4D NifTi or MINC volumetric image, and its corresponding mask (binary image)."
 
@@ -29,39 +32,48 @@ def getArgumentParser(parser = argparse.ArgumentParser(description = DESCRIPTION
 		help="[4D_image] [Mask] (default: %(default)s)", 
 		metavar=('*.nii.gz or *.mnc', '*.nii.gz or *.mnc'), 
 		default=['all_FA_skeletonised.nii.gz','mean_FA_skeleton_mask.nii.gz'])
+	inputtype = parser.add_mutually_exclusive_group(required=False)
+	inputtype.add_argument("-f", "--files", 
+		nargs='+', 
+		help="Neuroimage file location with wildcards. USAGE: -f {neuroimage} ...", 
+		metavar=('*.nii.gz or *.mnc'))
+	inputtype.add_argument("-l", "--filelist", 
+		nargs=1, 
+		help="Import neuroimages from text file. -l {1D text file}", 
+		metavar=('*.csv'))
+	parser.add_argument("-m", "--mask", 
+		nargs=1, 
+		help="Must be used with --files or --filelist. USAGE: -m {Mask}", 
+		metavar=('*.nii.gz or *.mnc'))
 	return parser
 
 def run(opts):
-	mask_name = opts.input[1]
+	if opts.mask:
+		mask_name = opts.mask[0]
+	else:
+		mask_name = opts.input[1]
 	if not os.path.isfile(mask_name):
-		print('Error %s not found. Please use -i option.' % mask_name)
+		print('Error: %s not found. Please use -i or -m option.' % mask_name)
 		quit()
-	img_mask = nib.load(mask_name)
+	img_mask = import_voxel_neuroimage(mask_name)
 	data_mask = img_mask.get_data()
 	affine_mask = img_mask.get_affine()
 	header_mask = img_mask.get_header()
+	mask_index = data_mask > 0.99
 
-	#check if minc file
-	img_all_name = opts.input[0]
-	_, file_ext = os.path.splitext(img_all_name)
-	if file_ext == '.gz':
-		_, file_ext = os.path.splitext(img_all_name)
-		if file_ext == '.mnc':
-			imgext = '.mnc'
-			img_all = nib.load(img_all_name)
-		else:
-			imgext = '.nii.gz'
-			os.system("zcat %s > temp_4d.nii" % img_all_name)
-			img_all = nib.load('temp_4d.nii')
-	elif file_ext == '.nii':
-		imgext = '.nii.gz' # default to zipped images
-		img_all = nib.load(img_all_name)
+	if opts.files:
+		nonzero_data = []
+		for image_path in opts.files:
+			nonzero_data.append(import_voxel_neuroimage(image_path, mask_index))
+		nonzero_data = np.array(nonzero_data).T
+	elif opts.filelist:
+		nonzero_data = []
+		for image_path in np.genfromtxt(opts.filelist[0], delimiter=','):
+			nonzero_data.append(import_voxel_neuroimage(image_path, mask_index))
+		nonzero_data = np.array(nonzero_data).T
 	else:
-		print('Error filetype for %s is not supported' % img_all_name)
-		quit()
+		nonzero_data = import_voxel_neuroimage(opts.input[0], mask_index)
 
-	data_all = img_all.get_data()
-	nonzero_data = data_all[data_mask>0.99]
 	# check mask
 	mean = np.mean(nonzero_data, axis=1)
 	if len(mean[mean==0]) != 0:
@@ -79,14 +91,11 @@ def run(opts):
 	np.save('python_temp/header_mask',header_mask)
 	np.save('python_temp/affine_mask',affine_mask)
 	np.save('python_temp/data_mask',data_mask)
-	np.save('python_temp/imgext',imgext)
 	num_voxel = nonzero_data.shape[0]
 	num_subjects = nonzero_data.shape[1]
 	np.save('python_temp/num_voxel',num_voxel)
 	np.save('python_temp/num_subjects',num_subjects)
 
-	os.system("rm temp_4d.nii")
-
 if __name__ == "__main__":
 	parser = getArgumentParser()
 	opts = parser.parse_args()