Add results postprocessing

Author: youennmereljourdan

postprocess.py

@@ -0,0 +1,47 @@
+import glob
+import os
+from argparse import ArgumentParser
+
+import nibabel as nib
+
+from postprocess.postprocess_service import PostprocessService
+
+
+def postproc():
+    postproc_srv = PostprocessService()
+
+    parser = ArgumentParser(description='Post processing of results')
+    parser.add_argument('--results', required=True, type=str, help='path to results')
+    args = parser.parse_args()
+
+    basedir = args.results
+    mean_path = os.path.join(basedir, 'mean_result.nii')
+    ds_path = os.path.join(basedir, 'dataset.csv')
+
+    ids = []
+    results = []
+    paths = glob.glob(os.path.join(basedir, '*/'), recursive=True)
+    for path in paths:
+        ids.append(os.path.basename(os.path.dirname(path)))
+        results.append(os.path.join(path, '_subject_id_01', 'result.nii'))
+
+    size = len(results)
+    if size < 1:
+        print(f'No results found in [{basedir}]')
+        return
+    else:
+        print(f'[{size}] results found in [{basedir}]')
+
+    print(f"Generating mean result image from [{size}] results...")
+    mean_nifti_image = postproc_srv.get_mean_image(results, 10)
+    nib.save(mean_nifti_image, mean_path)
+    print(f"Mean result image written to [{mean_path}]")
+
+    print(f"Building dataset from [{size}] results...")
+    dataframe = postproc_srv.get_dataframe(basedir, ids)
+    dataframe.to_csv(ds_path, index=False, sep=';')
+    print(f"Dataset CSV written to [{mean_path}]")
+
+
+if __name__ == '__main__':
+    postproc()

postprocess.sh

@@ -0,0 +1,16 @@
+#!/usr/bin/env bash
+
+#OAR -l walltime=2
+#OAR -O ./log/postprocess_log_%jobid%.stdout
+#OAR -E ./log/postprocess_log_%jobid%.stderr
+#OAR -q production
+
+TAG="fmri-confs-runner"
+
+BASE="/home/ymerel/empenn_group_storage/private/ymerel"
+RESULTS="$BASE/results/auditory"
+
+g5k-setup-docker -t
+docker pull ghcr.io/inria-empenn/fmri-confs-runner:latest
+docker run -u root -v "$RESULTS:/results" $TAG python postprocess.py --results /results
+

postprocess/__init__.py

@@ -0,0 +1,79 @@
+import numpy as np
+import pandas as pd
+from nibabel import load, Nifti1Image
+from nibabel.processing import resample_from_to
+from numpy import corrcoef, reshape
+from pandas import DataFrame
+from scipy.stats import spearmanr
+from typing import List
+
+
+class CorrelationService:
+    def compute_correlations(self, image, images: List[str]) -> DataFrame:
+        print(f"Computing correlations from [{image}] to [{len(images)}] images... ", end='')
+        dfs = []
+        for img in images:
+            corr = self.get_correlation_coefficient(image, img, 'spearman')
+            dfs.append(pd.DataFrame([[image, img, corr]], columns=['source', 'target', 'correlation']))
+        merged = pd.concat(dfs, ignore_index=True)
+        print("OK")
+        return merged.sort_values(by='correlation', ascending=False)
+
+    def get_correlation_coefficient(self,
+            file_1: str, file_2: str, method: str = 'pearson') -> float:
+        """ Return the correlation coefficient of two images.
+
+            Arguments :
+                - file_1, str - path to the first image
+                - file_2, str - path to the second image ; file_2 will be resampled on file_1
+                - method, str - either 'pearson', or 'spearman': the correlation method to use
+                - reslice_on_file_2, bool - set to :
+                    - True if you wish to reslice file_1 on file_2
+                    - False otherwise
+
+            Returns :
+                - _, float - the correlation coefficient of the two input images,
+                using the passed method
+        """
+
+        # Load images
+        image_1 = load(file_1)
+        image_2 = load(file_2)
+
+        # Set masking using NaN's
+        image_1 = self.mask_using_zeros(image_1)
+        image_2 = self.mask_using_zeros(image_2)
+
+        # Resample using nearest neighbours
+        image_2 = resample_from_to(image_2, image_1, order=0)
+
+        # Make 1D vectors from the images data
+        data_1 = reshape(image_1.get_fdata(), -1)
+        data_2 = reshape(image_2.get_fdata(), -1)
+
+        # Compute the correlation coefficient
+        if method == 'pearson':
+            return corrcoef(data_1, data_2)[0][1]
+        if method == 'spearman':
+            return spearmanr(data_1, data_2).correlation
+
+        raise AttributeError(f'Wrong correlation method provided: {method}.')
+
+    def mask_using_zeros(self, data_image: Nifti1Image) -> Nifti1Image:
+        """ Mask an image by replacing NaNs with zeros.
+
+            Arguments:
+                - data_image, nibabel.Nifti1Image : the image to mask
+
+            Returns:
+                - the masked image as nibabel.Nifti1Image
+        """
+
+        # Get data from the image
+        data = data_image.get_fdata()
+
+        # Replace NaNs by zeros
+        data[np.isnan(data)] = 0.0
+
+        # Return data as an image
+        return Nifti1Image(data, data_image.affine)

