ANTs SyGN Algorithm Experiment

Table of Contents

• 1  Summary
  • 1.1  Introduction
    • 1.1.1  Pseudo code for SyGN algorightm
  • 1.2  Objective
  • 1.3  Method
    • 1.3.1  algorithm comparison
    • 1.3.2  Image characterization
    • 1.3.3  Measurement
  • 1.4  Result
    • 1.4.1  Final converagence
    • 1.4.2  Fitting diagnosis with learning curve
    • 1.4.3  Verification of overfitting
  • 1.5  Conclusion
  • 1.6  Future work
  • 1.7  Usage warning
• 2  User-defined function
  • 2.1  generating_circle_dataset()
  • 2.2  Averaging_ANTsImage_list()
  • 2.3  plot_2DImage_deform()
  • 2.4  plot_2Dgrid_transform()
  • 2.5  plot_2DImage_histogram()
  • 2.6  SyN_groupwise_registration()
  • 2.7  run_experiment()
• 3  run_experiment
  • 3.1  SyNOnly
  • 3.2  SyNRA
  • 3.3  SyN
  • 3.4  Additional

Summary¶

Introduction¶

SyN is the state-of-the-art (SOTA) pair-wise image registration algorithm. With naturally extending it from pair-wise to group-wise, SyGN is able to co-registering a group of images to generate the group-specific template.

In the ants.registration() API, the typeofTransform argument provided 10 types of SyN option. They are either integrate with additional function or specialize for certain task. Importantly, we still lack of understanding on the capacity of SyGN on different type of data.

Pseudo code for SyGN algorightm¶

• Create the template by evenly-weighed summing all input images.
• For the i-th iterations:
  • For the j-th input image:
    • Pair-wise image registration on the input image to template,
    • accumulate transforms by weighted sum.
    • accumulate warped images by weighted-sum.
  • Update new template by backward transforming the accumulated warped image to template space with a gradient.
• Iterate over all j-images for i-iterations.
• Done.

Objective¶

The objective of this study is to testify the performance of the SyGN algoritm on dataset with different characteristic.

Method¶

algorithm comparison¶

We tested the fitting capacity of SyGN algorithm on different type of data. In this experiment, we using the following three options out of the total ten typeofTransform:

• "SyN": Symmetric normalization: Affine + deformable transformation, with mutual information as optimization metric.
• "SyNRA": Symmetric normalization: Rigid + Affine + deformable transformation, with mutual information as optimization metric.
• "SyNOnly": Symmetric normalization: no initial transformation, with mutual information as optimization metric. Assumes images are aligned by an inital transformation. Can be useful if you want to run an unmasked affine followed by masked deformable registration.

Image characterization¶

The comparison used the control experiment schema. Multiple set of data were generated with different characteristic. To keep the experiment simple, each dataset consist of four 2D image with 128 pixel by 128 pixel.

Three lower-level factors were testified:

• if_aligned: almost-aligned vs non-aligned
• if_inside_circle: without vs with inside circle
• if_equal_area: equal area vs unequal area

Two higher-level factors were testified:

• if_equal_intensity: intensity heterogeneity
• additional dataset: structural heterogeneity

Measurement¶

Based on the locally Euclidean property of Riemannian geometry, I use the Euclidean distance (L2-norm) to quantify the distance between the original and the warped image within each registration iteration:

$\text{L2-norm} = \sqrt{\sum_{\text{Comp}=1}^\text{n}{\text{Deform}^2}}$

where:

• $\text{Comp}$: the number of component. It is same as the dimention of the image.
• $\text{Deform}$: a defomation matrix for one component. It have same shape of the image.

Result¶

Final converagence¶

For Aligned dataset, template converged better with SyN/SyNRA than SyNOnly. For on-aligned datset, SyNOnly resulted better convergence in more homogeneous datset, and have lower variance. However, strong learners often cause overfitting.

Fitting diagnosis with learning curve¶

With SyNOnly, the template improved over iterations. With SyNRA and SyN, the template almost converged at the first iteration. Protential evidence for Overfitting!

Verification of overfitting¶

The affine transform in SyN and SyNRA introduced randomness and overfitting. The randomness stacked up within each iteration. The affine transform is supposed to align image, but it did more than that. Noticeably, affine transformation continuedly added unexpected warping over the registration.

Conclusion¶

• SyN and SyNRA are not appropriate for the structurally heterogeneous dataset, due to the affine transformation at every image-to-template registrations.
• SyNOnly is the better option, although it require many more iterations for the template to converge. It might require more than 10 iteration, although the documentation says “should be greater than 1 less than 10”.

Future work¶

Usage warning¶

• Issue: This notebook are supposed to created high volume of temporary files over the image generating function. I have tried to manage them by anchoring them at local folder: ./my_tmp/ and deleting them once they become useless. However, due to slow synchronization of the temporary files to servres/clouds, the deletion might fail and cause subsequent issue, such as saving new generated images and transformation.
  • Suggestion: make sure high speed connection to the servres/clouds to reduce the lagging.

In [12]:

import os
# make sure current TMPDIR was set before open jupyer kernel: 
# export TMPDIR=./my_tmp/
# some older bash: setenv TMPDIR './my_tmp' 
print(f'Current TMPDIR: {os.environ["TMPDIR"]}')
import shutil
import glob


import cv2 # image processing
from random import choices
import numpy as np # data structure
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt


import ants 
from ants import ants_image_io as iio
from ants import iMath

Current TMPDIR: ./my_tmp/

User-defined function¶

generating_circle_dataset()¶

In [3]:

