Merge branch 'main' of https://github.com/openneuropet/PET2BIDS

.github/workflows/phantoms.yaml

@@ -67,6 +67,18 @@ jobs:
             brew install dcm2niix
           fi
 
+      - name: Debug
+        uses: mxschmitt/action-tmate@v3
+        if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled }}
+        timeout-minutes: 15
+        with:
+          limit-access-to-actor: true
+
+      - name: Check dcm2niix is installed and on path
+        if: runner.os != 'Windows'
+        run:
+          dcm2niix -h 
+
       - name: Install dcm2niix windows
         if: runner.os == 'Windows'
         run: |

.gitignore

@@ -58,6 +58,12 @@ D:\\BIDS*
 *OpenNeuroPET-Phantoms*
 **PHANTOMS.zip
 ecat_validation/ECAT7_multiframe.v
+variables.env
+ecat_testing/matlab/
+ecat_testing/steps/
+ecat_testing/**/*.nii*
+*osem*
+*fbp*
+
+matlab/SiemensHRRTparameters.txt
 
-ecat_testing/
-SiemensHRRTparameters.txt

CITATION.cff

@@ -31,11 +31,17 @@ authors:
 - family-names: "Ganz-Benjaminsen"
   given-names: "Melanie"
   orcid: "https://orcid.org/0000-0002-9120-8098"
+- family-names: "Bilgel"
+  given-names: "Murat"
+  orcid: "https://orcid.org/0000-0001-5042-7422"
+- family-names: "Eierud"
+  given-names: "Cyrus"
+  orcid: "https://orcid.org/0000-0002-9942-676X"
 - family-names: "Pernet"
   given-names: "Cyril"
   orcid: "https://orcid.org/0000-0003-4010-4632"
 title: "PET2BIDS"
-version: 1.0.0
+version: 1.3.20240502
 url: "https://github.com/openneuropet/PET2BIDS"
 date-released: 2022-10-01
 license: MIT

JOSS_metadata/persons_info.txt

@@ -21,3 +21,5 @@
 0000-0002-9120-8098  Melanie Ganz-Benjaminse  Department of Computer Science, University of Copenhagen, Copenhagen, Denmark
 
 0000-0003-4010-4632  Cyril Pernet  Neurobiology Research Unit, Rigshospitalet, Copenhagen, Denmark
+
+0000-0001-5042-7422  Murat Bilgel  National Institute on Aging Intramural Research Program, Baltimore, MD, USA

contributors.md

@@ -3,7 +3,7 @@
 The following individuals have contributed to the PET2BIDS project (in alphabetical order). If you contributed and
 some icons needs to be added or your name is not listed, please add it.
 
-Murat Bilgel 💻 🐛 💡
+Murat Bilgel 💻 🐛 💡  
 Anthony Galassi 💻 📖 💬 🎨 💡 ⚠️  
 Melanie Ganz-Benjaminsen 🔍 💬 🤔 📋  
 Gabriel Gonzalez-Escamilla 💻 ⚠️ 🐛 👀  

ecat_testing/ints/README.md

