diff --git a/brkraw/lib/recoFunctions.py b/brkraw/lib/recoFunctions.py index c40e09b..d4e734c 100644 --- a/brkraw/lib/recoFunctions.py +++ b/brkraw/lib/recoFunctions.py @@ -68,7 +68,6 @@ def reco_phase_rotate(frame, Reco, actual_framenumber): # import variables RECO_rotate_all = get_value(Reco,'RECO_rotate') - if RECO_rotate_all.shape[1] > actual_framenumber: RECO_rotate = get_value(Reco,'RECO_rotate')[:, actual_framenumber] @@ -89,13 +88,12 @@ def reco_phase_rotate(frame, Reco, actual_framenumber): phase_matrix = np.ones_like(frame) for index in range(len(RECO_rotate)): f = np.arange(dims[index]) - #print(f) + if RECO_ft_mode in ['COMPLEX_FT', 'COMPLEX_FFT']: phase_vector = np.exp(1j*2*np.pi*RECO_rotate[index]*f) elif RECO_ft_mode in ['NO_FT', 'NO_FFT']: phase_vector = np.ones_like(f) elif RECO_ft_mode in ['COMPLEX_IFT', 'COMPLEX_IFFT']: - #print(RECO_rotate[index]) phase_vector = np.exp(1j*2*np.pi*(1-RECO_rotate[index])*f) else: raise ValueError('Your RECO_ft_mode is not supported') @@ -192,7 +190,7 @@ def reco_FT(frame, Reco, actual_framenumber): elif RECO_ft_mode in ['NO_FT', 'NO_FFT']: pass elif RECO_ft_mode in ['COMPLEX_IFT', 'COMPLEX_IFFT']: - frame = np.fft.fftn(frame) + frame = np.fft.ifftn(frame) #frame = sp.ifft(frame, axes=[0,1,2], center=False) else: raise ValueError('Your RECO_ft_mode is not supported') @@ -236,10 +234,10 @@ def reco_cutoff(frame, Reco, actual_framenumber): # Use function only if Reco.RECO_size is not equal to size(frame) dim_equal = True for i,j in zip(get_value(Reco,'RECO_size'), frame.shape): - dim_equal = (i==j) - - if not dim_equal: - + if i!=j: + dim_equal = False + + if not dim_equal: # Import variables RECO_offset = get_value(Reco,'RECO_offset')[:, actual_framenumber] RECO_size = get_value(Reco, 'RECO_size') @@ -247,9 +245,9 @@ def reco_cutoff(frame, Reco, actual_framenumber): # Cut the new part with RECO_size and RECO_offset pos_ges = [] for i in range(len(RECO_size)): - pos_ges.append(slice(RECO_offset[i], RECO_offset[i] + RECO_size[i])) newframe = frame[tuple(pos_ges)] + else: newframe = frame diff --git a/brkraw/lib/recoSigpy.py b/brkraw/lib/recoSigpy.py new file mode 100644 index 0000000..ac7fada --- /dev/null +++ b/brkraw/lib/recoSigpy.py @@ -0,0 +1,52 @@ +""" +Created on Sat Feb 08 2024 + DEVELOPED FOR BRUKER PARAVISION 360 datasets + Below code will work for 3D cartesian sequence + undersampled in the y-z phase view + Reconstructions are completed with SIGPY + + THIS IS A WORK IN PROGRESS due to a lack of CS + datasets from PV360.3.3 or higher + +@author: Tim Ho (UVA) +""" + +import sigpy as sp +import sigpy.mri as mr +import numpy as np + +from .utils import get_value + +def compressed_sensing_recon(data, acqp, meth, reco, lamda=0.01, method=None): + # Meta Variables + ky = get_value(meth,'PVM_EncGenSteps1') + kz = get_value(meth,'PVM_EncGenSteps2') + PVM_EncGenTotalSteps = get_value(meth, 'PVM_EncGenTotalSteps') + + kspaceShape = [1 for _ in range(7)] + RECO_ft_size = get_value(reco, 'RECO_ft_size') + NI = get_value(acqp, 'NI') + NR = get_value(acqp, 'NR') + kspaceShape[1:len(RECO_ft_size)+1] = RECO_ft_size + kspaceShape[0] = data.shape[1] + kspaceShape[5] = NI + kspaceShape[6] = NR + + # Resort Raw Sorted to K-SPACE + frame_sort = data.reshape((NR, PVM_EncGenTotalSteps, NI, kspaceShape[0], RECO_ft_size[0])) + + k_space = np.zeros(shape=kspaceShape, dtype=complex) + for index, (i,j) in enumerate(zip(ky,kz)): + k_space[:,:,int(i+max(ky)), j+max(kz)+1,0,:,:]= frame_sort[:,index,:,:,:].transpose(2,3,1,0) + + + N1, N2, N3, N4, N5, N6, N7 = k_space.shape + output = np.zeros(shape=(N2,N3,N4,N6,N7), dtype=complex) + for NR in range(N7): + for NI in range(N6): + k_space = np.squeeze(k_space[:,:,:,:,:,NI,NR]) + # Compressed Sensing + sens = mr.app.EspiritCalib(k_space).run() + output[:,:,:,NI,NR] = mr.app.L1WaveletRecon(k_space, sens, lamda=lamda).run() + + return output \ No newline at end of file diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py index c3860fa..62580e1 100644 --- a/brkraw/lib/recon.py +++ b/brkraw/lib/recon.py @@ -15,7 +15,7 @@ from copy import deepcopy -def recon(fid_binary, acqp, meth, reco, process = None, recoparts = 'default'): +def recon(fid_binary, acqp, meth, reco, process = 'image', recoparts = 'default'): """ Process FID -> Channel Sorted -> Frame-Sorted -> K-Sorted -> Image Parameters @@ -44,16 +44,9 @@ def recon(fid_binary, acqp, meth, reco, process = None, recoparts = 'default'): """ output = readBrukerRaw(fid_binary, acqp, meth) - if process == 'raw': return output - if get_value(meth,'PVM_EncCS') == 'Yes': - print(get_value(acqp, 'ACQ_scan_name' )) - print("Warning: Compressed Sensing NOT SUPPORTED...") - print("returning 'raw' sorting") - return output - # IF CORRECT SEQUENCES if 'rare' in get_value(acqp, 'ACQ_protocol_name').lower() or \ 'msme' in get_value(acqp, 'ACQ_protocol_name').lower() or \ @@ -62,6 +55,20 @@ def recon(fid_binary, acqp, meth, reco, process = None, recoparts = 'default'): 'mge' in get_value(acqp, 'ACQ_protocol_name').lower() or \ 'FLASH.ppg' == get_value(acqp, 'PULPROG'): + # Compressed Sensing + if get_value(meth,'PVM_EncCS') == 'Yes': + try: + import sigpy as sp + from . import recoSigpy + except ImportError: + raise ImportError('Sigpy Module Not Installed') + + print(get_value(acqp, 'ACQ_scan_name' )) + print("Warning: Compressed Sensing is not fully supported ...") + output = recoSigpy.compressed_sensing_recon(output, acqp, meth, reco) + return output + + # Full Cartesian Pipeline output = convertRawToFrame(output, acqp, meth) if process == 'frame': return output @@ -93,28 +100,80 @@ def readBrukerRaw(fid_binary, acqp, meth): Returns ------- - X : np.array + X : np.array [num_lines, channel, scan_size] """ - dt_code = np.dtype('float64') if get_value(acqp, 'ACQ_ScanPipeJobSettings')[0][1] == 'STORE_64bit_float' else np.dtype('int32') - fid = np.frombuffer(fid_binary, dt_code) - + from .reference import BYTEORDER, WORDTYPE + # META DATA NI = get_value(acqp, 'NI') NR = get_value(acqp, 'NR') - - ACQ_size = get_value(acqp, 'ACQ_jobs')[0][0] - numDataHighDim = np.prod(ACQ_size) - numSelectedRecievers = get_value(acqp, 'ACQ_ReceiverSelectPerChan').count('Yes') - nRecs = numSelectedRecievers - - jobScanSize = get_value(acqp, 'ACQ_jobs')[0][0] - dim1 = int(len(fid)/(jobScanSize*nRecs)) - - # Assume data is complex - X = fid[::2] + 1j*fid[1::2] - - # [num_lines, channel, scan_size] - X = np.reshape(X, [dim1, nRecs, int(jobScanSize/2)]) - + + dt_code = 'int32' + if get_value(acqp, 'ACQ_ScanPipeJobSettings') != None: + if get_value(acqp, 'ACQ_ScanPipeJobSettings')[0][1] == 'STORE_64bit_float': + dt_code = 'float64' + + if '32' in dt_code: + bits = 32 # Need to add a condition here + elif '64' in dt_code: + bits = 64 + DT_CODE = np.dtype(dt_code) + + BYTORDA = get_value(acqp, 'BYTORDA') + if BYTORDA == 'little': + DT_CODE = DT_CODE.newbyteorder('<') + elif BYTORDA == 'big': + DT_CODE = DT_CODE.newbyteorder('>') + + # Get FID FROM buffer + fid = np.frombuffer(fid_binary, DT_CODE) + + # Sort raw data + if '360' in get_value(acqp,'ACQ_sw_version'): + # METAdata for 360 + ACQ_size = get_value(acqp, 'ACQ_jobs')[0][0] + numDataHighDim = np.prod(ACQ_size) + nRecs = get_value(acqp, 'ACQ_ReceiverSelectPerChan').count('Yes') + + jobScanSize = get_value(acqp, 'ACQ_jobs')[0][0] + + # Assume data is complex + X = fid[::2] + 1j*fid[1::2] + + # [num_lines, channel, scan_size] + X = np.reshape(X, [-1, nRecs, int(jobScanSize/2)]) + + else: + # METAdata Versions Before 360 + nRecs = 1 # PV7 only save 1 channel + if get_value(acqp, 'ACQ_ReceiverSelect') != None: + nRecs = get_value(acqp, 'ACQ_ReceiverSelect').count('Yes') + + ACQ_size = get_value(acqp, 'ACQ_size' ) + if type(ACQ_size) == int: + ACQ_size = [ACQ_size] + + if get_value(acqp, 'GO_block_size') == 'Standard_KBlock_Format': + blocksize = int(np.ceil(ACQ_size[0]*nRecs*(bits/8)/1024)*1024/(bits/8)) + else: + blocksize = int(ACQ_size[0]*nRecs) + + # CHECK SIZE + if fid.size != blocksize*np.prod(ACQ_size[1:])*NI*NR: + raise Exception('readBrukerRaw 158: Error Size dont match') + + # Convert to Complex + fid = fid[::2] + 1j*fid[1::2] + fid = fid.reshape([-1,blocksize//2]) + + # Reshape Matrix [num_lines, channel, scan_size] + if blocksize != ACQ_size[0]*nRecs: + fid = fid[:,:ACQ_size[0]//2*nRecs] + fid = fid.reshape((-1,nRecs,ACQ_size[0]//2)) + X = fid.transpose(0,1,2) + + else: + X = fid.reshape((-1, nRecs, ACQ_size[0]//2)) + return X @@ -137,7 +196,7 @@ def convertRawToFrame(data, acqp, meth): results = data.copy() - # Turn Raw into frame + # Metadata NI = get_value(acqp, 'NI') NR = get_value(acqp, 'NR') @@ -155,6 +214,9 @@ def convertRawToFrame(data, acqp, meth): acqSizes = np.zeros(ACQ_dim) scanSize = get_value(acqp, 'ACQ_jobs')[0][0] + + if scanSize == 0: + scanSize = get_value(acqp,'ACQ_size')[0] isSpatialDim = [i == 'Spatial' for i in get_value(acqp, 'ACQ_dim_desc')] spatialDims = sum(isSpatialDim) @@ -167,15 +229,15 @@ def convertRawToFrame(data, acqp, meth): numresultsHighDim=np.prod(acqSizes[1:]) acqSizes[0] = scanSize - + if np.iscomplexobj(results): scanSize = int(acqSizes[0]/2) else: scanSize = acqSizes[0] - # Resort + # Start Resort based on if ACQ_dim>1: - # [..., num_lines, channel, scan_size] + # [num_readout, channel, scan_size] -> [channel, scan_size, num_readout] results = results.transpose((1,2,0)) results = results.reshape( @@ -200,16 +262,17 @@ def convertRawToFrame(data, acqp, meth): frame[:,:,:,ACQ_obj_order,:] = results else: + # ACQ = 1, Method is a Spectroscopy # Havent encountered this situation yet # Leaving code in just in case ''' results = np.reshape(results,(numSelectedRecievers, scanSize,1,NI,NR), order='F') return_out = np.zeros_like(results) return_out = np.transpose(results, (1, 2, 0, 3, 4)) - frame = return_out''' - raise 'Bug here 120' + frame = return_out + ''' + raise Exception('Bug here 120') - return frame @@ -228,9 +291,10 @@ def convertFrameToCKData(frame, acqp, meth): ------- data : np.array with size [scansize, ACQ_phase_factor, numDataHighDim/ACQ_phase_factor, numSelectedReceivers, NI, NR] - for simplified understand of the data structure + for simplified understanding of the data structure [x,y,z,_,n_channel,NI,NR] """ + # Metadata NI = get_value(acqp, 'NI') NR = get_value(acqp, 'NR') @@ -239,9 +303,11 @@ def convertFrameToCKData(frame, acqp, meth): ACQ_dim = get_value(acqp, 'ACQ_dim') numSelectedReceivers= frame.shape[2] - acqSizes = np.zeros(ACQ_dim) - + acqSizes = np.zeros(ACQ_dim) scanSize = get_value(acqp, 'ACQ_jobs')[0][0] + if scanSize == 0: + scanSize = get_value(acqp,'ACQ_size')[0] + acqSizes[0] = scanSize ACQ_size = acqSizes @@ -279,16 +345,13 @@ def convertFrameToCKData(frame, acqp, meth): # No zero-filling/interpolation available. PVM_EncZf = np.ones((ACQ_dim)) - # Resort - + # Start Resort Process frameData = frame.copy() # MGE with alternating k-space readout: Reverse every second scan. - if get_value(meth, 'EchoAcqMode') != None and get_value(meth,'EchoAcqMode') == 'allEchoes': frameData[:,:,:,1::2,:] = np.flipud(frame[:,:,:,1::2,:]) - # Calculate size of Cartesian k-space # Step 1: Anti-Aliasing ckSize = np.round(np.array(PVM_AntiAlias)*np.array(PVM_Matrix)) @@ -297,14 +360,11 @@ def convertFrameToCKData(frame, acqp, meth): reduceZf = 2*np.floor( (ckSize - ckSize/np.array(PVM_EncZf))/2 ) ckSize = ckSize - reduceZf - # # index of central k-space point (+1 for 1-based indexing in MATLAB) + # index of central k-space point (+1 for 1-based indexing in MATLAB) ckCenterIndex = np.floor(ckSize/2 + 0.25) + 1 - readStartIndex = int(ckSize[0]-scanSize + 1) - - # Reshape & store - # switch ACQ_dim + # Reshape & store based on dimension if ACQ_dim == 1: frameData = np.reshape(frameData,(scanSize, 1, 1, 1, numSelectedReceivers, NI, NR) , order='F') data = np.zeros((ckSize[0], 1, 1, 1, numSelectedReceivers, NI, NR), dtype=complex) @@ -339,21 +399,20 @@ def brkraw_Reco(kdata, reco, meth, recoparts = 'all'): elif recoparts == 'default': recoparts = ['quadrature', 'phase_rotate', 'zero_filling', 'FT', 'phase_corr_pi'] - # Other stuff - RECO_ft_mode = get_value(reco, 'RECO_ft_mode') - if '360' in meth.headers['title'.upper()]: - reco_ft_mode_new = [] - - for i in RECO_ft_mode: - if i == 'COMPLEX_FT' or i == 'COMPLEX_FFT': - reco_ft_mode_new.append('COMPLEX_IFT') - else: - reco_ft_mode_new.append('COMPLEX_FT') - - reco = set_value(reco, 'RECO_ft_mode', reco_ft_mode_new) - RECO_ft_mode = get_value(reco, 'RECO_ft_mode') - + # Metadata + RECO_ft_mode = get_value(reco, 'RECO_ft_mode') + # Adapt FT convention to acquisition version. + reco_ft_mode_new = [] + for i in RECO_ft_mode: + if i == 'COMPLEX_FT' or i == 'COMPLEX_FFT': + reco_ft_mode_new.append('COMPLEX_IFT') + else: + reco_ft_mode_new.append('COMPLEX_FT') + reco = set_value(reco, 'RECO_ft_mode', reco_ft_mode_new) + RECO_ft_mode = get_value(reco, 'RECO_ft_mode') + + # DIMS N1, N2, N3, N4, N5, N6, N7 = kdata.shape dims = kdata.shape[0:4] @@ -377,6 +436,7 @@ def brkraw_Reco(kdata, reco, meth, recoparts = 'all'): # --- START RECONSTRUCTION --- for recopart in recoparts: + if 'quadrature' in recopart: for NR in range(N7): for NI in range(N6): @@ -421,20 +481,17 @@ def brkraw_Reco(kdata, reco, meth, recoparts = 'all'): for chan in range(N5): reco_result[:,:,:,:,chan,NI,NR] = reco_phase_corr_pi(reco_result[:,:,:,:,chan,NI,NR], reco, map_index[(NI+1)*(NR+1)-1]) - ''' # There is a current bug with cutoff function if 'cutoff' in recopart: newdata_dims=[1, 1, 1, 1] reco_size = get_value(reco, 'RECO_size') newdata_dims[0:len(reco_size)] = reco_size newdata = np.zeros(shape=newdata_dims+[N5, N6, N7], dtype=np.complex128) - for NR in range(N7): for NI in range(N6): for chan in range(N5): newdata[:,:,:,:,chan,NI,NR] = reco_cutoff(reco_result[:,:,:,:,chan,NI,NR], reco, map_index[(NI+1)*(NR+1)-1]) reco_result=newdata - ''' if 'scale_phase_channels' in recopart: for NR in range(N7): @@ -470,6 +527,5 @@ def brkraw_Reco(kdata, reco, meth, recoparts = 'all'): for chan in range(N5): reco_result[:,:,:,:,chan,NI,NR] = reco_transposition(reco_result[:,:,:,:,chan,NI,NR], reco, map_index[(NI+1)*(NR+1)-1]) # --- End of RECONSTRUCTION Loop --- - - + return reco_result \ No newline at end of file diff --git a/brkraw/scripts/brkrecon.py b/brkraw/scripts/brkrecon.py index 05ffbbb..9105a1f 100644 --- a/brkraw/scripts/brkrecon.py +++ b/brkraw/scripts/brkrecon.py @@ -25,6 +25,8 @@ def main(): input_str = "input raw Bruker data" output_dir_str = "output directory name" + output_fnm_str = "output filename" + bids_opt = "create a JSON file contains metadata based on BIDS recommendation" subparsers = parser.add_subparsers(title='Sub-commands', description='To run this command, you must specify one of the functions listed' @@ -33,39 +35,125 @@ def main(): help='description', dest='function', metavar='command') - test = subparsers.add_parser("test", help='Run test mode') - test.add_argument("-i", "--input", help=input_str, type=str, default=None) - test.add_argument("-o", "--output", help=output_dir_str, type=str, default=None) + + nii = subparsers.add_parser("tonii", help='Convert all scans in a dataset to kspace, or complex image matrix') + nii.add_argument("input", help=input_str, type=str, default=None) + nii.add_argument("-b", "--bids", help=bids_opt, action='store_true') + nii.add_argument("-o", "--output", help=output_fnm_str, type=str, default=False) + nii.add_argument("-s", "--scanid", help="Scan ID, option to specify a particular scan to convert.", type=str) + nii.add_argument("-r", "--recoid", help="RECO ID (default=1), " + "option to specify a particular reconstruction id to convert", + type=int, default=1) + nii.add_argument("-t", "--subjecttype", help="override subject type in case the original setting was not properly set." + \ + "available options are (Biped, Quadruped, Phantom, Other, OtherAnimal)", type=str, default=None) + nii.add_argument("-p", "--position", help="override position information in case the original setting was not properly input." + \ + "the position variable can be defiend as _, " + \ + "available BodyParts are (Head, Foot, Tail) and sides are (Supine, Prone, Left, Right). (e.g. Head_Supine)", type=str, default=None) + + nii.add_argument("-f", "--formatting", help="FID processing methods" + \ + "available processing are (CKdata, image)", type=str, default='image') + + nii.add_argument("--ignore-slope", help='remove slope value from header', action='store_true') + nii.add_argument("--ignore-offset", help='remove offset value from header', action='store_true') + nii.add_argument("--ignore-rescale", help='remove slope and offset values from header', action='store_true') + nii.add_argument("--ignore-localizer", help='ignore the scan if it is localizer', action='store_true', default=True) args = parser.parse_args() - if args.function == 'test': - ipath = args.input - study = BrukerLoader(ipath) - print(ipath) + if args.function == 'tonii': + path = args.input + scan_id = args.scanid + reco_id = args.recoid + process = args.formatting + study = BrukerLoader(path) + slope, offset = set_rescale(args) + ignore_localizer = args.ignore_localizer if study.is_pvdataset: - - output = '{}_{}'.format(study._pvobj.subj_id,study._pvobj.study_id) - mkdir(output) - for id in study._avail.keys(): - fid_binary = study.get_fid(id) - acqp = study.get_acqp(id) - meth = study.get_method(id) - reco = study._pvobj.get_reco(id, 1) - process = 'image' - data = recon(fid_binary, acqp, meth, reco, process=process) + if args.output: + output = args.output + else: + output = '{}_{}'.format(study._pvobj.subj_id,study._pvobj.study_id) + if scan_id: + acqpars = study.get_acqp(int(scan_id)) + scanname = acqpars._parameters['ACQ_scan_name'] + scanname = scanname.replace(' ','-') + scan_id = int(scan_id) + reco_id = int(reco_id) - # Check if Image recontructed - if len(data.shape) > 3: - output_fname =f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}" - for c in range(data.shape[4]): - ni_img = nib.Nifti1Image(np.abs(np.squeeze(data)), affine=np.eye(4)) - nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) - print('NifTi file is generated... [{}]'.format(output_fname)) + if ignore_localizer and is_localizer(study, scan_id, reco_id): + print('Identified a localizer, the file will not be converted: ScanID:{}'.format(str(scan_id))) + else: + try: + recon2nifti(study, scan_id, reco_id, output, scanname, process) + except: + print('Conversion failed: ScanID:{}, RecoID:{}'.format(str(scan_id), str(reco_id))) + else: + for scan_id, recos in study._pvobj.avail_reco_id.items(): + acqpars = study.get_acqp(int(scan_id)) + scanname = acqpars._parameters['ACQ_scan_name'] + scanname = scanname.replace(' ','-') + if ignore_localizer and is_localizer(study, scan_id, recos[0]): + print('Identified a localizer, the file will not be converted: ScanID:{}'.format(str(scan_id))) + else: + try: + recon2nifti(study, scan_id, reco_id, output, scanname, process) + except Exception as e: + print('Conversion failed: ScanID:{}'.format(str(scan_id))) else: - print('{} is not PvDataset.'.format(ipath)) + print('{} is not PvDataset.'.format(path)) +def recon2nifti(pvobj, scan_id, reco_id, output, scanname, process): + output_fname = '{}-{}-{}-{}'.format(output, str(scan_id), reco_id, scanname) + visu_pars = pvobj._get_visu_pars(scan_id, 1) + method = pvobj._method[scan_id] + affine = pvobj._get_affine(visu_pars, method) + fid_binary = pvobj.get_fid(scan_id) + acqp = pvobj.get_acqp(scan_id) + reco = pvobj._pvobj.get_reco(scan_id, 1) + image = recon(fid_binary, acqp, method, reco, process = process,recoparts='default') + + if len(image.shape) > 7: + return + + #[x,y,z,_,n_channel,NI,NR] + image = image.transpose(0,1,2,5,4,6,3) + # MultiSlice Acq Correction + if '360' in acqp._parameters['ACQ_sw_version']: + if acqp._parameters['ACQ_dim'] == 2 and acqp._parameters['NSLICES'] > 1: + new_shape = list(image.shape) + new_shape[2] = acqp._parameters['NSLICES'] + new_shape[3] = int(new_shape[3]/acqp._parameters['NSLICES']) + image = image.reshape(new_shape) + image = image.transpose(1,0,2,3,4,5,6) + + else: + if acqp._parameters['ACQ_dim'] == 2 and acqp._parameters['NSLICES'] > 1: + new_shape = list(image.shape) + new_shape[2] = acqp._parameters['NSLICES'] + new_shape[3] = int(new_shape[3]/acqp._parameters['NSLICES']) + image = image.reshape(new_shape) + image = image.transpose(1,0,2,3,4,5,6) + + # [x, y, z, echo, channel, NR] + niiobj = nib.Nifti1Image(np.squeeze(np.abs(image)), affine) + niiobj = pvobj._set_nifti_header(niiobj, visu_pars, method, slope=False, offset=False) + niiobj.to_filename(output_fname+'-m'+'.nii.gz') + niiobj = nib.Nifti1Image(np.squeeze(np.angle(image)), affine) + niiobj = pvobj._set_nifti_header(niiobj, visu_pars, method, slope=False, offset=False) + niiobj.to_filename(output_fname+'-p'+'.nii.gz') + print('NifTi file is generated... [{}]'.format(output_fname)) + +def is_localizer(pvobj, scan_id, reco_id): + visu_pars = pvobj.get_visu_pars(scan_id, reco_id) + if 'VisuAcquisitionProtocol' in visu_pars.parameters: + ac_proc = visu_pars.parameters['VisuAcquisitionProtocol'] + if re.search('tripilot', ac_proc, re.IGNORECASE) or re.search('localizer', ac_proc, re.IGNORECASE): + return True + else: + return False + else: + return False if __name__ == '__main__': main() \ No newline at end of file diff --git a/brkraw/ui/main_win.py b/brkraw/ui/main_win.py index 5ee814f..97b8f66 100644 --- a/brkraw/ui/main_win.py +++ b/brkraw/ui/main_win.py @@ -35,6 +35,7 @@ def open_filediag(self): title = "Select file", filetypes = (("Zip compressed", "*.zip"), ("Paravision 6 format", "*.PVdatasets"), + ("Paravision 360 format", "*.PvDatasets") )) self._extend_layout() self._load_dataset() diff --git a/tests/compressed_sensing_example.py b/tests/compressed_sensing_example.py new file mode 100644 index 0000000..8ac522d --- /dev/null +++ b/tests/compressed_sensing_example.py @@ -0,0 +1,84 @@ +import os +import time + +import numpy as np +import matplotlib.pyplot as plt + +import brkraw as br +from brkraw.lib.parser import Parameter +from brkraw.lib.pvobj import PvDatasetDir +from brkraw.lib.utils import get_value, mkdir +import brkraw as br +from brkraw.lib.recon import recon + +import brkbart +from bart import bart + +import sigpy as sp +from sigpy import mri as mr +from sigpy import plot as pl + +import nibabel as nib + +PV_zipfile = "/home/jac/data/nmrsu_data/Tim_Wilson_Tim_207023_1_Default_20231202FeCl3_210800_360.3.4.PvDatasets" +data_loader = br.load(PV_zipfile) + +print(data_loader._avail.keys()) +ExpNum = 10 + +# Raw data processing for single job +fid_binary = data_loader.get_fid(ExpNum) +acqp = data_loader.get_acqp(ExpNum) +meth = data_loader.get_method(ExpNum) +reco = data_loader._pvobj.get_reco(ExpNum, 1) + +print(get_value(acqp, 'ACQ_protocol_name' ), ExpNum) +print('DIMS:', get_value(acqp, 'ACQ_dim' )) +# process OPTIONS: 'raw', 'frame', 'CKdata', 'image' +process = 'CKdata' + +# KSPACE DATA +data = recon(fid_binary, acqp, meth, reco, process=process, recoparts='all') + +# REFERENCE DATA +image = recon(fid_binary, acqp, meth, reco, process='image', + recoparts=['quadrature', 'phase_rotate', 'zero_filling', 'FT', 'phase_corr_pi', + 'cutoff', 'scale_phase_channels', 'sumOfSquares', 'transposition'] + ) +image = np.squeeze(image)[:,:,:,0] +data = data/np.max(np.abs(data)) # Normalize Data +print(data.shape) + +# SIMULATE POISSON UNDERSAMPLING +# NOTE coil dim is moved to first index for SIGPY standards +undersampled_data = np.zeros_like(data).transpose(4,0,1,2,3,5,6) +poisson = mr.poisson(data.shape[1:3],2) +for i in range(data.shape[4]): + for j in range(data.shape[5]): + for k in range(data.shape[6]): + undersampled_data[i,:,:,:,0,j,k] = data[:,:,:,0,i,j,k]*np.stack([poisson for _ in range(128)]) + +# Use only one set of data [coil,x,y,z,_,TE,TR] +undersampled_data = undersampled_data[:,:,:,:,0,0,0] +reconzf = sp.ifft(undersampled_data, axes=[1,2,3]) + + +# COMPRESSED SENSING ------------------- +sens = mr.app.EspiritCalib(undersampled_data).run() +#pl.ImagePlot(sens, z=0, title='Sensitivity Maps Estimated by ESPIRiT') +""" +lamda = 1 +img_sense = mr.app.SenseRecon(undersampled_data, sens, lamda=lamda).run() +""" +lamda = 0.0002 +img_L1 = mr.app.L1WaveletRecon(undersampled_data, sens, lamda=lamda).run() + +lamda = 0.001 +img_tv = mr.app.TotalVariationRecon(undersampled_data, sens, lamda=lamda).run() + + +#pl.ImagePlot(image, title='Fully Sampled') +#pl.ImagePlot(reconzf, title='Zero-fill Reconstruction') +#pl.ImagePlot(img_sense, title='SENSE Reconstruction') +#pl.ImagePlot(img_L1, title='L1Wavelet Reconstruction') +#pl.ImagePlot(img_tv, title='Total Variation Regularized Reconstruction') \ No newline at end of file diff --git a/tests/recon_cs.py b/tests/recon_cs.py new file mode 100644 index 0000000..3962f8b --- /dev/null +++ b/tests/recon_cs.py @@ -0,0 +1,32 @@ +import os +import time + +import numpy as np +import matplotlib.pyplot as plt + +import brkraw as br +from brkraw.lib.recon import * + +import sigpy as sp +from sigpy import mri as mr +from sigpy import plot as pl + +import nibabel as nib + + +# TEST WITH FUNCTION +folder = '/home/jac/data/data20230614' +MainDir = os.path.join(folder,'Price_Matt_7585_1_Default_CS_Trial_2_12718_360.3.4.PvDatasets') +print(MainDir) +i =8 +rawdata = br.load(os.path.join(MainDir)) +acqp = rawdata.get_acqp(i) +meth = rawdata.get_method(i) +reco = rawdata._pvobj.get_reco(i,1) +fid_binary = rawdata.get_fid(i) +data = recon(fid_binary, acqp, meth, reco, process='image') +print(data.shape) +pl.ImagePlot(np.squeeze(data)) + +data = recon(fid_binary, acqp, meth, reco, process='raw') +frame_test = convertRawToFrame(data, acqp, meth) diff --git a/tests/recon_testing.py b/tests/recon_testing_pv360.py similarity index 85% rename from tests/recon_testing.py rename to tests/recon_testing_pv360.py index 2f65e3b..c5d225e 100644 --- a/tests/recon_testing.py +++ b/tests/recon_testing_pv360.py @@ -35,7 +35,7 @@ reco = data_loader._pvobj.get_reco(ExpNum, 1) print(get_value(acqp, 'ACQ_protocol_name' ), ExpNum) - + # process OPTIONS: 'raw', 'frame', 'CKdata', 'image' process = 'image' @@ -46,13 +46,16 @@ data = data/np.max(np.abs(data)) # Normalize Data # Check if Image recontructed output = '{}_{}'.format(data_loader._pvobj.subj_id,data_loader._pvobj.study_id) - mkdir(output) + #mkdir(output) # Reconstructed Image Matrix is always 7-dimensional + #[x,y,z,_,n_channel,NI,NR] if len(data.shape) == 7: + print(data.shape) output_fname =f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}" + for c in range(data.shape[4]): - ni_img = nib.Nifti1Image(np.abs(np.squeeze(data)), affine=np.eye(4)) - nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) + ni_img = nib.Nifti1Image(np.abs(np.squeeze(data[:,:,:,:,c,:,:])), affine=np.eye(4)) + #nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) print('NifTi file is generated... [{}]'.format(output_fname)) \ No newline at end of file diff --git a/tests/recon_testing_pv6.py b/tests/recon_testing_pv6.py new file mode 100644 index 0000000..7fac76a --- /dev/null +++ b/tests/recon_testing_pv6.py @@ -0,0 +1,126 @@ +""" + DEVELOPED FOR BRUKER PARAVISION 6 datasets + Below code will work for cartesian sequence + GRE, MSME, RARE that were acquired with linear + phase encoding + IDK who (lucio) MAXYON + +@author: Tim Ho (UVA) +""" + +import os +import time + +import numpy as np +import matplotlib.pyplot as plt + +import brkraw as br +from brkraw.lib.parser import Parameter +from brkraw.lib.pvobj import PvDatasetDir +from brkraw.lib.utils import get_value, mkdir +from brkraw.lib.reference import BYTEORDER, WORDTYPE +import brkraw as br +from brkraw.lib.recon import * + +import nibabel as nib +import sigpy as sp + +PV_zipfile = '/home/jac/data/external_data/PV6/RARE_3D' +ExpNum = 35 +PV_zipfile = os.path.join(PV_zipfile,str(ExpNum)) +PV_zipfile = '/home/jac/data/external_data/PV6/RARE_3slice_packages/5' +# Raw data processing for single job +with open(os.path.join(PV_zipfile, 'method'), 'r') as f: + meth = Parameter(f.read().split('\n')) +with open(os.path.join(PV_zipfile, 'acqp'), 'r') as f: + acqp = Parameter(f.read().split('\n')) +print(os.path.join(PV_zipfile,'/pdata/1/reco')) +with open(os.path.join(PV_zipfile,'pdata/1','reco'), 'r') as f: + reco = Parameter(f.read().split('\n')) + +fid_path = os.path.join(PV_zipfile, 'fid') +print(fid_path) +if os.path.exists(fid_path): + with open(fid_path, 'rb') as f: + fid_binary = f.read() +else: + fid_path = os.path.join(PV_zipfile, 'rawdata.job0') + with open(fid_path, 'rb') as f: + fid_binary = f.read() +print(get_value(acqp, 'ACQ_sw_version')) +print(get_value(acqp, 'ACQ_protocol_name' ), ExpNum) + +# test functions +dt_code = np.dtype('int32') +bits = 32 # Follows dtype +dt_code = dt_code.newbyteorder('<') +fid = np.frombuffer(fid_binary, dt_code) +print(fid.shape) + +NI = get_value(acqp, 'NI') +NR = get_value(acqp, 'NR') +nRecs = 1 +if get_value(acqp,'ACQ_ReceiverSelect') != None: + nRecs = get_value(acqp, 'ACQ_ReceiverSelect').count('Yes') # THIS NEEDS TO BE EXAMINED BUT IDK HOW + +ACQ_size = get_value(acqp, 'ACQ_size' ) + +if get_value(acqp, 'GO_block_size') == 'Standard_KBlock_Format': + blocksize = int(np.ceil(ACQ_size[0]*nRecs*(bits/8)/1024)*1024/(bits/8)) +else: + blocksize = int(ACQ_size[0]*nRecs) + +# CHECK SIZE +print(fid.size, blocksize*np.prod(ACQ_size[1:])*NI*NR) +if fid.size != blocksize*np.prod(ACQ_size[1:])*NI*NR: + print('Error size dont match') + +# Reshape +fid = fid[::2] + 1j*fid[1::2] +fid = fid.reshape([-1,blocksize//2]) + +print(ACQ_size) + +# THIS REMOVES ZERO FILL (IDK THE PURPOSE FOR THIS) +if blocksize != ACQ_size[0]*nRecs: + print('a') + fid = fid[:,:ACQ_size[0]] + nRecs = 2 + fid = fid.reshape((-1,nRecs,ACQ_size[0]//2)) + fid = fid.transpose(0,1,2) + +print(fid.shape) +frame = convertRawToFrame(fid, acqp, meth) +CKdata = convertFrameToCKData(frame,acqp, meth) + +data = brkraw_Reco(CKdata, deepcopy(reco), meth, recoparts='all') + + +if len(data.shape) == 7: + print('data shape', data.shape) + output = '{}_{}'.format('pv6' , 'pv6') + mkdir(output) + output_fname =f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}" + + for c in range(data.shape[4]): + ni_img = nib.Nifti1Image(np.abs(np.squeeze(data[:,:,:,:,c,:,:])), affine=np.eye(4)) + nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) + print('NifTi file is generated... [{}]'.format(output_fname)) + +# Examine 2dseq file +with open(os.path.join(PV_zipfile, 'pdata/1/visu_pars'), 'r') as f: + visu_pars = Parameter(f.read().split('\n')) +dtype_code = np.dtype('{}{}'.format(BYTEORDER[get_value(visu_pars, 'VisuCoreByteOrder')], + WORDTYPE[get_value(visu_pars, 'VisuCoreWordType')])) + + +with open(os.path.join(PV_zipfile, 'pdata/1/2dseq'), 'rb') as f: + fid_binary = f.read() +_2dseq = np.frombuffer(fid_binary, dtype_code) +print(get_value(reco, 'RECO_size')[::-1]) +#plt.figure() +#plt.subplot(1,2,1) +#plt.imshow(_2dseq.reshape(get_value(reco, 'RECO_size')[::-1])[40,:,:]) +#plt.subplot(1,2,2) +#plt.imshow(np.abs(np.squeeze(data[:,:,:,0,0,0,0]))[:,:,40].T) +#plt.show() \ No newline at end of file diff --git a/tests/recon_testing_pv6_2.py b/tests/recon_testing_pv6_2.py new file mode 100644 index 0000000..16984a2 --- /dev/null +++ b/tests/recon_testing_pv6_2.py @@ -0,0 +1,130 @@ +""" + DEVELOPED FOR BRUKER PARAVISION 6 datasets + Below code will work for cartesian sequence + GRE, MSME, RARE that were acquired with linear + phase encoding + Tested with 2021_aswendt_gfap_pt_4wks + +@author: Tim Ho (UVA) +""" + +import os +import time + +import numpy as np +import matplotlib.pyplot as plt + +import brkraw as br +from brkraw.lib.parser import Parameter +from brkraw.lib.pvobj import PvDatasetDir +from brkraw.lib.utils import get_value, mkdir +import brkraw as br +from brkraw.lib.recon import * + +import nibabel as nib +import sigpy as sp + +root_path = '/home/jac/data/external_data/2021_aswendt_gfap_pt_4wks/MRI_raw_data/' +PV_dirs = [] +for root, dirs, files in os.walk(root_path, topdown=False): + for name in dirs: + if 'GV' in name: + PV_dirs.append(os.path.join(root, name)) + +for PV_zipfile in PV_dirs[:20]: + print(PV_zipfile) + data_loader = br.load(PV_zipfile) + for ExpNum in data_loader._avail.keys(): + + # Raw data processing for single job + fid_binary = data_loader.get_fid(ExpNum) + acqp = data_loader.get_acqp(ExpNum) + meth = data_loader.get_method(ExpNum) + reco = data_loader._pvobj.get_reco(ExpNum, 1) + #print(get_value(acqp, 'ACQ_sw_version')) + print(get_value(acqp, 'ACQ_protocol_name' ), ExpNum) + #print(get_value(acqp, 'ACQ_size' )) + try: + data = recon(fid_binary, acqp, meth, reco, process='image', recoparts='default') + #print(data.shape) + output = '{}_{}'.format(data_loader._pvobj.subj_id,data_loader._pvobj.study_id) + #mkdir(output) + if len(data.shape) == 7: + print(data.shape) + output_fname =f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}" + for c in range(data.shape[4]): + ni_img = nib.Nifti1Image(np.abs(np.squeeze(data[:,:,:,:,c,:,:])), affine=np.eye(4)) + #nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) + print('NifTi file is generated... [{}]'.format(output_fname)) + except Exception as e: + print(e, get_value(acqp, 'ACQ_protocol_name' ), ExpNum) + + """ + # test functions + dt_code = np.dtype('int32') + bits = 32 # Follows dtype + dt_code = dt_code.newbyteorder('<') + fid = np.frombuffer(fid_binary, dt_code) + + NI = get_value(acqp, 'NI') + NR = get_value(acqp, 'NR') + nRecs = 1 # THIS NEEDS TO BE EXAMINED BUT IDK HOW + ACQ_size = get_value(acqp, 'ACQ_size' ) + if type(ACQ_size) == int: + ACQ_size = [ACQ_size] + + if get_value(acqp, 'GO_block_size') == 'Standard_KBlock_Format': + blocksize = int(np.ceil(ACQ_size[0]*nRecs*(bits/8)/1024)*1024/(bits/8)) + else: + blocksize = int(ACQ_size[0]*nRecs) + + # CHECK SIZE + print(fid.size, blocksize*np.prod(ACQ_size[1:])*NI*NR) + if fid.size != blocksize*np.prod(ACQ_size[1:])*NI*NR: + print('Error Size dont match') + + # Reshape + fid = fid[::2] + 1j*fid[1::2] + fid = fid.reshape([-1,blocksize//2]) + + # THIS REMOVES ZERO FILL (IDK THE PURPOSE FOR THIS) + if blocksize != ACQ_size[0]*nRecs: + fid = fid[:,:ACQ_size[0]//2] + fid = fid.reshape((-1,ACQ_size[0]//2, nRecs)) + fid = fid.transpose(0,2,1) + + else: + #UNTESTED TIM FEB 12 2024 (IDK WHAT THIS DOES) + fid = fid.reshape((ACQ_size[0]//2, nRecs, -1)) + + frame = convertRawToFrame(fid, acqp, meth) + CKdata = convertFrameToCKData(frame,acqp, meth) + image = brkraw_Reco(CKdata, deepcopy(reco), meth, recoparts = 'default') + + print(fid.shape, frame.shape, CKdata.shape, image.shape) + """ +""" +frame = convertRawToFrame(fid, acqp, meth) +print(frame.shape) +#image = sp.ifft(np.squeeze(frame), axes=[0,1]) +CKdata = convertFrameToCKData(frame,acqp, meth) +print(CKdata.shape) +process = 'image' + +# test functions +data = brkraw_Reco(CKdata, deepcopy(reco), meth, recoparts = 'default') + +# ----------------------------------------------------------------- +data = recon(fid_binary, acqp, meth, reco, recoparts='default') +print(data.shape) + +output = '{}_{}'.format(data_loader._pvobj.subj_id,data_loader._pvobj.study_id) +mkdir(output) +if len(data.shape) == 7: + print(data.shape) + output_fname =f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}" + for c in range(data.shape[4]): + ni_img = nib.Nifti1Image(np.angle(np.squeeze(data[:,:,:,:,c,:,:])), affine=np.eye(4)) + nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) + print('NifTi file is generated... [{}]'.format(output_fname)) +""" \ No newline at end of file diff --git a/tests/recon_testing_pv6_3.py b/tests/recon_testing_pv6_3.py new file mode 100644 index 0000000..e61286a --- /dev/null +++ b/tests/recon_testing_pv6_3.py @@ -0,0 +1,109 @@ +""" + DEVELOPED FOR BRUKER PARAVISION 6 datasets + Below code will work for cartesian sequence + GRE, MSME, RARE that were acquired with linear + phase encoding + EALEXWater dataset + +@author: Tim Ho (UVA) +""" + +import os +import time + +import numpy as np +import matplotlib.pyplot as plt + +import brkraw as br +from brkraw.lib.parser import Parameter +from brkraw.lib.pvobj import PvDatasetDir +from brkraw.lib.utils import get_value, mkdir +import brkraw as br +from brkraw.lib.recon import * + +import nibabel as nib +import sigpy as sp + +PV_zipfile = '/home/jac/data/external_data/Test_Wat202402_15141272_1_Default_0206_Tomato_15141275_6.0.1.PvDatasets' +data_loader = br.load(PV_zipfile) + +for ExpNum in list(data_loader._avail.keys()): + # Raw data processing for single job + fid_binary = data_loader.get_fid(ExpNum) + acqp = data_loader.get_acqp(ExpNum) + meth = data_loader.get_method(ExpNum) + reco = data_loader._pvobj.get_reco(ExpNum, 1) + print(get_value(acqp, 'ACQ_sw_version')) + print(get_value(acqp, 'ACQ_protocol_name' ), ExpNum) + print(get_value(acqp, 'ACQ_size' )) + """ + # test functions + dt_code = np.dtype('int32') + bits = 32 # Follows dtype + dt_code = dt_code.newbyteorder('<') + fid = np.frombuffer(fid_binary, dt_code) + + NI = get_value(acqp, 'NI') + NR = get_value(acqp, 'NR') + nRecs = 4 # THIS NEEDS TO BE EXAMINED BUT IDK HOW + ACQ_size = get_value(acqp, 'ACQ_size' ) + + if get_value(acqp, 'GO_block_size') == 'Standard_KBlock_Format': + blocksize = int(np.ceil(ACQ_size[0]*nRecs*(bits/8)/1024)*1024/(bits/8)) + else: + blocksize = int(ACQ_size[0]*nRecs) + + # CHECK SIZE + print(fid.size, blocksize*np.prod(ACQ_size[1:])*NI*NR) + if fid.size != blocksize*np.prod(ACQ_size[1:])*NI*NR: + print('Error Size dont match') + + # Reshape + fid = fid[::2] + 1j*fid[1::2] + fid = fid.reshape([-1,blocksize//2]) + # THIS REMOVES ZERO FILL (IDK THE PURPOSE FOR THIS) + if blocksize != ACQ_size[0]*nRecs: + fid = fid[:,:ACQ_size[0]//2] + fid = fid.reshape((-1,ACQ_size[0]//2, nRecs)) + fid = fid.transpose(0,2,1) + + else: + print('hello') + #UNTESTED TIM FEB 12 2024 (IDK WHAT THIS DOES) + fid = fid.reshape((-1, nRecs, ACQ_size[0]//2)) + + + print(fid.shape) + frame = convertRawToFrame(fid, acqp, meth) + print(frame.shape) + #image = sp.ifft(np.squeeze(frame), axes=[0,1]) + CKdata = convertFrameToCKData(frame,acqp, meth) + print(CKdata.shape) + process = 'image' + + # test functions + data = brkraw_Reco(CKdata, deepcopy(reco), meth, recoparts = 'default') + plt.figure() + for c in range(data.shape[4]): + plt.subplot(1,4,c+1) + plt.imshow(np.abs(np.squeeze(data[:,:,:,:,c,0,:]))) + plt.show() + """ + # ----------------------------------------------------------------- + + data = recon(fid_binary, acqp, meth, reco, recoparts='all') + print(data.shape) + + if len(data.shape) == 7: + output = '{}_{}'.format(data_loader._pvobj.subj_id,data_loader._pvobj.study_id) + mkdir(output) + output_fname =f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}" + #plt.figure() + for c in range(data.shape[4]): + #plt.subplot(1,4,c+1) + #plt.imshow(np.abs(np.squeeze(data[:,:,7,:,c,0,:]))) + ni_img = nib.Nifti1Image(np.abs(np.squeeze(data[:,:,:,:,c,:,:])), affine=np.eye(4)) + nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) + #plt.show() + print('NifTi file is generated... [{}]'.format(output_fname)) + \ No newline at end of file diff --git a/tests/recon_testing_pv7.py b/tests/recon_testing_pv7.py new file mode 100644 index 0000000..3ba639d --- /dev/null +++ b/tests/recon_testing_pv7.py @@ -0,0 +1,99 @@ +""" + DEVELOPED FOR BRUKER PARAVISION 7 datasets + Below code will work for cartesian sequence + GRE, MSME, RARE that were acquired with linear + phase encoding + +@author: Tim Ho (UVA) +""" + +import os +import time + +import numpy as np +import matplotlib.pyplot as plt + +import brkraw as br +from brkraw.lib.parser import Parameter +from brkraw.lib.pvobj import PvDatasetDir +from brkraw.lib.utils import get_value, mkdir +import brkraw as br +from brkraw.lib.recon import * + +import nibabel as nib +import sigpy as sp + +PV_zipfile = '/home/jac/data/external_data/20231128_132106_hluna_piloFe_irm4_rata26_hluna_piloFe_irm4__1_1' +data_loader = br.load(PV_zipfile) + +ExpNum = 6 +for ExpNum in data_loader._avail.keys(): + # Raw data processing for single job + fid_binary = data_loader.get_fid(ExpNum) + acqp = data_loader.get_acqp(ExpNum) + meth = data_loader.get_method(ExpNum) + reco = data_loader._pvobj.get_reco(ExpNum, 1) + print(get_value(acqp, 'ACQ_sw_version')) + print(get_value(acqp, 'ACQ_protocol_name' ), ExpNum) + print(get_value(acqp, 'ACQ_size' )) + """ + # test functions + dt_code = np.dtype('int32') + bits = 32 # Follows dtype + dt_code = dt_code.newbyteorder('<') + fid = np.frombuffer(fid_binary, dt_code) + + NI = get_value(acqp, 'NI') + NR = get_value(acqp, 'NR') + nRecs = 1 # THIS NEEDS TO BE EXAMINED BUT IDK HOW + ACQ_size = get_value(acqp, 'ACQ_size' ) + + if get_value(acqp, 'GO_block_size') == 'Standard_KBlock_Format': + blocksize = int(np.ceil(ACQ_size[0]*nRecs*(bits/8)/1024)*1024/(bits/8)) + else: + blocksize = int(ACQ_size[0]*nRecs) + + # CHECK SIZE + print(fid.size, blocksize*np.prod(ACQ_size[1:])*NI*NR) + if fid.size != blocksize*np.prod(ACQ_size[1:])*NI*NR: + print('Error Size dont match') + + # Reshape + fid = fid[::2] + 1j*fid[1::2] + fid = fid.reshape([-1,blocksize//2]) + + # THIS REMOVES ZERO FILL (IDK THE PURPOSE FOR THIS) + if blocksize != ACQ_size[0]*nRecs: + fid = fid[:,:ACQ_size[0]//2] + fid = fid.reshape((-1,ACQ_size[0]//2, nRecs)) + fid = fid.transpose(0,2,1) + + else: + #UNTESTED TIM FEB 12 2024 (IDK WHAT THIS DOES) + fid = fid.reshape((ACQ_size[0]//2, nRecs, -1)) + + + print(fid.shape) + frame = convertRawToFrame(fid, acqp, meth) + print(frame.shape) + #image = sp.ifft(np.squeeze(frame), axes=[0,1]) + CKdata = convertFrameToCKData(frame,acqp, meth) + print(CKdata.shape) + process = 'image' + + # test functions + data = brkraw_Reco(CKdata, deepcopy(reco), meth, recoparts = 'default') + """ + # ----------------------------------------------------------------- + data = recon(fid_binary, acqp, meth, reco, recoparts='default') + print(data.shape) + + if len(data.shape) == 7: + output = '{}_{}'.format(data_loader._pvobj.subj_id,data_loader._pvobj.study_id) + mkdir(output) + output_fname =f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}" + for c in range(data.shape[4]): + ni_img = nib.Nifti1Image(np.angle(np.squeeze(data[:,:,:,:,c,:,:])), affine=np.eye(4)) + nib.save(ni_img, os.path.join(output,f"{acqp._parameters['ACQ_scan_name'].strip().replace(' ','_')}_C{c}.nii.gz")) + print('NifTi file is generated... [{}]'.format(output_fname)) + \ No newline at end of file