def generating_circle_dataset(sampling = 4, image_size = [128, 128], 
                              if_aligned = True,  if_inside_circle = True,
                              if_equal_area = True, if_equal_intensity = True, 
                              save_dir = None, show_image = False, if_verbose = False):
    """
    This is a image dataset generater, design for simulating data for evaluate algorithm. 
    It the mean values for radius, color and center were first generatd,
    then add the randomly generated residual to the means.
    --------
    ARGRUMENT
    sampling:
        int: number of sample image that are generated
    image_size: 
        array with 2 element
    if_aligned:
        If False, then not aligned;
        if True, then almost aligned,
    if_inside_circle: 
        If False, then no inside circle;
        if True, then has inside circle,
    if_equal_area:
        If False, then unequal area;
        if True, then equal area,
    if_equal_intensity:
        If False, then unequal intensity;
        if True, then equal intensity,
    save_dir:
        If None, then not save
        if path is given, then save to path
    show_image:
        if True, then show generated images,
        if False, otherwise.
    if_verbose:
        if True, then show message,
        if False, otherwise.

    """
    # # testing setup
    # sampling = 4
    # image_size = [128, 128]
    # if_aligned = False
    # if_inside_circle = True
    # if_equal_area = False
    # if_equal_intensity = False
    # show_image = True
    
    # intialize
    pi_x, pi_y = image_size
    circle_center = [int(pi_x/2), int(pi_y/2)]
    if if_verbose:
        print(f'image_size: {image_size}')
        print(f'circle_center: {circle_center}')

    # parameterization for mean
    mean_outer_radius = int(pi_x/4)
    mean_inner_radius = int(pi_x/8)
    mean_outer_color = int(255/3)
    mean_inner_color = int(255/3*2)
    mean_outer_radius_list = np.repeat(mean_outer_radius, sampling)
    mean_inner_radius_list = np.repeat(mean_inner_radius, sampling)
    mean_outer_color_list = np.repeat(mean_outer_color, sampling)
    mean_inner_color_list = np.repeat(mean_inner_color, sampling)
    if if_verbose:
        print(f'Mean color and radius')
        print(f' outer_radius_list: {mean_outer_radius_list}')
        print(f' mean_inner_radius_list: {mean_inner_radius_list}')
        print(f' mean_outer_color_list: {mean_outer_color_list}')
        print(f' mean_inner_color_list: {mean_inner_color_list}')

    # parameterization for residual
    radius_offset = int((mean_outer_radius - mean_inner_radius)/2)
    if if_equal_area:
        if if_verbose: print('Equal area:')
        outer_radius_list_res = np.repeat(choices(range(-radius_offset, radius_offset), k=1), sampling)
        inner_radius_list_res = np.repeat(choices(range(-radius_offset, radius_offset), k=1), sampling)
    else:
        if if_verbose: print('Unequal area')
        outer_radius_list_res = np.array(choices(range(-radius_offset, radius_offset), k=sampling))
        inner_radius_list_res = np.array(choices(range(-radius_offset, radius_offset), k=sampling))
    if if_verbose:
        print(f' outer_radius_list_res: {outer_radius_list_res}')
        print(f' inner_radius_list_res: {inner_radius_list_res}')

    color_offset = int((mean_inner_color - mean_outer_color)/2)
    if if_equal_intensity:
        if if_verbose: print('Equal intensity:')
        outer_color_list_res = np.repeat(choices(range(-color_offset, color_offset), k=1), sampling)
        inner_color_list_res = np.repeat(choices(range(-color_offset, color_offset), k=1), sampling)
    else:
        if if_verbose: print('Unequal intensity:')
        outer_color_list_res = np.array(choices(range(-color_offset, color_offset), k=sampling))
        inner_color_list_res = np.array(choices(range(-color_offset, color_offset), k=sampling))
    if if_verbose:
        print(f' outer_color_list_res: {outer_color_list_res}')
        print(f' inner_color_list_res: {inner_color_list_res}')

    if if_aligned:
        if if_verbose: print('Aligned')
        transition_x_list = choices(range(0,3), k=sampling)
        transition_y_list = choices(range(0,3), k=sampling)
    else:
        if if_verbose: print('Not aligned')
        transition_x_max_list = np.repeat(int(pi_x/2), sampling) - mean_outer_radius_list - outer_radius_list_res
        transition_y_max_list = np.repeat(int(pi_x/2), sampling) - mean_outer_radius_list - outer_radius_list_res
        transition_x_list = [choices(range(0, transition_x_max), k=1)[0] for transition_x_max in transition_x_max_list]
        transition_y_list = [choices(range(0, transition_y_max), k=1)[0] for transition_y_max in transition_y_max_list]
    if if_verbose:
        print(f' transition_x_list: {transition_x_list}')
        print(f' transition_y_list: {transition_y_list}')
    
    outer_radius_list = mean_outer_radius_list + outer_radius_list_res
    inner_radius_list = mean_inner_radius_list + inner_radius_list_res
    outer_color_list = mean_outer_color_list + outer_color_list_res
    inner_color_list = mean_inner_color_list + inner_color_list_res
    circle_center_x = np.repeat(circle_center[0], sampling) + transition_x_list
    circle_center_y = np.repeat(circle_center[1], sampling) + transition_y_list
    
    # collect image parameter
    image_para = {'outer_radius_list': outer_radius_list, 
          'inner_radius_list': inner_radius_list, 
          'outer_color_list': outer_color_list, 
          'inner_color_list': inner_color_list, 
          'circle_center_x': circle_center_x, 
          'circle_center_y': circle_center_y
         }
    if if_verbose:
        print('Adding residual to mean ')
        print(f' outer_radius_list: {outer_radius_list}')
        print(f' inner_radius_list: {inner_radius_list}')
        print(f' outer_color_list: {outer_color_list}')
        print(f' inner_color_list: {inner_color_list}')
        print(f' circle_center_x: {circle_center_x}')
        print(f' circle_center_y: {circle_center_y}')

    # create image
    my_circle_list = list()
    for s in range(sampling):
        # original image
        my_circle = np.zeros(shape=image_size, dtype=np.uint8 )
        if not if_inside_circle:
            my_circle = cv2.circle(my_circle, 
                                   center=[circle_center_x[s], circle_center_y[s]], 
                                   radius=outer_radius_list[s], 
                                   color=(int(outer_color_list[s])), 
                                   thickness= -1)
        else:
            my_circle = cv2.circle(my_circle, 
                                   center=[circle_center_x[s], circle_center_y[s]], 
                                   radius=outer_radius_list[s], 
                                   color=(int(outer_color_list[s])), 
                                   thickness= -1)
            my_circle = cv2.circle(my_circle, 
                                   center=[circle_center_x[s], circle_center_y[s]], 
                                   radius=inner_radius_list[s], 
                                   color=(int(inner_color_list[s])), 
                                   thickness= -1)
        my_circle_list.append(my_circle)

    # visualization
    if show_image:
        fig = plt.figure(figsize=(15, 15))
        columns = 4
        rows = 1
        for i in range(1, columns*rows +1):
            fig.add_subplot(rows, columns, i)
            plt.imshow(my_circle_list[i-1], cmap='gray', vmin=0, vmax=255)
        plt.show()
        
    # saving image
    if save_dir is not None:
        exp_dir= save_dir
        if not os.path.exists(exp_dir):
            if if_verbose: print(f'make path: {exp_dir}')
            os.mkdir(exp_dir)
        # create folder_dir and file, save file into folder_dir
        for i in range(sampling):
            folder_dir = exp_dir + 'image' + str(i+1).zfill(3) + '/'
            file = 'b_image' + str(i+1).zfill(3) + '.nii'
            if not os.path.exists(folder_dir):
                if if_verbose: print(f'make path: {folder_dir}')
                os.mkdir(folder_dir)
            if if_verbose: print(f'save: {folder_dir + file}')
            iio.image_write(ants.from_numpy(my_circle_list[i]), folder_dir + file)
            
    if if_verbose: print('\n\n')
    return my_circle_list, image_para

