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

def create_sample_data(self):
        """
        Create sample data files. The files created are:

        * analyze_image.hdr + .img: plain ANALYZE image
        * mgh_image.mgh: MGH image
        * nifti1_image.nii: NIfTI1 image
        * nifti2_image.nii:  NIfTI2 image
        * spm2_image.hdr + .img + .mat: SPM2 ANALYZE image
        * spm99_image.hdr + .img + .mat: SPM99 ANALYZE image
        """
        file_name = os.path.join(self.data_dir, "analyze_image")
        analyze_img = nib.AnalyzeImage(self.get_data(), np.eye(4))
        analyze_img.to_filename(file_name)
        file_name = os.path.join(self.data_dir, "mgh_image")
        mgh_img = nib.freesurfer.mghformat.MGHImage(self.get_data(),
                                                    np.eye(4))
        mgh_img.to_filename(file_name)
        file_name = os.path.join(self.data_dir, "nifti1_image")
        nifti1_img = nib.Nifti1Image(self.get_data(), self.get_affine())
        nifti1_img.to_filename(file_name)
        file_name = os.path.join(self.data_dir, "nifti2_image")
        nifti2_img = nib.Nifti2Image(self.get_data(), self.get_affine())
        nifti2_img.to_filename(file_name)
        file_name = os.path.join(self.data_dir, "spm2_image")
        spm2_img = nib.spm2analyze.Spm2AnalyzeImage(self.get_data(),
                                                    np.eye(4))
        spm2_img.to_filename(file_name)
        file_name = os.path.join(self.data_dir, "spm99_image")
        spm99_img = nib.spm99analyze.Spm99AnalyzeImage(self.get_data(),
                                                       np.eye(4))
        spm99_img.to_filename(file_name)

def MRISsampleVertexPatch(mris, mri_normals, mri, vno, wsize):
    if isinstance(mri_normals, nib.freesurfer.mghformat.MGHImage) == True:
        mri_normals = mri_normals.get_data().squeeze()
    if (len(mris) == 2):
        mris = mris[0]
    nz = mri_normals[vno,:]
    nz /= np.linalg.norm(nz)
    if (nz[2] < 0.9):
        nx = np.cross(nz, [0,0,1])
    else:
        nx = np.cross(nz, [0,1,0])
    nx /= np.linalg.norm(nx)
    ny = np.cross(nx, nz)
    ny /= np.linalg.norm(ny)
    return sample_patch(mri, mris[vno]+(nz), nx, ny, nz, wsize)

def run(opts):
	arg_surf = str(opts.input[0])
	thresh = float(opts.threshold[0])
	img_surf = nib.freesurfer.mghformat.load(arg_surf)
	data_full = img_surf.get_data()
	mask = data_full>thresh
	mask = np.array(mask*1).astype('>f4')
	nib.save(nib.freesurfer.mghformat.MGHImage(mask,img_surf.get_affine()),'bin_mask_%1.1f_%s' % (thresh,arg_surf))
	if opts.neg:
		data_full = data_full*-1
		mask = data_full>thresh
		mask = np.array(mask*1).astype('>f4')
		nib.save(nib.freesurfer.mghformat.MGHImage(mask,img_surf.get_affine()),'bin_mask_%1.1f_neg_%s' % (thresh,arg_surf))

def write_sobelZ(SobelZ,outdata_mask, affine, surface, hemi):
	outdata_mask[:,0,0] = SobelZ
	nib.save(nib.freesurfer.mghformat.MGHImage(outdata_mask,affine),"SobelZ_%s_%s_%s.mgh" % (medtype,surface,hemi))

img_data_lh = nib.freesurfer.mghformat.load("lh.all.%s.%s.mgh" % (surface,FWHM))
data_full_lh = img_data_lh.get_data()
data_lh = np.squeeze(data_full_lh)
affine_mask_lh = img_data_lh.get_affine()
n = data_lh.shape[1]
outdata_mask_lh = np.zeros_like(data_full_lh[:,:,:,1])
img_data_rh = nib.freesurfer.mghformat.load("rh.all.%s.%s.mgh" % (surface,FWHM))
data_full_rh = img_data_rh.get_data()
data_rh = np.squeeze(data_full_rh)
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)

