ENH: field prep

nmriprep/argprep/argprep.py

@@ -1,11 +1,12 @@
 import nibabel as nb
 import numpy as np
 
-from ..image import convert_nef_to_grey, read_tiff, save_slice
+from ..image import convert_nef_to_grey, save_slice
 from ..parser import get_argprep_parser
 from ..plotting import plot_curve, plot_mosaic, plot_single_slice
 from ..utils import find_files, inverse_rodbard, rodbard
 from .calibration import calibrate_standard
+from .fieldprep import find_fields
 
 
 def main():
@@ -14,22 +15,6 @@ def main():
     out_dir = src_dir.parent / 'preproc' if not args.output else args.output
     verbose = args.save_intermediate
 
-    # attempt to find flat field info
-    flatfield_correction = {}
-    if args.dark_field and args.flat_field:
-        flatfield_correction['dark'] = read_tiff(args.dark_field)
-        flatfield_correction['flat'] = read_tiff(args.flat_field)
-    else:
-        # try searching for a preprocessed darkfield
-        latest_ff = find_files(out_dir.glob('*flatfield.tiff'))
-        latest_df = find_files(out_dir.glob('*darkfield.tiff'))
-        if latest_ff and latest_df:
-            flatfield_correction['dark'] = read_tiff(latest_df[-1])
-            flatfield_correction['flat'] = read_tiff(latest_ff[-1])
-        else:
-            print('Skipping flat field correction...')
-            flatfield_correction = None
-
     # identify subjects for pipeline
     all_subject_ids = set(
         [
@@ -55,6 +40,18 @@ def main():
         fig_dir = sub_dir / 'figures'
         fig_dir.mkdir(exist_ok=True)
 
+        # attempt to find flat field info
+        flatfield_correction = {}
+        flatfield_correction['dark'] = find_fields(
+            args.dark_field, sub_dir.glob('*darkfield.tif*')
+        )
+        flatfield_correction['flat'] = find_fields(
+            args.flat_field, sub_dir.glob('*flatfield.tif*')
+        )
+        if any(v is None for v in flatfield_correction.values()):
+            print('Skipping flat field correction...')
+            flatfield_correction = None
+
         # calibrate standards to transform GV to radioactivity
         popt, std_rad, std_gv, std_stem = calibrate_standard(
             sub_id,

nmriprep/argprep/fieldprep.py

@@ -0,0 +1,53 @@
+from collections import defaultdict
+from pathlib import Path
+
+import numpy as np
+
+from ..image import convert_nef_to_grey, read_tiff, save_slice
+from ..parser import get_fieldprep_parser
+from ..utils import find_files, parse_kv
+
+
+def find_fields(user=None, search_map=None):
+    fieldpath = user if user is not None else find_files(search_map)[-1]
+    field = read_tiff(fieldpath) if fieldpath else None
+    return field
+
+
+def fieldprep():
+    args = get_fieldprep_parser().parse_args()
+
+    if args.dark_field:
+        raise NotImplementedError('Dark field processing coming soon...')
+
+    if args.flat_field:
+        ff_dir = args.flat_field
+        fnames = find_files(ff_dir.rglob('*flatfield*.nef'))
+        if len(fnames) < 1:
+            raise FileNotFoundError(f'No flat field files found in {ff_dir.absolute()}')
+
+        subdirs = defaultdict(list)
+        [subdirs[p.parent].append(p) for p in fnames]
+
+        for subdir in subdirs.keys():
+            sub_files = subdirs[subdir]
+            fname_parts = parse_kv(sub_files[0].stem)
+            out_stem = '_'.join(
+                f'{k}-{v}' for k, v in fname_parts.items() if 'flatfield' not in k
+            )
+            out_dir = (
+                Path(str(subdir).replace('sourcedata', 'preproc')).resolve()
+                if not args.output
+                else args.output
+            )
+            out_dir.mkdir(parents=True, exist_ok=True)
+
+            data = np.median(
+                np.stack(
+                    [convert_nef_to_grey(fname) for fname in sub_files],
+                    axis=2,
+                ),
+                axis=2,
+            )
+            save_slice(data, out_dir / f'{out_stem}_flatfield.tif')
+    return

nmriprep/parser.py

@@ -48,6 +48,25 @@ def get_argprep_parser():
     return parser
 
 
+def get_fieldprep_parser():
+    """Build parser object."""
+
+    parser = ArgumentParser(
+        description='Convert dark/flatfield .nef images for processing with argprep',
+        formatter_class=ArgumentDefaultsHelpFormatter,
+    )
+    parser.add_argument(
+        '--flat-field', help='Path to directory containing flat field image', type=Path
+    )
+    parser.add_argument(
+        '--dark-field', help='Path to directory containing dark field image', type=Path
+    )
+    parser.add_argument(
+        '--output', help='Optionally specify output directory', type=Path
+    )
+    return parser
+
+
 def get_roiextract_parser():
     """Build parser object."""
 

pyproject.toml

@@ -36,6 +36,7 @@ Homepage = "https://github.com/brainkcl/nmriprep"
 
 [project.scripts]
 argprep = "nmriprep.argprep.argprep:main"
+fieldprep = "nmriprep.argprep.fieldprep:fieldprep"
 roi_extract = "nmriprep.measure:roi_extract"
 
 [tool.hatch.build]