Averaging_ANTsImage_list()¶

In [13]:

def Averaging_ANTsImage_list(image_list, weights, filename):
    '''
    Averaging ANTsImage_list to result an initial_template
    '''
    initial_template = image_list[0] * 0
    for i in range(len(image_list)):
        temp = image_list[i] * weights[i]
        temp = ants.resample_image_to_target(temp, initial_template)
        initial_template = initial_template + temp
    iio.image_write(initial_template, filename)
    return(initial_template)

plot_2DImage_deform()¶

In [14]:

def plot_2DImage_deform(w, cmap = 'hsv', title = None):
    '''
    plot foward and inverse transforms, as well as making comparison
    ----
    argument: w is the outpout of ants.registration()
    '''
    fwdtransforms = iio.image_read(w["fwdtransforms"][0])
    invtransforms = iio.image_read(w["invtransforms"][1])

    vmin = min(list([fwdtransforms.min(), invtransforms.min()]))
    vmax = max(list([fwdtransforms.max(), invtransforms.max()]))

    fwdtransforms_x = fwdtransforms.numpy()[:,:,0]
    fwdtransforms_y = fwdtransforms.numpy()[:,:,1]
    invtransforms_x = invtransforms.numpy()[:,:,0]
    invtransforms_y = invtransforms.numpy()[:,:,1]

    fig, ((ax1, ax3, ax5),(ax2, ax4, ax6)) = plt.subplots(2,3, figsize=(18,11))
    fig.suptitle(title, fontsize=20)

    ax1.set_title('fwdtransforms_x')
    im1 = ax1.imshow(fwdtransforms_x, 
                     cmap= cmap, vmin=vmin, vmax=vmax)
    ax2.set_title('fwdtransforms_y')
    im2 = ax2.imshow(fwdtransforms_y, 
                     cmap= cmap, vmin=vmin, vmax=vmax)

    ax3.set_title('-1 * invtransforms_x')
    im3 = ax3.imshow(invtransforms_x * (-1.0), 
                     cmap= cmap, vmin=vmin, vmax=vmax)
    ax4.set_title('-1 * invtransforms_y ')
    im4 = ax4.imshow(invtransforms_y * (-1.0), 
                     cmap= cmap, vmin=vmin, vmax=vmax)

    ax5.set_title('abs_diff_Transform_x = \n abs(fwdtransforms_x - invtransforms_x)')
    im5 = ax5.imshow(abs(fwdtransforms_x + invtransforms_x), 
                     cmap= cmap, vmin=vmin, vmax=vmax)
    ax6.set_title('abs_diff_Transform_y = \n abs(fwdtransforms_y -  invtransforms_y)')
    im6 = ax6.imshow(abs(fwdtransforms_y +  invtransforms_y), 
                     cmap= cmap, vmin=vmin, vmax=vmax)

    fig.subplots_adjust(right=0.93)
    cbar_ax = fig.add_axes([0.95, 0.15, 0.02, 0.7])
    fig.colorbar(im1, cax=cbar_ax)
    plt.show()
    

plot_2Dgrid_transform()¶

In [15]:

def plot_2Dgrid_transform(fixed, moving, w, weighted_warpedmovout, 
                          tranform = 'fwdtransforms', 
                          vmin =0, vmax =255, cmap = 'gray',
                          line_spacing = 5, title = None):
    '''
    plot fixed, moving, new_imgrid
    ----
    argument: w is the outpout of ants.registration()
    '''
    imgrid = moving.numpy()
    imgrid[0::line_spacing,:]=133
    imgrid[:, 0::line_spacing]=133
    imgrid = ants.from_numpy(imgrid)
    my_warpedmovout = ants.apply_transforms(fixed, imgrid, w[tranform])
    
    # ploting
    fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(18,6))
    fig.suptitle(title, fontsize=20)

    ax1.set_title('Original oving image')
    im1 = ax1.imshow(moving.numpy(),
                     cmap= cmap, vmin=vmin, vmax=vmax)
    ax2.set_title('warped moving image')
    im2 = ax2.imshow(my_warpedmovout.numpy(),
                     cmap= cmap, vmin=vmin, vmax=vmax)
    ax3.set_title('weighted warped moving image')
    im3 = ax3.imshow(weighted_warpedmovout.numpy(),
                     cmap= cmap, vmin=vmin, vmax=vmax)

    fig.subplots_adjust(right=0.93)
    cbar_ax = fig.add_axes([0.95, 0.15, 0.02, 0.7])
    fig.colorbar(im1, cax=cbar_ax)
    plt.show()
    

plot_2DImage_histogram()¶

In [16]:

def plot_2DImage_histogram(diff, title=None):
    
    diff = diff.numpy()
    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12,5))
    fig.suptitle(title, fontsize=20)
    
    ax1.set_title('Image differencing')
    ax1.imshow(diff, cmap='gray', vmin=0, vmax=255)
    
    ax2.set_title('Histogram of the Image differencing')
    diff = diff.flatten()
    diff = diff[abs(diff) > 0.1]
    ax2.hist(diff, 30, density=False, facecolor='grey', alpha=1)
    
    plt.show()
    

SyN_groupwise_registration()¶

In [17]:

def SyN_groupwise_registration(dataDirectory, 
                               initial_template = None,
                               iterations=1,
                               gradient_step=0.2,
                               blending_weight=0.75,
                               weights=None,
                               type_of_transform = "SyNOnly",
                               if_plot = False,
                               if_verbose = False,
                               if_save = False
                              ):
    '''
    This function performs groupwise registration with ANTs SyN to obtain group-specific template. 
    It is modified from the original ANTs.build_template() function
    ----------
    Arguments:
        dataDirectory: 
        initial_template: 
            If None (default), generate initial_template by averaging all image. 
            If ANTsImage type, then use it as initial_template.
        iterations:
            int type. 
            Number of iterations
        gradient_step:
            float type. 
            The percentage of backward transform from template (average of warped images) back to the moving image.
        blending_weight:
            float type
        weights:
            If None (default), generate even weights for all input images.
            If list of number of all input images, then use the list as weights
        type_of_transform:
            str type. "SyNOnly" (default). the option of type_of_transform is same as the type_of_transform in ANTs.registration()
        if_plot:
            boolean type.
            If False (default), then no plot output
            If True, then plot output
        if_verbose:
            boolean type.
            If False (default), then no procedure message output
            If True, then procedure message output
        if_save:
            boolean type.
            If False (default), then no saving new template
            If True, then saving new template
    ----------
    Return:
        xavg:
            final template as in ANTsImage type
        L2norm_result:
            Trainning output as in pd.dataframe type. This data saved during each image registration and iteration
            The dataframe consist of four column:
                L2norm_fwdtransforms_list: forward transform of moving images to template, with 2 components
                L2norm_wavg_list: accumulate weight: accumulate of weighted forward transform within the same iteration
                iterations_index: counter for iterations
                img_index: counter for image registration
    '''
    
    # batch load .nii files
    ImageFiles = glob.glob(dataDirectory + "**/b_*.nii" , recursive=True)
    ImageFiles.sort()
    image_list = list()
    for i in range(len(ImageFiles)):
      image_list.append(ants.image_read(ImageFiles[i], dimension = 2))

    # print(f'Found {len(ImageFiles)} matching .nifti file(s):')
    # for file in ImageFiles:
    #     print(file)
    
    # check weight, create if not given
    if weights is None:
        weights = np.repeat(1.0 / len(image_list), len(image_list))
    weights = [x / sum(weights) for x in weights]

    # check initial_template, create if not given
    if initial_template is None:
        initial_template_name = dataDirectory + 'template_0.nii'
        if os.path.exists(initial_template_name):
            if if_verbose: print(f'using existing initial_template: {initial_template_name}')
            initial_template = ants.image_read(dataDirectory + 'template_0.nii')
        else: 
            if if_verbose: print(f'creating  initial_template: {initial_template_name}')
            initial_template = Averaging_ANTsImage_list(image_list = image_list, 
                                                        weights = weights, 
                                                        filename= initial_template_name)
            
    # empty list for collecting data
    L2norm_fwdtransforms_list = list()
    L2norm_wavg_list = list()
    iterations_index = list()
    img_index = list()

    xavg = initial_template.clone()
    for i in range(iterations):  
        if if_plot:
            print(f'Visualizing: template_{i}')
            plt.imshow(xavg.numpy(), cmap='gray', vmin=0, vmax=255)
            plt.show()
        for k in range(len(image_list)):
            iterations_index.append(i)
            img_index.append(k)
            w1 = ants.registration(xavg, image_list[k], type_of_transform=type_of_transform, verbose = False)
            fwdtransforms = iio.image_read(w1["fwdtransforms"][0])
            if k == 0:
                wavg = fwdtransforms * weights[k]
                xavgNew = w1["warpedmovout"] * weights[k]
            else:
                wavg = wavg + fwdtransforms * weights[k]
                xavgNew = xavgNew + w1["warpedmovout"] * weights[k]
            
            # collecting data per image registration and per iteration
            L2norm_fwdtransforms = np.sqrt(np.square(fwdtransforms.numpy()[:,:,0]) + np.square(fwdtransforms.numpy()[:,:,1])).sum()
            L2norm_wavg = np.sqrt(np.square(wavg.numpy()[:,:,0]) + np.square(wavg.numpy()[:,:,1])).sum()
            L2norm_fwdtransforms_list.append(L2norm_fwdtransforms)
            L2norm_wavg_list.append(L2norm_wavg)
            
            if if_verbose:
                print(f'fwdtransforms.abs().sum() of (iterations_{i+1}/{iterations}, image_{k+1}/{len(image_list)}): {L2norm_fwdtransforms}')   
                print(f'wavg.abs().sum() of (iterations_{i+1}/{iterations}, image_{k+1}/{len(image_list)}): {L2norm_wavg}')   
            if if_plot:
                plot_2DImage_deform(w1, cmap = 'hsv', 
                                    title = 'Image_'+ str(k+1) + ' deformation')
                plot_2Dgrid_transform(xavg, image_list[k], w1, xavgNew, 
                                      title= 'warpping image_'+ str(k+1) + ' with grid')
        if if_verbose: print(f'backward transforming the fully warpped image to scaled template with gradient = {-gradient_step}')
        wscl = (-1.0) * gradient_step
        wavg = wavg * wscl
        wavgfn = dataDirectory + 'wavg_' + str(i+1) + '.nii.gz'
        iio.image_write(wavg, wavgfn)
        xavg = ants.apply_transforms(xavgNew, xavgNew, wavgfn)
        if blending_weight is not None:
            xavg = xavg * blending_weight + iMath(xavg, "Sharpen") * (1.0 - blending_weight)

        # save iteration result
        if if_save:
            xavg_name = dataDirectory + 'template_' + str(i+1) + '.nii'
            print(f'New template saving as: {xavg_name}')
            iio.image_write(xavg, xavg_name)
        # visualize iteration result
        if if_plot & if_save: 
            plot_2DImage_histogram(diff = xavg- xavgNew, 
                                   title = 'Template differencing: fully_warped and stepback_warped')

    if if_plot:
        print(f'Visualizing final template: template_{i+1}')
        plt.imshow(xavg.numpy(), cmap='gray', vmin=0, vmax=255)
        plt.show()

    # packaging collected data into dataframe
    d = {'L2norm_fwdtransforms_list': L2norm_fwdtransforms_list,
         'L2norm_wavg_list': L2norm_wavg_list,
         'iterations_index': iterations_index,
         'img_index': img_index
        }
    L2norm_result = pd.DataFrame(d)

    if if_plot:
        L2norm_result_long = pd.melt(L2norm_result, 
                                     id_vars=['iterations_index', 'img_index'],
                                     value_vars=['L2norm_fwdtransforms_list', 'L2norm_wavg_list'])
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
        fig.suptitle('Trainning report')

        sns.lineplot(data = L2norm_result, 
                     x = 'iterations_index', y = 'L2norm_fwdtransforms_list',
                     hue = 'img_index', marker='.', ax = ax1)

        sns.lineplot(data = L2norm_result, 
                     x = 'img_index', y = 'L2norm_wavg_list',
                     hue = 'iterations_index', marker='.', ax = ax2)

        plt.show()
        
    return xavg, L2norm_result

run_experiment()¶

In [18]:

def run_experiment(tmp_experiment = './tmp_experiment_dataset/',
                   sampling = 4, 
                   type_of_transform = 'SyNOnly', 
                   iterations = 10
                  ):
    '''
    Run experiment base on Orthogonal Experimental Design. 
    This function incoporate the tmp file management that remove all tmp once the experiment done.
    This function consist of two man sub-functions:
        generating_circle_dataset() generate 16 dataset
        SyN_groupwise_registration() create group-wise dataset.
    
    '''

    # design matrix
    if_aligned = np.tile(np.repeat([True, False], 2),4)
    if_inside_circle = np.tile(np.repeat([False, True], 4),2)
    if_equal_area = np.tile([True, False], 8)
    if_equal_intensity = np.repeat([True, False], 8)
    save_dir = [tmp_experiment + 'Dataset_' + str(i).zfill(3) + '/' for i in np.arange(1,1+16)]
    d = {'if_aligned': if_aligned,
         'if_inside_circle': if_inside_circle,
         'if_equal_area': if_equal_area,
         'if_equal_intensity': if_equal_intensity,
         'save_dir': save_dir
        }
    dm = pd.DataFrame(data=d)

    # Generating image
    img_para_batch = list()
    for i in range(16):
        _, img_para_list = generating_circle_dataset(sampling = sampling,
                                                    if_aligned = dm.if_aligned[i], 
                                                    if_inside_circle = dm.if_inside_circle[i],
                                                    if_equal_area = dm.if_equal_area[i],
                                                    if_equal_intensity = dm.if_equal_intensity[i],
                                                    save_dir = dm.save_dir[i], 
#                                                         save_dir = None, 
                                                    show_image = False,
                                                    if_verbose = False)
        img_para_batch.append(img_para_list)

    # registration and collecting result
    L2norm_result_df = pd.DataFrame()
    for d in dm.save_dir:
        print(f'SyN groupwise registration for {d}')
        _, L2norm_result = SyN_groupwise_registration(dataDirectory = d,
                                                      initial_template = None,
                                                      iterations=iterations,
                                                      gradient_step=0.2,
                                                      blending_weight=0.75,
                                                      weights=None,
                                                      type_of_transform = type_of_transform,
                                                      if_plot = False,
                                                      if_verbose = False,
                                                      if_save = False)
        L2norm_result['save_dir'] = d
        L2norm_result_df = pd.concat([L2norm_result_df, L2norm_result])
        !find './my_tmp/' -type f -iname 'tmp*.mat' -maxdepth 1 | xargs rm
        !find './my_tmp/' -type f -iname '*.nii.gz' -maxdepth 1 | xargs rm

    # join design matrix and result
    L2norm_result_df['iterations_img_counter'] = L2norm_result_df.iterations_index * sampling + L2norm_result_df.img_index + 1
    L2norm_result_dm = L2norm_result_df.set_index('save_dir').\
        join(dm.set_index('save_dir')).\
        reset_index()
    return L2norm_result_dm

run_experiment¶

In [19]:

# experiment setup
repetition = 15 # repetition of the experiment
sampling = 100 # number of sample for each repetition
iterations = 10 # with in each group-wise registration
    
tmp_experiment = './tmp_experiment_dataset/'
if not os.path.exists(tmp_experiment):
    print(f'make path: {tmp_experiment}')
    os.mkdir(tmp_experiment)
    
experiment_result = './experiment_result/'
if not os.path.exists(experiment_result):
    print(f'make path: {experiment_result}')
    os.mkdir(experiment_result)
    
tmp_experiment_rep_list = [tmp_experiment + 'repetition_' + str(r).zfill(3) + '/' for r in range(repetition)]

SyNOnly¶

In [30]:

type_of_transform = 'SyNOnly'

In [1]:

L2norm_result_dm_rep = pd.DataFrame()
for r in range(repetition):
    
    if os.path.exists(tmp_experiment_rep_list[r]):
        shutil.rmtree(tmp_experiment_rep_list[r])
    print(f'repetition={r}')

    if not os.path.exists(tmp_experiment_rep_list[r]):
        print(f'make path: {tmp_experiment_rep_list[r]}')
        os.mkdir(tmp_experiment_rep_list[r])
    L2norm_result_dm = run_experiment(tmp_experiment = tmp_experiment_rep_list[r],
                                      sampling = sampling,
                                      type_of_transform = type_of_transform,
                                      iterations = iterations
                                     )
    
    L2norm_result_dm['repetition'] = r
    L2norm_result_dm_rep = pd.concat([L2norm_result_dm_rep, L2norm_result_dm])
    
    print(f'remove patyh: {tmp_experiment_rep_list[r]}')
    shutil.rmtree(tmp_experiment_rep_list[r])
