How to use the nibabel.freesurfer.mghformat.load function in nibabel

def run(opts):

	allsurf = opts.input[0]
	#load surface data
	img_data = nib.freesurfer.mghformat.load(allsurf)
	img_data_full = img_data.get_data()

	img_data = np.squeeze(img_data_full)
	n = img_data.shape[1]

	if opts.label:
		label = opts.label[0]
		lable_name = os.path.basename(opts.label[0].split('.mgh',1)[0])
		img_label = nib.freesurfer.mghformat.load(label)
		data_label = img_label.get_data()
		data_label = np.squeeze(data_label)

		if opts.range:
			lowerLabel = int(opts.range[0])
			upperLabel = int(opts.range[1]) + 1
			num_label = upperLabel - lowerLabel
			print('Extracting from %d to %d label for a total of %d mean values' % (lowerLabel,upperLabel, num_label))
			outMeanValues = np.zeros([n,num_label])
			counter=0
			for i in range(lowerLabel,upperLabel):
				outMeanValues[:,counter]=img_data[data_label==(i)].mean(axis=0)
				counter+=1
		else:
			num_label = int(data_label.max())
			print('Extracting %d mean label values' % num_label)

print "Usage: %s [1D Subgrouping Variable] [surface (area or thickness)]" % (str(sys.argv[0]))
	print "  "
	print "Careful. Creates a subgroub based on missing data from an analysis."
	print "The subgrouping variable should be a 1D text file list with missing"
	print "variables coded as a string (e.g. NA or NaN)."
	print "--- You must be in same directory as python_temp ---"
	print "*******************************************************************"
else:
	cmdargs = str(sys.argv)
	arg_subgroupvariable = str(sys.argv[1])
	surftype = str(sys.argv[2])
	if len(sys.argv)==3:
		fwhm = str('03B')
	else:
		fwhm = str(sys.argv[3])
	img_surf_lh = nib.freesurfer.mghformat.load('lh.all.%s.%s.mgh' % (surftype,fwhm))
	img_surf_rh = nib.freesurfer.mghformat.load('rh.all.%s.%s.mgh' % (surftype,fwhm))
	subgroupvariable = np.genfromtxt(arg_subgroupvariable, delimiter=',')
	masking_variable=np.isfinite(subgroupvariable)
	data_full_lh = img_surf_lh.get_data()
	data_lh = np.squeeze(data_full_lh)
	if data_lh.shape[1] > len(subgroupvariable):
		print "Error. Number of subjects doesn't equal length of subgrouping variable"
		exit()
	affine_mask_lh = img_surf_lh.get_affine()
	outdata_mask_lh = data_full_lh[:,:,:,1]
	subdata_lh=data_lh[:,masking_variable]


	data_full_rh = img_surf_rh.get_data()
	data_rh = np.squeeze(data_full_rh)
	affine_mask_rh = img_surf_rh.get_affine()

def run(opts):
	surface = str(opts.input[0])
	num_cores = int(opts.input[1])
	surf_name = surface.split('.mgh',1)[0]
	if len(surface.split('.00.',1))==1:
		print("Please input unsmoothed surface. e.g., ?h.all.???.00.mgh")
		exit()
	surf_gen = surface.split('.00.',1)[0]
	hemi = surface.split('.',1)[0]
	if not os.path.exists('python_temp_area'):
		os.mkdir('python_temp_area')
	img_surf = nib.freesurfer.mghformat.load(surface)
	data_full = img_surf.get_data()
	data_short = np.squeeze(data_full)
	affine_mask = img_surf.get_affine()
	num_vertex = data_short.shape[0]
	low_values_indices = data_short < 0
	data_short[low_values_indices] = 0
	bc_data = boxcox(data_short, num_vertex, num_cores)
	np.save('python_temp_area/bc_data',bc_data)
	out_bc_data=data_full
	out_bc_data[:,0,0,:] = bc_data
	nib.save(nib.freesurfer.mghformat.MGHImage(out_bc_data,affine_mask),"%s.boxcox.mgh" % surf_name)
	if opts.nosmoothing:
		exit()
	else:
		print("Smoothing %s.boxcox.mgh" %  surf_name) 
		os.system = os.popen("$FREESURFER_HOME/bin/mri_surf2surf --hemi %s --s fsaverage --sval %s.boxcox.mgh --fwhm 3 --cortex --tval %s.boxcox.03B.mgh" % (hemi, surf_name, surf_gen))