#create masks
if opts.fmri:
	maskthresh = opts.fmri
	bin_mask_lh = np.logical_or(mean_lh > maskthresh, mean_lh < (-1*maskthresh))

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]
	num_subjects = data_short.shape[1]
	outdata_mask = data_full[:,:,:,1]
	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)
	print "Smoothing %s.boxcox.mgh" %  surf_name 
	os.system = os.popen("mri_surf2surf --hemi %s --s fsaverage --sval %s.boxcox.mgh --fwhm 3 --cortex --tval %s.03B.boxcox.mgh" % (hemi, surf_name, surf_gen))

def savemgh_v2(image_array, index, imagename, affine=None):
	if not imagename.endswith('mgh'):
		imagename += '.mgh'
	outdata = image_array.astype(np.float32, order = "C")
	if image_array.ndim == 1:
		imgout = np.zeros((index.shape[0],index.shape[1],index.shape[2]))
		imgout[index]=outdata
	elif image_array.shape[1] > 1:
		imgout = np.zeros((index.shape[0],index.shape[1],index.shape[2],image_array.shape[1]))
		imgout[index]=outdata
	else:
		imgout = np.zeros((index.shape[0],index.shape[1],index.shape[2]))
		imgout[index]=outdata[:,0]
	nib.save(nib.freesurfer.mghformat.MGHImage(imgout.astype(np.float32, order = "C"),affine=affine),imagename)

""" Unconditionally set affine into the header """
        hdr = self._header
        shape = np.array(self._dataobj.shape[:3])

        # for more information, go through save_mgh.m in FreeSurfer dist
        voxelsize = voxel_sizes(self._affine)
        Mdc = self._affine[:3, :3] / voxelsize
        c_ras = self._affine.dot(np.hstack((shape / 2.0, [1])))[:3]

        # Assign after we've had a chance to raise exceptions
        hdr['delta'] = voxelsize
        hdr['Mdc'] = Mdc.T
        hdr['Pxyz_c'] = c_ras


load = MGHImage.load
save = MGHImage.instance_to_filename

Transformation matrix

    copy_header: boolean, optional
        Indicated if the header of the reference image should be used to
        create the new image

    Returns
    -------

    new_img: image
        An image which has the same type as the reference image.
    """
    if data.dtype == bool:
        default_dtype = np.int8
        if (LooseVersion(nibabel.__version__) >= LooseVersion('1.2.0') and
                isinstance(ref_img, nibabel.freesurfer.mghformat.MGHImage)):
            default_dtype = np.uint8
        data = as_ndarray(data, dtype=default_dtype)
    header = None
    if copy_header:
        header = copy.copy(ref_img.get_header())
        header['scl_slope'] = 0.
        header['scl_inter'] = 0.
        header['glmax'] = 0.
        header['cal_max'] = np.max(data) if data.size > 0 else 0.
        header['cal_max'] = np.min(data) if data.size > 0 else 0.
    return ref_img.__class__(data, affine, header=header)

img_data_lh = nib.freesurfer.mghformat.load("lh.all.%s.%s.mgh" % (surface,FWHM))
	data_full_lh = img_data_lh.get_data()
	data_lh = np.squeeze(data_full_lh)
	affine_mask_lh = img_data_lh.get_affine()
	n = data_lh.shape[1] # num_subjects
	outdata_mask_lh = np.zeros_like(data_full_lh[:,:,:,1])
	img_data_rh = nib.freesurfer.mghformat.load("rh.all.%s.%s.mgh" % (surface,FWHM))
	data_full_rh = img_data_rh.get_data()
	data_rh = np.squeeze(data_full_rh)
	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")
		if opts.inputsurfs: