changed order of gibbs unringing

Author: Ian-Goodall-Halliwell

micaflow/cli.py

@@ -797,7 +797,8 @@ def format_help(self):
         "--output", required=True, help="Output path for the denoised image (.nii.gz)"
     )
     denoise_parser.add_argument("--b0-denoise", action='store_true', help="Denoise b0 volumes separately (default: False)")
-
+    denoise_parser.add_argument("--gibbs", action='store_true', help="Apply Gibbs ringing correction (default: False)")
+    denoise_parser.add_argument("--threads", type=int, help="Number of threads to use (default: 1)")
 
     # Motion Correction command
     motion_corr_parser = subparsers.add_parser(

micaflow/resources/Snakefile

@@ -537,14 +537,16 @@ if RUN_DWI:
             bvec = BVEC_FILE
         output:
             denoised = f"{TEMP_DIR}/{SUBJECT}_{SESSION}_denoised_DWI.nii.gz"
-        threads: LIGHT_THREADS
+        threads: HEAVY_THREADS
         shell:
             """
             micaflow denoise \
                 --input {input.moving} \
                 --bval {input.bval} \
                 --bvec {input.bvec} \
-                --output {output.denoised}
+                --output {output.denoised} \
+                --gibbs \
+                --threads {threads}
             """
     rule dwi_b0_extraction:
         input:
@@ -606,7 +608,7 @@ if RUN_DWI:
             output:
                 corrected = f"{TEMP_DIR}/{SUBJECT}_{SESSION}_denoised_bias-corrected_DWI.nii.gz",
                 b0_corrected = f"{TEMP_DIR}/{SUBJECT}_{SESSION}_biascorrected-b0.nii.gz"
-            threads: HEAVY_THREADS
+            threads: LIGHT_THREADS
             shell:
                 """
                 micaflow bias_correction \
@@ -615,7 +617,6 @@ if RUN_DWI:
                     --b0-output {output.b0_corrected} \
                     --output {output.corrected} \
                     --shell-dimension {SHELL_DIMENSION} \
-                    --gibbs \
                     --threads {threads}
                 """
         rule b0_synthseg:
@@ -945,7 +946,7 @@ if RUN_DWI:
             output:
                 corrected = f"{TEMP_DIR}/{SUBJECT}_{SESSION}_denoised_bias-corrected_DWI.nii.gz",
                 b0_corrected = f"{TEMP_DIR}/{SUBJECT}_{SESSION}_biascorrected-b0.nii.gz"
-            threads: HEAVY_THREADS
+            threads: LIGHT_THREADS
             shell:
                 """
                 micaflow bias_correction \
@@ -955,7 +956,6 @@ if RUN_DWI:
                     --output {output.corrected} \
                     --mask {input.mask} \
                     --shell-dimension {SHELL_DIMENSION} \
-                    --gibbs \
                     --threads {threads}
                 """
 

micaflow/scripts/bias_correction.py

@@ -435,6 +435,24 @@ def print_help_message():
       {YELLOW}--output{RESET} FLAIR_corrected.nii.gz \\
       {YELLOW}--mode{RESET} 3d
     
+    {BLUE}# Example 5: Custom shell dimension{RESET}
+    micaflow bias_correction \\
+      {YELLOW}--input{RESET} DWI.nii.gz \\
+      {YELLOW}--output{RESET} DWI_corrected.nii.gz \\
+      {YELLOW}--b0{RESET} b0.nii.gz \\
+      {YELLOW}--b0-output{RESET} b0_corrected.nii.gz \\
+      {YELLOW}--shell-dimension{RESET} 4
+    
+    {BLUE}# Example 6: With Gibbs ringing removal{RESET}
+    micaflow bias_correction \\
+      {YELLOW}--input{RESET} DWI.nii.gz \\
+      {YELLOW}--output{RESET} DWI_corrected.nii.gz \\
+      {YELLOW}--mask{RESET} brain_mask.nii.gz \\
+      {YELLOW}--b0{RESET} b0.nii.gz \\
+      {YELLOW}--b0-output{RESET} b0_corrected.nii.gz \\
+      {YELLOW}--mode{RESET} 4d \\
+      {YELLOW}--gibbs{RESET}
+    
     {CYAN}{BOLD}────────── WHAT IS BIAS FIELD? ────────────────────────{RESET}
     
     {GREEN}Bias field (intensity non-uniformity):{RESET}
@@ -628,7 +646,7 @@ def bias_field_correction_3d(image_path, output_path, mask_path=None, gibbs=Fals
 
 
 def bias_field_correction_4d(image_path, mask_path=None, output_path=None, 
-                             b0_path=None, b0_corrected_path=None, shell_dimension=3, gibbs=False, threads=1):
+                             b0_path=None, b0_corrected_path=None, shell_dimension=3):
     """
     Apply N4 bias field correction to a 4D diffusion image.
     
@@ -654,8 +672,6 @@ def bias_field_correction_4d(image_path, mask_path=None, output_path=None,
     shell_dimension : int, default=3
         Dimension along which diffusion volumes are organized (0-indexed).
         For standard NIfTI: dimension 3 (4th dimension).
-    gibbs : bool, optional
-        If True, apply Gibbs ringing removal before bias correction.
     
     Returns
     -------
@@ -721,24 +737,13 @@ def bias_field_correction_4d(image_path, mask_path=None, output_path=None,
     print(f"{CYAN}Loading 4D diffusion image...{RESET}")
     img = ants.image_read(image_path)
 
-    if gibbs:
-        print(f"{CYAN}Running Gibbs ringing removal on 4D data...{RESET}")
-        arr = img.numpy()
-        gibbs_removal(arr, slice_axis=2, n_points=3, inplace=True, num_processes=threads)
-        img = img.new_image_like(arr)
-
     img_data = img.numpy()
     print(f"  Image shape: {img_data.shape}")
     print(f"  Shell dimension: {shell_dimension}")
 
     if b0_path:
         print(f"{CYAN}Loading b=0 image...{RESET}")
         b0_img = ants.image_read(b0_path)
-        if gibbs:
-            print(f"{CYAN}Running Gibbs ringing removal on b=0 image...{RESET}")
-            arr_b0 = b0_img.numpy()
-            gibbs_removal(arr_b0, slice_axis=2, n_points=3, inplace=True, num_processes=threads)
-            b0_img = b0_img.new_image_like(arr_b0)
     else:
         b0_img = None
 
@@ -1054,7 +1059,7 @@ def run_bias_field_correction(image_path, output_path, mask_path=None, mode="aut
         if mode == "4d":
             return bias_field_correction_4d(
                 image_path, mask_path, output_path, 
-                b0_path, b0_corrected_path, shell_dimension, gibbs, threads
+                b0_path, b0_corrected_path, shell_dimension
             )
         else:  # 3d
             return bias_field_correction_3d(image_path, output_path, mask_path, gibbs, threads)

micaflow/scripts/denoise.py

@@ -46,6 +46,14 @@
     --output <path/to/denoised_dwi.nii.gz> \\
     --b0-denoise
 
+# With Gibbs ringing removal
+micaflow denoise \\
+    --input <path/to/dwi.nii.gz> \\
+    --bval <path/to/dwi.bval> \\
+    --bvec <path/to/dwi.bvec> \\
+    --output <path/to/denoised_dwi.nii.gz> \\
+    --gibbs
+
 Python API Usage:
 ----------------
 >>> from micaflow.scripts.denoise import run_denoise
@@ -67,6 +75,15 @@
 ...     output="denoised_dwi.nii.gz",
 ...     b0_denoising=True
 ... )
+>>>
+>>> # With Gibbs ringing removal
+>>> output = run_denoise(
+...     moving="raw_dwi.nii.gz",
+...     moving_bval="dwi.bval", 
+...     moving_bvec="dwi.bvec",
+...     output="denoised_dwi.nii.gz",
+...     gibbs=True
+... )
 
 Pipeline Integration:
 --------------------
@@ -118,6 +135,7 @@
 import os
 import numpy as np
 from dipy.denoise.patch2self import patch2self
+from dipy.denoise.gibbs import gibbs_removal
 from dipy.io.gradients import read_bvals_bvecs
 from colorama import init, Fore, Style
 
@@ -158,6 +176,8 @@ def print_help_message():
     {CYAN}{BOLD}─────────────────── OPTIONAL ARGUMENTS ───────────────────{RESET}
       {YELLOW}--b0-denoise{RESET}: Denoise b0 volumes separately (default: False)
                    {MAGENTA}Experimental - not recommended for most cases{RESET}
+      {YELLOW}--gibbs{RESET}     : Apply Gibbs ringing removal after denoising
+      {YELLOW}--threads{RESET}   : Number of threads for Gibbs removal (default: 1)
     
     {CYAN}{BOLD}──────────────────── EXAMPLE USAGE ──────────────────────{RESET}
     
@@ -176,6 +196,14 @@ def print_help_message():
       {YELLOW}--output{RESET} denoised_dwi.nii.gz \\
       {YELLOW}--b0-denoise{RESET}
     
+    {BLUE}# With Gibbs ringing removal{RESET}
+    micaflow denoise \\
+      {YELLOW}--input{RESET} raw_dwi.nii.gz \\
+      {YELLOW}--bval{RESET} dwi.bval \\
+      {YELLOW}--bvec{RESET} dwi.bvec \\
+      {YELLOW}--output{RESET} denoised_dwi.nii.gz \\
+      {YELLOW}--gibbs{RESET}
+    
     {CYAN}{BOLD}────────────────── PATCH2SELF ALGORITHM ──────────────────{RESET}
     
     {GREEN}How it works:{RESET}
@@ -226,7 +254,7 @@ def print_help_message():
     print(help_text)
 
 
-def run_denoise(moving, moving_bval, moving_bvec, output, b0_denoising=False):
+def run_denoise(moving, moving_bval, moving_bvec, output, b0_denoising=False, gibbs=False, threads=1):
     """
     Denoise diffusion-weighted images using the Patch2Self algorithm.
     
@@ -251,6 +279,10 @@ def run_denoise(moving, moving_bval, moving_bvec, output, b0_denoising=False):
     b0_denoising : bool, optional
         If True, denoise b0 volumes separately. Default: False.
         Standard practice is to exclude b0 volumes from denoising.
+    gibbs : bool, optional
+        If True, apply Gibbs ringing removal after denoising. Default: False.
+    threads : int, optional
+        Number of threads to use for Gibbs ringing removal. Default: 1.
         
     Returns
     -------
@@ -293,6 +325,15 @@ def run_denoise(moving, moving_bval, moving_bvec, output, b0_denoising=False):
     ...     output="denoised_dwi.nii.gz",
     ...     b0_denoising=True
     ... )