def run(opts):
	arg_subgroupvariable = str(opts.input[0])
	surftype = str(opts.input[1])
	fwhm = str(opts.fwhm[0])
	img_surf_lh = nib.freesurfer.mghformat.load('lh.all.%s.%s.mgh' % (surftype,fwhm))
	img_surf_rh = nib.freesurfer.mghformat.load('rh.all.%s.%s.mgh' % (surftype,fwhm))
	subgroupvariable = np.genfromtxt(arg_subgroupvariable, delimiter=',')
	masking_variable=np.isfinite(subgroupvariable)
	data_full_lh = img_surf_lh.get_data()
	data_lh = np.squeeze(data_full_lh)
	if data_lh.shape[1] > len(subgroupvariable):
		print("Error. Number of subjects doesn't equal length of subgrouping variable")
		exit()
	affine_mask_lh = img_surf_lh.get_affine()
	subdata_lh=data_lh[:,masking_variable]

	data_full_rh = img_surf_rh.get_data()
	data_rh = np.squeeze(data_full_rh)
	affine_mask_rh = img_surf_rh.get_affine()
	subdata_rh=data_rh[:,masking_variable]

print("[%d] : %s\t%s" % (counter, roi, astr))
		else:
			print("[%d] : %s\t%s\t\tCONTAINS %d MISSING VARIABLES!" % (counter, roi, astr, num_missing))
	print("\nChecking image lengths [# subjects = %d]" % num_subjects)
	for sf in folders:
		temp_img = ("%s/lh.all.%s.%s.mgh" % (sf,surface,FWHM))
		if not os.path.isfile(temp_img):
			print ("Error: %s not found." % temp_img)
		else:
			temp_num_img = nib.freesurfer.mghformat.load(temp_img).shape[-1]
			if temp_num_img == num_subjects:
				print("%s ...OK" % temp_img)
			else:
				print("Error: Length of %s [%d] does not match number of subjects[%d]" % (temp_img, temp_num_img, num_subjects))
			temp_img = ("%s/rh.all.%s.%s.mgh" % (sf,surface,FWHM))
			temp_num_img = nib.freesurfer.mghformat.load(temp_img).shape[-1]
			if temp_num_img == num_subjects:
				print("%s ...OK" % temp_img)
			else:
				print("Error: Length of %s [%d] does not match number of subjects[%d]" % (temp_img, temp_num_img, num_subjects))

index_lh = np.array(np.genfromtxt("temp.lh.verticeslist", delimiter=',').astype(np.int))
	index_rh = np.array(np.genfromtxt("temp.rh.verticeslist", delimiter=',').astype(np.int))
	bin_mask_lh = np.zeros_like(mean_lh)
	bin_mask_lh[index_rh]=1
	bin_mask_lh[mean_lh==0]=0
	bin_mask_lh = bin_mask_lh.astype(bool)
	bin_mask_rh = np.zeros_like(mean_rh)
	bin_mask_rh[index_rh]=1
	bin_mask_lh[mean_lh==0]=0
	bin_mask_rh = bin_mask_rh.astype(bool)
elif opts.binmask:
	print("Loading masks")
	img_binmgh_lh = nib.freesurfer.mghformat.load(opts.binmask[0])
	binmgh_lh = img_binmgh_lh.get_data()
	binmgh_lh = np.squeeze(binmgh_lh)
	img_binmgh_rh = nib.freesurfer.mghformat.load(opts.binmask[0])
	binmgh_rh = img_binmgh_rh.get_data()
	binmgh_rh = np.squeeze(binmgh_rh)
	bin_mask_lh = binmgh_lh>1
	bin_mask_rh = binmgh_rh>1