L2norm_result_dm_rep = L2norm_result_dm_rep.sort_values(by = ['repetition',
                                                              'save_dir',
                                                              'iterations_index', 
                                                              'img_index']).\
    reset_index(drop = True)
L2norm_result_dm_rep.to_csv(experiment_result + type_of_transform + '.csv', index=False)

In [31]:

L2norm_result_dm_rep = pd.read_csv(experiment_result + type_of_transform + '.csv')

L2norm_result_dm_subset = L2norm_result_dm_rep[ L2norm_result_dm_rep['if_aligned']==True].\
    reset_index(drop = True)

g = sns.FacetGrid(L2norm_result_dm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', hue = 'if_inside_circle',
                  row_order=[True, False], col_order=[True, False], hue_order = [False, True],
                  height = 5 )
g.map_dataframe(sns.lineplot, x='iterations_img_counter', y = 'L2norm_wavg_list', marker ='.', ci = 68)
g.add_legend()
g.set(ylim=(0, 17000))
g.fig.subplots_adjust(top=0.9)
g.fig.suptitle('Group-wise registration report: SyNOnly, on aligned dataset')
# g.set(yscale = 'log')

L2norm_result_dm_subset = L2norm_result_dm_rep[ L2norm_result_dm_rep['if_aligned']==False].\
    reset_index(drop = True)

g = sns.FacetGrid(L2norm_result_dm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', hue = 'if_inside_circle',
                  row_order=[True, False], col_order=[True, False], hue_order = [False, True],
                  height = 5 )
g.map_dataframe(sns.lineplot, x='iterations_img_counter', y = 'L2norm_wavg_list', marker ='.', ci = 68)
g.add_legend()
g.set(ylim=(0, 17000))
g.fig.subplots_adjust(top=0.9)
g.fig.suptitle('Group-wise registration report: SyNOnly, on non-aligned dataset')
# g.set(yscale = 'log')

Out[31]:

Text(0.5, 0.98, 'Group-wise registration report: SyNOnly, on non-aligned dataset')

SyNRA¶

In [32]:

type_of_transform = 'SyNRA'

In [ ]:

L2norm_result_dm_rep = pd.DataFrame()
for r in range(repetition):
    
    if os.path.exists(tmp_experiment_rep_list[r]):
        shutil.rmtree(tmp_experiment_rep_list[r])
    print(f'repetition={r}')

    if not os.path.exists(tmp_experiment_rep_list[r]):
        print(f'make path: {tmp_experiment_rep_list[r]}')
        os.mkdir(tmp_experiment_rep_list[r])
    L2norm_result_dm = run_experiment(tmp_experiment = tmp_experiment_rep_list[r],
                                      sampling = sampling,
                                      type_of_transform = type_of_transform,
                                      iterations = iterations
                                     )
    
    L2norm_result_dm['repetition'] = r
    L2norm_result_dm_rep = pd.concat([L2norm_result_dm_rep, L2norm_result_dm])
    
    print(f'remove patyh: {tmp_experiment_rep_list[r]}')
    shutil.rmtree(tmp_experiment_rep_list[r])
L2norm_result_dm_rep = L2norm_result_dm_rep.sort_values(by = ['repetition',
                                                              'save_dir',
                                                              'iterations_index', 
                                                              'img_index']).\
    reset_index(drop = True)
L2norm_result_dm_rep.to_csv(experiment_result + type_of_transform + '.csv', index=False)

In [33]:

L2norm_result_dm_rep = pd.read_csv(experiment_result + type_of_transform + '.csv')

L2norm_result_dm_subset = L2norm_result_dm_rep[ L2norm_result_dm_rep['if_aligned']==True].\
    reset_index(drop = True)

g = sns.FacetGrid(L2norm_result_dm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', hue = 'if_inside_circle',
                  row_order=[True, False], col_order=[True, False], hue_order = [False, True],
                  height = 5 )
g.map_dataframe(sns.lineplot, x='iterations_img_counter', y = 'L2norm_wavg_list', marker ='.', ci = 68)
g.add_legend()
g.set(ylim=(0, 17000))
g.fig.subplots_adjust(top=0.9)
g.fig.suptitle('Group-wise registration report: SyNRA, on aligned dataset')
# g.set(yscale = 'log')

L2norm_result_dm_subset = L2norm_result_dm_rep[ L2norm_result_dm_rep['if_aligned']==False].\
    reset_index(drop = True)

g = sns.FacetGrid(L2norm_result_dm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', hue = 'if_inside_circle',
                  row_order=[True, False], col_order=[True, False], hue_order = [False, True],
                  height = 5 )
g.map_dataframe(sns.lineplot, x='iterations_img_counter', y = 'L2norm_wavg_list', marker ='.', ci = 68)
g.add_legend()
g.set(ylim=(0, 17000))
g.fig.subplots_adjust(top=0.9)
g.fig.suptitle('Group-wise registration report: SyNRA, on non-aligned dataset')
# g.set(yscale = 'log')

Out[33]:

Text(0.5, 0.98, 'Group-wise registration report: SyNRA, on non-aligned dataset')

SyN¶

In [34]:

type_of_transform = 'SyN'

In [16]:

L2norm_result_dm_rep = pd.DataFrame()
for r in range(repetition):
    
    if os.path.exists(tmp_experiment_rep_list[r]):
        shutil.rmtree(tmp_experiment_rep_list[r])
    print(f'repetition={r}')

    if not os.path.exists(tmp_experiment_rep_list[r]):
        print(f'make path: {tmp_experiment_rep_list[r]}')
        os.mkdir(tmp_experiment_rep_list[r])
    L2norm_result_dm = run_experiment(tmp_experiment = tmp_experiment_rep_list[r],
                                      sampling = sampling,
                                      type_of_transform = type_of_transform,
                                      iterations = iterations
                                     )
    
    L2norm_result_dm['repetition'] = r
    L2norm_result_dm_rep = pd.concat([L2norm_result_dm_rep, L2norm_result_dm])
    
    print(f'remove patyh: {tmp_experiment_rep_list[r]}')
    shutil.rmtree(tmp_experiment_rep_list[r])