@@ -0,0 +1,60 @@
+# matlab_numpy_ints
+
+This is a very simple example of testing matlab vs python reading and writing of ecat files scaling up to a more
+complex set of testing that reads real ecat data and writes it to niftis.
+
+## Dead Simple Testing / Proofing 
+
+Then simplest test cases that use "imaging" data made of small 2 x 1 arrays of integer data can found in:
+
+`ints.m` and `ints.py`
+
+## Outline of Testing on Real Data (e.g. steps required to read and write ecat into Nifti)
+
+A part of this debugging process is to capture the state of inputs and outputs through each iteration of the python
+and matlab code. These outputs are cataloged and saved in this directory (`ecat_testing/ints/`) within subdirectory 
+`steps/`. The naming, order, and description of each step are laid out in the table below for easy reference. At each
+of these steps the outputs are saved in the `steps/` directory with the filestem (and the python or matlab code used to 
+make them) as described. These temporary files are written if and only if the environment variable `ECAT_SAVE_STEPS=1` 
+is set.
+
+| Step | Description                                           | Matlab        | Python         | FileStem                     |
+|------|-------------------------------------------------------|---------------|----------------|------------------------------|
+| 1    | Read main header                                      | `read_ECAT7.m` | `read_ecat.py` | `1_read_mh_ecat*`            |
+| 2    | Read subheaders                                       | `read_EACT7.m` | `read_ecat.py` | `2_read_sh_ecat*`            |
+| 3    | Determine file/data type                              | `read_ECAT7.m` | `read_ecat.py` | `3_determine_data_type*`     |
+| 4    | Read image data                                       | `read_ECAT7.m` | `read_ecat.py` | `4_read_img_ecat*`           |
+| 5    | scale if calibrated                                   | `read_ECAT7.m` | `read_ecat.py` | `5_scale_img_ecat*`          |
+| 6    | Pass Data to ecat2nii                                 | `ecat2nii.m`  | `ecat2nii.py`  | `6_ecat2nii*`                |
+| 7    | Flip ECAT data into Nifti space                       | `ecat2nii.m`  | `ecat2nii.py`  | `7_flip_ecat2nii*`           |
+| 8    | Rescale to 16 bits                                    | `ecat2nii.m`  | `ecat2nii.py`  | `8_rescale_to_16_ecat2nii*`  |
+| 9    | Calibration Units Scaling                             | `ecat2nii.m`  | `ecat2nii.py`  | `9_scal_cal_units_ecat2nii*` |
+| 10   | Save to Nifti                                         | `ecat2nii.m`  | `ecat2nii.py`  | `10_save_nii_ecat2nii*`      |
+| 11   | Check to see if values are recorded as they should be | `ecat2nii.m`   | `ecat2nii.py`   | `11_read_saved_nii*`          |
+
+
+1. Read main header: this is the first step in reading the ecat file, the main header contains information about the
+   file, the subheaders, and the image data. These will be saved as jsons for comparison.
+2. Read subheaders: the subheaders contain information about the image data, these will be saved as jsons for comparison
+   as well.
+3. Determine file/data type: this step is to determine the type of data in the image data, this will be saved a json or 
+   a text file with a single line specifying the datatype and endianness, but probably a json.
+4. Read image data: this is the final step in reading the ecat file, the image data is read in and will be examined at 
+   3 different time points (if available). E.g. the first frame, the middle frame, and the final frame. Only a single 2D
+   slice will be saved from each of the time points, and it too will be taken from the "middle" of its 3D volume. We're
+   only attempting to compare whether python and matlab have done a decent job of reading the data as recorded in.
+5. Repeat step 4 but scale the data if it should be scaled
+6. Save objects for comparison (as best as one can) before they are passed to the ecat2nii function. This will include
+   the mainheader, subheaders, and image data.
+7. Return the transformed data to the nifti space from ecat. This follows the 3 flip dimension steps performed across 
+   the 3D image data. This output will use the same frames as step 4 and 5.
+8. Rescale the data to 16 bits: this should only occur if the data is 12bit as is sometimes the case with ecat data. As 
+   a note to self, attempting these steps in Python will start to lead to wildly different values when compared to 
+   matlab. It's most likely not necessary to do this step as the data is handled in numpy, but this writer won't promise
+   to eat his hat if it turns out to be necessary.
+9. Calibration Units: (Search for 'calibration_units == 1' to locate this in code) Here we can potentially alter the 
+   data again by scaling, rounding, and converting from int to float. As in steps 4 and 5 we will save the first, 
+   middle, and last frames of the data as 2D slices in the middle of their respective 3D volumes for comparison.
+10. Save to nifti: the data is saved to a nifti file and the output is saved for comparison as .nii files named in line
+   with the FileStem column above.
+11. Additional Steps: TBD

ecat_testing/ints/check_ecat_read.py

@@ -0,0 +1,14 @@
+import numpy
+import nibabel
+import dotenv
+import os
+
+
+
+# load the environment variables
+dotenv.load_dotenv('variables.env')
+
+# get the path to the input ecat file from the environment variables
+ecat_file = os.getenv('WELL_BEHAVED_ECAT_FILE')
+wustl_fbp_ecat=os.getenv('WUSTL_FBP_ECAT')
+wustl_osem_ecat=os.getenv('WUSTL_OSEM_ECAT')

ecat_testing/ints/ints.m

@@ -0,0 +1,55 @@
+env = loadenv('variables.env');
+ecat_file=env('WELL_BEHAVED_ECAT_FILE');
+wustl_fbp_ecat=env('WUSTL_FBP_ECAT');
+wustl_osem_ecat=env('WUSTL_OSEM_ECAT');
+
+% convert the following ints to bytes and write them to a file
+ints_to_write = [-32768, 32767]
+
+
+% use the same format as the python version of this code
+destination_file = 'matlab_sample_bytes.bytes'
+
+% write all the ints as bytes to the file
+fid = fopen(destination_file, 'w');
+disp("writing ints_to_write to sample file");
+disp(ints_to_write);
+fwrite(fid, ints_to_write, 'int16', 'b');
+%for i = 1:length(ints_to_write)
+%    fwrite(fid, typecast(ints_to_write(i), 'int16'), 'int16', 'b');
+%end
+fclose(fid);
+
+% read out the bytes as written to the sample file
+fid = fopen('matlab_sample_bytes.bytes', 'r');
+bytes = fread(fid);
+disp("bytes written from matlab in matlab_sample_bytes.byte and read with fread (no arguments)");
+disp(bytes);
+fclose(fid);
+
+% read in bytes from python, at least I have a good idea of what that is written as
+fid = fopen('python_sample_bytes.bytes', 'r');
+python_bytes = fread(fid);
+disp("bytes written from python in python_sample_bytes.byte and read with fread (no arguments)");
+disp(python_bytes);
+
+various_data_types = {'int16=>int16', 'int16', 'uint16'};
+
+disp("now we open the matlab file and read the bytes using multiple arguments for fread");
+
+% read in the bytes as int16
+
+calibration_factor = 0.7203709074867248;
+
+for t = 1:length(various_data_types)
+    % oddly enough varying the second argument to fread doesn't seem to change the output
+    fid = fopen(destination_file, 'r', 'ieee-be');
+    various_bytes = fread(fid, various_data_types{t});
+    disp(['datatype used for reading in bytes: ', various_data_types{t}]);
+    disp(various_bytes);
+    disp(various_data_types{t});
+    fclose(fid);
+
+    % scale the data to see what happens
+    scaled_bytes = calibration_factor * various_bytes
+end