postprocess/correlation_service.py

@@ -0,0 +1,57 @@
+import os
+
+import nibabel as nib
+import numpy as np
+import pandas as pd
+from typing import List
+
+from postprocess.correlation_service import CorrelationService
+
+
+class PostprocessService:
+    corr_srv = CorrelationService()
+
+    def get_dataframe(self, path, ids: List[str]) -> pd.DataFrame:
+        dataframes = []
+        ref_contrast = os.path.join(path, 'ref', '_subject_id_01', 'result.nii')
+        mean = os.path.join(path, 'mean_result.nii')
+        for conf_id in ids:
+            contrast = os.path.join(path, conf_id, '_subject_id_01', 'result.nii')
+            config = os.path.join(path, conf_id, 'config.csv')
+            df = pd.read_csv(config, delimiter=';').astype(bool)
+            df['id'] = conf_id
+            df['from_ref'] = self.corr_srv.get_correlation_coefficient(contrast, ref_contrast, 'spearman')
+            df['from_mean'] = self.corr_srv.get_correlation_coefficient(contrast, mean, 'spearman')
+            dataframes.append(df)
+
+        return pd.concat(dataframes, ignore_index=True)
+
+    def get_mean_image(self, inputs: list, batch_size: int) -> nib.Nifti1Image:
+        total_sum = None
+        count = 0
+
+        for i in range(0, len(inputs), batch_size):
+            batch_paths = inputs[i:i + batch_size]
+            batch_images = [nib.load(path).get_fdata() for path in batch_paths]
+
+            # Stack the batch images into a single numpy array
+            batch_array = np.stack(batch_images, axis=0)
+
+            # Calculate the sum of the batch
+            batch_sum = np.sum(batch_array, axis=0)
+
+            # Accumulate the sum and count
+            if total_sum is None:
+                total_sum = batch_sum
+            else:
+                total_sum += batch_sum
+
+            count += len(batch_paths)
+
+        # Calculate the mean image
+        mean_image = total_sum / count
+
+        # Create a new NIfTI image with the mean data
+        mean_nifti = nib.Nifti1Image(mean_image, affine=nib.load(inputs[0]).affine)
+
+        return mean_nifti

postprocess/postprocess_service.py

@@ -2,8 +2,8 @@
 
 #OAR --array 10
 #OAR -l walltime=4
-#OAR -O ./log/log_%jobid%.stdout
-#OAR -E ./log/log_%jobid%.stderr
+#OAR -O ./log/run_config_log_%jobid%.stdout
+#OAR -E ./log/run_config_log_%jobid%.stderr
 #OAR -q production
 
 TAG="fmri-confs-runner"
@@ -25,7 +25,7 @@ fi
 echo "Running configuration is [$CONF]"
 
 g5k-setup-docker -t
-docker build . -t $TAG
+docker pull ghcr.io/inria-empenn/fmri-confs-runner:latest
 if [ "$OAR_ARRAY_INDEX" -eq 1 ]; then
     # write ref config only for the first job
     docker run -u root -v "$DATA:/data" -v "$RESULTS:/results" -v "$WORK:/work" -v "$CONFIGS:/configs" $TAG python run.py --configs "/configs/$CONF" --data /data/data_desc.json --ref /configs/config_ref.csv