Script for generating figures and exploring results for the purposes of my Multi-Velocity LDDMM final project for Comp 775 at UNC (Fall 2012).

experiments.py

import sh
import gen
from matplotlib import pyplot as plt
import matplotlib.cm as cm
import numpy as np
import SimpleITK as sitk
import sys

lddmm_command = "/Users/cpatrick/calatk-build/CalaTK-Build/bin/LDDMM"
apply_map_command = "/Users/cpatrick/calatk-build/CalaTK-Build/bin/applyMap"

lddmm = sh.Command(lddmm_command)
apply_map = sh.Command(apply_map_command)


def main():
    name = sys.argv[1]
    thismodule = sys.modules[__name__]
    cur_function = getattr(thismodule, name)
    cur_function()


def exp1():
    fixed_filename = 'fixed.mha'
    moving_filename = 'moving.mha'
    output_filename = 'lddmm_moving.mha'
    deformation_filename = 'lddmm_deformation.mha'
    config_filename = 'lddmm_sample.json'
    gen.make_square_image_2d_fixed(fixed_filename)
    gen.make_square_image_2d_moving(moving_filename)

    output = lddmm(fixed_filename, moving_filename, deformation_filename,
        '--config', config_filename, '--solver', 'relaxation')
    print output
    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title('Output after Deformation')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data)

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input)

    plt.figure()
    plt.title('Fixed')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data)

    plt.figure()
    plt.title('Moving')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data)
    plt.show()


def exp2():
    img = np.zeros((128, 128), dtype=np.uint8)
    img[20:80, 20:80] = 128
    frequency_image = np.fft.rfft2(img)
    out_image = np.fft.irfft2(frequency_image)
    plt.figure()
    plt.title('After FFT')
    plt.imshow(frequency_image)
    plt.show()
    plt.figure()
    plt.title('After FFT -> IFFT')
    plt.imshow(out_image)
    plt.show()


def exp3_bullseye():
    fixed_filename = 'bullseye2DTimePoint0.mha'
    moving_filename = 'bullseye2DTimePoint1.mha'
    output_filename = 'lddmm_bullseye.mha'
    deformation_filename = 'lddmm_bullseye_deformation.mha'
    config_filename = 'lddmm_sample.json'

    output = lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'relaxation')
    print output
    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title('Output after Deformation')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data)

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input)

    plt.figure()
    plt.title('Fixed')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data)

    plt.figure()
    plt.title('Moving')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data)
    plt.show()


def exp4():
    """Brad suggested running with a growing circle."""
    fixed_filename = 'circle.mha'
    moving_filename = 'blank.mha'
    output_filename = 'lddmm_circle.mha'
    deformation_filename = 'lddmm_circle_deformation.mha'
    config_filename = 'test_lddmm_config.json'

    img = np.zeros((32, 32), dtype=np.float)
    cx, cy = 16, 16  # The center of circle
    radius = 5
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, fixed_filename)

    cx, cy = 16, 16  # The center of circle
    radius = 10
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, moving_filename)

    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'relaxation', '--config', config_filename,
        '--logLevel', '4', _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title("LDDMM Generated Deformation Field")
    data = sitk.GetArrayFromImage(sitk.ReadImage(deformation_filename))
    plt.quiver(range(0, 32), range(0, 32), data[0], data[1], pivot='mid', color='r')

    plt.figure()
    plt.title('LDDMM Output after Deformation')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_output.eps')

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_difference.eps')

    plt.figure()
    plt.title('Fixed')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_source.eps')

    plt.figure()
    plt.title('Moving')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_target.eps')
    plt.show()


