update python to go along with matlab changes in 9beee53

ecat_testing/read_matlab_nii.py

@@ -0,0 +1,42 @@
+import nibabel
+import pathlib
+import numpy
+
+steps_dir = pathlib.Path(__file__).parent / 'steps'
+
+def first_middle_last_frames_to_text(four_d_array_like_object, output_folder, step_name='_step_name_'):
+    frames = [0, four_d_array_like_object.shape[-1] // 2, four_d_array_like_object.shape[-1] - 1]
+    frames_to_record = []
+    for f in frames:
+        frames_to_record.append(four_d_array_like_object[:, :, :, f])
+
+    # now collect a single 2d slice from the "middle" of the 3d frames in frames_to_record
+    for index, frame in enumerate(frames_to_record):
+        numpy.savetxt(output_folder / f"{step_name}_frame_{frames[index]}.tsv",
+                      frames_to_record[index][:, :, frames_to_record[index].shape[2] // 2],
+                      delimiter="\t", fmt='%s')
+
+path_to_matlab_nii = pathlib.Path('matlab_nii.nii')
+path_to_python_nii = pathlib.Path('python_nii.nii.gz')
+
+
+python_nii = nibabel.load(str(path_to_python_nii))
+matlab_nii = nibabel.load(str(path_to_matlab_nii))
+
+python_data = python_nii.get_fdata()
+matlab_data = matlab_nii.get_fdata()
+
+# compare the two arrays in each
+print(numpy.allclose(python_data, matlab_data, rtol=0.5))
+
+
+# subtract the two arrays
+diff = python_data - matlab_data
+print(f"difference max and min: {diff.max()}, {diff.min()}")
+print(f"mean difference: {diff.mean()}")
+
+# save diff as nii
+diff_nii = python_nii.__class__(diff, python_nii.affine, python_nii.header)
+nibabel.save(diff_nii, steps_dir / '12_diff_between_written_niis.nii.gz')
+
+first_middle_last_frames_to_text(diff, steps_dir, step_name='13_diff_between_written_niis')

pypet2bids/pypet2bids/ecat2nii.py

@@ -46,7 +46,7 @@ def ecat2nii(ecat_main_header=None,
     """
     Converts an ECAT file into a nifti and a sidecar json, used in conjunction with
     :func:`pypet2bids.read_ecat.read_ecat`
-
+    
     :param ecat_main_header: the main header of an ECAT file
     :param ecat_subheaders: the subheaders for each frame of the ECAT file
     :param ecat_pixel_data: the imaging/pixel data from the ECAT file
@@ -75,33 +75,18 @@ def ecat2nii(ecat_main_header=None,
     nifti_file_w_out_extension = os.path.splitext(str(pathlib.Path(nifti_file).name))[0]
 
     # if already read nifti file skip re-reading
-    if (
-        ecat_main_header is None
-        and ecat_subheaders is None
-        and ecat_pixel_data is None
-        and ecat_file
-    ):
+    if ecat_main_header is None and ecat_subheaders is None and ecat_pixel_data is None and ecat_file:
         # collect ecat_file
         main_header, sub_headers, data = read_ecat(ecat_file=ecat_file)
-    elif (
-        ecat_file is None
-        and type(ecat_main_header) is dict
-        and type(ecat_subheaders) is list
-        and type(ecat_pixel_data) is numpy.ndarray
-    ):
-        main_header, sub_headers, data = (
-            ecat_main_header,
-            ecat_subheaders,
-            ecat_pixel_data,
-        )
+    elif ecat_file is None and type(ecat_main_header) is dict and type(ecat_subheaders) is list and type(
+            ecat_pixel_data) is numpy.ndarray:
+        main_header, sub_headers, data = ecat_main_header, ecat_subheaders, ecat_pixel_data
     else:
-        raise Exception(
-            "Must pass in filepath for ECAT file or "
-            "(ecat_main_header, ecat_subheaders, and ecat_pixel data "
-            f"got ecat_file={ecat_file}, type(ecat_main_header)={type(ecat_main_header)}, "
-            f"type(ecat_subheaders)={type(ecat_subheaders)}, "
-            f"type(ecat_pixel_data)={type(ecat_pixel_data)} instead."
-        )
+        raise Exception("Must pass in filepath for ECAT file or "
+                        "(ecat_main_header, ecat_subheaders, and ecat_pixel data "
+                        f"got ecat_file={ecat_file}, type(ecat_main_header)={type(ecat_main_header)}, "
+                        f"type(ecat_subheaders)={type(ecat_subheaders)}, "
+                        f"type(ecat_pixel_data)={type(ecat_pixel_data)} instead.")
 
     # debug step #6 view data as passed to ecat2nii method
     if ecat_save_steps == '1':
@@ -128,20 +113,18 @@ def ecat2nii(ecat_main_header=None,
     pixel_data_type = data.dtype
 
     # check for TimeZero supplied via kwargs
-    if kwargs.get("TimeZero", None):
-        TimeZero = kwargs["TimeZero"]
+    if kwargs.get('TimeZero', None):
+        TimeZero = kwargs['TimeZero']
     else:
         logger.warning("Metadata TimeZero is missing -- set to ScanStart or empty to use the scanning time as "
               "injection time")
 
     # get image shape
     img_shape = data.shape
-    shape_from_headers = (
-        sub_headers[0]["X_DIMENSION"],
-        sub_headers[0]["Y_DIMENSION"],
-        sub_headers[0]["Z_DIMENSION"],
-        main_header["NUM_FRAMES"],
-    )
+    shape_from_headers = (sub_headers[0]['X_DIMENSION'],
+                          sub_headers[0]['Y_DIMENSION'],
+                          sub_headers[0]['Z_DIMENSION'],
+                          main_header['NUM_FRAMES'])
 
     # make sure number of data elements matches frame number
     single_frame = False
@@ -150,8 +133,7 @@ def ecat2nii(ecat_main_header=None,
     if img_shape != shape_from_headers and not single_frame:
         raise Exception(
             f"Mis-match between expected X,Y,Z, and Num. Frames dimensions ({shape_from_headers} obtained from headers"
-            f"and shape of imaging data ({img_shape}"
-        )
+            f"and shape of imaging data ({img_shape}")
 
     # format data into acceptable shape for nibabel, by first creating empty matrix
     img_temp = numpy.zeros(shape=(sub_headers[0]['X_DIMENSION'],
@@ -167,27 +149,17 @@ def ecat2nii(ecat_main_header=None,
     prompts, randoms = [], []
 
     # load frame data into img temp
-    for index in reversed(
-        range(img_shape[3])
-    ):  # Don't throw stones working from existing matlab code
+    for index in reversed(range(img_shape[3])):  # Don't throw stones working from existing matlab code
         print(f"Loading frame {index + 1}")
         img_temp[:, :, :, index] = numpy.flip(numpy.flip(numpy.flip(
             data[:, :, :, index] * sub_headers[index]['SCALE_FACTOR'], 1), 2), 0)
         start.append(sub_headers[index]['FRAME_START_TIME'] * 60)  # scale to per minute
         delta.append(sub_headers[index]['FRAME_DURATION'] * 60)  # scale to per minute
 
-        if main_header.get("SW_VERSION", 0) >= 73:
+        if main_header.get('SW_VERSION', 0) >= 73:
             # scale both to per minute
-            prompts.append(
-                sub_headers[index]["PROMPT_RATE"]
-                * sub_headers[index]["FRAME_DURATION"]
-                * 60
-            )
-            randoms.append(
-                sub_headers[index]["RANDOM_RATE"]
-                * sub_headers[index]["FRAME_DURATION"]
-                * 60
-            )
+            prompts.append(sub_headers[index]['PROMPT_RATE'] * sub_headers[index]['FRAME_DURATION'] * 60)
+            randoms.append(sub_headers[index]['RANDOM_RATE'] * sub_headers[index]['FRAME_DURATION'] * 60)
         else:
             # this field is not available in ecat 7.2
             prompts.append(0)
@@ -242,29 +214,20 @@ def ecat2nii(ecat_main_header=None,
         final_image = numpy.round(final_image) * sca
 
     qoffset_x = -1 * (
-        (
-            (sub_headers[0]["X_DIMENSION"] * sub_headers[0]["X_PIXEL_SIZE"] * 10 / 2)
-            - sub_headers[0]["X_PIXEL_SIZE"] * 5
-        )
-    )
+        ((sub_headers[0]['X_DIMENSION'] * sub_headers[0]['X_PIXEL_SIZE'] * 10 / 2) - sub_headers[0][
+            'X_PIXEL_SIZE'] * 5))
 
     qoffset_y = -1 * (
-        (
-            (sub_headers[0]["Y_DIMENSION"] * sub_headers[0]["Y_PIXEL_SIZE"] * 10 / 2)
-            - sub_headers[0]["Y_PIXEL_SIZE"] * 5
-        )
-    )
+        ((sub_headers[0]['Y_DIMENSION'] * sub_headers[0]['Y_PIXEL_SIZE'] * 10 / 2) - sub_headers[0][
+            'Y_PIXEL_SIZE'] * 5))
 
     qoffset_z = -1 * (
-        (
-            (sub_headers[0]["Z_DIMENSION"] * sub_headers[0]["Z_PIXEL_SIZE"] * 10 / 2)
-            - sub_headers[0]["Z_PIXEL_SIZE"] * 5
-        )
-    )
+        ((sub_headers[0]['Z_DIMENSION'] * sub_headers[0]['Z_PIXEL_SIZE'] * 10 / 2) - sub_headers[0][
+            'Z_PIXEL_SIZE'] * 5))
 
     # build affine if it's not included in function call
     if not affine:
-        mat = np.diag([sub_headers[0]['X_PIXEL_SIZE'],
+        mat = numpy.diag([sub_headers[0]['X_PIXEL_SIZE'],
                        sub_headers[0]['Y_PIXEL_SIZE'],
                        sub_headers[0]['Z_PIXEL_SIZE']]) * 10
         affine = nibabel.affines.from_matvec(mat, [qoffset_x, qoffset_y, qoffset_z])