ecat_testing/ints/ints.py

@@ -0,0 +1,97 @@
+import numpy
+import nibabel
+import dotenv
+import os
+
+# load the environment variables
+dotenv.load_dotenv('variables.env')
+
+# get the path to the input ecat file from the environment variables
+ecat_file = os.getenv('WELL_BEHAVED_ECAT_FILE')
+wustl_fbp_ecat=os.getenv('WUSTL_FBP_ECAT')
+wustl_osem_ecat=os.getenv('WUSTL_OSEM_ECAT')
+
+
+def create_sample_bytes(destination_file='python_sample_bytes.bytes'):
+    """
+    Create a file with two 16-bit integers, one positive and one negative.
+
+    :param destination_file: _description_, defaults to 'sample_bytes.bytes'
+    :type destination_file: str, optional
+    """
+    bytes_to_write = [-32768, 32767]
+    with open(destination_file, 'wb') as file:
+        byte_order = 'big'
+        signed = True
+        for b in bytes_to_write:
+            print(f"Writing {b} to {destination_file}, byte order: {byte_order}, signed: {signed}")
+            file.write(b.to_bytes(2, byteorder=byte_order, signed=True))
+
+
+def read_bytes(file_name='python_sample_bytes.bytes'):
+    """
+    Open a byte file and read in the bytes without any extra fluff.
+
+    :param file_name: the file to read bytes from, defaults to 'sample_bytes.bytes'
+    :type file_name: str, optional
+    :return: the bytes read from the file
+    :rtype: bytes
+    """
+    with open(file_name, 'rb') as file:
+        bytes_read = file.read()
+    print(f"Read these bytes from {file_name}: {bytes_read}")
+    return bytes_read
+
+
+# create the bytes to be read
+create_sample_bytes()
+# read the bytes
+read_nums = read_bytes()
+
+# these are some of the datatypes we wish to test with numpy
+# i2 = integer 16 signed, >i2 = integer 16 big endian signed, H = integer 16 unsigned
+dtypes = ['i2', '>i2', 'H']
+
+# create a dictionary to hold the arrays we create with numpy
+arrays_by_dtype = {}
+
+# iterate through the datatypes and create arrays with numpy
+for d in dtypes:
+    numpy_int16 = numpy.frombuffer(read_nums, dtype=d)
+    print(f"Reading bytes with numpy given datatype: {d}\nArray is listing this as : {numpy_int16.dtype} {numpy_int16}")
+    arrays_by_dtype[d] = numpy_int16
+
+print(f"Arrays by dtype: {arrays_by_dtype}")
+
+# next we go through the same steps that the ecat converter does but in miniature since a 1x2 array is easier for us
+# hairless apes to comprehend than a near as enough to make no difference N dimensional array
+
+# scale it the calibration factor we have been dealing with in these wustl ecats is 0.7203709074867248
+calibration_factor = 0.7203709074867248
+
+scaled_arrays = {}
+
+for k, v in arrays_by_dtype.items():
+    # try recasting the  scaled array after the multiplication
+    scaled_arrays[k] = v * calibration_factor
+
+print(f"These are the arrays after being scaled by {calibration_factor}: {scaled_arrays}")
+
+# write these out to nifti's
+for k, v in scaled_arrays.items():
+    nifti = nibabel.Nifti1Image(v, affine=numpy.eye(4))
+    input_data_type = v.dtype
+    nibabel_data_type = nifti.get_data_dtype()
+    print(f"Input data type: {input_data_type}, Nibabel data type: {nibabel_data_type}")
+    nibabel.save(nifti, f"nibabel_{k}.nii.gz")
+    print(f"Saved array to numpy_{k}.nii.gz")
+    print(f"loading that array results in the following: {nibabel.load(f'nibabel_{k}.nii.gz').get_fdata()}")
+
+# what happens if we don't scale the arrays and write them to nifti
+for k, v in arrays_by_dtype.items():
+    input_data_type = v.dtype
+    nibabel_data_type = nifti.get_data_dtype()
+    print(f"Input data type: {input_data_type}, Nibabel data type: {nibabel_data_type}")
+    nibabel.save(nifti, f"nibabel_{k}_unscaled.nii.gz")
+    print(f"Saved array to numpy_{k}_unscaled.nii.gz")
+    print(f"loading that array results in the following: {nibabel.load(f'nibabel_{k}_unscaled.nii.gz').get_fdata()}")