<jupyter_start><jupyter_text>Data PreprocessingDevelopping a pipeline to normalize and clean the MRI images.   Data Import<jupyter_code># export
# default_exp preprocessing
import pandas as pd
from pathlib import Path
import SimpleITK as sitk
from faimed3d.all import *<jupyter_output><empty_output><jupyter_text>Data is stored on a secure network drive<jupyter_code>DATA_DIR = Path('/media/ScaleOut/vahldiek/MRI/SIJ')<jupyter_output><empty_output><jupyter_text>The preferred structure of the dataframe changed during the project. So I need to adapt it a bit.<jupyter_code>df = pd.read_csv(DATA_DIR/'dataset_training_477.csv')
df = df[df.has_TIRM_and_T1 == 1]
df['path_TIRM'] = [fn.replace('../data/', '') for fn in df.path_TIRM]
df['path_T1'] = [fn.replace('../data/', '') for fn in df.path_T1]<jupyter_output><empty_output><jupyter_text>ExplorationUsing the DicomExplorer from `faimed3d` multiple images are viewed and analyzed. It seems, the pixel/voxel values have a strongly skewed distribution with most pixel-values beeing on the lower end of the histogram and a long tail of high intensities.<jupyter_code>subset = pd.concat([df[df.cohort == cohort].sample(5) for cohort in df.cohort.unique()]) # viewing a representative subset
fn = DATA_DIR/'TRAINING/01012/T1' # large patient
fn = DATA_DIR/'TRAINING/01093/T1/' # act. changes, struct changes, SPA and ASAS compatible
#fn = DATA_DIR/'TRAINING/01084/TIRM/' # no act but struct
image = TensorDicom3D.create(fn)
image = Resample3D((20, 224, 224), (3.5, 1, 1))(image)
DicomExplorer(image, figsize = (5,5)).show()<jupyter_output>Statistics:
  Mean px value:     500.3769226074219 
  Std of px values:  473.8222961425781 
  Min px value:      0.0 
  Max px value:      2323.333251953125 
  Median px value:   329.7845153808594 
  Skewness:          1.2070366144180298 
  Kurtosis:          0.37535691261291504 

Tensor properties 
  Tensor shape:      (20, 224, 224)
  Tensor dtype:      torch.float32<jupyter_text>The pixel values show a heavy right tailed distribution and different maximum intensitites.  Most images contain parts of the MRI table, the lower image areas are black, which explains part of the skewed distributions.  The image areas with very high intensities are probably from CSF.  This observation is true for all the 20 sampled images I viewed.  Clipping the lower range of the pixels seems to be no problem,as it only turns the very hypointense background to black. Clipping very high pixel intensities could also be save, but as active inflammation is also hyperintens, one needs to be carefull.   Preprocessing steps Field Bias Correction<jupyter_code># export
def field_bias_correction(image, numberFittingLevels = 4, numberOfIteration = [50]):
    corrector = sitk.N4BiasFieldCorrectionImageFilter()
    corrector.SetMaximumNumberOfIterations(numberOfIteration * numberFittingLevels)
    bias_corrected = corrector.Execute(
        sitk.Cast(image.as_sitk(), 
                  sitk.sitkFloat32))
    bias_corrected = TensorDicom3D.create(bias_corrected)
    bias_corrected._metadata = image._metadata
    return bias_corrected
# slow
bias_corrected_image = field_bias_correction(image)
DicomExplorer(bias_corrected_image, figsize = (5,5)).show()<jupyter_output>Statistics:
  Mean px value:     481.5676574707031 
  Std of px values:  428.35113525390625 
  Min px value:      0.0 
  Max px value:      1889.69287109375 
  Median px value:   331.3974914550781 
  Skewness:          1.0523874759674072 
  Kurtosis:          -0.09217572212219238 

Tensor properties 
  Tensor shape:      (20, 224, 224)
  Tensor dtype:      torch.float32<jupyter_text>Denoising<jupyter_code># export
def denoising(image):
    denoised_image = sitk.CurvatureFlow(
        image1=image.as_sitk(),
        timeStep=0.125,
        numberOfIterations=3)
    denoised_image = TensorDicom3D.create(denoised_image)
    denoised_image._metadata = image._metadata
    return denoised_image