else:
	bin_mask_lh = mean_lh>0
	bin_mask_rh = mean_rh>0
data_lh = data_lh[bin_mask_lh]
num_vertex_lh = data_lh.shape[0]
data_rh = data_rh[bin_mask_rh]
num_vertex_rh = data_rh.shape[0]
num_vertex = num_vertex_lh + num_vertex_rh
all_vertex = data_full_lh.shape[0]

if opts.input: 
#load variables

data_lh = np.squeeze(data_full_lh)
				affine_mask_lh = img_data_lh.affine
				outdata_mask_lh = np.zeros_like(data_full_lh[:,:,:,1])
				mask_lh = data_lh.mean(1)!=0
				img_data_rh = nib.freesurfer.mghformat.load("%s/rh.all.%s.%s.mgh" % (sfolder, surface,FWHM))
				data_full_rh = img_data_rh.get_data()
				data_rh = np.squeeze(data_full_rh)
				affine_mask_rh = img_data_rh.affine
				outdata_mask_rh = np.zeros_like(data_full_rh[:,:,:,1])
				mask_rh = data_rh.mean(1)!=0
				data.append(np.hstack((data_lh[mask_lh].T,data_rh[mask_rh].T)))
				num_vertex_lh = data_lh[mask_lh].shape[0]
				all_vertex = data_full_lh.shape[0]
			else:
				data_lh = np.squeeze(nib.freesurfer.mghformat.load("%s/lh.all.%s.%s.mgh" % (sfolder, surface,FWHM)).get_data())
				data_rh = np.squeeze(nib.freesurfer.mghformat.load("%s/rh.all.%s.%s.mgh" % (sfolder, surface,FWHM)).get_data())
				data.append(np.hstack((data_lh[mask_lh].T,data_rh[mask_rh].T)))
				data_lh = data_rh = []
		data = np.array(data, dtype=np.float32, order = 'c')
		nonzero = np.zeros_like(data)
		nonzero[:] = np.copy(data)

		#TFCE
		if opts.vertextriangularmesh:
			# 3 Neighbour vertex connectity
			print("Creating adjacency set")
			if opts.vertexsrf:
				v_lh, faces_lh = nib.freesurfer.read_geometry(opts.vertexsrf[0])
				v_rh, faces_rh = nib.freesurfer.read_geometry(opts.vertexsrf[1])
			else:
				v_lh, faces_lh = nib.freesurfer.read_geometry("%s/fsaverage/surf/lh.sphere" % os.environ["SUBJECTS_DIR"])
				v_rh, faces_rh = nib.freesurfer.read_geometry("%s/fsaverage/surf/rh.sphere" % os.environ["SUBJECTS_DIR"])

num_perm = int(sys.argv[6])
	num_cores = int(sys.argv[7])
	FWHM = '03B'

	if not os.path.exists("python_temp_med_%s" % surface):
		os.mkdir("python_temp_med_%s" % surface)

	pred_x = np.genfromtxt(arg_predictor, delimiter=",")
	covars = np.genfromtxt(arg_covars, delimiter=",")
	depend_y = np.genfromtxt(arg_depend, delimiter=",")

	np.save("python_temp_med_%s/pred_x" % surface,pred_x)
	np.save("python_temp_med_%s/covars" % surface,covars)
	np.save("python_temp_med_%s/depend_y" % surface,depend_y)

	img_lh = nib.freesurfer.mghformat.load("lh.all.%s.%s.mgh" % (surface,FWHM))
	data_full_lh = img_lh.get_data()
	data_lh = np.squeeze(data_full_lh)
	affine_mask_lh = img_lh.get_affine()
	num_vertex = data_lh.shape[0]
	num_subjects = data_lh.shape[1]
	outdata_mask_lh = data_full_lh[:,:,:,1]

	img_rh = nib.freesurfer.mghformat.load("rh.all.%s.%s.mgh" % (surface,FWHM))
	data_full_rh = img_rh.get_data()
	data_rh = np.squeeze(data_full_rh)
	affine_mask_rh = img_rh.get_affine()
	outdata_mask_rh = data_full_rh[:,:,:,1]

	np.save("python_temp_med_%s/num_subjects" % surface,num_subjects)
	np.save("python_temp_med_%s/num_vertex" % surface,num_vertex)
	np.save("python_temp_med_%s/affine_mask_lh" % surface,affine_mask_lh)

