divide the intensity of a set of medical images by the percentile of the dataset

normalize-percentile.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Normalize the intensity of a set of images by 
finding the median 99th percentile of all images
and dividing each image by that value
Author: Jacob Reinhold
"""

import os
import re
import sys
from argparse import ArgumentParser
from pathlib import Path

import numpy as np
import torch
import torchio as tio


def main() -> int:
    parser = ArgumentParser(description="Normalize images for image processing")
    parser.add_argument("top_dir", type=Path)
    parser.add_argument("out_dir", type=Path)
    parser.add_argument("-i", "--include", type=str, default=[".*"], nargs="+")
    parser.add_argument("-e", "--exclude", type=str, default=[None], nargs="+")
    parser.add_argument("-p", "--percentile", type=float, default=99.5)
    parser.add_argument("-v", "--verbose", action="store_true")
    parser.add_argument("-d", "--dry-run", action="store_true")
    args = parser.parse_args()
    include_progs = [re.compile(inc) for inc in args.include]
    exclude_progs = [re.compile(exc) for exc in args.exclude if exc is not None]
    percentiles = []
    for root, dirs, files in os.walk(args.top_dir):
        for _fn in files:
            if all(inc.search(_fn) is None for inc in include_progs) or any(
                exc.search(_fn) is not None for exc in exclude_progs
            ):
                continue
            fn = Path(root) / _fn
            if args.verbose:
                print(f"Computing percentile: {str(fn)}")
            image = tio.ScalarImage(fn)
            if not args.dry_run:
                p = np.percentile(image.numpy(), args.percentile)
                percentiles.append(p)
            else:
                percentiles.append(np.random.randn())
    percentile = np.median(percentiles)
    std = np.std(percentiles)
    if args.verbose:
        print(f"Median: {percentile}; Std: {std}")
    for root, dirs, files in os.walk(args.top_dir):
        for _fn in files:
            if all(inc.search(_fn) is None for inc in include_progs) or any(
                exc.search(_fn) is not None for exc in exclude_progs
            ):
                continue
            fn = Path(root) / _fn
            if args.verbose:
                print(f"Normalizing: {str(fn)}")
            image = tio.ScalarImage(fn)
            out_fn = str(fn).replace(str(fn.parents[len(fn.parents) - 2]) + "/", "")
            out_path = args.out_dir / out_fn
            if args.verbose:
                print(f"Saving normalized image: {str(out_path)}")
            if not args.dry_run:
                out_path.parent.mkdir(parents=True, exist_ok=True)
                normalized = image.numpy() / percentile
                image.set_data(torch.from_numpy(normalized))
                image.save(out_path)
    return 0


if __name__ == "__main__":
    sys.exit(main())

If you use this script in an academic paper, please cite the paper:

    @inproceedings{reinhold2019evaluating,
      title={Evaluating the impact of intensity normalization on {MR} image synthesis},
      author={Reinhold, Jacob C and Dewey, Blake E and Carass, Aaron and Prince, Jerry L},
      booktitle={Medical Imaging 2019: Image Processing},
      volume={10949},
      pages={109493H},
      year={2019},
      organization={International Society for Optics and Photonics}}