def exp5():
    """This is a small circle moving a few pixels."""
    fixed_filename = 'exp5_circle.mha'
    moving_filename = 'exp5_circle_moved.mha'
    output_filename = 'exp5_lddmm_circle.mha'
    deformation_filename = 'exp5_lddmm_circle_deformation.mha'
    config_filename = 'test_lddmm_config.json'

    img = np.zeros((32, 32), dtype=np.float)
    cx, cy = 16, 16  # The center of circle
    radius = 5
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, fixed_filename)

    img *= 0
    cx, cy = 21, 21  # The center of circle
    radius = 5
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, moving_filename)

    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'relaxation', '--config', config_filename,
        '--logLevel', '4', _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title('LDDMM Generated Deformed Output')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)

    plt.figure()
    plt.title("LDDMM Generated Deformation Field")
    data = sitk.GetArrayFromImage(sitk.ReadImage(deformation_filename))
    plt.quiver(data[0], data[1], pivot='mid', color='r')

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Fixed')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Moving')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.show()


def exp6():
    """Running LDDMM on Danielle's gradients"""
    fixed_filename = 'fixed_slide_2d.mha'
    moving_filename = 'moving_slide_2d.mha'
    output_filename = 'exp6_lddmm_slide_2d.mha'
    deformation_filename = 'exp6_lddmm_slide_2d_deformation.mha'
    config_filename = 'test_lddmm_config.json'
    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'relaxation', '--config', config_filename,
        '--logLevel', '4', _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title("LDDMM Generated Deformation Field")
    data = sitk.GetArrayFromImage(sitk.ReadImage(deformation_filename))
    plt.quiver(data[0][2:-2][2:-2], data[1][2:-2][2:-2], pivot='mid', color='r')

    plt.figure()
    plt.title('LDDMM Deformed Output')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('slide_lddmm_output.eps')

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)
    plt.savefig('slide_lddmm_difference.eps')

    plt.figure()
    plt.title('Source Image')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('slide_lddmm_source.eps')

    plt.figure()
    plt.title('Target Image')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('slide_lddmm_target.eps')
    plt.show()


def exp7():
    """Running LDDMM on Danielle's gradients"""
    fixed_filename = 'fixed_slide_2d.mha'
    moving_filename = 'moving_slide_2d.mha'
    output_filename = 'exp7_lddmm_slide_2d.mha'
    deformation_filename = 'exp7_lddmm_slide_2d_deformation.mha'
    config_filename = 'test_lddmm_config.json'
    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'multirelaxation', '--config', config_filename,
        '--logLevel', '4', _err_to_out=True, _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title('Output after Deformation')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Fixed')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Moving')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.show()


def monkey_lddmm():
    """Running LDDMM on monkey brains"""
    fixed_filename = 'monkey_01_02.mhd'
    moving_filename = 'monkey_01_03.mhd'
    output_filename = 'monkey_lddmm_output.mha'
    deformation_filename = 'monkey_lddmm_deformation.mha'
    config_filename = 'test_lddmm_config.json'
    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'relaxation', '--config', config_filename,
        '--logLevel', '4', _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title("LDDMM Generated Deformation Field")
    data = sitk.GetArrayFromImage(sitk.ReadImage(deformation_filename))
    plt.quiver(data[0][2:-2][2:-2], data[1][2:-2][2:-2], pivot='mid', color='r')

    plt.figure()
    plt.title('LDDMM Deformed Output')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('monkey_lddmm_output.eps')

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)
    plt.savefig('monkey_lddmm_difference.eps')

    plt.figure()
    plt.title('Source Image')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('monkey_lddmm_source.eps')

    plt.figure()
    plt.title('Target Image')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('monkey_lddmm_target.eps')
    plt.show()


def multi_circle():
    """Brad suggested running with a growing circle."""
    fixed_filename = 'circle.mha'
    moving_filename = 'blank.mha'
    output_filename = 'multi_circle.mha'
    deformation_filename = 'multi_circle_deformation.mha'
    config_filename = 'test_lddmm_config.json'

    img = np.zeros((32, 32), dtype=np.float)
    cx, cy = 16, 16  # The center of circle
    radius = 5
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, fixed_filename)

    cx, cy = 16, 16  # The center of circle
    radius = 10
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, moving_filename)

    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'multirelaxation', '--config', config_filename,
        '--logLevel', '4', _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title("Multi-LDDMM Generated Deformation Field")
    data = sitk.GetArrayFromImage(sitk.ReadImage(deformation_filename))
    plt.quiver(range(0, 32), range(0, 32), data[0], data[1], pivot='mid', color='r')

    plt.figure()
    plt.title('Multi-LDDMM Output after Deformation')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_multi_output.eps')

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)
    plt.savefig('cirle_multi_difference.eps')

    plt.figure()
    plt.title('Fixed')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_source.eps')

    plt.figure()
    plt.title('Moving')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_target.eps')
    plt.show()


