ANTs SyGN Algorithm Experiment
Weiquan Luo
2022-01-05
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.
- For the j-th input image:
- 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:
L2-norm=√∑nComp=1Deform2
where:
- Comp: the number of component. It is same as the dimention of the image.
- 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¶
SyNandSyNRAare not appropriate for the structurally heterogeneous dataset, due to the affine transformation at every image-to-template registrations.SyNOnlyis 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
