""" Created on 10/08/2017 @author: Niklas Pallast Neuroimaging & Neuroengineering Department of Neurology University Hospital Cologne """ from __future__ import print_function import os import sys import numpy as np import nibabel as nib import nibabel.nifti1 as nii import pv_parseBruker_md_np as pB import P2_IDLt2_mapping as mapT2 class Bruker2Nifti: def __init__(self, study, expno, procno, rawfolder, procfolder, ftype='NIFTI_GZ'): self.study = study self.expno = str(expno) self.procno = str(procno) self.rawfolder = rawfolder self.procfolder = procfolder self.ftype = ftype def read_2dseq(self, map_raw=False, pv6=False, sc=1.0): study = self.study expno = self.expno procno = self.procno rawfolder = self.rawfolder self.acqp = pB.parsePV(os.path.join(rawfolder, study, expno, 'acqp')) self.method = pB.parsePV(os.path.join(rawfolder, study, expno, 'method')) self.subject = pB.parsePV(os.path.join(rawfolder, study, 'subject')) # get header information datadir = os.path.join(rawfolder, study, expno, 'pdata', procno) #self.d3proc = pB.parsePV(os.path.join(datadir, 'd3proc')) # removed for PV6 self.visu_pars = pB.parsePV(os.path.join(datadir, 'visu_pars')) hdr = pB.getNiftiHeader(self.visu_pars, sc=sc) #print("hdr:", hdr) if hdr is None or not isinstance(hdr[12], str): return # read '2dseq' file f_id = open(os.path.join(datadir, '2dseq'), 'rb') data = np.fromfile(f_id, dtype=np.dtype(hdr[12])).reshape(hdr[1], hdr[2], hdr[3], hdr[4], order='F') f_id.close() # map to raw data range (PV6) if map_raw: visu_core_data_slope = np.array(map(float, self.visu_pars['VisuCoreDataSlope'].split()), dtype=np.float32) visu_core_data_offs = np.array(map(float, self.visu_pars['VisuCoreDataOffs'].split()), dtype=np.float32) visu_core_data_shape = list(data.shape) visu_core_data_shape[:2] = (1, 1) if pv6: data = data / visu_core_data_slope.reshape(visu_core_data_shape) else: data = data * visu_core_data_slope.reshape(visu_core_data_shape) data = data + visu_core_data_offs.reshape(visu_core_data_shape) # NIfTI image #data = np.flip(data, axis=(1,2)) #flip z-axis nim = nii.Nifti1Image(data, None) # NIfTI header #header = nim.header header = nim._header #print("header:"); print(header) header['pixdim'] = [0.0, hdr[5], hdr[6], hdr[7], hdr[8], 0.0, 0.0, 0.0] #nim.setXYZUnit('mm') header.set_xyzt_units(xyz='mm', t=None) #nim.header = header #header = nim.get_header() #print("header:"); print(header) # write header in xml structure #xml = pB.getXML(datadir + "/") xml = pB.getXML(os.path.join(datadir, 'visu_pars')) #print("xml:"); print(xml) # add protocol information (method, acqp, visu_pars, d3proc) to Nifti's header extensions #nim.extensions += ('comment', xml) #extension = nii.Nifti1Extension('comment', xml) self.hdr = hdr self.nim = nim self.xml = xml def save_nifti(self, subfolder=''): procfolder = os.path.join(self.procfolder, self.study) if not os.path.isdir(procfolder): os.mkdir(procfolder) if "Localizer" in self.acqp['ACQ_protocol_name']: procfolder = os.path.join(self.procfolder, self.study, "Localizer") elif "DTI" in self.acqp['ACQ_protocol_name'] or "Diffusion" in self.acqp['ACQ_protocol_name']: procfolder = os.path.join(self.procfolder, self.study, "DTI") elif "fMRI" in self.acqp['ACQ_protocol_name']: procfolder = os.path.join(self.procfolder, self.study, "fMRI") elif "Turbo" in self.acqp['ACQ_protocol_name']: procfolder = os.path.join(self.procfolder, self.study, "T2w") elif "MSME" in self.acqp['ACQ_protocol_name']: procfolder = os.path.join(self.procfolder, self.study, "T2map") else: procfolder = os.path.join(self.procfolder, self.study, "Others") if not os.path.isdir(procfolder): os.mkdir(procfolder) if self.ftype == 'NIFTI_GZ': ext = 'nii.gz' elif self.ftype == 'NIFTI': ext = 'nii' elif self.ftype == 'ANALYZE': ext = 'img' else: ext = 'nii.gz' fname = '.'.join([self.study, self.expno, self.procno, ext]) # write Nifti file print(os.path.join(procfolder, fname)) if not hasattr(self, 'nim'): return nib.save(self.nim, os.path.join(procfolder, fname)) return os.path.join(procfolder, fname) def save_table(self, subfolder= ''): procfolder = os.path.join(self.procfolder, self.study) if not os.path.isdir(procfolder): os.mkdir(procfolder) procfolder = os.path.join(self.procfolder, self.study, subfolder) if not os.path.isdir(procfolder): os.mkdir(procfolder) #dw_bval_each = float(self.method['PVM_DwBvalEach']) if 'PVM_DwEffBval' in self.method: dw_eff_bval = np.array(list(map(float, self.method['PVM_DwEffBval'].split())), dtype=np.float32) #print("dw_bval_each:", dw_bval_each) #print("dw_eff_bval:"); print(dw_eff_bval) if 'PVM_DwAoImages' in self.method: dw_ao_images = int(self.method['PVM_DwAoImages']) if 'PVM_DwNDiffDir' in self.method: dw_n_diff_dir = int(self.method['PVM_DwNDiffDir']) #print("dw_ao_images:", dw_ao_images) #print("dw_n_diff_dir:", dw_n_diff_dir) if 'PVM_DwDir' in self.method: dw_dir = np.array(list(map(float, self.method['PVM_DwDir'].split())), dtype=np.float32) dw_dir = dw_dir.reshape((dw_n_diff_dir, 3)) nd = dw_ao_images + dw_n_diff_dir bvals = np.zeros(nd, dtype=np.float32) dwdir = np.zeros((nd, 3), dtype=np.float32) bvals[dw_ao_images:] = dw_eff_bval[dw_ao_images:] dwdir[dw_ao_images:] = dw_dir fname = '.'.join([self.study, self.expno, self.procno, 'btable', 'txt']) print(os.path.join(procfolder, fname)) # Open btable file to write binary (windows format) #fid = open(os.path.join(procfolder, fname), 'wb') - py 2.6 fid = open(os.path.join(procfolder, fname),mode='w',buffering=-1) for i in range(nd): fid.write("%.4f" % (bvals[i],) + " %.8f %.8f %.8f" % tuple(dwdir[i])) #print("%.4f" % (bvals[i],) + " %.8f %.8f %.8f" % tuple(dwdir[i]), end="\r\n", file=fid) - py 2.6 # Close file fid.close() fname = '.'.join([self.study, self.expno, self.procno, 'bvals', 'txt']) print(os.path.join(procfolder, fname)) # Open bvals file to write binary (unix format) fid = open(os.path.join(procfolder, fname), mode='w', buffering=-1) #fid = open(os.path.join(procfolder, fname), 'wb') - py 2.6 fid.write(" ".join("%.4f" % (bvals[i],) for i in range(nd))) #print(" ".join("%.4f" % (bvals[i],) for i in range(nd)), end=chr(10), file=fid) - py 2.6 # Close bvals file fid.close() fname = '.'.join([self.study, self.expno, self.procno, 'bvecs', 'txt']) print(os.path.join(procfolder, fname)) # Open bvecs file to write binary (unix format) fid = open(os.path.join(procfolder, fname), mode='w', buffering=-1) #fid = open(os.path.join(procfolder, fname), 'wb') - py 2.6 for k in range(3): fid.write(" ".join("%.8f" % (dwdir[i,k],) for i in range(nd))) #print(" ".join("%.8f" % (dwdir[i,k],) for i in range(nd)), end=chr(10), file=fid)- py 2.6 # Close bvecs file fid.close() if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description='Convert ParaVision to NIfTI') requiredNamed = parser.add_argument_group('Required named arguments') requiredNamed.add_argument('-i','--input_folder', help='raw data folder') # parser.add_argument('-o','--output_folder', help='output data folder') # parser.add_argument('study', help='study name') # parser.add_argument('expno', help='experiment number') # parser.add_argument('procno', help='processed (reconstructed) images number') parser.add_argument('-f','--model', help='T2_2p (default) : Two parameter T2 decay S(t) = S0 * exp(-t/T2)\n' 'T2_3p : Three parameter T2 decay S(t) = S0 * exp(-t/T2) + C' , nargs='?', const='T2_2p', type=str, default='T2_2p') parser.add_argument('-u','--upLim', help='upper limit of TE - default: 100', nargs='?', const=100, type=int, default=100) parser.add_argument('-s','--snrLim', help='upper limit of SNR - default: 1.5', nargs='?', const=1.5, type=float, default=1.5) parser.add_argument('-k','--snrMethod', help='Brummer ,Chang, Sijbers', nargs='?', const='Brummer', type=str, default='Brummer') parser.add_argument('-m', '--map_raw', action='store_true', help='get the real values') parser.add_argument('-p', '--pv6', action='store_true', help='ParaVision 6') parser.add_argument('-t', '--table', action='store_true', help='save b-values and diffusion directions') args = parser.parse_args() input_folder = None # raw data folder if args.input_folder is not None: input_folder = args.input_folder if not os.path.isdir(input_folder): sys.exit("Error: '%s' is not an existing directory." % (input_folder,)) list = os.listdir(input_folder) listOfScans = [s for s in list if s.isdigit()] if len(listOfScans) == 0: sys.exit("Error: '%s' contains no numbered scans." % (input_folder,)) print('Start to process '+str(len(listOfScans))+' scans...') procno ='1' study=input_folder.split('/')[len(input_folder.split('/'))-1] print(study) img = [] for expno in np.sort(listOfScans): path = os.path.join(input_folder, expno, 'pdata', procno) if not os.path.isdir(path): sys.exit("Error: '%s' is not an existing directory." % (path,)) if os.path.exists(os.path.join(path,'2dseq')): img = Bruker2Nifti(study, expno, procno, os.path.split(input_folder)[0], input_folder, ftype='NIFTI_GZ') img.read_2dseq(map_raw=args.map_raw, pv6=args.pv6) resPath = img.save_nifti() if resPath is None: continue if 'VisuAcqEchoTime' in img.visu_pars: echoTime = img.visu_pars['VisuAcqEchoTime'] echoTime = np.fromstring(echoTime, dtype=float, sep=' ') if len(echoTime) > 3: mapT2.getT2mapping(resPath,args.model,args.upLim,args.snrLim,args.snrMethod,echoTime) if resPath is not None: pathlog = os.path.dirname(os.path.dirname(resPath)) pathlog = os.path.join(pathlog, 'data.log') logfile = open(pathlog, 'w') logfile.write(img.subject['coilname']) logfile.close()