def multi_circle_2():
    """Brad suggested running with a growing circle."""
    fixed_filename = 'circle.mha'
    moving_filename = 'blank.mha'
    output_filename = 'multi_circle.mha'
    deformation_filename = 'multi_circle_deformation.mha'
    config_filename = 'test_lddmm_config.json'

    img = np.zeros((32, 32), dtype=np.float)
    cx, cy = 16, 16  # The center of circle
    radius = 5
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, fixed_filename)

    cx, cy = 16, 16  # The center of circle
    radius = 10
    y, x = np.ogrid[-radius: radius, -radius: radius]
    index = x ** 2 + y ** 2 <= radius ** 2
    img[cy - radius:cy + radius, cx - radius:cx + radius][index] = 255
    itk_img = sitk.GetImageFromArray(img)
    sitk.WriteImage(itk_img, moving_filename)

    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'multirelaxation',  # '--config', config_filename,
        '--logLevel', '4', _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title("Multi-LDDMM Generated Deformation Field")
    data = sitk.GetArrayFromImage(sitk.ReadImage(deformation_filename))
    plt.quiver(range(0, 32), range(0, 32), data[0], data[1], pivot='mid', color='r')

    plt.figure()
    plt.title('Multi-LDDMM Output after Deformation')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_multi_output.eps')

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)
    plt.savefig('cirle_multi_difference.eps')

    plt.figure()
    plt.title('Fixed')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_source.eps')

    plt.figure()
    plt.title('Moving')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('cirle_lddmm_target.eps')
    plt.show()


def monkey_multi():
    """Running Multi-LDDMM on monkey brains"""
    fixed_filename = 'monkey_01_02.mhd'
    moving_filename = 'monkey_01_03.mhd'
    output_filename = 'monkey_multi_output.mha'
    deformation_filename = 'monkey_multi_deformation.mha'
    config_filename = 'test_lddmm_config.json'
    for line in lddmm(fixed_filename, moving_filename, deformation_filename,
        '--solver', 'multirelaxation',  # '--config', config_filename,
        '--logLevel', '4', _iter=True):
        print line

    apply_map(deformation_filename, fixed_filename, output_filename)

    plt.figure()
    plt.title("LDDMM Generated Deformation Field")
    data = sitk.GetArrayFromImage(sitk.ReadImage(deformation_filename))
    plt.quiver(data[0][2:-2][2:-2], data[1][2:-2][2:-2], pivot='mid', color='r')

    plt.figure()
    plt.title('Multi-LDDMM Deformed Output')
    data = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('monkey_multi_output.eps')

    plt.figure()
    plt.title('Difference Between Input and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_input = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(generated_output - given_input, cmap=cm.Greys_r)

    plt.figure()
    plt.title('Difference Between Moving and Deformation')
    generated_output = sitk.GetArrayFromImage(sitk.ReadImage(output_filename))
    given_moving = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(generated_output - given_moving, cmap=cm.Greys_r)
    plt.savefig('monkey_multi_difference.eps')

    plt.figure()
    plt.title('Source Image')
    data = sitk.GetArrayFromImage(sitk.ReadImage(fixed_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('monkey_lddmm_source.eps')

    plt.figure()
    plt.title('Target Image')
    data = sitk.GetArrayFromImage(sitk.ReadImage(moving_filename))
    plt.imshow(data, cmap=cm.Greys_r)
    plt.savefig('monkey_lddmm_target.eps')
    plt.show()

if __name__ == '__main__':
    main()