+    >>> 
+    >>> # With Gibbs ringing removal
+    >>> output = run_denoise(
+    ...     moving="raw_dwi.nii.gz",
+    ...     moving_bval="dwi.bval",
+    ...     moving_bvec="dwi.bvec",
+    ...     output="denoised_dwi.nii.gz",
+    ...     gibbs=True
+    ... )
     """
     # Validate input files exist
     for filepath, name in [(moving, "Input DWI"), 
@@ -334,7 +375,18 @@ def run_denoise(moving, moving_bval, moving_bvec, output, b0_denoising=False):
         b0_threshold=50,
         b0_denoising=b0_denoising,
     )
+    env = os.environ.copy()
+    env["OPENBLAS_NUM_THREADS"] = "1"
+    env["OMP_NUM_THREADS"] = "1"
+    env["MKL_NUM_THREADS"] = "1"
+    env["VECLIB_MAXIMUM_THREADS"] = "1"
+    env["NUMEXPR_NUM_THREADS"] = "1"
+    os.environ.update(env)
 
+    if gibbs:
+        print(f"\n{CYAN}Applying Gibbs ringing removal...{RESET}")
+        gibbs_removal(denoised, slice_axis=2, n_points=3, inplace=True, num_processes=threads)
+
     print(f"\n{CYAN}Saving denoised image...{RESET}")
     nib.save(nib.Nifti1Image(denoised, moving_image.affine), output)
 
@@ -380,6 +432,17 @@ def run_denoise(moving, moving_bval, moving_bvec, output, b0_denoising=False):
         action='store_true', 
         help="Denoise b0 volumes separately (experimental, default: False)."
     )
+    parser.add_argument(
+        "--gibbs", 
+        action='store_true', 
+        help="Apply Gibbs ringing removal after denoising."
+    )
+    parser.add_argument(
+        "--threads", 
+        type=int, 
+        default=1, 
+        help="Number of threads for Gibbs removal (default: 1)."
+    )
 
     args = parser.parse_args()
 
@@ -389,7 +452,9 @@ def run_denoise(moving, moving_bval, moving_bvec, output, b0_denoising=False):
             args.bval, 
             args.bvec, 
             args.output, 
-            args.b0_denoise
+            args.b0_denoise,
+            args.gibbs,
+            args.threads
         )
 
         print(f"\n{GREEN}{BOLD}Denoising successfully completed!{RESET}")