|
| 1 | +import h5py |
| 2 | +import numpy as np |
| 3 | +import xml.etree.ElementTree as etree |
| 4 | +from typing import Dict, List, Optional, Sequence, Tuple |
| 5 | + |
| 6 | + |
| 7 | +def load_file_T2(fname: str) -> Tuple: |
| 8 | + """ |
| 9 | + Load T2 fastmri file. |
| 10 | + |
| 11 | + Parameters: |
| 12 | + ----------- |
| 13 | + fname : str |
| 14 | + Path to the h5 fastmri file. |
| 15 | + |
| 16 | + Returns: |
| 17 | + -------- |
| 18 | + Tuple |
| 19 | + A tuple containing the kspace, calibration_data, hdr, im_recon, and attributes of the file. |
| 20 | + """ |
| 21 | + |
| 22 | + with h5py.File(fname, "r") as hf: |
| 23 | + kspace = hf["kspace"][:] |
| 24 | + calibration_data = hf["calibration_data"][:] |
| 25 | + hdr = hf["ismrmrd_header"][()] |
| 26 | + im_recon = hf["reconstruction_rss"][:] |
| 27 | + atts = dict() |
| 28 | + atts['max'] = hf.attrs['max'] |
| 29 | + atts['norm'] = hf.attrs['norm'] |
| 30 | + atts['patient_id'] = hf.attrs['patient_id'] |
| 31 | + atts['acquisition'] = hf.attrs['acquisition'] |
| 32 | + |
| 33 | + return kspace, calibration_data, hdr, im_recon, atts |
| 34 | + |
| 35 | + |
| 36 | +def load_file_dwi(fname: str) -> Tuple: |
| 37 | + """ |
| 38 | + Load DWI fastmri file. |
| 39 | + |
| 40 | + Parameters: |
| 41 | + ----------- |
| 42 | + fname : str |
| 43 | + Path to the h5 fastmri file. |
| 44 | + |
| 45 | + Returns: |
| 46 | + -------- |
| 47 | + Tuple |
| 48 | + A tuple containing the kspace, calibration_data, hdr, and coil sensitivity maps. |
| 49 | + """ |
| 50 | + |
| 51 | + with h5py.File(fname, 'r') as f: |
| 52 | + kspace = f['kspace'][:] |
| 53 | + calibration = f['calibration_data'][:] |
| 54 | + coil_sens_maps = f['coil_sens_maps'][:] |
| 55 | + #phase_corr = f['phase_correction'][:] |
| 56 | + |
| 57 | + ismrmrd_header = f['ismrmrd_header'][()] |
| 58 | + hdr = get_regridding_params(ismrmrd_header) |
| 59 | + |
| 60 | + return kspace, calibration, coil_sens_maps, hdr |
| 61 | + |
| 62 | + |
| 63 | +def get_padding(hdr: str) -> float: |
| 64 | + """ |
| 65 | + Extract the padding value from an XML header string. |
| 66 | +
|
| 67 | + Parameters: |
| 68 | + ----------- |
| 69 | + hdr : str |
| 70 | + The XML header string. |
| 71 | +
|
| 72 | + Returns: |
| 73 | + -------- |
| 74 | + float |
| 75 | + The padding value calculated as (x - max_enc)/2, where x is the readout dimension and |
| 76 | + max_enc is the maximum phase-encoding dimension. |
| 77 | + """ |
| 78 | + et_root = etree.fromstring(hdr) |
| 79 | + lims = ["encoding", "encodingLimits", "kspace_encoding_step_1"] |
| 80 | + enc_limits_max = int(et_query(et_root, lims + ["maximum"])) + 1 |
| 81 | + enc = ["encoding", "encodedSpace", "matrixSize"] |
| 82 | + enc_x = int(et_query(et_root, enc + ["x"])) |
| 83 | + padding = (enc_x - enc_limits_max)/2 |
| 84 | + |
| 85 | + return padding |
| 86 | + |
| 87 | + |
| 88 | +def et_query(root: etree.Element, qlist: Sequence[str], namespace: str = "http://www.ismrm.org/ISMRMRD") -> str: |
| 89 | + """ |
| 90 | + ElementTree query function. |
| 91 | + |
| 92 | + This function queries an XML document using ElementTree. |
| 93 | + |
| 94 | + Parameters: |
| 95 | + ----------- |
| 96 | + root : Element |
| 97 | + Root of the XML document to search through. |
| 98 | + qlist : Sequence of str |
| 99 | + A sequence of strings for nested searches, e.g., ["Encoding", "matrixSize"]. |
| 100 | + namespace : str, optional |
| 101 | + XML namespace to prepend query. |
| 102 | + |
| 103 | + Returns: |
| 104 | + -------- |
| 105 | + str |
| 106 | + The retrieved data as a string. |
| 107 | + """ |
| 108 | + s = "." |
| 109 | + prefix = "ismrmrd_namespace" |
| 110 | + |
| 111 | + ns = {prefix: namespace} |
| 112 | + |
| 113 | + for el in qlist: |
| 114 | + s = s + f"//{prefix}:{el}" |
| 115 | + |
| 116 | + value = root.find(s, ns) |
| 117 | + if value is None: |
| 118 | + raise RuntimeError("Element not found") |
| 119 | + |
| 120 | + return str(value.text) |
| 121 | + |
| 122 | + |
| 123 | +def zero_pad_kspace_hdr(hdr: str, unpadded_kspace: np.ndarray) -> np.ndarray: |
| 124 | + """ |
| 125 | + Perform zero-padding on k-space data to have the same number of |
| 126 | + points in the x- and y-directions. |
| 127 | +
|
| 128 | + Parameters |
| 129 | + ---------- |
| 130 | + hdr : str |
| 131 | + The XML header string. |
| 132 | + unpadded_kspace : array-like of shape (sl, ro , coils, pe) |
| 133 | + The k-space data to be padded. |
| 134 | +
|
| 135 | + Returns |
| 136 | + ------- |
| 137 | + padded_kspace : ndarray of shape (sl, ro_padded, coils, pe_padded) |
| 138 | + The zero-padded k-space data, where ro_padded and pe_padded are |
| 139 | + the dimensions of the readout and phase-encoding directions after |
| 140 | + padding. |
| 141 | +
|
| 142 | + Notes |
| 143 | + ----- |
| 144 | + The padding value is calculated using the `get_padding` function, which |
| 145 | + extracts the padding value from the XML header string. If the difference |
| 146 | + between the readout dimension and the maximum phase-encoding dimension |
| 147 | + is not divisible by 2, the padding is applied asymmetrically, with one |
| 148 | + side having an additional zero-padding. |
| 149 | +
|
| 150 | + """ |
| 151 | + padding = get_padding(hdr) |
| 152 | + if padding%2 != 0: |
| 153 | + padding_left = int(np.floor(padding)) |
| 154 | + padding_right = int(np.ceil(padding)) |
| 155 | + else: |
| 156 | + padding_left = int(padding) |
| 157 | + padding_right = int(padding) |
| 158 | + padded_kspace = np.pad(unpadded_kspace, ((0,0),(0,0),(0,0), (padding_left,padding_right))) |
| 159 | + |
| 160 | + return padded_kspace |
| 161 | + |
| 162 | + |
| 163 | +def get_regridding_params(hdr: str) -> Dict: |
| 164 | + """ |
| 165 | + Extracts regridding parameters from header XML string. |
| 166 | +
|
| 167 | + Parameters |
| 168 | + ---------- |
| 169 | + hdr : str |
| 170 | + Header XML string. |
| 171 | +
|
| 172 | + Returns |
| 173 | + ------- |
| 174 | + dict |
| 175 | + A dictionary containing the extracted parameters. |
| 176 | +
|
| 177 | + """ |
| 178 | + res = { |
| 179 | + 'rampUpTime': None, |
| 180 | + 'rampDownTime': None, |
| 181 | + 'flatTopTime': None, |
| 182 | + 'acqDelayTime': None, |
| 183 | + 'echoSpacing': None |
| 184 | + } |
| 185 | + |
| 186 | + et_root = etree.fromstring(hdr) |
| 187 | + namespace = {'ns': "http://www.ismrm.org/ISMRMRD"} |
| 188 | + |
| 189 | + for node in et_root.findall('ns:encoding/ns:trajectoryDescription/ns:userParameterLong', namespace): |
| 190 | + if node[0].text in res.keys(): |
| 191 | + res[node[0].text] = float(node[1].text) |
| 192 | + |
| 193 | + return res |
| 194 | + |
| 195 | + |
| 196 | +def save_recon(outp_dict: Dict[str, any], output_path: str) -> None: |
| 197 | + """ |
| 198 | + Save reconstruction results to an HDF5 file. |
| 199 | +
|
| 200 | + Parameters |
| 201 | + ---------- |
| 202 | + outp_dict : dict |
| 203 | + A dictionary containing the reconstructed images, with the image names as keys. |
| 204 | + output_path : str |
| 205 | + The file path to save the reconstructed images. |
| 206 | +
|
| 207 | + Returns |
| 208 | + ------- |
| 209 | + None |
| 210 | + """ |
| 211 | + |
| 212 | + hf = h5py.File(output_path, "w") |
| 213 | + for key, outp in outp_dict.items(): |
| 214 | + hf.create_dataset(key, data=outp) |
| 215 | + hf.close() |
0 commit comments