From 05a50c9aec613d80f37c341b194dbc3fd32f11a8 Mon Sep 17 00:00:00 2001 From: Tim Date: Sat, 20 Apr 2024 14:03:18 -0400 Subject: [PATCH 1/9] recon tests not updated to new structure --- tests/compressed_sensing_example.py | 84 ------------------ tests/recon_cs.py | 32 ------- tests/recon_testing_pv360.py | 61 ------------- tests/recon_testing_pv6.py | 126 --------------------------- tests/recon_testing_pv6_2.py | 130 ---------------------------- tests/recon_testing_pv6_3.py | 109 ----------------------- tests/recon_testing_pv7.py | 99 --------------------- 7 files changed, 641 deletions(-) delete mode 100644 tests/compressed_sensing_example.py delete mode 100644 tests/recon_cs.py delete mode 100644 tests/recon_testing_pv360.py delete mode 100644 tests/recon_testing_pv6.py delete mode 100644 tests/recon_testing_pv6_2.py delete mode 100644 tests/recon_testing_pv6_3.py delete mode 100644 tests/recon_testing_pv7.py diff --git a/tests/compressed_sensing_example.py b/tests/compressed_sensing_example.py deleted file mode 100644 index 8ac522d..0000000 --- a/tests/compressed_sensing_example.py +++ /dev/null @@ -1,84 +0,0 @@ -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 deleted file mode 100644 index 3962f8b..0000000 --- a/tests/recon_cs.py +++ /dev/null @@ -1,32 +0,0 @@ -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_pv360.py b/tests/recon_testing_pv360.py deleted file mode 100644 index c5d225e..0000000 --- a/tests/recon_testing_pv360.py +++ /dev/null @@ -1,61 +0,0 @@ -""" - 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 - #[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[:,:,:,:,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 deleted file mode 100644 index 7fac76a..0000000 --- a/tests/recon_testing_pv6.py +++ /dev/null @@ -1,126 +0,0 @@ -""" - 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 deleted file mode 100644 index 16984a2..0000000 --- a/tests/recon_testing_pv6_2.py +++ /dev/null @@ -1,130 +0,0 @@ -""" - 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 deleted file mode 100644 index e61286a..0000000 --- a/tests/recon_testing_pv6_3.py +++ /dev/null @@ -1,109 +0,0 @@ -""" - 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 deleted file mode 100644 index 3ba639d..0000000 --- a/tests/recon_testing_pv7.py +++ /dev/null @@ -1,99 +0,0 @@ -""" - 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 From fe7faae5652a56908ba5d243ce89edcb34bd68b2 Mon Sep 17 00:00:00 2001 From: Tim Date: Sat, 20 Apr 2024 22:12:59 -0400 Subject: [PATCH 2/9] Remake Recon Plugin --- brkraw/lib/recoFunctions.py | 305 ++------------- brkraw/lib/recon.py | 725 +++++++++++------------------------- 2 files changed, 247 insertions(+), 783 deletions(-) diff --git a/brkraw/lib/recoFunctions.py b/brkraw/lib/recoFunctions.py index d4e734c..654291c 100644 --- a/brkraw/lib/recoFunctions.py +++ b/brkraw/lib/recoFunctions.py @@ -3,302 +3,75 @@ 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') +def phase_rotate(frame, RECO_rotate, framenumber): - if RECO_rotate_all.shape[1] > actual_framenumber: - RECO_rotate = get_value(Reco,'RECO_rotate')[:, actual_framenumber] + if RECO_rotate.shape[1] > framenumber: + RECO_rotate = RECO_rotate[:, framenumber] + RECO_rotate -= 0.5 else: - RECO_rotate = get_value(Reco,'RECO_rotate')[:,0] + RECO_rotate = 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]] + dims = [frame.shape[0], frame.shape[1], frame.shape[2]] - # start process - phase_matrix = np.ones_like(frame) + # Create Shift matrix in KSPACE + phase_matrix = np.ones(shape=dims,dtype=complex) for index in range(len(RECO_rotate)): f = np.arange(dims[index]) - - 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']: - phase_vector = np.exp(1j*2*np.pi*(1-RECO_rotate[index])*f) - else: - raise ValueError('Your RECO_ft_mode is not supported') - + phase_vector = np.exp(1j*2*np.pi*(1-RECO_rotate[index])*f) if index == 0: - phase_matrix *= np.tile(phase_vector[:,np.newaxis,np.newaxis,np.newaxis], [1, dims[1], dims[2], dims[3]]) + phase_matrix *= np.tile(phase_vector[:,np.newaxis,np.newaxis], [1, dims[1], dims[2]]) elif index == 1: - phase_matrix *= np.tile(phase_vector[np.newaxis,:,np.newaxis,np.newaxis], [dims[0], 1, dims[2], dims[3]]) + phase_matrix *= np.tile(phase_vector[np.newaxis,:,np.newaxis], [dims[0], 1, dims[2]]) 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') + tmp = np.zeros((1,1,dims[2]), dtype=complex) + tmp[0,0,:] = phase_vector + phase_matrix *= np.tile(tmp, [dims[0], dims[1], 1]) + return phase_matrix + +# Replace with zero padding +def zero_filling(frame, RECO_ft_size, signal_position=np.array([0.5,0.5,0.5])): # 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'))]) - + not_Equal = any([(i != j) for i,j in zip(frame.shape,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]) + dims = (frame.shape[0], frame.shape[1], frame.shape[2]) # start process + newframe = np.zeros(RECO_ft_size, dtype=complex) + startpos = np.zeros(len(RECO_ft_size), dtype=int) + pos_ges = [None] * 3 - # 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 - + 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]] = 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.ifftn(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) - +def phase_corr(frame): + checkerboard = np.ones(shape=frame.shape[:3]) # 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): - 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') - - # 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 + checkerboard[::2,::2,::2] = -1 + checkerboard[1::2,1::2,::2] = -1 + checkerboard[::2,1::2,1::2] = -1 + checkerboard[1::2,::2,1::2] = -1 + checkerboard *= -1 + return checkerboard \ No newline at end of file diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py index 62580e1..43cab16 100644 --- a/brkraw/lib/recon.py +++ b/brkraw/lib/recon.py @@ -1,531 +1,222 @@ # -*- 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) + DEVELOPED FOR BRUKER PARAVISION 6, 7, and 360 datasets + Below are methods implemented for FID sorting into KSPACE + and simple image reconstruction. + MOST sequences are assumed to be cartesian + + This package was primarily created as an + intermediate for other image reconstruction + software such as MIRT and BART + +# NOTES +# EPI is not fully working yet +# T1_IG_FLASH_flc has an issue with size +# NEED to test compress sense for RARE """ -from .utils import get_value, set_value -from .recoFunctions import * +from .recoFunctions import phase_rotate, phase_corr, zero_filling +from ..api.brkobj.scan import ScanObj import numpy as np -from copy import deepcopy - - -def recon(fid_binary, acqp, meth, reco, process = 'image', 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 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'): - - # 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 +import warnings + +SUPPORTED_PROTOCOLS = ['rare','localizer' ,'gre', 'msme', + 'mge','dess', 'fisp', 'flash'] + +class Reconstruction: + def __init__(self, scanobj:'ScanObj', reco_id:'int'=1) -> None: + self.acqp = scanobj.acqp + self.method = scanobj.method + self.fid = scanobj.get_fid() + self.CS = True if self.method.get('PVM_EncCS')=='Yes' else False + self.NI = self.acqp['NI'] + self.NR = self.acqp['NR'] + self.NRecs = 1 + self.reco_id = reco_id + self.info = scanobj.get_info(self.reco_id) + self.protocol = self.info.protocol + self.reco = scanobj.get_reco(self.reco_id).reco + self.supported_protocol = any([True for i in SUPPORTED_PROTOCOLS if i in self.protocol['protocol_name'].lower()]) + + # 1) Convert Buffer to a np array + def sort_fid(self): + """ Sorts FID into a 3D np matrix [num_readouts, channel, scan_size] + + Returns + ------- + X : np.array [num_lines, channel, scan_size] + """ + # META DATA + dt_code = 'int32' + if self.acqp.get('ACQ_ScanPipeJobSettings') != None: + if self.acqp['ACQ_ScanPipeJobSettings'][0][1] == 'STORE_64bit_float': + dt_code = 'float64' + + bits = 64 if '64' in dt_code else 32 + DT_CODE = np.dtype(dt_code) + + BYTORDA = self.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(self.fid.read(), DT_CODE) + # Check Version and Sort fid data + if '360' in self.protocol['sw_version']: + # METAdata for 360 + self.NRecs = self.acqp['ACQ_ReceiverSelectPerChan'].count('Yes') + scanSize = self.acqp['ACQ_jobs'][0][0] + X = fid[::2] + 1j*fid[1::2] - # Full Cartesian Pipeline - 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 [num_lines, channel, scan_size] - """ - from .reference import BYTEORDER, WORDTYPE - # META DATA - NI = get_value(acqp, 'NI') - NR = get_value(acqp, 'NR') - - 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) + else: + # METAdata Versions Before 360 + self.NRecs = self.acqp['ACQ_ReceiverSelect'].count('Yes') if self.acqp['ACQ_ReceiverSelect'] != None else 1 + ACQ_size = self.acqp['ACQ_size'] if isinstance(self.acqp['ACQ_size'],int) else self.acqp['ACQ_size'] + scanSize = ACQ_size[0] + if self.acqp['GO_block_size'] == 'Standard_KBlock_Format': + blocksize = int(np.ceil(ACQ_size[0]*self.NRecs*(bits/8)/1024)*1024/(bits/8)) + else: + blocksize = int(ACQ_size[0]*self.NRecs) - BYTORDA = get_value(acqp, 'BYTORDA') - if BYTORDA == 'little': - DT_CODE = DT_CODE.newbyteorder('<') - elif BYTORDA == 'big': - DT_CODE = DT_CODE.newbyteorder('>') + # CHECK SIZE + if fid.size != blocksize*np.prod(ACQ_size[1:])*self.self.NI*self.self.NR: + raise ValueError('Error FID size dont match') - # Get FID FROM buffer - fid = np.frombuffer(fid_binary, DT_CODE) + # Convert to Complex + X = fid[::2] + 1j*fid[1::2] + X = X.reshape([-1,blocksize//2]) - # 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] + # Reshape Matrix [num_lines, channel, scan_size] + if blocksize != scanSize*self.NRecs: + X = X[:,:scanSize//2*self.NRecs] # [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)) + X = X.reshape((-1, self.NRecs, scanSize//2)) + + return X + + # 2) Convert to KSPACE + def sort_kspace(self, fid = None): + """ + FID = [num_lines, channel, scan_size] + KSPACE = [kx,ky,kz,NRec,NI,NR] + """ + if fid == None: + fid = self.sort_fid() + if not self.supported_protocol: + warnings.warn("SEQUENCE PROTOCOL {} NOT SUPPORTED YET...\nreturning readout sorted".format(self.acqp.get('ACQ_scan_name' ))) + return fid + + N_total, _, Nreadout = fid.shape + + # Meta data + dims = self.acqp.get('ACQ_dim') + obj_order = [self.acqp.get('ACQ_obj_order')] if isinstance(self.acqp.get('ACQ_obj_order'), int) else self.acqp.get( 'ACQ_obj_order') + phase_factor = self.method.get('PVM_RareFactor') if 'rare' in self.method.get('Method').lower() else self.acqp.get('ACQ_phase_factor') + + PVM_Matrix = self.method.get('PVM_Matrix') + kSize = np.round(np.array(self.method.get('PVM_AntiAlias'))*np.array(PVM_Matrix)) + zerofill = 2*np.floor( (kSize - kSize/np.array(self.method.get('PVM_EncZf')))/2 ) + kSize = kSize - zerofill + center = np.floor(kSize/2) + + EncMatrix = self.method.get('PVM_EncMatrix') + NPE = self.method.get('PVM_EncGenTotalSteps') if self.CS else np.prod(EncMatrix[1:]) + phase_encode2 = [0] + if dims == 3: + phase_encode2 = self.method.get('PVM_EncSteps2') + phase_encode2 = (phase_encode2 + center[2]).astype(int) + phase_encode1 = self.method.get('PVM_EncSteps1') + phase_encode1 = (phase_encode1 + center[1]).astype(int) + + if self.method.get('PVM_IsEpiScan') == 'Yes': + fid = fid.reshape(N_total,self.NRecs,int(Nreadout/kSize[0]),int(kSize[0])).transpose(2,0,1,3) + fid = fid.reshape(-1,self.NRecs,int(kSize[0])) + fid[1::2,:,:] = fid[1::2,:,::-1] + Nreadout = int(kSize[0]) + readStart = int(kSize[0]-Nreadout) + + assert np.prod(fid.shape) == (Nreadout*NPE*self.NI*self.NRecs*self.NR), 'Method calculated size does not match size of fid' + + temp = np.zeros([int(kSize[0]), int(kSize[1]),int(kSize[2]) if dims == 3 else 1, self.NRecs, self.NI, self.NR], dtype=complex) + if self.CS: + warnings.warn('Compressed Sensing has only been tested on undersampled GRE sequences') + phase_index1 = (self.method.get('PVM_EncGenSteps1') + center[1]).astype(int) + phase_index2 = (self.method.get('PVM_EncGenSteps2') + center[2]).astype(int) + fid = fid.reshape((self.NR,NPE,self.NI,self.NRecs,Nreadout)).transpose(4,1,3,2,0) + for index, (i,j) in enumerate(zip(phase_index1,phase_index2)): + temp[readStart:,i,j,:,:,:] = fid[:,index,:,:,:] + fid = np.zeros_like(temp) + fid[:,:,:,:,obj_order,:] = temp 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 - - -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() - - # Metadata - 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] - - 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) - - 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] - - # Start Resort based on - if ACQ_dim>1: - # [num_readout, channel, scan_size] -> [channel, scan_size, num_readout] - 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: - # 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 Exception('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 understanding of the data structure - [x,y,z,_,n_channel,NI,NR] - """ - # Metadata - 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] - if scanSize == 0: - scanSize = get_value(acqp,'ACQ_size')[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)) + fid = fid.reshape((self.NR,-1,self.NI,phase_factor,self.NRecs,Nreadout)).transpose(0,2,4,1,3,5) + fid = fid.reshape((self.NR,self.NI,self.NRecs,NPE,Nreadout)).transpose((4,3,2,1,0)) + fid = fid.reshape(Nreadout, int(EncMatrix[1]), int(EncMatrix[2]) if dims == 3 else 1, self.NRecs, self.NI, self.NR, order = 'F') + temp[readStart:,phase_encode1,:,:,:,:] = fid[:,:,phase_encode2,:,:,:] + fid = np.zeros_like(temp) + fid[:,:,:,:,obj_order,:] = temp + + if self.method.get('EchoAcqMode') == 'allEchoes': + fid[:,:,:,:,1::2,:] = fid[::-1,:,:,:,1::2,:] + + return fid + + def process_kspace(self): + kspace = self.sort_kspace() + if len(kspace.shape) != 6: + return kspace + # Shift Object + map_index= np.reshape(np.arange(0,kspace.shape[4]*kspace.shape[5]), (kspace.shape[5], kspace.shape[4]) ).flatten() + for NR in range(self.NR): + for NI in range(self.NI): + kspace[:,:,:,:,NI,NR] *= np.tile(phase_rotate(kspace[:,:,:,:,NI,NR], + self.reco.get('RECO_rotate'), + map_index[(NI+1)*(NR+1)-1])[:,:,:,np.newaxis], + [1,1,1,self.NRecs]) + + # Zeropad KSPACE + newdata_dims=[1, 1, 1] + RECO_ft_size = self.reco.get('RECO_ft_size') + newdata_dims[0:len(RECO_ft_size)] = RECO_ft_size + newdata = np.zeros(shape=newdata_dims+[self.NRecs, self.NI, self.NR], dtype=complex) + for NR in range(self.NR): + for NI in range(self.NI): + for chan in range(self.NRecs): + newdata[:,:,:,chan,NI,NR] = zero_filling(kspace[:,:,:,chan,NI,NR], RECO_ft_size).reshape(newdata[:,:,:,chan,NI,NR].shape) + + return newdata - PVM_EncZf=get_value(meth, 'PVM_EncZf') - - if PVM_EncZf == None: - # No zero-filling/interpolation available. - PVM_EncZf = np.ones((ACQ_dim)) - - # 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)) - # 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 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) - 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'] - # 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] - 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]) - - 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 + # 4) CONVERT TO IMAGE SPACE if FULLY SAMPLED CARTESIAN + def reconstruct(self, kspace=None, rms=True): + if kspace == None: + kspace = self.process_kspace() + if len(kspace.shape) != 6: + return kspace # sorted fid + if self.CS: + return kspace # zero padded kspace + + # Always FT and correct Phase + image = np.fft.fftshift(np.fft.ifftn(kspace, axes=(0,1,2))) + image *= np.tile(phase_corr(image)[:,:,:,np.newaxis,np.newaxis,np.newaxis], + [1,1,1,self.NRecs,self.NI,self.NR]) + if rms: + image = np.sqrt(np.mean(np.square(np.abs(image)), axis=3)) + return image + +def reconstruction(scanobj,process='image', **kwargs): + # Ensure Scans are Image Based + ACQ_dim_desc = [scanobj.acqp.get('ACQ_dim_desc')] if isinstance(scanobj.acqp.get('ACQ_dim_desc'), str) else scanobj.acqp.get('ACQ_dim_desc') + if 'Spectroscopic' in ACQ_dim_desc: + warnings.warn('Scan is spectroscopic') + process = 'readout' + + # Reconstruction Processing + recoObj = Reconstruction(scanobj) + if process == 'readout': + return recoObj.sort_fid() + elif process == 'kspace': + return recoObj.process_kspace() + return recoObj.reconstruct(rms=kwargs['rms'] if 'rms' in kwargs.keys() else True) \ No newline at end of file From 2bf7f750b763a5e3e946283578bb2e87bf304e91 Mon Sep 17 00:00:00 2001 From: Tim Date: Sat, 20 Apr 2024 23:05:00 -0400 Subject: [PATCH 3/9] fix size check --- brkraw/lib/recon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py index 43cab16..3ea18db 100644 --- a/brkraw/lib/recon.py +++ b/brkraw/lib/recon.py @@ -82,7 +82,7 @@ def sort_fid(self): blocksize = int(ACQ_size[0]*self.NRecs) # CHECK SIZE - if fid.size != blocksize*np.prod(ACQ_size[1:])*self.self.NI*self.self.NR: + if fid.size != blocksize*np.prod(ACQ_size[1:])*self.NI*self.NR: raise ValueError('Error FID size dont match') # Convert to Complex From 64d5555f9fe208817efc7fd1469cd26f98e90b18 Mon Sep 17 00:00:00 2001 From: Tim Date: Sat, 20 Apr 2024 23:16:11 -0400 Subject: [PATCH 4/9] parameter guards --- brkraw/lib/recon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py index 3ea18db..8247ec6 100644 --- a/brkraw/lib/recon.py +++ b/brkraw/lib/recon.py @@ -73,7 +73,7 @@ def sort_fid(self): else: # METAdata Versions Before 360 - self.NRecs = self.acqp['ACQ_ReceiverSelect'].count('Yes') if self.acqp['ACQ_ReceiverSelect'] != None else 1 + self.NRecs = self.acqp.get('ACQ_ReceiverSelect').count('Yes') if self.acqp.get('ACQ_ReceiverSelect') != None else 1 ACQ_size = self.acqp['ACQ_size'] if isinstance(self.acqp['ACQ_size'],int) else self.acqp['ACQ_size'] scanSize = ACQ_size[0] if self.acqp['GO_block_size'] == 'Standard_KBlock_Format': From 75dd1b8e5f4d24d9aa047cf73b18d2ba2a66e6e9 Mon Sep 17 00:00:00 2001 From: Tim Date: Sun, 21 Apr 2024 00:05:07 -0400 Subject: [PATCH 5/9] limit fftshift --- brkraw/lib/recon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py index 8247ec6..144f72c 100644 --- a/brkraw/lib/recon.py +++ b/brkraw/lib/recon.py @@ -199,7 +199,7 @@ def reconstruct(self, kspace=None, rms=True): return kspace # zero padded kspace # Always FT and correct Phase - image = np.fft.fftshift(np.fft.ifftn(kspace, axes=(0,1,2))) + image = np.fft.fftshift(np.fft.ifftn(kspace, axes=(0,1,2)), axes=(0,1,2)) image *= np.tile(phase_corr(image)[:,:,:,np.newaxis,np.newaxis,np.newaxis], [1,1,1,self.NRecs,self.NI,self.NR]) if rms: From 872834f0862afaafb42e5f84cf8c9999a8c19b8c Mon Sep 17 00:00:00 2001 From: Tim Date: Sun, 21 Apr 2024 13:15:27 -0400 Subject: [PATCH 6/9] sigpy not necessary --- brkraw/lib/recoSigpy.py | 52 ----------------------------------------- 1 file changed, 52 deletions(-) delete mode 100644 brkraw/lib/recoSigpy.py diff --git a/brkraw/lib/recoSigpy.py b/brkraw/lib/recoSigpy.py deleted file mode 100644 index ac7fada..0000000 --- a/brkraw/lib/recoSigpy.py +++ /dev/null @@ -1,52 +0,0 @@ -""" -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 From 3ca27534e867f7f5d649793d6791840732b4aaea Mon Sep 17 00:00:00 2001 From: Tim Date: Mon, 22 Apr 2024 00:51:54 -0400 Subject: [PATCH 7/9] minor edits and add recon test file --- brkraw/lib/recon.py | 35 ++++++++++++++++--------------- tests/recon_api_test.py | 46 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 17 deletions(-) create mode 100644 tests/recon_api_test.py diff --git a/brkraw/lib/recon.py b/brkraw/lib/recon.py index 144f72c..9efec11 100644 --- a/brkraw/lib/recon.py +++ b/brkraw/lib/recon.py @@ -24,6 +24,21 @@ SUPPORTED_PROTOCOLS = ['rare','localizer' ,'gre', 'msme', 'mge','dess', 'fisp', 'flash'] +def reconstruction(scanobj,process='image', **kwargs): + # Ensure Scans are Image Based + ACQ_dim_desc = [scanobj.acqp.get('ACQ_dim_desc')] if isinstance(scanobj.acqp.get('ACQ_dim_desc'), str) else scanobj.acqp.get('ACQ_dim_desc') + if 'Spectroscopic' in ACQ_dim_desc: + warnings.warn('Scan is spectroscopic') + process = 'readout' + + # Reconstruction Processing + recoObj = Reconstruction(scanobj) + if process == 'readout': + return recoObj.sort_fid() + elif process == 'kspace': + return recoObj.process_kspace() + return recoObj.reconstruct(rms=kwargs['rms'] if 'rms' in kwargs.keys() else True) + class Reconstruction: def __init__(self, scanobj:'ScanObj', reco_id:'int'=1) -> None: self.acqp = scanobj.acqp @@ -177,9 +192,9 @@ def process_kspace(self): map_index[(NI+1)*(NR+1)-1])[:,:,:,np.newaxis], [1,1,1,self.NRecs]) - # Zeropad KSPACE - newdata_dims=[1, 1, 1] + # Zeropad KSPACE RECO_ft_size = self.reco.get('RECO_ft_size') + newdata_dims=[1, 1, 1] newdata_dims[0:len(RECO_ft_size)] = RECO_ft_size newdata = np.zeros(shape=newdata_dims+[self.NRecs, self.NI, self.NR], dtype=complex) for NR in range(self.NR): @@ -205,18 +220,4 @@ def reconstruct(self, kspace=None, rms=True): if rms: image = np.sqrt(np.mean(np.square(np.abs(image)), axis=3)) return image - -def reconstruction(scanobj,process='image', **kwargs): - # Ensure Scans are Image Based - ACQ_dim_desc = [scanobj.acqp.get('ACQ_dim_desc')] if isinstance(scanobj.acqp.get('ACQ_dim_desc'), str) else scanobj.acqp.get('ACQ_dim_desc') - if 'Spectroscopic' in ACQ_dim_desc: - warnings.warn('Scan is spectroscopic') - process = 'readout' - - # Reconstruction Processing - recoObj = Reconstruction(scanobj) - if process == 'readout': - return recoObj.sort_fid() - elif process == 'kspace': - return recoObj.process_kspace() - return recoObj.reconstruct(rms=kwargs['rms'] if 'rms' in kwargs.keys() else True) \ No newline at end of file + \ No newline at end of file diff --git a/tests/recon_api_test.py b/tests/recon_api_test.py new file mode 100644 index 0000000..169cde4 --- /dev/null +++ b/tests/recon_api_test.py @@ -0,0 +1,46 @@ +from brkraw.app import tonifti as tonii +from brkraw.lib.recon import reconstruction + +import numpy as np +import matplotlib.pyplot as plt + +import nibabel as nib + +fiddata_path = 'YOUR.PvDatasets' +studyobj = tonii.brkraw.BrkrawToNifti(fiddata_path) +scan_id = 10 # REPLACE WITH YOUR SCAN NUMBER +reco_id = 1 +scanobj = studyobj.get_scan(scan_id) + +reconobj = np.abs(np.squeeze(reconstruction(scanobj, rms=True))) +affine = studyobj.get_affine(scan_id, reco_id) +dataobj = studyobj.get_dataobj(scan_id, reco_id) + +print(reconobj.shape) +plt.subplot(1,2,1) +if len(reconobj.shape) == 2: + plt.imshow(reconobj[:,:].T) + plt.title(f"{scan_id},{scanobj.acqp['ACQ_scan_name']}") +elif len(reconobj.shape) == 3: + plt.imshow(reconobj[:,:,reconobj.shape[2]//2].T) + plt.title(f"{scan_id},{scanobj.acqp['ACQ_scan_name']}") +elif len(reconobj.shape) == 4: + plt.imshow(reconobj[:,:,reconobj.shape[2]//2,0].T) + plt.title(f"{scan_id},{scanobj.acqp['ACQ_scan_name']}") +plt.subplot(1,2,2) +if len(reconobj.shape) == 2: + plt.imshow(dataobj[:,:].T) + plt.title(f"{scan_id},dataobj") +elif len(reconobj.shape) == 3: + plt.imshow(dataobj[:,:,dataobj.shape[2]//2].T) + plt.title(f"{scan_id},dataobj") +elif len(reconobj.shape) == 4: + plt.imshow(dataobj[:,:,dataobj.shape[2]//2,0].T) + plt.title(f"{scan_id},dataobj") +plt.show() + +print(reconobj.shape, dataobj.shape) +assert np.prod(dataobj.shape) == np.prod(reconobj.shape), "Shape mismatched" + +niiobj = nib.Nifti1Image(reconobj/np.max(np.abs(reconobj)), affine) +niiobj.to_filename('reconfile.nii.gz') From 136e9bcfeeef16d4be58708b4788c87da5d766f2 Mon Sep 17 00:00:00 2001 From: Tim Date: Mon, 22 Apr 2024 01:36:09 -0400 Subject: [PATCH 8/9] brkrecon removed --- brkraw/scripts/brkrecon.py | 159 ------------------------------------- pyproject.toml | 3 +- 2 files changed, 1 insertion(+), 161 deletions(-) delete mode 100644 brkraw/scripts/brkrecon.py diff --git a/brkraw/scripts/brkrecon.py b/brkraw/scripts/brkrecon.py deleted file mode 100644 index 9105a1f..0000000 --- a/brkraw/scripts/brkrecon.py +++ /dev/null @@ -1,159 +0,0 @@ -""" -@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" - 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' - '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') - - 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 == '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: - 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) - - 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(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/pyproject.toml b/pyproject.toml index c1a68cb..1b1ece0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,5 +46,4 @@ style = "pep440" [project.scripts] brkraw = "brkraw.scripts.brkraw:main" -brk-backup = 'brkraw.scripts.brk_backup:main' -brkrecon = "brkraw.scripts.brkrecon:main" \ No newline at end of file +brk-backup = 'brkraw.scripts.brk_backup:main' \ No newline at end of file From 6157d31f8f711f7add27d6354538c3eada586548 Mon Sep 17 00:00:00 2001 From: Tim Date: Mon, 22 Apr 2024 11:51:33 -0400 Subject: [PATCH 9/9] minor change in zero fill --- brkraw/lib/recoFunctions.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/brkraw/lib/recoFunctions.py b/brkraw/lib/recoFunctions.py index 654291c..96f6c22 100644 --- a/brkraw/lib/recoFunctions.py +++ b/brkraw/lib/recoFunctions.py @@ -5,8 +5,8 @@ in recon.py """ -from .utils import get_value import numpy as np +import warnings def phase_rotate(frame, RECO_rotate, framenumber): @@ -41,7 +41,8 @@ def zero_filling(frame, RECO_ft_size, signal_position=np.array([0.5,0.5,0.5])): not_Equal = any([(i != j) for i,j in zip(frame.shape,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') + warnings.warn('Signal needs to be between 0 and 1\nDefaulting to 0.5') + signal_position=np.array([0.5,0.5,0.5]) # calculate additional variables dims = (frame.shape[0], frame.shape[1], frame.shape[2])