L2norm_result_dm_rep = L2norm_result_dm_rep.sort_values(by = ['repetition',
                                                              'save_dir',
                                                              'iterations_index', 
                                                              'img_index']).\
    reset_index(drop = True)
L2norm_result_dm_rep.to_csv(experiment_result + type_of_transform + '.csv', index=False)

In [35]:

L2norm_result_dm_rep = pd.read_csv(experiment_result + type_of_transform + '.csv')

L2norm_result_dm_subset = L2norm_result_dm_rep[ L2norm_result_dm_rep['if_aligned']==True].\
    reset_index(drop = True)
g = sns.FacetGrid(L2norm_result_dm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', hue = 'if_inside_circle',
                  row_order=[True, False], col_order=[True, False], hue_order = [False, True],
                  height = 5 )
g.map_dataframe(sns.lineplot, x='iterations_img_counter', y = 'L2norm_wavg_list', marker ='.', ci = 68)
g.set(ylim=(0, 17000))
g.add_legend()
g.fig.subplots_adjust(top=0.9)
g.fig.suptitle('Group-wise registration report: SyN, on aligned dataset')
# g.set(yscale = 'log')

L2norm_result_dm_subset = L2norm_result_dm_rep[ L2norm_result_dm_rep['if_aligned']==False].\
    reset_index(drop = True)

g = sns.FacetGrid(L2norm_result_dm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', hue = 'if_inside_circle',
                  row_order=[True, False], col_order=[True, False], hue_order = [False, True],
                  height = 5 )
g.map_dataframe(sns.lineplot, x='iterations_img_counter', y = 'L2norm_wavg_list', marker ='.', ci = 68)
g.add_legend()
g.set(ylim=(0, 17000))
g.fig.subplots_adjust(top=0.9)
g.fig.suptitle('Group-wise registration report: SyN, on non-aligned dataset')
# g.set(yscale = 'log')

Out[35]:

Text(0.5, 0.98, 'Group-wise registration report: SyN, on non-aligned dataset')

In [81]:

SyNOnly_L2norm_result_dm_rep = pd.read_csv(experiment_result + 'SyNOnly' + '.csv')
last_SyNOnly = SyNOnly_L2norm_result_dm_rep[(SyNOnly_L2norm_result_dm_rep['iterations_index']==9) & (SyNOnly_L2norm_result_dm_rep['img_index']==99)]
last_SyNOnly = last_SyNOnly[['if_aligned', 'if_inside_circle', 'if_equal_area', 'if_equal_intensity', 'L2norm_wavg_list']]
last_SyNOnly['typeofTransform'] = 'SyNOnly'

SyNRA_L2norm_result_dm_rep = pd.read_csv(experiment_result + 'SyNRA' + '.csv')
last_SyNRA = SyNRA_L2norm_result_dm_rep[(SyNRA_L2norm_result_dm_rep['iterations_index']==9) & (SyNRA_L2norm_result_dm_rep['img_index']==99)]
last_SyNRA = last_SyNRA[['if_aligned', 'if_inside_circle', 'if_equal_area', 'if_equal_intensity', 'L2norm_wavg_list']]
last_SyNRA['typeofTransform'] = 'SyNRA'

SyN_L2norm_result_dm_rep = pd.read_csv(experiment_result + 'SyN' + '.csv')
last_SyN = SyN_L2norm_result_dm_rep[(SyN_L2norm_result_dm_rep['iterations_index']==9) & (SyN_L2norm_result_dm_rep['img_index']==99)]
last_SyN = last_SyN[['if_aligned', 'if_inside_circle', 'if_equal_area', 'if_equal_intensity', 'L2norm_wavg_list']]
last_SyN['typeofTransform'] = 'SyN'

last_L2norm = pd.concat([last_SyNOnly, last_SyNRA, last_SyN]).reset_index(drop = True)
last_L2norm

Out[81]:

	if_aligned	if_inside_circle	if_equal_area	if_equal_intensity	L2norm_wavg_list	typeofTransform
0	True	False	True	True	1101.862793	SyNOnly
1	True	False	False	True	1664.642212	SyNOnly
2	False	False	True	True	950.564331	SyNOnly
3	False	False	False	True	1963.993164	SyNOnly
4	True	True	True	True	1750.734375	SyNOnly
...	...	...	...	...	...	...
715	False	False	False	False	2278.352295	SyN
716	True	True	True	False	1484.059082	SyN
717	True	True	False	False	2062.238281	SyN
718	False	True	True	False	2869.317871	SyN
719	False	True	False	False	2435.089844	SyN

720 rows × 6 columns

In [107]:

last_L2norm_subset = last_L2norm[last_L2norm['if_aligned']==True].reset_index(drop = True)
g = sns.FacetGrid(last_L2norm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', 
                  row_order=[True, False], col_order=[True, False],
                  height = 5 )
g.map_dataframe(sns.boxplot, hue = 'typeofTransform', y = 'L2norm_wavg_list', x = 'if_inside_circle',)
g.add_legend()
g.fig.subplots_adjust(top=0.9)
g.set(ylim=(0, 4500))
g.fig.suptitle('Group-wise registration: final loss on aligned dataset')

last_L2norm_subset = last_L2norm[last_L2norm['if_aligned']==False].reset_index(drop = True)
g = sns.FacetGrid(last_L2norm_subset, 
                  row = 'if_equal_intensity', col ='if_equal_area', 
                  row_order=[True, False], col_order=[True, False],
                  height = 5 )
g.map_dataframe(sns.boxplot, hue = 'typeofTransform', y = 'L2norm_wavg_list', x = 'if_inside_circle',)
g.add_legend()
g.fig.subplots_adjust(top=0.9)
g.set(ylim=(0, 4500))
g.fig.suptitle('Group-wise registration: final loss on non-aligned dataset')

Out[107]:

Text(0.5, 0.98, 'Group-wise registration: final loss on non-aligned dataset')

Additional¶

In [20]:

contours = np.array([[48,48], [48,80], [80,80], [80,48]])
my_rectangle = np.zeros((128,128))
my_rectangle = cv2.fillPoly(my_rectangle, pts = [contours], color =(255))

contours = np.array([[60, 60], [48,64], [60, 68], [64,80], [68, 68], [80,64], [68, 60], [64, 48]])
my_star = np.zeros((128,128))
my_star = cv2.fillPoly(my_star, pts = [contours], color =(255))

contours = np.array([[52, 52], [48,64], [52, 76], [64,80], [76, 76], [80,64], [76, 52], [64, 48]])
my_octangle = np.zeros((128,128))
my_octangle = cv2.fillPoly(my_octangle, pts = [contours], color =(255))

my_circle = np.zeros(shape=[128, 128], dtype=np.uint8 )
my_circle = cv2.circle(my_circle, center=(64, 64), radius=16, color=(255), thickness= -1)

myshape = [my_rectangle, my_star, my_octangle, my_circle]

# visualization
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4, figsize=(18,5))

im1 = ax1.imshow(my_rectangle, cmap='gray', vmin=0, vmax=255)
im2 = ax2.imshow(my_star, cmap='gray', vmin=0, vmax=255)
im3 = ax3.imshow(my_octangle, cmap='gray', vmin=0, vmax=255)
im4 = ax4.imshow(my_circle, cmap='gray', vmin=0, vmax=255)


fig.subplots_adjust(right=0.93)
cbar_ax = fig.add_axes([0.95, 0.15, 0.02, 0.7])
fig.colorbar(im1, cax=cbar_ax)
plt.show()


tmp_experiment = './tmp_experiment_dataset/'
exp_dir= tmp_experiment + 'Dataset_' + 'additional/'
if not os.path.exists(exp_dir):
    print(f'make path: {exp_dir}')
    os.mkdir(exp_dir)

# create folder_dir and file, save file into folder_dir
for i in range(4):
    folder_dir = exp_dir + 'polygon' + str(i).zfill(3) + '/'
    file = 'b_polygon' + str(i).zfill(3) + '.nii'
    if not os.path.exists(folder_dir):
        print(f'make path: {folder_dir}')
        os.mkdir(folder_dir)
    print(f'save: {folder_dir + file}')
#     nib.save(nib.Nifti1Image(circlexxx[i], np.eye(4)), folder_dir + file) 
    iio.image_write(ants.from_numpy(myshape[i]), folder_dir + file)

make path: ./tmp_experiment_dataset/Dataset_additional/polygon000/
save: ./tmp_experiment_dataset/Dataset_additional/polygon000/b_polygon000.nii
make path: ./tmp_experiment_dataset/Dataset_additional/polygon001/
save: ./tmp_experiment_dataset/Dataset_additional/polygon001/b_polygon001.nii
make path: ./tmp_experiment_dataset/Dataset_additional/polygon002/
save: ./tmp_experiment_dataset/Dataset_additional/polygon002/b_polygon002.nii
make path: ./tmp_experiment_dataset/Dataset_additional/polygon003/
save: ./tmp_experiment_dataset/Dataset_additional/polygon003/b_polygon003.nii

In [ ]:

repetition = 15 # repetition of the experiment
iterations = 10 # with in each group-wise registration

type_of_transform = ['SyNOnly','SyNRA', 'SyN']

L2norm_result_rep = pd.DataFrame()
for t in type_of_transform:
#     print(t)
    for r in range(repetition):
#         print(r)
        _, L2norm_result = SyN_groupwise_registration(dataDirectory = exp_dir,
                                                      initial_template = None,
                                                      iterations = iterations,
                                                      gradient_step=0.2,
                                                      blending_weight=0.75,
                                                      weights=None,
                                                      type_of_transform = t,
                                                      if_plot = False,
                                                      if_verbose = False,
                                                      if_save = False)

        L2norm_result['repetition'] = r
        L2norm_result['typeofTransform'] = t
        L2norm_result_rep = pd.concat([L2norm_result_rep, L2norm_result])

In [154]:

L2norm_result_rep['iterations_img_counter'] = L2norm_result_rep.iterations_index * 4 + L2norm_result_rep.img_index + 1
L2norm_result_rep = L2norm_result_rep.sort_values(by = ['typeofTransform', 'iterations_img_counter', 'repetition']).\
    reset_index(drop = True)
L2norm_result_rep.to_csv(experiment_result + t + '_additional' + '.csv', index=False)
L2norm_result_rep

Out[154]:

	L2norm_fwdtransforms_list	L2norm_wavg_list	iterations_index	img_index	repetition	typeofTransform	iterations_img_counter
0	2224.970703	556.242676	0	0	0	SyN	1
1	2203.076904	550.769226	0	0	1	SyN	1
2	2696.454102	674.113525	0	0	2	SyN	1
3	3369.958984	842.489746	0	0	3	SyN	1
4	3349.531006	837.382751	0	0	4	SyN	1
...	...	...	...	...	...	...	...
1795	3188.124268	3026.107422	9	3	10	SyNRA	40
1796	10124.568359	6922.455078	9	3	11	SyNRA	40
1797	3530.222656	5200.533691	9	3	12	SyNRA	40
1798	4601.507812	4964.696289	9	3	13	SyNRA	40
1799	4414.956055	6745.846191	9	3	14	SyNRA	40

1800 rows × 7 columns

In [156]:

sns.lineplot(data = L2norm_result_rep, 
             x = 'iterations_img_counter', 
             y = 'L2norm_wavg_list', 
             hue = 'typeofTransform' )

Out[156]:

<AxesSubplot:xlabel='iterations_img_counter', ylabel='L2norm_wavg_list'>

In [158]:

_, L2norm_result = SyN_groupwise_registration(dataDirectory = exp_dir,
                                              initial_template = None,
                                              iterations = 5,
                                              gradient_step=0.2,
                                              blending_weight=0.75,
                                              weights=None,
                                              type_of_transform = 'SyNRA',
                                              if_plot = True,
                                              if_verbose = False,
                                              if_save = False)

Visualizing: template_0

Visualizing: template_1

Visualizing: template_2

Visualizing: template_3

Visualizing: template_4

Visualizing final template: template_5