# slow
denoised_image = denoising(image)
DicomExplorer(denoised_image, figsize = (5,5)).show()<jupyter_output><empty_output><jupyter_text>Rescaling<jupyter_code>def rescaling(image, lwr=0, upr=255):
    resacleFilter = sitk.RescaleIntensityImageFilter()
    resacleFilter.SetOutputMaximum(upr)
    resacleFilter.SetOutputMinimum(lwr)
    scaled_image = resacleFilter.Execute(image.as_sitk())
    return TensorDicom3D.create(scaled_image)
# slow
scaled_image = rescaling(image)
DicomExplorer(scaled_image, figsize = (5,5)).show()<jupyter_output><empty_output><jupyter_text>Clipping> Clip the tail of the image<jupyter_code># export
def clip_tail(image, lwr_quant=0, upr_quant=0.99):
    lwr = image.quantile(lwr_quant)
    upr = image.quantile(upr_quant)
    return image.clip(lwr, upr)
# slow
clipped_image = clip_tail(image)
DicomExplorer(clipped_image, figsize = (5,5)).show()<jupyter_output><empty_output><jupyter_text>Histogram Matching<jupyter_code>spacing = TensorDicom3D.create(DATA_DIR/list(subset.path_T1)[0]).get_spacing()
origin = TensorDicom3D.create(DATA_DIR/list(subset.path_T1)[0]).get_origin()
direction = TensorDicom3D.create(DATA_DIR/list(subset.path_T1)[0]).get_direction()
# slow
batch = torch.stack([Resample3D((20, 224, 224), (3.5, 1, 1))(TensorDicom3D.create(DATA_DIR/fn)) for fn in subset.path_T1])
DicomExplorer(batch.mean(0)).show()
# slow
reference_image = batch.mean(0)
reference_image.set_spacing(spacing)
reference_image.set_direction(direction)
reference_image.set_origin(origin)
# export
def hist_matching(image, reference):
    hist_matching = sitk.HistogramMatchingImageFilter()
    hist_matching.ThresholdAtMeanIntensityOn()
    hist_matched_image = hist_matching.Execute(image.as_sitk(), reference_image.as_sitk())
    return TensorDicom3D.create(hist_matched_image)
# slow
hist_matched_image = hist_matching(image, reference_image)
DicomExplorer(hist_matched_image).show()
# export
class Pipeline():
    def __init__(self, functions:(list, tuple)):
        assert isinstance(functions, (list, tuple)), 'Functions need to be in a list or tuple'
        self.functions = functions
        
    def append(self, f):
        self.functions.append(f)
    
    def __call__(self, x):
        for f in self.functions: x = f(x)
        return x
pipe = Pipeline([Resample3D((20, 224, 224), (3.5, 1, 1)), clip_tail, denoising, field_bias_correction])
# disabled
clipped_and_denoised_image = denoising(clipped_image)
clipped_and_bias_corrected_image = field_bias_correction(clipped_image)
clipped_and_denoised_and_bias_corrected = pipe(image)
# disabled
clipped_and_denoised_and_bias_corrected_and_hist_matched = hist_matching(clipped_and_denoised_and_bias_corrected, reference_image)
# disabled
images = [image, clipped_image, denoised_image, 
          bias_corrected_image, clipped_and_denoised_image, clipped_and_bias_corrected_image, 
          clipped_and_denoised_and_bias_corrected, clipped_and_denoised_and_bias_corrected_and_hist_matched]

labels = ['original', 'clipped_image', 'denoised_image', 
          'bias_corrected_image', 'clipped_and_denoised_image', 
          'clipped_and_bias_corrected_image', 'clipped_and_denoised_and_bias_corrected', 'clipped_and_denoised_and_bias_corrected_and_hist_matched']
# disabled
DicomExplorer(clip_tail(clipped_and_denoised_and_bias_corrected, upr_quant=0.995),  cmap='Greys_r').show()
# disabled
ListViewer(images, labels, cmap='Greys_r').show()
# export
preprocessing_pipeline = Pipeline([Resample3D((20, 224, 224), (3.5, 1, 1)), clip_tail, field_bias_correction])<jupyter_output><empty_output><jupyter_text>It is now possible to loop through all images in the dataset, apply the preprocessing pipeline and then save the preprovessed image as NIfTI.<jupyter_code># hide
from nbdev.export import *
notebook2script()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09540363074546157, 'fraction_numerical': 0.045834942706321616, 'mean_word_length': 5.467943380516236, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '559152', 'license': 'Apache-2.0', 'n_tokens_mistral': 2929, 'n_tokens_neox': 2638, 'n_words': 579, 'path': 'nbs/01_preprocessing.ipynb', 'repo_name': 'kbressem/spa', 'size': '7767'}