diff --git a/brkraw/lib/recoFunctions.py b/brkraw/lib/recoFunctions.py new file mode 100644 index 0000000..c40e09b --- /dev/null +++ b/brkraw/lib/recoFunctions.py @@ -0,0 +1,306 @@ +# -*- coding: utf-8 -*- +""" + DEVELOPED FOR BRUKER PARAVISION 360 datasets + Functions below are made to support functions + in recon.py + +@author: Tim Ho (UVA) +""" + +from .utils import get_value +import numpy as np + +def reco_qopts(frame, Reco, actual_framenumber): + + # import variables + RECO_qopts = get_value(Reco, 'RECO_qopts') + + # claculate additional parameters + dims = [frame.shape[0], frame.shape[1], frame.shape[2], frame.shape[3]] + + # check if the qneg-Matrix is necessary: + use_qneg = False + if (RECO_qopts.count('QUAD_NEGATION') + RECO_qopts.count('CONJ_AND_QNEG')) >= 1: + use_qneg = True + qneg = np.ones(frame.shape) # Matrix containing QUAD_NEGATION multiplication matrix + + # start process + for i in range(len(RECO_qopts)): + if RECO_qopts[i] == 'COMPLEX_CONJUGATE': + frame = np.conj(frame) + elif RECO_qopts[i] == 'QUAD_NEGATION': + if i == 0: + qneg = qneg * np.tile([[1, -1]], [np.ceil(dims[0]/2), dims[1], dims[2], dims[3]]) + elif i == 1: + qneg = qneg * np.tile([[1], [-1]], [dims[0], np.ceil(dims[1]/2), dims[2], dims[3]]) + elif i == 2: + tmp = np.zeros([1, 1, dims[2], 2]) + tmp[0, 0, :, :] = [[1, -1]] + qneg = qneg * np.tile(tmp, [dims[0], dims[1], np.ceil(dims[2]/2), dims[3]]) + elif i == 3: + tmp = np.zeros([1, 1, 1, dims[3], 2]) + tmp[0, 0, 0, :, :] = [[1, -1]] + qneg = qneg * np.tile(tmp, [dims[0], dims[1], dims[2], np.ceil(dims[3]/2)]) + elif RECO_qopts[i] == 'CONJ_AND_QNEG': + frame = np.conj(frame) + if i == 0: + qneg = qneg * np.tile([[1, -1]], [np.ceil(dims[0]/2), dims[1], dims[2], dims[3]]) + elif i == 1: + qneg = qneg * np.tile([[1], [-1]], [dims[0], np.ceil(dims[1]/2), dims[2], dims[3]]) + elif i == 2: + tmp = np.zeros([1, 1, dims[2], 2]) + tmp[0, 0, :, :] = [[1, -1]] + qneg = qneg * np.tile(tmp, [dims[0], dims[1], np.ceil(dims[2]/2), dims[3]]) + elif i == 3: + tmp = np.zeros([1, 1, 1, dims[3], 2]) + tmp[0, 0, 0, :, :] = [[1, -1]] + qneg = qneg * np.tile(tmp, [dims[0], dims[1], dims[2], np.ceil(dims[3]/2)]) + + if use_qneg: + if qneg.shape != frame.shape: + qneg = qneg[0:dims[0], 0:dims[1], 0:dims[2], 0:dims[3]] + frame = frame * qneg + + return frame + + +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] + else: + RECO_rotate = get_value(Reco,'RECO_rotate')[:,0] + + if isinstance( get_value(Reco,'RECO_ft_mode'), list): + if any(x != get_value(Reco,'RECO_ft_mode')[0] for x in get_value(Reco,'RECO_ft_mode')): + raise ValueError('It''s not allowed to use different transfomations on different Dimensions: ' + Reco['RECO_ft_mode']) + RECO_ft_mode = get_value(Reco,'RECO_ft_mode')[0] + else: + RECO_ft_mode = get_value(Reco,'RECO_ft_mode') + + # calculate additional variables + dims = [frame.shape[0], frame.shape[1], frame.shape[2], frame.shape[3]] + + # start process + 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') + + if index == 0: + phase_matrix *= np.tile(phase_vector[:,np.newaxis,np.newaxis,np.newaxis], [1, dims[1], dims[2], dims[3]]) + elif index == 1: + phase_matrix *= np.tile(phase_vector[np.newaxis,:,np.newaxis,np.newaxis], [dims[0], 1, dims[2], dims[3]]) + elif index == 2: + tmp = np.zeros((1,1,dims[2],1), dtype=complex) + tmp[0,0,:,0] = phase_vector + phase_matrix *= np.tile(tmp, [dims[0], dims[1], 1, dims[3]]) + elif index == 3: + tmp = np.zeros((1,1,1,dims[3]), dtype=complex) + tmp[0,0,0,:] = phase_vector + phase_matrix *= np.tile(tmp, [dims[0], dims[1], dims[2], 1]) + + frame *= phase_matrix + return frame + + +def reco_zero_filling(frame, Reco, actual_framenumber, signal_position): + # check input + RECO_ft_mode = get_value(Reco,'RECO_ft_mode') + + # Check if Reco.RECO_ft_size is not equal to size(frame) + not_Equal = any([(i != j) for i,j in zip(frame.shape,get_value(Reco, 'RECO_ft_size'))]) + + if not_Equal: + if any(signal_position > 1) or any(signal_position < 0): + raise ValueError('signal_position has to be a vector between 0 and 1') + + RECO_ft_size = get_value(Reco,'RECO_ft_size') + + # check if ft_size is correct: + for i in range(len(RECO_ft_size)): + if RECO_ft_size[i] < frame.shape[i]: + raise ValueError('RECO_ft_size has to be bigger than the size of your data-matrix') + + # calculate additional variables + dims = (frame.shape[0], frame.shape[1], frame.shape[2], frame.shape[3]) + + # start process + + # Dimensions of frame and RECO_ft_size doesn't match? -> zero filling + if not_Equal: + newframe = np.zeros(RECO_ft_size, dtype=complex) + startpos = np.zeros(len(RECO_ft_size), dtype=int) + pos_ges = [None] * 4 + + for i in range(len(RECO_ft_size)): + diff = RECO_ft_size[i] - frame.shape[i] + 1 + startpos[i] = int(np.floor(diff * signal_position[i] + 1)) + if startpos[i] > RECO_ft_size[i]: + startpos[i] = RECO_ft_size[i] + pos_ges[i] = slice(startpos[i] - 1, startpos[i] - 1 + dims[i]) + + newframe[pos_ges[0], pos_ges[1], pos_ges[2], pos_ges[3]] = frame + else: + newframe = frame + + del startpos, pos_ges + + else: + newframe = frame + + return newframe + + +def reco_FT(frame, Reco, actual_framenumber): + """ + Perform Fourier Transform on the input frame according to the specified RECO_ft_mode in the Reco dictionary. + + Args: + frame: ndarray + Input frame to perform Fourier Transform on + Reco: dict + Dictionary containing the specified Fourier Transform mode (RECO_ft_mode) + actual_framenumber: int + Index of the current frame + + Returns: + frame: ndarray + Output frame after Fourier Transform has been applied + """ + + # Import variables + RECO_ft_mode = get_value(Reco,'RECO_ft_mode')[0] + + # Start process + if RECO_ft_mode in ['COMPLEX_FT', 'COMPLEX_FFT']: + frame = np.fft.fftn(frame) + #frame = sp.fft(frame, axes=[0,1,2], center=False) + 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 = sp.ifft(frame, axes=[0,1,2], center=False) + else: + raise ValueError('Your RECO_ft_mode is not supported') + + return frame + + +def reco_phase_corr_pi(frame, Reco, actual_framenumber): + # start process + checkerboard = np.ones(shape=frame.shape) + + # Use NumPy broadcasting to alternate the signs + checkerboard[::2,::2,::2,0] = -1 + checkerboard[1::2,1::2,::2,0] = -1 + + checkerboard[::2,1::2,1::2,0] = -1 + checkerboard[1::2,::2,1::2,0] = -1 + + frame = frame * checkerboard * -1 + + return frame + + +def reco_cutoff(frame, Reco, actual_framenumber): + """ + Crops the input frame according to the specified RECO_size in the Reco dictionary. + + Args: + frame: ndarray + Input frame to crop + Reco: dict + Dictionary containing the specified crop size (RECO_size) and offset (RECO_offset) + actual_framenumber: int + Index of the current frame + + Returns: + newframe: ndarray + Cropped output frame + """ + + # 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: + + # Import variables + RECO_offset = get_value(Reco,'RECO_offset')[:, actual_framenumber] + RECO_size = get_value(Reco, 'RECO_size') + + # 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 + + return newframe + + +def reco_scale_phase_channels(frame, Reco, channel): + # check input + reco_scale = get_value(Reco,'RecoScaleChan') + if not isinstance(reco_scale, list): + reco_scale = [reco_scale] + if channel <= len(reco_scale) and reco_scale != None: + scale = reco_scale[int(channel)] + else: + scale = 1.0 + + reco_phase = get_value(Reco,'RecoPhaseChan') + + if not isinstance(reco_phase, list): + reco_phase = [reco_phase] + if channel <= len(reco_phase) and reco_phase != None: + phase = reco_phase[int(channel)] + else: + phase = 0.0 + + spFactor = scale * np.exp(1j * phase * np.pi / 180.0) + # multiply each pixel by common scale and phase factor + frame = spFactor * frame + return frame + + +def reco_sumofsquares(frame, Reco): + out = np.sqrt( np.sum(np.square(np.abs(frame)), axis=4, keepdims=True) ) + return out + + +def reco_transposition(frame, Reco, actual_framenumber): + # Import variables + RECO_transposition = get_value(Reco,'RECO_transposition')[actual_framenumber - 1] + + # Calculate additional variables + dims = [frame.shape[i] for i in range(4)] + + # Start process + if RECO_transposition > 0: + ch_dim1 = (RECO_transposition % 4) + ch_dim2 = RECO_transposition - 1 + new_order = list(range(4)) + new_order[int(ch_dim1)] = ch_dim2 + new_order[int(ch_dim2)] = ch_dim1 + frame = np.transpose(frame, new_order) + frame = np.reshape(frame, dims) + + return frame \ No newline at end of file diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py new file mode 100644 index 0000000..c3860fa --- /dev/null +++ b/brkraw/lib/recon.py @@ -0,0 +1,475 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat Jan 20 10:06:38 2024 + DEVELOPED FOR BRUKER PARAVISION 360 datasets + Below code will work for cartesian sequence + GRE, MSME, RARE that were acquired with linear + phase encoding + +@author: Tim Ho (UVA) +""" + +from .utils import get_value, set_value +from .recoFunctions import * +import numpy as np +from copy import deepcopy + + +def recon(fid_binary, acqp, meth, reco, process = None, recoparts = 'default'): + """ Process FID -> Channel Sorted -> Frame-Sorted -> K-Sorted -> Image + + Parameters + ---------- + fid_binary : bytestring + + acqp : dict (brkraw Parameter structure) + + meth : dict (brkraw Parameter structure) + + reco : dict (brkraw Parameter structure) + + process: 'raw', 'frame', 'CKdata', 'image' + + recoparts: 'default' or list() + List Options: ['quadrature', 'phase_rotate', 'zero_filling', 'FT', 'phase_corr_pi', + 'cutoff', 'scale_phase_channels', 'sumOfSquares', 'transposition'] + + 'default': + recoparts = ['quadrature', 'phase_rotate', 'zero_filling', + 'FT', 'phase_corr_pi'] + + Returns + ------- + output : np.array + + """ + 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 \ + 'localizer' in get_value(acqp, 'ACQ_protocol_name').lower() or \ + 'gre' in get_value(acqp, 'ACQ_protocol_name').lower() or \ + 'mge' in get_value(acqp, 'ACQ_protocol_name').lower() or \ + 'FLASH.ppg' == get_value(acqp, 'PULPROG'): + + output = convertRawToFrame(output, acqp, meth) + if process == 'frame': + return output + + output = convertFrameToCKData(output, acqp, meth) + if process == 'CKdata': + return output + + output = brkraw_Reco(output, deepcopy(reco), meth, recoparts = recoparts) + + else: + print("Warning: SEQUENCE PROTOCOL {} NOT SUPPORTED...".format(get_value(acqp, 'ACQ_scan_name' ))) + print("returning 'raw' sorting") + + return output + + +def readBrukerRaw(fid_binary, acqp, meth): + """ Sorts FID into a 3D np matrix [num_readouts, channel, scan_size] + + Parameters + ---------- + fid_binary : bytestring + + acqp : dict (brkraw Parameter structure) + + meth : dict (brkraw Parameter structure) + + + Returns + ------- + X : np.array + """ + 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) + + 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)]) + + return X + + +def convertRawToFrame(data, acqp, meth): + """ Prelinary raw fid sorting + + Parameters + ---------- + data : np.array with size [num_readouts, channel, scan_size] + + acqp : dict (brkraw Parameter structure) + + meth : dict (brkraw Parameter structure) + + Returns + ------- + frame : np.array with size + [scansize, ACQ_phase_factor, numDataHighDim/ACQ_phase_factor, numSelectedReceivers, NI, NR] + """ + + results = data.copy() + + # Turn Raw into frame + NI = get_value(acqp, 'NI') + NR = get_value(acqp, 'NR') + + if 'rare' in get_value(meth,'Method').lower(): + ACQ_phase_factor = get_value(meth,'PVM_RareFactor') + else: + ACQ_phase_factor = get_value(acqp,'ACQ_phase_factor') + + ACQ_obj_order = get_value(acqp, 'ACQ_obj_order') + if isinstance(ACQ_obj_order, int): + ACQ_obj_order = [ACQ_obj_order] + + ACQ_dim = get_value(acqp, 'ACQ_dim') + numSelectedRecievers= results.shape[-2] + + acqSizes = np.zeros(ACQ_dim) + scanSize = get_value(acqp, 'ACQ_jobs')[0][0] + + isSpatialDim = [i == 'Spatial' for i in get_value(acqp, 'ACQ_dim_desc')] + spatialDims = sum(isSpatialDim) + + if isSpatialDim[0]: + if spatialDims == 3: + acqSizes[1:] = get_value(meth, 'PVM_EncMatrix')[1:] + elif spatialDims == 2: + acqSizes[1] = get_value(meth, 'PVM_EncMatrix')[1] + + numresultsHighDim=np.prod(acqSizes[1:]) + acqSizes[0] = scanSize + + if np.iscomplexobj(results): + scanSize = int(acqSizes[0]/2) + else: + scanSize = acqSizes[0] + + # Resort + if ACQ_dim>1: + # [..., num_lines, channel, scan_size] + results = results.transpose((1,2,0)) + + results = results.reshape( + int(numSelectedRecievers), + int(scanSize), + int(ACQ_phase_factor), + int(NI), + int(numresultsHighDim/ACQ_phase_factor), + int(NR), order='F') + + # reorder to [scansize, ACQ_phase_factor, numDataHighDim/ACQ_phase_factor, numSelectedReceivers, NI, NR] + results = np.transpose(results, (1, 2, 4, 0, 3, 5)) + + results = results.reshape( + int(scanSize), + int(numresultsHighDim), + int(numSelectedRecievers), + int(NI), + int(NR), order='F') + + frame = np.zeros_like(results) + frame[:,:,:,ACQ_obj_order,:] = results + + else: + # 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' + + + return frame + + +def convertFrameToCKData(frame, acqp, meth): + """Frame to unprocessed KSPACE + + Parameters + ---------- + frame : np.array with size [num_readouts, channel, scan_size] + + acqp : dict (brkraw Parameter structure) + + meth : dict (brkraw Parameter structure) + + Returns + ------- + data : np.array with size + [scansize, ACQ_phase_factor, numDataHighDim/ACQ_phase_factor, numSelectedReceivers, NI, NR] + for simplified understand of the data structure + [x,y,z,_,n_channel,NI,NR] + """ + NI = get_value(acqp, 'NI') + NR = get_value(acqp, 'NR') + + ACQ_phase_factor = get_value(acqp,'ACQ_phase_factor') + ACQ_obj_order = get_value(acqp, 'ACQ_obj_order') + ACQ_dim = get_value(acqp, 'ACQ_dim') + numSelectedReceivers= frame.shape[2] + + acqSizes = np.zeros(ACQ_dim) + + scanSize = get_value(acqp, 'ACQ_jobs')[0][0] + acqSizes[0] = scanSize + ACQ_size = acqSizes + + isSpatialDim = [i == 'Spatial' for i in get_value(acqp, 'ACQ_dim_desc')] + spatialDims = sum(isSpatialDim) + + if isSpatialDim[0]: + if spatialDims == 3: + acqSizes[1:] = get_value(meth, 'PVM_EncMatrix')[1:] + elif spatialDims == 2: + acqSizes[1] = get_value(meth, 'PVM_EncMatrix')[1] + numDataHighDim=np.prod(acqSizes[1:]) + + if np.iscomplexobj(frame): + scanSize = int(acqSizes[0]/2) + else: + scanSize = acqSizes[0] + + if ACQ_dim==3: + + PVM_EncSteps2=get_value(meth, 'PVM_EncSteps2') + assert PVM_EncSteps2 != None + + PVM_Matrix = get_value(meth, 'PVM_Matrix') + PVM_EncSteps1 = get_value(meth,'PVM_EncSteps1') + + PVM_AntiAlias = get_value(meth, 'PVM_AntiAlias') + if PVM_AntiAlias == None: + # No anti-aliasing available. + PVM_AntiAlias = np.ones((ACQ_dim)) + + PVM_EncZf=get_value(meth, 'PVM_EncZf') + + if PVM_EncZf == None: + # No zero-filling/interpolation available. + PVM_EncZf = np.ones((ACQ_dim)) + + # Resort + + 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)) + # Step 2: Interpolation + + 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) + ckCenterIndex = np.floor(ckSize/2 + 0.25) + 1 + + readStartIndex = int(ckSize[0]-scanSize + 1) + + + # Reshape & store + # switch ACQ_dim + 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) + data[readStartIndex-1:,0,0,0,:,:,:] = frameData + + elif ACQ_dim == 2: + frameData=np.reshape(frameData,(scanSize, int(ACQ_size[1]), 1, 1, numSelectedReceivers, NI, NR) , order='F') + data=np.zeros([int(ckSize[0]), int(ckSize[1]), 1, 1, numSelectedReceivers, NI, NR], dtype=complex) + encSteps1indices = (PVM_EncSteps1 + ckCenterIndex[1] - 1).astype(int) + data[readStartIndex-1:,encSteps1indices,0,0,:,:,:] = frameData.reshape(data[readStartIndex-1:,encSteps1indices,0,0,:,:,:].shape, order='F') + + elif ACQ_dim == 3: + frameData = np.reshape(frameData,(scanSize, int(ACQ_size[1]), int(ACQ_size[2]), 1, numSelectedReceivers, NI, NR) , order='F') + data=np.zeros([int(ckSize[0]), int(ckSize[1]), int(ckSize[2]), 1, numSelectedReceivers, NI, NR], dtype=complex) + encSteps1indices = (PVM_EncSteps1 + ckCenterIndex[1] - 1).astype(int) + encSteps2indices = (PVM_EncSteps2 + ckCenterIndex[2] - 1).astype(int) + + data[readStartIndex-1:,list(encSteps1indices),:,:,:,:,:] = frameData[:,:,list(encSteps2indices),:,:,:,:] + + else: + raise 'Unknown ACQ_dim with useMethod' + + return data + + +def brkraw_Reco(kdata, reco, meth, recoparts = 'all'): + reco_result = kdata.copy() + + if recoparts == 'all': + recoparts = ['quadrature', 'phase_rotate', 'zero_filling', 'FT', 'phase_corr_pi', + 'cutoff', 'scale_phase_channels', 'sumOfSquares', 'transposition'] + 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') + + # Adapt FT convention to acquisition version. + N1, N2, N3, N4, N5, N6, N7 = kdata.shape + + dims = kdata.shape[0:4] + for i in range(4): + if dims[i]>1: + dimnumber = (i+1) + + NINR=kdata.shape[5]*kdata.shape[6] + signal_position=np.ones(shape=(dimnumber,1))*0.5 + + same_transposition = True + RECO_transposition = get_value(reco, 'RECO_transposition') + + if not isinstance(RECO_transposition, list): + RECO_transposition = [RECO_transposition] + for i in RECO_transposition: + if i != RECO_transposition[0]: + same_transposition = False + + map_index= np.reshape( np.arange(0,kdata.shape[5]*kdata.shape[6]), (kdata.shape[6], kdata.shape[5]) ).flatten() + + # --- START RECONSTRUCTION --- + for recopart in recoparts: + if 'quadrature' in recopart: + for NR in range(N7): + for NI in range(N6): + for channel in range(N5): + reco_result[:,:,:,:,channel,NI,NR] = reco_qopts(kdata[:,:,:,:,channel,NI,NR], reco, map_index[(NI+1)*(NR+1)-1]) + + + if 'phase_rotate' in recopart: + for NR in range(N7): + for NI in range(N6): + for channel in range(N5): + reco_result[:,:,:,:,channel,NI,NR] = reco_phase_rotate(kdata[:,:,:,:,channel,NI,NR], reco, map_index[(NI+1)*(NR+1)-1]) + + + if 'zero_filling' in recopart: + RECO_ft_size = get_value(reco,'RECO_ft_size') + + newdata_dims=[1, 1, 1, 1] + + newdata_dims[0:len(RECO_ft_size)] = RECO_ft_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_zero_filling(reco_result[:,:,:,:,chan,NI,NR], reco, map_index[(NI+1)*(NR+1)-1], signal_position).reshape(newdata[:,:,:,:,chan,NI,NR].shape) + + reco_result=newdata + + + if 'FT' in recopart: + RECO_ft_mode = get_value(reco,'RECO_ft_mode')[0] + for NR in range(N7): + for NI in range(N6): + for chan in range(N5): + reco_result[:,:,:,:,chan,NI,NR] = reco_FT(reco_result[:,:,:,:,chan,NI,NR], reco, map_index[(NI+1)*(NR+1)-1]) + + + if 'phase_corr_pi' in recopart: + for NR in range(N7): + for NI in range(N6): + 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): + for NI in range(N6): + for chan in range(N5): + reco_result[:,:,:,:,chan,NI,NR] = reco_scale_phase_channels(reco_result[:,:,:,:,chan,NI,NR], reco, chan) + + + if 'sumOfSquares' in recopart: + reco_result = np.sqrt( np.sum(np.square(np.abs(reco_result)), axis=4, keepdims=True)) + reco_result = reco_result[:,:,:,:,:1,:,:] + reco_result = np.real(reco_result) + + N5 = 1 # Update N5 + + if 'transposition' in recopart: + if same_transposition: + # import variables + RECO_transposition = RECO_transposition[0] + # calculate additional variables: + + # start process + if RECO_transposition > 0: + ch_dim1 = (RECO_transposition % len(kdata.shape)) + ch_dim2 = RECO_transposition - 1 + new_order = [0, 1, 2, 3] + new_order[int(ch_dim1)] = int(ch_dim2) + new_order[int(ch_dim2)] = int(ch_dim1) + reco_result = reco_result.transpose(new_order + [4, 5, 6]) + else: + for NR in range(N7): + for NI in range(N6): + 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/lib/utils.py b/brkraw/lib/utils.py index ef96b30..14e12a2 100644 --- a/brkraw/lib/utils.py +++ b/brkraw/lib/utils.py @@ -154,6 +154,13 @@ def get_value(pars, key): else: return pars.parameters[key] +def set_value(pars, key, value): + if key not in pars.parameters.keys(): + print('Created a key') + pars.parameters[key] = value + else: + pars.parameters[key] = value + return pars def is_all_element_same(listobj): if listobj is None: diff --git a/brkraw/scripts/brkrecon.py b/brkraw/scripts/brkrecon.py new file mode 100644 index 0000000..05ffbbb --- /dev/null +++ b/brkraw/scripts/brkrecon.py @@ -0,0 +1,71 @@ +""" +@author: Timothy Ho (UVA) +""" +# -*- coding: utf-8 -*- +from operator import index +from ..lib.errors import * +from .. import BrukerLoader, __version__ +from ..lib.utils import set_rescale, save_meta_files, mkdir +from ..lib.recon import recon +import argparse +import os, re +import sys + +import numpy as np +import nibabel as nib + +_supporting_bids_ver = '1.2.2' + + +def main(): + parser = argparse.ArgumentParser(prog='brkrecon', + description="Brkrecon command-line interface") + parser.add_argument("-v", "--version", action='version', version='%(prog)s v{}'.format(__version__)) + + + input_str = "input raw Bruker data" + output_dir_str = "output directory name" + + subparsers = parser.add_subparsers(title='Sub-commands', + description='To run this command, you must specify one of the functions listed' + 'below next to the command. For more information on each function, ' + 'use -h next to the function name to call help document.', + 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) + + args = parser.parse_args() + + if args.function == 'test': + ipath = args.input + study = BrukerLoader(ipath) + print(ipath) + + 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) + + # 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)) + else: + print('{} is not PvDataset.'.format(ipath)) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 1b1ece0..c1a68cb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,4 +46,5 @@ style = "pep440" [project.scripts] brkraw = "brkraw.scripts.brkraw:main" -brk-backup = 'brkraw.scripts.brk_backup:main' \ No newline at end of file +brk-backup = 'brkraw.scripts.brk_backup:main' +brkrecon = "brkraw.scripts.brkrecon:main" \ No newline at end of file diff --git a/tests/recon_testing.py b/tests/recon_testing.py new file mode 100644 index 0000000..2f65e3b --- /dev/null +++ b/tests/recon_testing.py @@ -0,0 +1,58 @@ +""" + DEVELOPED FOR BRUKER PARAVISION 360 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 recon + +import nibabel as nib + + +#PV_zipfile = "/home/jac/data/PV36034results20230630/Wilson_Tim_26171_1_Default_RAREvfl_26174_360.3.4.PvDatasets" +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) + +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_protocol_name' ), ExpNum) + + # process OPTIONS: 'raw', 'frame', 'CKdata', 'image' + process = 'image' + + # test functions + start_time = time.time() + data = recon(fid_binary, acqp, meth, reco, process=process) + print("{} convert {} seconds".format (process, time.time()-start_time)) + 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) + + # Reconstructed Image Matrix is always 7-dimensional + if len(data.shape) == 7: + 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)) + \ No newline at end of file