if len(a) > 10:
			astr = '[n>10]'
		else:
			astr = ','.join(a.astype(np.str))
			astr = '['+astr+']'
		if num_missing == 0:
			print("[%d] : %s\t%s" % (counter, roi, astr))
		else:
			print("[%d] : %s\t%s\t\tCONTAINS %d MISSING VARIABLES!" % (counter, roi, astr, num_missing))
	print("\nChecking image lengths [# subjects = %d]" % num_subjects)
	for sf in folders:
		temp_img = ("%s/lh.all.%s.%s.mgh" % (sf,surface,FWHM))
		if not os.path.isfile(temp_img):
			print ("Error: %s not found." % temp_img)
		else:
			temp_num_img = nib.freesurfer.mghformat.load(temp_img).shape[-1]
			if temp_num_img == num_subjects:
				print("%s ...OK" % temp_img)
			else:
				print("Error: Length of %s [%d] does not match number of subjects[%d]" % (temp_img, temp_num_img, num_subjects))
			temp_img = ("%s/rh.all.%s.%s.mgh" % (sf,surface,FWHM))
			temp_num_img = nib.freesurfer.mghformat.load(temp_img).shape[-1]
			if temp_num_img == num_subjects:
				print("%s ...OK" % temp_img)
			else:
				print("Error: Length of %s [%d] does not match number of subjects[%d]" % (temp_img, temp_num_img, num_subjects))

affine_mask_rh = img_data_rh.get_affine()
outdata_mask_rh = np.zeros_like(data_full_rh[:,:,:,1])
if not os.path.exists("lh.mean.%s.%s.mgh" % (surface,FWHM)):
	mean_lh = np.sum(data_lh,axis=1)/data_lh.shape[1]
	outmean_lh = np.zeros_like(data_full_lh[:,:,:,1])
	outmean_lh[:,0,0] = mean_lh
	nib.save(nib.freesurfer.mghformat.MGHImage(outmean_lh,affine_mask_lh),"lh.mean.%s.%s.mgh" % (surface,FWHM))
	mean_rh = np.sum(data_rh,axis=1)/data_rh.shape[1]
	outmean_rh = np.zeros_like(data_full_rh[:,:,:,1])
	outmean_rh[:,0,0] = mean_rh
	nib.save(nib.freesurfer.mghformat.MGHImage(outmean_rh,affine_mask_rh),"rh.mean.%s.%s.mgh" % (surface,FWHM))
else:
	img_mean_lh = nib.freesurfer.mghformat.load("lh.mean.%s.%s.mgh" % (surface,FWHM))
	mean_full_lh = img_mean_lh.get_data()
	mean_lh = np.squeeze(mean_full_lh)
	img_mean_rh = nib.freesurfer.mghformat.load("rh.mean.%s.%s.mgh" % (surface,FWHM))
	mean_full_rh = img_mean_rh.get_data()
	mean_rh = np.squeeze(mean_full_rh)

#TFCE
if opts.triangularmesh:
	print "Creating adjacency set"
	# 3 Neighbour vertex connectity
	v_lh, faces_lh = nib.freesurfer.read_geometry("%s/fsaverage/surf/lh.sphere" % os.environ["SUBJECTS_DIR"])
	v_rh, faces_rh = nib.freesurfer.read_geometry("%s/fsaverage/surf/rh.sphere" % os.environ["SUBJECTS_DIR"])
	adjac_lh = create_adjac(v_lh,faces_lh)
	adjac_rh = create_adjac(v_rh,faces_rh)
elif opts.adjfiles:
	print "Loading prior adjacency set"
	arg_adjac_lh = opts.adjfiles[0]
	arg_adjac_rh = opts.adjfiles[1]
	adjac_lh = np.load(arg_adjac_lh)